Mathematica 9 is now available
Student Support Forum
-----
Student Support Forum: 'Numerically solve coupled equations' topicStudent Support Forum > General > Archives > "Numerically solve coupled equations"

< Previous Comment | Next Comment >Help | Reply To Comment | Reply To Topic
Author Comment/Response
Bill Simpson
10/18/12 00:48am

Modifying your code slightly

In[1]:= Δ=4;
Ux=1;
Uy=1+Δ/3;
Uz=2 Δ;
e1=1/2;
e2=12;
f1=1/2;
f2=1-f1;
Solve[f1 ((e1-eff1)/(1+Da*(e1-eff1)))+f2 ((e2-eff1)/(1+Da*(e2-eff1)))==0,eff1]

Out[9]= {
{eff1 -> (2 + 25*Da - Sqrt[4 + 529*Da^2])/(4*Da)},
{eff1 -> (2 + 25*Da + Sqrt[4 + 529*Da^2])/(4*Da)}
}

Likewise for eff2,3.

Using these sidesteps the numerical process of trying to find a solution on every iteration.

Next I rearrange your code to try to sidestep the non-numeric problem.

In[20]:= fa[eff1_/; NumericQ[eff1], eff2_/; NumericQ[eff2], eff3_/; NumericQ[eff3]] := 1/(4*π)*NIntegrate[(Sin[θ]^3*Cos[ϕ]^2)/(Ux^2*((Sin[θ]^2*Cos[ϕ]^2*eff1)/Ux^2 + (Sin[θ]^2*Sin[ϕ]^2*eff2)/Uy^2 + (Cos[ϕ]^2*eff3)/Uz^2)), {ϕ, 0, 2 π}, {θ, 0, π}];
fb[eff1_/;NumericQ[eff1], eff2_/; NumericQ[eff2], eff3_/; NumericQ[eff3]] := 1/(4*π)*NIntegrate[(Sin[θ]^3*Sin[ϕ]^2)/(Uy^2*((Sin[θ]^2*Cos[ϕ]^2*eff1)/Ux^2 + (Sin[θ]^2*Sin[ϕ]^2*eff2)/Uy^2 + (Cos[ϕ]^2*eff3)/Uz^2)), {ϕ, 0, 2 π}, {θ, 0, π}];
fc[eff1_/; NumericQ[eff1], eff2_/; NumericQ[eff2], eff3_/; NumericQ[eff3]] := 1/(4*π)*NIntegrate[(Sin[θ]*Cos[ϕ]^2)/(Uz^2*((Sin[θ]^2*Cos[ϕ]^2*eff1)/Ux^2 + (Sin[θ]^2*Sin[ϕ]^2*eff2)/Uy^2 + (Cos[ϕ]^2*eff3)/Uz^2)), {ϕ, 0, 2 π}, {θ, 0, π}]; NMinimize[
Da = fa[eff1, eff2, eff3];
Db = fb[eff1, eff2, eff3];
Dc = fc[eff1, eff2, eff3];
neweff1 = (2+25*Da-Sqrt[4+529*Da^2])/(4*Da);
neweff2 = (2+25*Db-Sqrt[4+529*Db^2])/(4*Db);
neweff3 = (2+25*Dc-Sqrt[4+529*Dc^2])/(4*Dc);
Norm[neweff1-eff1]+Norm[neweff2-eff2] + Norm[neweff3-eff3], {eff1, eff2, eff3}
]

That has substituted the Rho expression inside each NIntegrate, included conditions so that NMinimize can't attempt a symbolic evaluation, which would make NIntegrate fail, and tries to minimize the change in eff1,2,3 as you go round your circular dependencies.

That is slow and it complains that your integrals may be oscillatory, which is not surprising, but if I haven't made a mistake then it may have dealt with your non-numeric issue.

That still leaves the issue of potential negative real parts, but perhaps you can see if you can get any solution at all out of this and then look into which of the two conjugate solutions to use for solving each of eff1,2,3.

URL: ,

Subject (listing for 'Numerically solve coupled equations')
Author Date Posted
Numerically solve coupled equations Daniel 10/17/12 4:23pm
Re: Numerically solve coupled equations Bill Simpson 10/18/12 00:48am
Re: Numerically solve coupled equations Daniel 10/23/12 01:29am
< Previous Comment | Next Comment >Help | Reply To Comment | Reply To Topic