MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: NIntegrate and speed

  • To: mathgroup at smc.vnet.net
  • Subject: [mg116823] Re: NIntegrate and speed
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Tue, 1 Mar 2011 05:23:59 -0500 (EST)

Mea culpa... yes, you did. I saw it in another post after writing mine  
(Bob Hanlon?), but didn't notice it in yours.

It's Short Attention Span Theater, here.

Bobby

On Mon, 28 Feb 2011 11:43:52 -0600, Daniel Lichtblau <danl at wolfram.com>  
wrote:

> 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
>
> Didn't I do that? Under item (2) somewhere?
>
> Daniel
>
>
>> 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: [mg116780] 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


  • Prev by Date: Re: Vector Runge-Kutta ODE solver with compilation?
  • Next by Date: problem to define a complex quantity
  • Previous by thread: Re: NIntegrate and speed
  • Next by thread: Re: NIntegrate and speed