Finding roots of a nonlinear equations

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

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.

Sign in to comment.

Answers (1)

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

Asked:

on 29 Dec 2012

Community Treasure Hunt

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

Start Hunting!