{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "" 0 21 "" 0 1 0 0 0 1 0 0 0 0 2 0 0 0 0 1 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 
2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 69 "Compute the L-V cu
rve for water as predicted by the Van der Waals EoS" }}{PARA 0 "" 0 "
" {TEXT -1 53 "by minimization of the Helmoltz Free Energy A(V,T,n):" 
}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 21 "Van der
 Waals EoS is:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 26 "pvdw := nRT/(V-b) - a/V^2;" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 85 "We are concerned with isotherma
l minimization and our EoS predicts the properties of " }}{PARA 0 "" 
0 "" {TEXT -1 35 "both phase, therefore we can write " }}{PARA 0 "" 0 
"" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 17 "A = A0 - int(p,V)" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "where A0 \+
is an arbitrary constant that we can set to 0" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "A := -int(pvdw,V);
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 73 "Next \+
we need the Van der Waals constants, we can evaluate these by making" 
}}{PARA 0 "" 0 "" {TEXT -1 85 "use of the fact that at the critical po
int the first and second derivatives of p with" }}{PARA 0 "" 0 "" 
{TEXT -1 82 "respect to V must vanish. And noting that the critical po
int is at 374 C, 212 bar," }}{PARA 0 "" 0 "" {TEXT -1 36 "with R = 8.3
14 J/mol, and Vcrit ~ 3." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 30 "avdw := solve(diff(pvdw,V),a):" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 52 "bvdwn := solve(subs (a=avdw, V=3,diff(pvdw,V,V
)),b);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "avdwn := subs(V=3,nRT=(37
4+273)*8.314,b=bvdwn,avdw);" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 71 "set up the minimization using pseudocompounds t
o represent A(V,n), note" }}{PARA 0 "" 0 "" {TEXT -1 40 "we specify th
e temperature the next line" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 26 "nRTK := (324 + 273)*8.314;" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 46 "pvdwn := subs(a=avdwn,b=bvdwn,nRT=nRTK, pvdw);
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "Anum := subs(a=avdwn,b=bvdwn,nR
T=nRTK, A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "plot(pvdwn,V=1.1..16
);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "plot(Anum,V=1.1...20);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "dVol := 00.2: pi := 0: Vol := 1.1:j
 := 'j':" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "while Vol < 20 do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "   pi := pi + 1:" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 17 "   nh2o[pi] := 1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 19 "   vh2o[pi] := Vol;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "   ah2
o[pi] := subs(V=Vol,Anum);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "   Vo
l := Vol + dVol;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "end do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "Atot := sum(ah2o[j]*alpha[j],j=1..p
i):" }{MPLTEXT 0 21 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "ntot := \+
sum(nh2o[j]*alpha[j],j=1..pi):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "v
tot := sum(vh2o[j]*alpha[j],j=1..pi):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 36 "constraints := \{vtot = 3, ntot = 1\}:" }{MPLTEXT 0 21 0 "" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "with(simplex):" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 20 "do the minimization:" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "ans
 := minimize( Atot, constraints, NONNEGATIVE ):" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 85 "filter through the result
s to get the non-zero alphas, the index points to the volume" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "i := 0: j
 := 0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "while i < pi do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "   i := i + 1:" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 28 "   if (rhs(ans[i]) > 0) then" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 17 "      j := j + 1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
24 "      rans[j] := ans[i];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "   \+
end if:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "end do:" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 16 "rans[1];rans[2];" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 88 "use the volumes identifie
d by the indexes on the alpha's to back calculate the pressure:" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "i1 \+
:= 5; i2 := 26;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "subs(V=vh2o[i1+1
],pvdwn);subs(V=vh2o[i2],pvdwn);" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 45 "next back calculate by solving A = n*mu -
 p*V" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 13 "with(linalg):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "avec := vect
or(2,[ah2o[i1],ah2o[i2]]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "mat :
= matrix(2,2,[-vh2o[i1],nh2o[i1],-vh2o[i2],nh2o[i2]]):" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 19 "linsolve(mat,avec);" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 87 "and lastly compare the roots of
 the vdw at the pressure obtained by minimization to the" }}{PARA 0 "
" 0 "" {TEXT -1 38 "volumes of the stable pseudocompounds:" }}{PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "solve(pvdwn
-483);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "vh2o[i1],vh2o[i2];" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 16 "refine estimate:" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "vmn
1 :=vh2o[i1-1];vmx1 :=vh2o[i1+1];dv1 := (vmx1-vmn1)/20;" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 58 "vmn2 :=vh2o[i2-1];vmx2 :=vh2o[i2+1];dv2 := (vm
x2-vmn2)/20;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "pi := 0: Vol := vmn
1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "while Vol < vmx1 do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "   pi := pi + 1:" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 17 "   nh2o[pi] := 1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 19 "   vh2o[pi] := Vol;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "   ah2
o[pi] := subs(V=Vol,Anum);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "   Vo
l := Vol + dv1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "end do:" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 12 "Vol := vmn2:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "while Vol < vmx2 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 16 "   pi := pi + 1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "   nh2o[p
i] := 1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "   vh2o[pi] := Vol;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "   ah2o[pi] := subs(V=Vol,Anum);" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "   Vol := Vol + dv2;" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 7 "end do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "j
 := 'j';" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "Atot := sum(ah2o[j]*alp
ha[j],j=1..pi):" }{MPLTEXT 0 21 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
38 "ntot := sum(nh2o[j]*alpha[j],j=1..pi):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 38 "vtot := sum(vh2o[j]*alpha[j],j=1..pi):" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 36 "constraints := \{vtot = 3, ntot = 1\}:" }}
{PARA 0 "" 0 "" {MPLTEXT 0 21 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 20 "do the minimization:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "ans := minimize( Atot, constraints,
 NONNEGATIVE ):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 15 "i := 0: j := 0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 16 "while i < pi do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "   i := i
 + 1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "   if (rhs(ans[i]) > 0) th
en" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "      j := j + 1;" }}{PARA 0 
"> " 0 "" {MPLTEXT 1 0 24 "      rans[j] := ans[i];" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 10 "   end if:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "en
d do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "rans[1];rans[2];" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 83 "re-assign i1 and i2 to th
e indices on alpha here, and repeat refinement if desired:" }}{PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "i1 := 12; i
2 := 30;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "subs(V=vh2o[i1],pvdwn);
subs(V=vh2o[i2],pvdwn);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "
" }}}}{MARK "2 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 
0 1 2 33 1 1 }