Re: NIntegrate and speed
- To: mathgroup at smc.vnet.net
- Subject: [mg116850] Re: NIntegrate and speed
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Wed, 2 Mar 2011 04:35:24 -0500 (EST)
- References: <ikihpt$7tj$1@smc.vnet.net> <4D6CD43D.8030408@umbc.edu>
- Reply-to: drmajorbob at yahoo.com
But actually... the symbolic timings are extremely variable, depending on
whether Mathematica has cached the Integrate result.
NIntegrate is fast, but Integrate can take 20 seconds.
Bobby
On Tue, 01 Mar 2011 14:10:19 -0600, DrMajorBob <btreat1 at austin.rr.com>
wrote:
> I may have posted the wrong code, but the only difference is using 3/2
> rather than 1.5. I see no reason that should make a huge difference, but
> it obviously does!
>
> The same change is no help in the double integral.
>
> These are timings on my machine (3-yr-old iMac 24/2.8/4GB/1TB/SD).
>
> 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]]/(1 + (k*rd)^2)^(3/2)*ff[k], {k, 0, \[Infinity]}]}]
>
> {0.712826, {\[Pi] BesselJ[0, 8000 t], 0.000424068}}
>
> (*Symbolic evaluation of ff ala Bobby*)
> 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]}]}]
>
> {20.7246, {\[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)^(3/2), {k,
> 0, \[Infinity]}]}]
>
> {34.8137, {ff[t], 0.000424067}}
>
> Bobby
>
> On Tue, 01 Mar 2011 05:10:53 -0600, Kevin J. McCann
> <Kevin.McCann at umbc.edu> wrote:
>
>> 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 mma8 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
>>>>
>>>
>>>
>
>
--
DrMajorBob at yahoo.com