Re: Beginner--Help on using FindRoot to solve the system of equations
- To: mathgroup at smc.vnet.net
- Subject: [mg66222] Re: Beginner--Help on using FindRoot to solve the system of equations
- From: "ben" <benjamin.friedrich at gmail.com>
- Date: Thu, 4 May 2006 05:21:19 -0400 (EDT)
- References: <e39kql$cud$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi Xiaochu,
2d-findroot is numerically subtle,
i don't know why exactly it breaks down in your example,
but there are good textbooks on this.
But what you can do is first solve analytically for Dtilde1,
substitute this into the second equation, and
apply 1d-findroot for Dtilde2.
Code is below.
Bye
Ben
(* specify dR *)
dR:=1;
(* solve analytically for Dtilde1, substitute this into second equation
*)
eqn=mu1-mu2/.Solve[P1==P2,Dtilde1][[1]];
Plot[Re[eqn],{Dtilde2,0,10}];
(* solve this *)
sol1=FindRoot[{0\[Equal]eqn},{Dtilde2,2}]
sol2=FindRoot[{0\[Equal]eqn},{Dtilde2,4}]//Chop
(* get corresponding Dtilde1 *)
(Dtilde1/.Solve[P1\[Equal]P2,Dtilde1][[1]])/.sol1
(Dtilde1/.Solve[P1\[Equal]P2,Dtilde1][[1]])/.sol2
xiaochu at che.utexas.edu schrieb:
> This is a code to solve vapor-liquid equilibrium by van der Waals Equation of State.
>
> I don't why this code is not working, please help, thanks very much!
>
>
> (* Name of EOS *)
>
> EOSName = "Van der Waals";
>
> Ttilde = .;
>
> Dtilde = .;
>
> Z = 3/(3-Dtilde)-9*Dtilde/8/Ttilde;
>
> (* Related variables *)
>
> Psi1 = (-9*Dtilde)/(8*Ttilde) - Log[3-Dtilde];
>
> lnB = (1-Z)-Psi1;
>
> mu = -lnB+Log[Dtilde];
>
> P = Dr Ttilde Z;
>
> DtildeG=10^-14;
>
> DtildeL=2.91;
>
> Dtilde1=.;
>
> Dtilde2=.;
>
> mu1=mu/.Dtilde->Dtilde1;
>
> mu2=mu/.Dtilde->Dtilde2;
>
> P1=P/.Dtilde->Dtilde1;
>
> P2=P/.Dtilde->Dtilde2;
>
> Ttilde = 0.1;
>
> Result=
> FindRoot[
> {P1==P2,mu1==mu2},
> {Dtilde1,DtildeG},
> {Dtilde2,DtildeL},
> MaxIterations->1000,
> WorkingPrecision->16];
>
> Link to the forum page for this post:
> http://www.mathematica-users.org/webMathematica/wiki/wiki.jsp?pageName=Special:Forum_ViewTopic&pid=10101#p10101
> Posted through http://www.mathematica-users.org [[postId=10101]]