How do I fix this so the Iteration will converge?

function [v,err,count] = Clapeyrons(P)

% Calculate molar volume as predicted by the Ideal Gas Law equation
% using the Newton-Raphson algorithm with initial estimate
% determined by ideal gas law.

% Inputs:
R = 0.082057; % Ideal gas constant in L atm / K mol
T = 293; % Temperature in K
Tc = 416.90; % Critical temperature of Cl2 in K
Pc = 78.72918; % Critical Pressure of Cl2 in atm
a = ((R^2)*(Tc^(5/2)))/(9*Pc*(2^(1/3)-1));
b = (R*Tc*(2^(1/3)-1))/(3*Pc);
% P = 1; % Pressure in atm

% Outputs:
% v equals molar volume in L/mol as predicted by 
% equation for gas considered at temperature T and pressure P.
% err equals modulus of function evaluated at approximate root.
% count is number of iterations taken by Newton-Raphson algorithm.

% Version 1: created 24/04/19. Author: Paul Curran
% Version 2: created 06/02/20. Author: Savana Stewart

if (~isscalar(P)) || (~isreal(P)) || P <= 0
    error('Input argument P must be positive real scalar.')
end

Iteration_limit = 20; % maximum number of iterations permitted
Tolerance=10^-7; % maximum acceptable value for modulus of

A = (a*P)/((R^2)*(T^(5/2)));
B = (b*P)/(R*T);
v = R * T / P; 
Z = (P*v)/(R*T);
C = A-B-B^2;
f1 = [1 -1 C -A*B]; % f = (Z^3) - (Z^2) + (A-B-(B^2))*Z - (A*B)
f = polyval(f1,Z);
v = Z*R*T/P; 

for count = 1:Iteration_limit + 1
    % Terminate with error message if iteration limit exceeded:
    if count == Iteration_limit + 1
        error('Iteration limit reached. Iteration did not converge.') 
    end

    % Terminate iteration if function is sufficiently small at current
    % estimate
    if abs(f)<Tolerance
        break
    end
    %f_dash = 3*Z^2 - 2*Z + (A - B - (B^2));   % Evaluate derivative of f
    f_dash = [0 3 -2 C];
    Z = Z - (f1./f_dash);   % Newton-Raphson iteration
    
end
err = abs(f); % Error is magnitude of f(v) at final root estimate
end

%% 
% COPY THE FOLLOWING INTO THE COMMAND WINDOW
% P = [1 1.5 2 2.5 3 5 10 15 25 50 100];
% v = zeros(1,length(P));
% for i=1:length(P)
% [u,err,count]=Clapeyrons(P(i));
% v(i)=u;
% end 
% 
% P1=[1:0.1:100];
% V1=R*T./P1;
% plot(P1,V1)
% hold on
% plot(P,v,'or')

I only changed the format a little and it now says 'Iteration limit reached. Iteration did not converge.' as it says to do when the iteration does not converge. How can I fix it?

Answer 1

[Geoff Hayes]

Savana - I haven't run the above code and so haven't stepped through it with the debugger, but if you were to do so, what would you see? The for loop is exited early only when the following condition is satisfied

    % Terminate iteration if function is sufficiently small at current
    % estimate
    if abs(f)<Tolerance
        break
    end

but nowhere in the body of this loop does f ever change and so it makes sense that you go through each iteration of the loop and see the "Iteration limit reached. Iteration did not converge" message. I suspect that you will want to re-compute f somewhere in this loop with the new Z. But I don't think that will be enough since you will want to compute the derivative too. I think that your code should look more like the following (assuming that f1 and f_dash are the coefficients for the polynomials of the function and its derivative)

    Z = Z - polyval(f1,Z)/polyval(f_dash,Z); 

and you could see how close Z is said to zero on each iteration.