newton raphson methon and symbolic functions

Hello, I have to write a code that calculates how much the spherical tank must be filled to contain 30m3 if R=3m using Newton Raphson method. The volume is calculated according to the expression 𝑉=𝜋*ℎ^2* [3*𝑅−ℎ]/ 3. I have to print the aproximations and relative error. I have tried writing the code multipule times, but Ialways have errors with derivative of a function and/or solving the equation for the exact answer (functions diff and solve). This is my code so far:

R=3;
h=0:0.1:6;
V=pi.*h.^2.*(3*R-h)/3;
plot(h,V,'k');

V_desired=30;
x=3;
max_iterations=100;
tolerance=1e-6;

answ=solve('pi*h^2*(3*3-h/3==V_desired', h);
A(1,1)=('iteration');
A(1,2)=('aproximation');
A(1,3)=('relative error');

for iteration=1:max_iteration
    H=volume(x, R);
    HD=dvolume(x,R);
    x=x-H/HD;
    error=(answ-x)/answ*100;

    A(iteration+1, 1)=iteration;
    A(iteration+1, 2)=x;
    A(iteration+1, 3)=error;

    if error<tolerance
        break
    end
end

disp(A);

This is a file for function:

function vol=volume(h,R)
vol=pi*h^2*(3*R-h)/3;
end

And this one is for derivative of that function:

functiondvol=dvolumen(h,R)
dvol=diff(volume(h,R), h)
end

Any answer would help, thank you!

Answer 1

[Torsten]

Better differentiate your function with respect to h in advance and pass the result to "dvolume" so that you only need to substitute the actual h value.

Further I suggest using a while-loop instead of a for-loop and to add the computation of the error in volume(h), not only in h.

R=3;
h=0:0.1:6;
V=pi.*h.^2.*(3*R-h)/3;
plot(h,V,'k');

V_desired=30;
xold=3;
max_iterations=100;
tolerance=1e-6;

%answ=solve('pi*h^2*(3*R-h)/3==V_desired', h);
%A(1,1)=('iteration');
%A(1,2)=('aproximation');
%A(1,3)=('relative error');
A(1,1) = 0;
A(1,2) = xold;
A(1,3) = NaN;

for iteration=1:max_iterations
    H=volume(xold, R, V_desired);
    HD=dvolume(xold,R, V_desired);
    x=xold-H/HD;
    error=abs((xold-x)/xold)*100;

    A(iteration+1, 1)=iteration;
    A(iteration+1, 2)=x;
    A(iteration+1, 3)=error;

    xold = x;
    
    if error<tolerance
        break
    end
end

disp(A);

         0    3.0000       NaN
    1.0000    2.0610   31.2989
    2.0000    2.0270    1.6492
    3.0000    2.0269    0.0067
    4.0000    2.0269    0.0000

This is a file for function:

function vol=volume(h,R, V_desired)
vol=pi*h^2*(3*R-h)/3-V_desired;
end

And this one is for derivative of that function:

function dvol=dvolume(h,R, V_desired)
%dvol=diff(volume(h,R), h)
syms hsym 
dvol = diff(pi*hsym^2*(3*R-hsym)/3-V_desired,hsym);
dvol = double(subs(dvol,hsym,h));
end

Answer 2

[Hassaan]

    % Define the radius of the tank and the desired volume
    R = 3;
    v_desired = 30;

    % Initial guess for h
    h = 1;

    % Define tolerance and maximum number of iterations
    tolerance = 1e-6;
    max_iterations = 100;

    % Newton-Raphson iteration
    for iteration = 1:max_iterations
        current_volume = volume(h, R);
        current_derivative = dvolume(h, R);
        h_new = h - (current_volume - v_desired) / current_derivative;

        % Calculate relative error
        error = abs((h_new - h) / h_new);

        % Display current iteration, h value, and error
        fprintf('Iteration %d: h = %.6f, Error = %.6f\n', iteration, h_new, error);

        % Check if the error is within the tolerance
        if error < tolerance
            break;
        end

        % Update h for the next iteration
        h = h_new;
    end

Iteration 1: h = 2.376526, Error = 0.579218
Iteration 2: h = 2.037410, Error = 0.166445
Iteration 3: h = 2.026919, Error = 0.005176
Iteration 4: h = 2.026906, Error = 0.000007
Iteration 5: h = 2.026906, Error = 0.000000

    % Display the final height
    fprintf('\nFinal height h = %.6f after %d iterations\n', h, iteration);

Final height h = 2.026906 after 5 iterations

% Volume function
function vol = volume(h, R)
    vol = pi * h^2 * (3 * R - h) / 3;
end

% Corrected derivative function
function dvol = dvolume(h, R)
    syms hs % Define a symbolic variable for h
    vol = pi * hs^2 * (3 * R - hs) / 3; % Define the volume as a symbolic expression
    dvolSym = diff(vol, hs); % Take the symbolic derivative with respect to hs
    dvol = subs(dvolSym, hs, h); % Substitute the symbolic variable with the numeric value of h
    dvol = double(dvol); % Convert symbolic result to numeric
end

------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Answer 3

[Sulaymon Eshkabilov]

Here is another alt. solution:

R=3;
h=0:0.1:6;
V=pi.*h.^2.*(3*R-h)/3;
plot(h,V,'k');
IDX = find(ceil(V)==30);
hold on
plot(h(IDX), V(IDX), 'rd', 'MarkerSize', 13, 'MarkerFaceColor','m')
legend('V', 'Solution point of h')

V_desired=30;
x=0.5;     % Initial guess value
max_iterations=100;
tolerance=1e-6;
syms hh
EQN = pi*hh^2*(3*R-hh)/3==V_desired;
answ=solve(EQN, hh);
Answer = answ(2);
Answer = double(Answer);
disp(['Iters   ', '  Appr  ', ' Error '])

Iters     Appr   Error 

for iteration=1:max_iterations
    Error=(Answer-x)/Answer;
    A(iteration, :)=[iteration, x, abs(Error)];
    fprintf('%d %15.6f %15.7f \n', A(iteration,:));
   [Vol, dVol]=VOLUME(x);
    x=x-Vol/dVol;
    if abs(Error)<tolerance
        disp('Iteration is halted becasue error tolerance is reached')
        break
    end
end

1        0.500000       0.7533186 
2        3.714896       0.8327916 
3        1.975809       0.0252093 
4        2.027236       0.0001632 
5        2.026906       0.0000000 

Iteration is halted becasue error tolerance is reached

function [Vol, dVol]=VOLUME(h)
R = 3;
V_desired=30;
Vol=pi*h^2*(3*R-h)/3-V_desired;
syms hh 
dV=diff(pi*hh^2*(3*R-hh)/3-30-V_desired, hh);
dVol=double(subs(dV, hh,h));
end