Finding roots of a nonlinear equations
Show older comments
I am trying to find the spinodal points for methane. I set up the function 'PREOS', which have two roots. All parameters are constant, T(temperature) and a(f(T)) are constants also. V is variable who's value I am interested in finding at which 'PREOS' is zero. The parameters and the method that I am using to find the roots are given below. I am getting the same root, 11.1038, at different temperatures like 100,90. And the second root I can not find, I have tried different initial values guess. Any idea about finding all the real roots of nonlinear function on matlab?
function y=PREOS(V,R,T,a,b) y=[-R*T/(V-b)+(2*a*(V+b))/(V^2+2*V*b-b^2)];
% PR EOS for calculating spinodal points of Methane
% Methane Critial Properties, K, cm3 and bar w=0.012; Tc=190.6; Pc=45.99; R=83.14; % Input Temperature T=120; Tr=T/Tc; %EOS Parameters psi=0.45724; omega=0.07780; b=omega*R*Tc/Pc; alpha=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2; a=psi*alpha*(R^2)*(Tc^2)/Pc; x=fzero(@(V) PREOS(V,R,T,a,b),10)
1 Comment
Roger Stafford
on 29 Dec 2012
Why use 'fzero' which requires an estimate when you can use the direct solution to a quadratic equation with no estimate needed? The two roots are easily obtained in terms of a square root involving a, b, R, and T in a comparatively simple expression. If the numerator of either root involves a small difference between two large quantities, you can rationalize the numerator to get better accuracy.
Answers (1)
Roger Stafford
on 29 Dec 2012
For your equation
-R*T/(V-b)+(2*a*(V+b))/(V^2+2*V*b-b^2) = 0
the two roots are quite simple:
s = (R*T+sqrt(2*R^2*T^2-4*a*R*T+4*a^2))/(2*a-R*T);
V1 = b*s;
V2 = -b/s;
For improved accuracy I have rationalized the second root's numerator under the assumption that R and T are both positive. Note that these roots will always be of opposite sign.
Categories
Find more on Physics in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!