1

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.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
Menzel Wang
  • 95
  • 5
  • 2
    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 – Michael E2 Aug 11 '15 at 03:29
  • 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! – Michael E2 Aug 11 '15 at 03:30
  • 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. – Michael E2 Aug 11 '15 at 03:32

1 Answers1

1

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.

obsolesced
  • 501
  • 4
  • 12