secant method in matlab

Hello,

I had to write a code that finds how much the spherical tank must be filled to contaiIn 30m3 if R=3m using secant method. I also had to find the exact solution and errors that are calculated by the formula error=(exact-aprox)/exact*100.

This is my code

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

V_desired=30;
x1=4;
x2=3;
number_iterations=50;

syms hh
ans=solve(pi*hh^2*(3*R-hh)/3==V_trazeni,hh);
answ=ans(2);
answ=double(answ);

for iterations=1:number_iterations
    fx1=volume(x1,R,V_desired);
    fx2=volume(x2,R,V_desired);

    x3=(x1*fx2-x2*fx1)/(fx2-fx1);

    error=abs((rj-x3)/rj)*100;

    matrix(iteracija,1)=iterations;
    matrix(iteracija,2)=x3;
    matrix(iteracija,3)=error;

    x1=x2;
    x2=x3;
end

disp(matrix);

function vol = volume(h,R,V_desired)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
vol=pi*h^2*(3*R-h)/3-V_desired;
end

But I always get this as a result

    1.0000    2.0249    0.0000
    2.0000    2.0268    0.0000
    3.0000    2.0269    0.0000
    4.0000    2.0269    0.0000
    5.0000    2.0269    0.0000
    6.0000    2.0269    0.0000
    7.0000       NaN    0.0000
    8.0000       NaN    0.0000
    9.0000       NaN    0.0000
   10.0000       NaN    0.0000
   11.0000       NaN    0.0000
   12.0000       NaN    0.0000
   13.0000       NaN    0.0000
   14.0000       NaN    0.0000
   15.0000       NaN    0.0000
   16.0000       NaN    0.0000
   17.0000       NaN    0.0000
   18.0000       NaN    0.0000
   19.0000       NaN    0.0000
   20.0000       NaN    0.0000
   21.0000       NaN    0.0000
   22.0000       NaN    0.0000
   23.0000       NaN    0.0000
   24.0000       NaN    0.0000
   25.0000       NaN    0.0000
   26.0000       NaN    0.0000
   27.0000       NaN    0.0000
   28.0000       NaN    0.0000
   29.0000       NaN    0.0000
   30.0000       NaN    0.0000
   31.0000       NaN    0.0000
   32.0000       NaN    0.0000
   33.0000       NaN    0.0000
   34.0000       NaN    0.0000
   35.0000       NaN    0.0000
   36.0000       NaN    0.0000
   37.0000       NaN    0.0000
   38.0000       NaN    0.0000
   39.0000       NaN    0.0000
   40.0000       NaN    0.0000
   41.0000       NaN    0.0000
   42.0000       NaN    0.0000
   43.0000       NaN    0.0000
   44.0000       NaN    0.0000
   45.0000       NaN    0.0000
   46.0000       NaN    0.0000
   47.0000       NaN    0.0000
   48.0000       NaN    0.0000
   49.0000       NaN    0.0000
   50.0000       NaN    0.0000
   

What is the problem here and how do I fix it?

Thank you!

Answer 1

[Hassaan]

R = 3;
V_desired = 30; % Desired volume
x1 = 4; % Initial guess
x2 = 3; % Second initial guess
number_iterations = 50;

syms hh
% Solve for exact solution
eqn = pi*hh^2*(3*R-hh)/3 == V_desired;
ans = solve(eqn, hh);
answ = double(ans(2)); % Choose the appropriate root
rj = answ; % Exact solution

matrix = zeros(number_iterations, 3); % Preallocate matrix

for iterations = 1:number_iterations
    fx1 = volume(x1, R, V_desired);
    fx2 = volume(x2, R, V_desired);
    
    if (fx2 - fx1) == 0
        break; % Avoid division by zero
    end
    
    x3 = (x1*fx2 - x2*fx1) / (fx2 - fx1);
    error = abs((rj - x3)/rj)*100;
    
    matrix(iterations, 1) = iterations;
    matrix(iterations, 2) = x3;
    matrix(iterations, 3) = error;
    
    % Update for next iteration
    x1 = x2;
    x2 = x3;
    
    % Break if the solution has effectively converged
    if error < 1e-6
        matrix = matrix(1:iterations, :); % Trim the matrix
        break;
    end
end

disp(matrix);

    1.0000    2.0249    0.0980
    2.0000    2.0268    0.0071
    3.0000    2.0269    0.0000
    4.0000    2.0269    0.0000

function vol = volume(h, R, V_desired)
    % Calculate volume difference from desired
    vol = pi*h^2*(3*R-h)/3 - V_desired;
end

Remember, the choice of initial guesses (x1 and x2) significantly affects the convergence of the secant method, so make sure they are reasonable based on the physical problem. If the method still fails to converge, consider using different initial guesses or a more robust root-finding method.

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

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.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering