I've tried the answers in similar posts but they don't seem to work. As per title, I need to double integrate a complicated quickly oscillatory function. I've checked and there are no poles, the function is well behaved and falls to 0 quickly. The function has three variables: r, td and k. I want to find a plot in terms of k and integrate r and td. I need to integrate td from 0 to infinity and r from td to infinity (thus the Boole[r > td] and both limits set from 0 to infinity). I've tried different integration methods like QuasiMonteCarlo, which yields some result but with a ton of error, specially for big k, or LevinRule which is the most natural but yields an error and an absurd result (a super big number to the power of a super big number). The error is:

XXX is a Levin function of differential order 72 which exceeds value of \option "MaxOrder" -> 50. Treating XXX as a non-Levin function

Where XXX is a long expression related to my input (but weirdly changed in some places)

The code I'm using (for the time being for a given k) is

F0[td_, r_] := 2 (r^2 - td^2)^2 (r^2 + 6 r + 12);

F1[td_, r_] := 
  2 (r^2 - td^2) (-r^2 (r^3 + 4 r^2 + 12 r + 24) + 
     td^2 (r^3 + 12 r^2 + 60 r + 120));

F2[td_, r_] := 
  1/2 (r^4 (r^4 + 4 r^3 + 20 r^2 + 72 r + 144) - 
     2 td^2 r^2 (r^4 + 12 r^3 + 84 r^2 + 360 r + 720) + 
     td^4 (r^4 + 20 r^3 + 180 r^2 + 840 r + 1680));

Itdr[td_, r_] := 
  Exp[td/2] + Exp[-td/2] + (td^2 - r^2 - 4 r)/(4 r) Exp[-r/2];

Integrand[td_, r_] := 
  k^3/(12 \[Pi]) (Exp[-r/2] Cos[k td])/(
   r^3 Itdr[td, r]) (SphericalBesselJ[0, k r] F0[td, r] + 
     SphericalBesselJ[1, k r]/(k r)  F1[td, r] + 
     SphericalBesselJ[2, k r]/(k r)^2 F2[td, r]);

 k = 0.1;
limit = Infinity;
 Boole[rd > td] Integrand[td, rd], {rd, 0, limit}, {td, 0, limit}, Method -> "LevinRule"]

I've also tried changing the limits to be finite (after all, the function drops quickly) but this doesn't work specially well. Any idea on what I should try next? Any help is greatly appreciated.

Edit: some plots of the function. For 'small' k=1/10 the function has this form (for different values of td)

k=1/10, td=0 k=1/10, td=10 k=1/10, td = 100

You can see that it converges nice and easy. For higher values of k the oscillation is super fast. For k=100:

k=100, td =0 k=100, td=10

It still converges, but the fast oscillation makes it hard to find a reasonable integration

  • $\begingroup$ Could you include a plot? $\endgroup$ May 29 '20 at 15:40
  • $\begingroup$ Sure, I have edited the initial question $\endgroup$
    – Salva
    Jun 2 '20 at 9:04

I tried the following:

  1. Change the definition of your integrand function to use Set rather than SetDelayed, i.e. Integrand[td_, r_] = ...; there is no need for that static expression to be continuously re-evaluated.
  2. Rather than with Boole, write your limits of integration explicitly: NIntegrate[Integrand[td, rd], {td, 0, limit}, {rd, td, limit}].
  3. Try the default automatically selected methods before LevinRule which may be slower.

With those in mind, and using your definitions for the rest of the function:

Integrand[td_, r_] = k^3/(12 Pi) (Exp[-r/2] Cos[k td])/(r^3 Itdr[td, r]) (SphericalBesselJ[0, k r] F0[td, r] + SphericalBesselJ[1, k r]/(k r) F1[td, r] + SphericalBesselJ[2, k r]/(k r)^2 F2[td, r]);

k = 1/10;

  {limit, NIntegrate[Integrand[td, rd], {td, 0, limit}, {rd, td, limit}]}, 
  {limit, {1, 5, 100, 200, 300, 600, 1000}}

As you mentioned, the integrand appears pretty well-behaved and converges quickly, but higher values of limit make the calculation painfully slow. Perhaps you don't need them though:

  PlotRange -> All, Joined -> True,
  Mesh -> All, MeshStyle -> Directive[PointSize[0.02], Red]

plot of value of integral vs. increasing limit of integration

It is typically preferable to avoid names starting with capital letters for user-defined quantities (e.g. Integrand) to clearly distinguish them from the built-ins, which are all capitalized.

  • $\begingroup$ Thanks a lot for your answer! Why the F functions should be defined with " := " while Integrand should be defined with "=" ? I agree that since the integral converges quickly I can avoid large values of limit. However since I need to do a plot over k and the higher the value of k the slower the convergence I'll have to do this carefully. I just tried your suggestions and it's taking a while. How much time did it took for you to run your code? And last question, if I increase k the performance gets worse in terms of time and precission. Any suggestion that can help? $\endgroup$
    – Salva
    Jun 2 '20 at 8:59
  • $\begingroup$ @Salva In your case it doesn’t really matter whether you define F with Set or SetDelayed. Ultimately the fact that Integrando is defined with Set means that it will be called only once anyway. I just didn’t want to change the code too much with minute changes that would not affect the performance anyway. $\endgroup$
    – MarcoB
    Jun 2 '20 at 13:41
  • $\begingroup$ @Salva It may have taken a couple of minutes to generate the plot I showed. It’s definitely slow for high values of k. I am not sure if there’s anything to be done to improve performance for high k. As you said, your function becomes particularly difficult to integrate due to its oscillatory nature at high k. You should try a few different integration methods and see if you chance upon one that goes faster. $\endgroup$
    – MarcoB
    Jun 2 '20 at 13:46
  • $\begingroup$ @Salva Since the integral value is pretty small for high k, I’d also suggest setting a finite AccuracyGoal in NIntegrate (say, AccuracyGoal -> 6). The default setting of PrecisionGoal may cause NIntegrate to have a hard time deciding that a solution has converged when the value of the integral is (close to) zero. Setting an explicit accuracy requirement allows it to call the solution “good enough” and move on. Worth a try. $\endgroup$
    – MarcoB
    Jun 2 '20 at 13:48
  • $\begingroup$ I'll try with AccuracyGoal. I had already tried with several methods but do not yield satisfactory results. For example, montercarlo methods give a big error while levinrule, which I expected to be the best method here, spit the error shown above. I have been able to produce reasonable results playing with the limits up to k = 1000, which is already good enough. However I'm struggling with k=10^4. The numerical error is higher than the result. But I'll play around with accuracygoal and if I cannot make it work I'll call it a day and use results up to k=1000. Once again, thanks a lot!! $\endgroup$
    – Salva
    Jun 2 '20 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.