Student Support Forum > General > > "Numerically solve coupled equations"

 Post Reply: Name: Email Address: Please send email when my message is replied to. Url (optional): Subject: Message: view original message? Attachment (optional): Please answer this: 5+4 =

 Original Message (ID '69005') By Bill Simpson: 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.