ABSTRACT (original article): We describe a simple method for tracking solutions of nonlinear equations f(u,α) = 0 through turning points (also known as limit or saddle-node bifurcation points). Our implementation makes use of symbolic software such as Mathematica to derive an exact system of nonlinear ODE equations to follow the solution path, using a parameterization closely related to arc length. We illustrate our method with examples taken from the engineering literature, including examples that involve nonlinear boundary value problems that have been discretized by finite difference methods. Since the code requirement to implement the method is modest, we believe the method is ideal for demonstrating continuation methods in the classroom. CITATION (original article): Higgins, B. G., & Binous, H. (2010). A Simple Method for Tracking Turning Points in Parameter Space. JOURNAL OF CHEMICAL ENGINEERING OF JAPAN, 43(12), 1035–1042. DOI: https://doi.org/10.1252/jcej.10we122
We present a systematic method for the determination of all molar volumes of the van der Waals equation of state. The proposed method is based on Homotopy. The reduced form of this cubic equation of state is given by: where , and are the reduced temperature, pressure and volume.Without loss of generality, let us choose =0.5 and =0.1 and determine the corresponding three values of . To do so, we consider the following Homotopy function: where . Clearly, when , one recovers the ideal gas law and the problem has an obvious solution since it is a linear equation! Homotopy consists in starting with the obvious solution for t=0 and changing t until , which will yield the needed molar volumes.
f()=+--
v
P
3
2
v
v
1
3
8
3
T
T
P
v
T
P
v
H,t=tg()+(1-t)f()
v
v
v
f()=-
v
P
v
8
3
T
t=0
t=1
In[]:=
P=0.1;T=0.5;
f[v_]:=(P+3/v^2)(v-1/3)-8/3T(*vanderWaalscubicequationofstate*)
NSolve[f[v]==0,v](*allthreerootsusingNSolve*)
Out[]=
{{v11.0287},{v2.23166},{v0.406301}}
g[v_]:=Pv-8/3T(*idealgaslaw*)
h[v_,t_]:=tf[v]+(1-t)g[v](*Homotopyfunction*)
MyGrad[fon_,var_List]:=D[fon,#]&/@var(*Jacobianmatrix*)
In[]:=
J=MyGrad[h[x,t],{x,t}]
Out[]=
0.1(1-t)+t0.1-+,0.+0.1+-+x-0.1x
6-+x
1
3
3
x
3
2
x
3
2
x
1
3
s is a continuation parameter
In[]:=
Js=J/.xx[s]/.tt[s]
Out[]=
0.1(1-t[s])+t[s]0.1-+,0.+0.1+-+x[s]-0.1x[s]
6-+x[s]
1
3
3
x[s]
3
2
x[s]
3
2
x[s]
1
3
In[]:=
Eq[1]=D[x[s],{s,1}]First[Drop[Js,1]]
Out[]=
′
x
3
2
x[s]
1
3
In[]:=
Eq[2]=D[t[s],{s,1}]-First[Drop[Js,{2}]]
Out[]=
′
t
6-+x[s]
1
3
3
x[s]
3
2
x[s]
WhenEvent keeps track of the molar volumes obtained for t=1
In[]:=
{sol,points}=Flatten[Reap@Quiet@NDSolve[{Eq[1],Eq[2],t[0]0,x[0]8/3T/P,WhenEvent[t[s]==1,{Sow[{t[s],Log[x[s]]}]}]},{x[s],t[s]},{s,-1000,1000}],1]
Out[]=
x[s]InterpolatingFunction
[s],t[s]InterpolatingFunction
[s],{{1.,2.4005},{1.,0.802745},{1.,-0.900661}}
In[]:=
data=Table[Evaluate@{t[s],x[s]}/.sol,{s,-35.5,50,0.01}];
All three molar volumes
In[]:=
points
Out[]=
{{1.,2.4005},{1.,0.802745},{1.,-0.900661}}
As expected we find three roots when t=1; log scale is used for clarity
In[]:=
plt2=Show[ListLogPlot[data,AspectRatio->1,AspectRatio->1,Frame->True,GridLines->Automatic,PlotStyle->Darker@Red,Joined->True],PlotRange->{{-3,4},Automatic},Epilog->{PointSize[0.025],Darker@Green,Point[points[[1]]],Point[points[[2]]],Point[points[[3]]],Dashed,Darker@Blue,Line[{{1,-20},{1,20}}]},FrameLabel->{Style["t",16,Italic],Style["molar volume",16]},ImageSize->450{1,1},PlotLabel->Row[{Style["molar volumes :",14,Darker@Red]," ",Style[E^points[[1,2]],Darker@Magenta,14]," ",Style[E^points[[2,2]],Darker@Green,14]," ",Style[E^points[[3,2]],Darker@Blue,14]}]]
Out[]=
Further Information
Further Information
For details about the procedure see Ref.1 as well as the link to the lecture by the late Prof. Swan at MIT: https://youtu.be/muFAQx5dUdU
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
All roots of the Van der Waals equation of state using homotopy
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, June 25, 2025
https://community.wolfram.com/groups/-/m/t/3484884
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, June 25, 2025
https://community.wolfram.com/groups/-/m/t/3484884