Re: NIntegrate and speed
- To: mathgroup at smc.vnet.net
- Subject: [mg116826] Re: NIntegrate and speed
- From: "Kevin J. McCann" <Kevin.McCann at umbc.edu>
- Date: Wed, 2 Mar 2011 04:31:03 -0500 (EST)
- References: <ikihpt$7tj$1@smc.vnet.net>
A couple of comments. (1) Bobby, your times are amazingly fast. I am
running a 3.6 GHz quad-core i5 processor with 4G of memory and my times
are quite a bit longer than yours; (2) my timings with ff done
numerically are about twice as as fast as the symbolic:
(* Symbolic evaluation of ff ala Bobby *)
R = 8000; Z = 1; rd = 3500;
Timing[Clear[ff];
ff[t_] = Integrate[Cos[t R Sin[\[Theta]]], {\[Theta], 0, \[Pi]},
Assumptions -> {t >= 0}]; {ff[t],
i2 = NIntegrate[(
Exp[-k Abs[Z]] ff[k])/(1 + (k rd)^2)^1.5, {k, 0, \[Infinity]}]}]
{22.907, {\[Pi] BesselJ[0, 8000 t], 0.000424068}}
(* Numerical evaluation of ff ala Kevin *)
R = 8000; Z = 1; rd = 3500;
Clear[ff];
ff[t_?NumericQ] :=
NIntegrate[Cos[t R Sin[\[Theta]]], {\[Theta], 0, \[Pi]}];
Timing[{ff[t],
i2 = NIntegrate[(
Exp[-k Abs[Z]] ff[k])/(1 + (k rd)^2)^1.5, {k, 0, \[Infinity]}]}]
{13.484, {ff[t], 0.000424067}}
I guess it depends a lot on the computer. Incidentally, I am running
64-bit Mathematica 8 under 64-bit Vista.
Kevin
On 3/1/2011 5:29 AM, DrMajorBob wrote:
> The inner integration is much faster done symbolically:
>
> Timing[Clear[ff];
> ff[t_] =
> Integrate[Cos[t*R*Sin[\[Theta]]], {\[Theta], 0, \[Pi]},
> Assumptions -> {t>= 0}];
> {ff[t], i2 =
> NIntegrate[
> Exp[-k*Abs[Z]]/(1 + (k*rd)^2)^1.5*ff[k], {k, 0, \[Infinity]}]}
> ]
>
> {0.725109, {\[Pi] BesselJ[0, 8000 t], 0.000424068}}
>
> Bobby
>
> On Mon, 28 Feb 2011 03:59:54 -0600, Daniel Lichtblau<danl at wolfram.com>
> wrote:
>
>>
>>
>> ----- Original Message -----
>>> From: "Marco Masi"<marco.masi at ymail.com>
>>> To: mathgroup at smc.vnet.net
>>> Sent: Sunday, February 27, 2011 3:35:46 AM
>>> Subject: NIntegrate and speed
>>> 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
>>
>> You presented it as an iterated integral. Mathematically that is fine
>> but from a language semantics viewpoint you now have a small problem. It
>> is that the outer integral cannot correctly do symbolic analysis of its
>> integrand but it may try to do so anyway. In essence, the outer
>> integrand "looks" to be nonnumerical until actual values areplugged in
>> for the outer variable of integration.
>>
>> There are (at least) two ways to work around this. One is to recast as a
>> double (as opposed to iterated) integral.
>>
>> Timing[i1 =
>> NIntegrate[
>> Exp[-k*Abs[Z]]/(1 + (k*rd)^2)^(3/2)*
>> Cos[k*R*Sin[\[Theta]]], {\[Theta], 0, \[Pi]}, {k,
>> 0, \[Infinity]}]]
>> {39.733, 0.0004240679194556497}
>>
>> An alternative is to define the inner function as a black box that only
>> evaluates for numeric input. In that situation the outer NIntegrate will
>> not attempt to get cute with its integrand.
>>
>> ff[t_?NumericQ] :=
>> NIntegrate[Cos[t* R*Sin[\[Theta]]], {\[Theta], 0, \[Pi]}]
>>
>> In[90]:= Timing[
>> i2 = NIntegrate[
>> Exp[-k* Abs[Z]]/(1 + (k* rd)^2)^1.5 *ff[k], {k, 0, \[Infinity]}]]
>> Out[90]= {26.63, 0.0004240673399701612}
>>
>>
>>> 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?
>>
>> Use the correct symbolic form of the inner integral. It involves BesselJ
>> rather than BesselJZero.
>>
>> In[91]:= ff2[t_] =
>> Integrate[Cos[t* Sin[\[Theta]]], {\[Theta], 0, \[Pi]},
>> Assumptions -> Element[t, Reals]]
>> Out[91]= \[Pi] BesselJ[0, Abs[t]]
>>
>> In[92]:= Timing[
>> i3 = NIntegrate[
>> Exp[-k Abs[Z]]/(1 + (k *rd)^2)^(3/2)* ff2[k*R], {k,
>> 0, \[Infinity]}]]
>> Out[92]= {0.7019999999999982, 0.0004240679192434893}
>>
>> Not surprisingly this is much faster, and will help to get you past the
>> speed bumps you allude to below.
>>
>>>
>>> 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.
>>
>> Compile will not help because most of the time will be spent in
>> NIntegrate code called from the virtual machine of the run time library
>> (that latter if you compile to C). Evaluating in parallel should help.
>> Also there might be option settings that allow NIntegrate to handle this
>> faster than by default but without significant degradation in quality of
>> results. Here is a set of timings using a few different methods, and
>> have PrecisionGoal set fairly low (three digits).
>>
>> In[109]:= Table[
>> Timing[NIntegrate[
>> Exp[-k Abs[Z]]/(1 + (k *rd)^2)^(3/2)* \[Pi] BesselJ[0,
>> Abs[k*R]], {k, 0, \[Infinity]}, PrecisionGoal -> 3,
>> Method -> meth]], {meth, {Automatic, "DoubleExponential",
>> "Trapezoidal", "RiemannRule"}}]
>>
>> During evaluation of In[109]:= NIntegrate::ncvb: NIntegrate failed to
>> converge to prescribed accuracy after 9 recursive bisections in k near
>> {k} = {0.0002724458978988764}. NIntegrate obtained
>> 0.00042483953211734914` and 0.000012161444876769691` for the integral
>> and error estimates.>>
>>
>> Out[109]= {{0.6709999999999923,
>> 0.0004240678889181539}, {0.0150000000000432,
>> 0.0004240644189596502}, {0.03100000000000591,
>> 0.0004240644189596502}, {0.04699999999996862,
>> 0.0004248395321173491}}
>>
>> I rather suspect there are more option tweaks that could make this
>> faster still without appreciable degradation in quality of results.
>>
>> Daniel Lichtblau
>> Wolfram Research
>>
>
>