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
evaluated: http://www.cs.cmu.edu/~adamchik/articles/integr/mier.nb
(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