Re: an integral containing BesselJ

• To: mathgroup at smc.vnet.net
• Subject: [mg67501] Re: an integral containing BesselJ
• From: ab_def at prontomail.com
• Date: Wed, 28 Jun 2006 03:52:39 -0400 (EDT)
• References: <e784gb\$fph\$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```dimmechan at yahoo.com wrote:
> In Mathematica 5.2 i took the result:
>
> In[1]:=
> \!\(Integrate[\(Log[x]\/\@\(1 + x\^2\)\) BesselJ[0, x], {x, 0, 8}] //
> Timing\)
>
> Out[1]=
> {46.171 Second,0}
>
> which is a special case of:
>
> In[2]:=
> \!\(Integrate[\(Log[a\ x]\/\@\(1 + x\^2\)\) BesselJ[0,
>              x], {x, 0, 8}] // Timing\)
>
> Out[2]=
> \!\({32.18800000000001`\ Second, BesselI[0, 1\/2]\ BesselK[0, 1\/2]\
> Log[a]}\)
>
> However performing the numerical integration gave the result:
>
> In[4]:=
> \!\(NIntegrate[\(Log[x]\/\@\(1 +
>           x\^2\)\) BesselJ[0, x], {x, 0, 8}, Method \[Rule]
>         Oscillatory] // Timing\)
>
> Out[4]=
> {0.75 Second,-0.997939}
>
> I also perfrorm the symbolic integration in Mathematica 4.0 and took a
> lengthy result, which agrees numerically with the obtained result from
> the numerical integration:
>
> In[1]:=
> \!\(Integrate[\(Log[a\ x]\/\@\(1 + x\^2\)\) BesselJ[0, x], {x,
>         0, \[Infinity]}] // Timing\)
>
> Out[1]=
> \!\(\*
>   RowBox[{"{",
>     RowBox[{\(18.594`\ Second\), ",",
>       RowBox[{\(1\/12\ \((\(-1\) +
>               HypergeometricPFQ[{1\/2}, {1, 1},
>                 1\/4])\)\ \((6\ EulerGamma\^2 - \[Pi]\^2 -
>               6\ EulerGamma\ Log[4])\)\),
>         "+", \(1\/12\ \((6\ EulerGamma\^2 - \[Pi]\^2 -
>               6\ EulerGamma\ Log[4] +
>               12\ BesselI[0, 1\/2]\ BesselK[0, 1\/2]\ Log[a])\)\), "+",
>
>         RowBox[{\(1\/\(2\ \@\[Pi]\)\),
>           RowBox[{"(",
>             RowBox[{\(Log[4]\), " ",
>               RowBox[{"(",
>
>                 RowBox[{\(1\/2\ \@\[Pi]\ \((\(-1\) +
>                         HypergeometricPFQ[{1\/2}, {1, 1}, 1\/4])\)\
> Log[4]\),
>                   "+",
>                   RowBox[{
>                     UnderoverscriptBox["\[Sum]", \(K\$14174 = 1\),
>                       InterpretationBox["\[Infinity]",
>                         DirectedInfinity[
>                         1]]], \(\(2\^\(\(-1\) - 2\ K\$14174\)\ Gamma[
>                             1\/2 + K\$14174]\ PolyGamma[0,
>                             1\/2 + K\$14174]\)\/Gamma[1 +
> K\$14174]\^3\)}], "+",
>
>                   RowBox[{
>                     UnderoverscriptBox["\[Sum]", \(K\$14174 = 1\),
>                       InterpretationBox["\[Infinity]",
>                         DirectedInfinity[
>                         1]]], \(-\(\(2\^\(\(-1\) - 2\ K\$14174\)\ Gamma[
>                               1\/2 + K\$14174]\ PolyGamma[0,
>                               1 + K\$14174]\)\/Gamma[1 +
> K\$14174]\^3\)\)}]}],
>                 ")"}]}], ")"}]}], "+",
>         FractionBox[
>           RowBox[{\((\(-2\)\ EulerGamma + Log[4])\), " ",
>             RowBox[{"(",
>
>               RowBox[{\(EulerGamma\ \@\[Pi]\ \((\(-1\) +
>                       HypergeometricPFQ[{1\/2}, {1, 1}, 1\/4])\)\),
> "+",
>                 RowBox[{
>                   UnderoverscriptBox["\[Sum]", \(K\$14447 = 1\),
>                     InterpretationBox["\[Infinity]",
>                       DirectedInfinity[
>                       1]]], \(\(4\^\(-K\$14447\)\ Gamma[
>                           1\/2 + K\$14447]\ PolyGamma[0,
>                           1 + K\$14447]\)\/Gamma[1 + K\$14447]\^3\)}]}],
>               ")"}]}], \(2\ \@\[Pi]\)], "+",
>         RowBox[{\(1\/\(3\ \@\[Pi]\)\),
>           RowBox[{"(",
>             RowBox[{\(1\/8\ \[Pi]\^\(5/2\)\ \((\(-1\) +
>                     HypergeometricPFQ[{1\/2}, {1, 1}, 1\/4])\)\),
>               "-", \(3\/8\ \@\[Pi]\ \((\(-1\) +
>                     HypergeometricPFQ[{1\/2}, {1, 1}, 1\/4])\)\
> Log[4]\^2\),
>               "+",
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(-\(\(3\ 2\^\(\(-3\) - 2\ K\$14601\)\ Gamma[
>                           1\/2 + K\$14601]\ Log[16]\ PolyGamma[0,
>                           1\/2 + K\$14601]\)\/Gamma[1 +
> K\$14601]\^3\)\)}], "+",
>
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(-\(\(3\ 2\^\(\(-3\) - 2\ K\$14601\)\ Gamma[
>                           1\/2 +
>                             K\$14601]\ PolyGamma[0, 1\/2 + \
> K\$14601]\^2\)\/Gamma[1 + K\$14601]\^3\)\)}], "+",
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(\(3\ 2\^\(\(-3\) - 2\ K\$14601\)\ Gamma[
>                         1\/2 + K\$14601]\ Log[16]\ PolyGamma[0,
>                         1 + K\$14601]\)\/Gamma[1 + K\$14601]\^3\)}], "+",
>
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(\(3\ 2\^\(\(-2\) - 2\ K\$14601\)\ Gamma[
>                         1\/2 + K\$14601]\ PolyGamma[0,
>                         1\/2 + K\$14601]\ PolyGamma[0,
>                         1 + K\$14601]\)\/Gamma[1 + K\$14601]\^3\)}], "+",
>
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(-\(\(3\ 2\^\(\(-3\) - 2\ K\$14601\)\ Gamma[
>                           1\/2 +
>                             K\$14601]\ PolyGamma[0, 1 +
> K\$14601]\^2\)\/Gamma[1 \
> + K\$14601]\^3\)\)}], "+",
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(-\(\(3\ 2\^\(\(-3\) - 2\ K\$14601\)\ Gamma[
>                           1\/2 + K\$14601]\ PolyGamma[1,
>                           1\/2 + K\$14601]\)\/Gamma[1 +
> K\$14601]\^3\)\)}], "+",
>
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$14601 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(\(3\ 2\^\(\(-3\) - 2\ K\$14601\)\ Gamma[
>                         1\/2 + K\$14601]\ PolyGamma[1,
>                         1 + K\$14601]\)\/Gamma[1 + K\$14601]\^3\)}]}],
> ")"}]}],
>         "+",
>         RowBox[{\(1\/\(3\ \@\[Pi]\)\),
>           RowBox[{"(",
>
>             RowBox[{\(3\/2\ EulerGamma\^2\ \@\[Pi]\ \((\(-1\) +
>                     HypergeometricPFQ[{1\/2}, {1, 1}, 1\/4])\)\),
>               "+", \(1\/8\ \[Pi]\^\(5/2\)\ \((\(-1\) +
>                     HypergeometricPFQ[{1\/2}, {1, 1}, 1\/4])\)\), "+",
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$15304 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(\(3\ 4\^\(-K\$15304\)\ EulerGamma\ Gamma[
>                         1\/2 + K\$15304]\ PolyGamma[0,
>                         1 + K\$15304]\)\/Gamma[1 + K\$15304]\^3\)}], "+",
>
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$15304 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(\(3\ 2\^\(\(-1\) - 2\ K\$15304\)\ Gamma[
>                         1\/2 +
>                           K\$15304]\ PolyGamma[0, 1 +
> K\$15304]\^2\)\/Gamma[1 + \
> K\$15304]\^3\)}], "+",
>               RowBox[{
>                 UnderoverscriptBox["\[Sum]", \(K\$15304 = 1\),
>                   InterpretationBox["\[Infinity]",
>                     DirectedInfinity[
>                     1]]], \(-\(\(3\ 2\^\(\(-2\) - 2\ K\$15304\)\ Gamma[
>                           1\/2 + K\$15304]\ PolyGamma[1,
>                           1 + K\$15304]\)\/Gamma[1 +
> K\$15304]\^3\)\)}]}],
>             ")"}]}]}]}], "}"}]\)
>
> In[2]:=
> N[%[[2]]]
>
> Out[2]=
> 0.0580341\[InvisibleSpace] + 0.0833333 (-12.6717 + 11.7973 Log[a])
>
> In[3]:=
> % /. a -> 1
>
> Out[3]=
> -0.997939
>
> Can anyone help me?

Here's a link to a paper that explains how such integrals can be
(also have a look at issac90.pdf). First we find the transforms of the
factors:

In[1]:= Integrate[Log[1/x]/Sqrt[1 + 1/x^2]*x^(s - 1), {x, 0, Infinity}]

In[2]:= Integrate[BesselJ[0, x]*x*x^(s - 1), {x, 0, Infinity}]

In[3]:= fs = Simplify[%%*%, -1 < s < 0]

Out[3]= (2^(-2 + s)*Gamma[-(s/2)]*Gamma[(1 + s)/2]^2*(PolyGamma[0,
-(s/2)] - PolyGamma[0, (1 + s)/2]))/(Sqrt[Pi]*Gamma[1/2 - s/2])

The original integral is equal to 1/(2*Pi*I) times the integral of fs
along the line Re[s] == -1/2:

In[4]:= 1/(2*Pi*I)*NIntegrate[fs, {s, -I*Infinity, -1/2, I*Infinity}]

Out[4]= -0.9979393746397847 - 2.6504622345529306*^-17*I

We can represent fs as a sum of the derivatives of gamma function
ratios:

In[5]:= fs == D[2^(-2 + s)*Gamma[-s/2 + a]*Gamma[(1 + s)/2]^2/
(Sqrt[Pi]*Gamma[(1 - s)/2]) -
1/2*2^(-2 + s)*Gamma[-s/2]*Gamma[(1 + s)/2 + a]^2/
(Sqrt[Pi]*Gamma[(1 - s)/2]),
a] /. a -> 0 // Simplify

Out[5]= True

Since the integral of the ratio of gamma functions is a Meijer
function, the answer is expressible in terms of the derivatives of
MeijerG:

In[6]:= D[1/(4*Sqrt[Pi])*MeijerG[{{1/2 - a}, {}}, {{0, 0}, {0}}, 1/4] -
1/(8*Sqrt[Pi])*MeijerG[{{1/2}, {}}, {{a, a}, {0}}, 1/4],
a] /. a -> 0 // Simplify

Out[6]= (1/(8*Sqrt[Pi]))*(-Derivative[{{0}, {}}, {{0, 1}, {0}},
0][MeijerG][{{1/2}, {}}, {{0, 0}, {0}}, 1/4] - Derivative[{{0}, {}},
{{1, 0}, {0}}, 0][MeijerG][{{1/2}, {}}, {{0, 0}, {0}}, 1/4] -
2*Derivative[{{1}, {}}, {{0, 0}, {0}}, 0][MeijerG][{{1/2}, {}}, {{0,
0}, {0}}, 1/4])

This is probably not very useful, because Mathematica cannot evaluate
N[%6], and FunctionExpand[%6] gives Indeterminate.

Using the same technique of differentiation with respect to a
parameter, we can express the integral in terms of the derivatives of
HypergeometricPFQ:

In[7]:= int = (1/24)*(-2*BesselI[0, 1/2]^2*(Pi^2 +
6*EulerGamma*(-EulerGamma + Log[4])) + (9*EulerGamma -
3*Log[16])*Derivative[{0}, {0, 1}, 0][HypergeometricPFQ][{1/2}, {1, 1},
1/4] + 3*Log[16]*Derivative[{1}, {0, 0}, 0][HypergeometricPFQ][{1/2},
{1, 1}, 1/4] + 3*Derivative[{1}, {0, 1}, 0][HypergeometricPFQ][{1/2},
{1, 1}, 1/4] - 3*Derivative[{2}, {0, 0}, 0][HypergeometricPFQ][{1/2},
{1, 1}, 1/4] - 15*EulerGamma*Derivative[{0, 1}, {0, 0, 0},
0][HypergeometricPFQ][{1/2, 1}, {1, 1, 1}, 1/4] - 12*Derivative[{0, 1},
{0, 0, 1}, 0][HypergeometricPFQ][{1/2, 1}, {1, 1, 1}, 1/4] -
3*Derivative[{0, 2}, {0, 0, 0}, 0][HypergeometricPFQ][{1/2, 1}, {1, 1,
1}, 1/4] + 9*Derivative[{1, 1}, {0, 0, 0}, 0][HypergeometricPFQ][{1/2,
1}, {1, 1, 1}, 1/4]);

Mathematica seems to get the approximations for mixed partial
derivatives wrong:

In[8]:= N[Derivative[{1, 1}, {0, 0, 0}, 0][
HypergeometricPFQ][{1/2, 1}, {1, 1, 1}, 1/4],
20]

Out[8]= 4.523987519373309478408412473406`20.*^40

This looks like a bug, because even the simplest finite difference
formula can give a good approximation:

In[9]:= Block[{\$MaxExtraPrecision = 100},
int /.
Derivative[{si: (0)..., 1, si2: (0)..., 1, si3: (0)...},
{(0)..}, 0][HypergeometricPFQ][a_, b_, z_] :>
((f[h, k] - f[h, -k] - f[-h, k] + f[-h, -k])/(4*h*k) /.
f -> Function @@ {HypergeometricPFQ[
a + {si, #, si2, #2, si3}, b, z]}) /.
{h -> 10^-25, k -> 10^-25} //
N[#, 50]&
]

Out[9]=
-0.9979393746360745614717131506926624475792306906555576417821764083894`50.

Generalizing to Log[a*x] is easy:

Integrate[Log[a*x]/Sqrt[1 + x^2]*BesselJ[0, x], {x, 0, Infinity}] ==

Integrate[(Log[x] + Log[a])/Sqrt[1 + x^2]*BesselJ[0, x],

{x, 0, Infinity}] ==

int + BesselI[0, 1/2]*BesselK[0, 1/2]*Log[a]

Maxim Rytin
m.r at inbox.ru

```

• Prev by Date: Re: solving an equation with sums
• Next by Date: Re: solving an equation with sums
• Previous by thread: an integral containing BesselJ
• Next by thread: Scale ContourPlot axes