Re: NIntegrate and speed
- To: mathgroup at smc.vnet.net
- Subject: [mg116787] Re: NIntegrate and speed
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Mon, 28 Feb 2011 04:59:54 -0500 (EST)
----- 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