Re: NIntegrate HypergeometricU
- To: mathgroup at smc.vnet.net
- Subject: [mg84867] Re: NIntegrate HypergeometricU
- From: "David W.Cantrell" <DWCantrell at sigmaxi.net>
- Date: Wed, 16 Jan 2008 23:10:54 -0500 (EST)
- References: <fmkf3a$9fg$1@smc.vnet.net>
"subo_dha at yahoo.com" <subodha.g at gmail.com> wrote:
> I am trying to NIntegrate the function HypergeometricU[a,b,z] with
> a>=1, b>=1,z=R/sin[x]^2 where integration is from 0 to Pi/2. For lower
> values of R this gives answers and for Higher values of R this is
> giving errors .
>
> \!\(NIntegrate::"singd" \(\(:\)\(\ \)\) "NIntegrate's singularity
> handling
> has failed at point \
> \!\({x}\)=\!\({8.145014137541`5.15492657977389*^-15}\) for the
> specified \
> precision goal. Try using larger values for any of $MaxExtraPrecision
> or the \
> options WorkingPrecision, or SingularityDepth and MaxRecursion"\)
>
> Even though I change the working precisio,SingularityDepth and
> MaxRecursion it is not helping me much. Please advice me how to solve
> this problem.
I suggest an alternative approach. Although Mathematica cannnot calculate a
symbolic result for your definite integral with respect to x, it succeeds
if we help it by making the simple substitution u = Sin[x]. Your integral
then becomes
In[8]:= Integrate[HypergeometricU[a, b, r/u^2]/Sqrt[1 - u^2], {u, 0, 1}]
Out[8]= If[Re[a] > -(1/2),
(1/2)*Sqrt[Pi]*((Sqrt[Pi]*Gamma[1 - b])/Gamma[1 + a - b] - (2*Gamma[1/2 +
a]*Gamma[1/2 - b]*HypergeometricPFQ[{1/2, 1/2 + a}, {3/2, 1/2 + b},
r])/(Sqrt[1/r]*Gamma[a]*Gamma[1 + a - b]) + ((1/r)^(-1 + b)*Gamma[-(1/2) +
b]*HypergeometricPFQ[{1 - b, 1 + a - b}, {3/2 - b, 2 - b}, r])/((-1 +
b)*Gamma[a])), Integrate[HypergeometricU[a, b, r/u^2]/Sqrt[1 - u^2], {u, 0,
1}, Assumptions -> Re[a] <= -(1/2)]]
In[9]:= sol = FullSimplify[%, a >= 1 && b >= 1 && r > 0]
Out[9]= (Pi*(2*Gamma[a]*Gamma[1 - b] - 4^b*r^(1 - b)*Gamma[2 - 2*b]*Gamma[1
+ a - b]*Gamma[-(1/2) + b]*HypergeometricPFQRegularized[{1 - b, 1 + a - b},
{3/2 - b, 2 - b}, r] - 2*Pi*Sqrt[r]*Gamma[1/2 +
a]*HypergeometricPFQRegularized[{1/2, 1/2 + a}, {3/2, 1/2 + b},
r]*Sec[b*Pi]))/(4*Gamma[a]*Gamma[1 + a - b])
I suspect that that result is basically correct. However, at least two
cautions should be noted.
Caution 1: Numerical evaluation of sol may require care. For example,
In[15]:= N[sol /. {a -> 2, b -> 5/3, r -> 100}]
Out[15]= -7.06858
is misleading. Instead, use
In[16]:= N[sol /. {a -> 2, b -> 5/3, r -> 100}, 12]
Out[16]= 0.0000576344250174
Caution 2: For some values of the parameters (integers and integers+1/2, I
think), sol is Indeterminate, even though your integral is well defined.
For example,
In[17]:= sol /. {a -> 2, b -> 3, r -> 2}
During evaluation of In[17]:= \[Infinity]::indet: Indeterminate \
expression 0/4 3 Sqrt[\[Pi]] ComplexInfinity ComplexInfinity \
encountered. >>
Out[17]= Indeterminate
although
In[18]:= NIntegrate[
HypergeometricU[a, b, r/Sin[x]^2] /. {a -> 2, b -> 3, r -> 2}, {x,0,Pi/2}]
Out[18]= 0.147262
You might wish to compare that with
In[19]:= N[sol /. {a -> 2, b -> 2999/1000, r -> 2}, 12]
Out[19]= 0.147180371332
In[20]:= N[sol /. {a -> 2, b -> 3001/1000, r -> 2}, 12]
Out[20]= 0.147343999426
I suspect that Mathematica might not be able to simplify
Limit[sol /. {a -> 2, r -> 2}, b -> 3], but I'm not sure
because I stopped waiting for a result.
David W. Cantrell