How do I fix this so the Iteration will converge?

12 次查看（过去 30 天）
Quincey Jones
Quincey Jones 2020-2-13
回答： Geoff Hayes 2020-2-13
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?

请先登录，再回答此问题。

回答（1 个）

Geoff Hayes
Geoff Hayes 2020-2-13
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.

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by