Re: NIntegrate and speed
- To: mathgroup at smc.vnet.net
- Subject: [mg116785] Re: NIntegrate and speed
- From: Andrew Moylan <amoylan at wolfram.com>
- Date: Mon, 28 Feb 2011 04:59:32 -0500 (EST)
You're right, the inner integral is a Bessel function:
In[1]:= Integrate[Cos[k R Sin[\[Theta]]], {\[Theta], 0, \[Pi]},
Assumptions -> {R > 0, k > 0}]
Out[1]= \[Pi] BesselJ[0, k R]
So, you can use a one-dimensional NIntegrate:
In[2]:= R = 8000; Z = 1;
rd = 3500;
In[4]:= NIntegrate[
Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 (\[Pi] BesselJ[0, k R]), {k,
0, \[Infinity]}] // Timing
Out[4]= {0.474709, 0.000424068}
You can try to speed up the result by choosing a different method.
In this case, NIntegrate selects "ExtrapolatingOscillatory" as the method:
In[5]:= NIntegrate[
Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 (\[Pi] BesselJ[0, k R]), {k,
0, \[Infinity]}, Method -> "ExtrapolatingOscillatory"] // Timing
Out[5]= {0.434378, 0.000424068}
The new "LevinRule" general oscillatory method is actually faster for this integral:
In[6]:= NIntegrate[
Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 (\[Pi] BesselJ[0, k R]), {k,
0, \[Infinity]}, Method -> "LevinRule"] // Timing
Out[6]= {0.135486, 0.000424068}
To do many NIntegrates in a loop, Parallelize would be the first thing to try:
In[7]:= Clear[R];
Parallelize[
Table[NIntegrate[
Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 (\[Pi] BesselJ[0, k R]), {k,
0, \[Infinity]}, Method -> "LevinRule"], {R, 6000, 9000, 1000}]]
Out[8]= {0.00052313, 0.000470404, 0.000424068, 0.000383576}
I hope this helps,
Andrew Moylan
Wolfram Research
On Feb 27, 2011, at 8:35 PM, Marco Masi wrote:
> I have the following problems with NIntegrate.
>
> 1) I would like to make the following double numerical integral converge without errors
>
> R = 8000; Z = 1; rd = 3500;
> NIntegrate[Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 (NIntegrate[Cos[k R Sin[\[Theta]]], {\[Theta], 0, \[Pi]}]), {k, 0, \[Infinity]}]
>
> It tells non numerical values present and I don't understand why, since it evaluates finally a numerical value? 0.000424067
>
> 2) Isn't the second integrand a cylindrical Bessel function of order 0? So, I expected that
> NIntegrate[Exp[-k Abs[Z]]/(1 + (k rd)^2)^1.5 BesselJZero[0, k R], {k, 0, \[Infinity]}] doing the same job. But it fails to converge and gives 0.00185584- i4.96939*10^-18. Trying with WorkingPrecision didn't make things better. How can this be fixed?
>
> 3) The above Nintegrals will go into a loop and should be evaluated as fast as possible. How? With Compile, CompilationTarget -> "C", Paralleization, etc.?
>
> Any suggestions?
>
> Marco.
>