1
$\begingroup$

I'm trying to solve two nonlinear algebra equations for two unknown parameters, U and Tf. Since some terms in these equations contain integration, and the integration also contains U and Tf. The main problem I have encountered is that Mathematica function NIntegrate cannot get numerical results, even in FindRoot.

My code is as following: 

(*Define functions*) DI1[x_] := E^((Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] (-R - x))/( 2 U) + ((K - U^2) (R + x) - 2 U Log[R + x])/(2 U)) HypergeometricU[( K Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] + U^2 Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4])/( K^2 + 2 K U^2 + 4 CC K U^2 + U^4), 0, ( Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] (R + x))/U] DI2[x_] := E^((Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] (-R - x))/( 2 U) + ((K - U^2) (R + x) - 2 U Log[R + x])/(2 U)) LaguerreL[-(K Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] + U^2 Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4])/(K^2 + 2 K U^2 + 4 CC K U^2 + U^4), -1, ( Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] (R + x))/U] I1[x_] := NIntegrate[DI1[s], {s, x, ∞}] I2[x_] := NIntegrate[DI2[s], {s, 0, x}] (*Define these parameters for simplicity*) a = Sqrt[K^2 + 2 K U^2 + 4 CC K U^2 + U^4] C1 := (K Q (K + U^2) Gamma[(K + U^2)/a])/(a U) C2 := (K Q (K + U^2) Gamma[(K + U^2)/a] DI1[0])/(a U DI2[0]) - ( K Tf (2 I1[0] - DI1[0]))/(U (I1[0] - DI1[0]) DI2[0]) C3 := (K Q (K + U^2) Gamma[(K + U^2)/a] I1[0])/(a U) - (Tf I1[0] )/( I1[0] - DI1[0] ) (*Define equations*) equ1[U_, Tf_] := -C1*(-U/K DI1'[0] + DI1[0]) + C2*(-U/K DI2'[0] + DI2[0]) - (-Tf DI1[0])/(I1[0] - U/K DI1[0]) == 1/Le Exp[-U*Le*R]/R^2 1/ Integrate[s^-2 Exp[-U*Le*s], {s, R, ∞}] equ2[U_, Tf_] := 1/Le Exp[-U*Le*R]/R^2 1/ Integrate[ s^-2 Exp[-U*Le*s], {s, R, ∞}] == (σ + (1 - σ) Tf)^2 Exp[ Z/2 (Tf - 1)/(σ + (1 - σ) Tf)] (*FInd roots of equ1 and equ2*) Block[{K = 1, CC = 2, Q = 1, σ = 0.15, Z = 10, Le = 1, R = 1, Tf = 0.8}, FindRoot[equ1[U, Tf], {U, 1}]]

(Edit) I think the main problem is the definition of I1[x_] := NIntegrate[DI1[s], {s, x, ∞}]. While FindRoot function is running, NIntegrate encounters non-numerical results. But if I change the definition from NIntegrate to Integrate, this integration cannot be calculated out.

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful $\endgroup$ Commented Aug 11, 2015 at 3:29
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ Commented Aug 11, 2015 at 3:30
  • $\begingroup$ It seems the standard greeting (above) had not been extended to you. In particular, don't forget to accept the answer, if any, that solves your problems. $\endgroup$ Commented Aug 11, 2015 at 3:32

1 Answer 1

Reset to default
1
$\begingroup$

I encountered a problem copying and pasting your code, as HypergeometricU and LaguerreL begin on new lines, but I believe you meant for them to be part of the definitions for DI1 and DI2. Fixing that, the problem of non-numerical integrand comes from undefined global U, which can be set inside Block, and FindRoot attempting to evaluate its argument symbolically, which can be prevented by restricting its argument's definition to numeric values:

block[u_,expr_]:=Block[{K=1,CC=2,Q=1,\[Sigma]=0.15,Z=10,Le=1,R=1,Tf=0.8,U=u},expr] SetAttributes[block,HoldRest] equ1[U_?NumericQ]:=block[U,equ1[U,Tf]]

However, preventing FindRoot from manipulating its argument symbolically removes its ability to handle equations, as evaluating the equations numerically returns only true or false, so redefine the equations as functions:

equ1[U_, Tf_] := -C1*(-U/K DI1'[0] + DI1[0]) + C2*(-U/K DI2'[0] + DI2[0]) - (-Tf DI1[0])/(I1[0] - U/K DI1[0]) - 1/Le Exp[-U*Le*R]/R^2 1/ Integrate[s^-2 Exp[-U*Le*s], {s, R, \[Infinity]}] equ2[U_, Tf_] := 1/Le Exp[-U*Le*R]/R^2 1/ Integrate[s^-2 Exp[-U*Le*s], {s, R, \[Infinity]}] - (\[Sigma] + (1 - \[Sigma]) Tf)^2 Exp[Z/2 (Tf - 1)/(\[Sigma] + (1 - \[Sigma]) Tf)]

Now calling FindRoot[equ1[u],{u,1}] returns a host of errors from divergent I1[0] and Integrate[s^-2 Exp[-U*Le*s], {s, R, \[Infinity]}], with the latter requiring U>0 since Le=1. For I1[0], we have

block[U, Integrate[DI1[s], {s, 0, Infinity}]] // FullSimplify

$$\int_0^{\infty } \frac{e^{-\frac{(s+1) \left(U^2+\sqrt{U^4+10 U^2+1}-1\right)}{2 U}} HypergeometricU\left(\frac{U^2+1}{\sqrt{U^4+10 U^2+1}},0,\frac{(s+1) \sqrt{U^4+10 U^2+1}}{U}\right)}{s+1} \, ds$$

which again requires U>0 to converge.

One way to restrict the domain of FindRoot to U>0 is

Abs[U] /. FindRoot[equ1[Abs[U]], {U, 1}]

0.209635

There is no real root for equ2:

block[U, equ2[U, Tf]] // FullSimplify

ConditionalExpression[-0.206495 + 1/(1 - E^U U Gamma[0, U]), Re[U] > 0]

which is positive for all U>0.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.