1
$\begingroup$

I am trying to perform numerically the following integral

$$\int_0^8\text{d}x\,\text{Re}\left[\frac{e^{-\frac{a^2}{2}-\frac{x^2}{2}} x^4 \sin (b x)\left(e^{-i c x} \text{erfc}\left(\frac{-c +i x }{\sqrt{2}}\right)+e^{i c x} \text{erfc}\left(\frac{c +i x }{\sqrt{2}}\right)\right)}{2 \left(4 d^2 x^2+9\right)^6}\right] $$

For $a=b=c=5,\,d=0.0002$ the form of the integrand is the following

Plot[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]])
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}],{x,0,10}, PlotPoints -> 50000, PlotRange -> All]

enter image description here

However, different ways of choosing the integration region give me (very) different values for the integral

NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) 
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x,0,8}, MaxRecursion -> 22,AccuracyGoal -> 12]

-5.2007*10^-12

NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]])
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x, 0, 4, 8}, MaxRecursion -> 22, AccuracyGoal -> 12]

7.25566*10^-12

NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) 
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x, 0, 2, 4, 6, 8}, MaxRecursion -> 22, AccuracyGoal -> 12]

9.74422*10^-16

NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]])
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x, 0, 1, 2, 3, 4, 5, 6, 7, 8}, MaxRecursion -> 22, AccuracyGoal -> 12]

9.88541*10^-16

Total[Table[NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]])
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x, a - 1, a}, MaxRecursion -> 22, AccuracyGoal -> 12], {a, 1, 8}]]

1.42922*10^-8

Total[Table[NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]])
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x, 2 (a - 1), 2 a}, MaxRecursion -> 22, AccuracyGoal -> 12], {a, 1, 4}]]

-5.15013*10^-8

Total[Table[  NIntegrate[Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c -1/2 a^2)
(Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]])
/. {b -> 5, d -> 0.001/5, c -> 5, a -> 5}], {x, 4 (a - 1), 4 a}, MaxRecursion -> 22,AccuracyGoal -> 12], {a, 1, 2}]]

2.14488*10^-9

The plot of the integrand looks far too easy to have so different results using different techniques. What may be going wrong? Is there an efficient way of computing the right value?

$\endgroup$
5
$\begingroup$

The problem is that you are not using large enough precision goal for an oscillatory integrand with very small absolute values over the integration range.

Look at the log plot -- there are more oscillations than the ones visible with Plot:

enter image description here

The remedy is to use higher precision goal. Try this:

NIntegrate[
 Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/
     2 E^(-(1/2) x^2 - I x c - 1/2 a^2) (Erfc[(I x - c)/Sqrt[2]] + 
      E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) /. {b -> 5, d -> 0.001/5, 
    c -> 5, a -> 5}], {x, 0, 4, 8}, MaxRecursion -> 22, 
 PrecisionGoal -> 12, 
 Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
   "SymbolicProcessing" -> 0}]

(* 9.88564*10^-16 *)

The use of the method option is not crucial, it just speeds things up.

What may be going wrong?

Let us look into the sampling points used in the integration. Below I use an internal package for this, but EvaluationMonitor,Sow, and Reap can be used instead.

Needs["Integration`NIntegrateUtilities`"]

Your first integral, with the simple range specification, finishes with applying the integration rule once:

NIntegrateSamplingPoints@
 NIntegrate[
  Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/
      2 E^(-(1/2) x^2 - I x c - 1/2 a^2) (Erfc[(I x - c)/Sqrt[2]] + 
       E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) /. {b -> 5, d -> 0.001/5, 
     c -> 5, a -> 5}], {x, 0, 8}, MaxRecursion -> 22, 
  AccuracyGoal -> 12]

enter image description here

Your second integral finishes with applying the integration rule twice:

NIntegrateSamplingPoints@
 NIntegrate[
  Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/
      2 E^(-(1/2) x^2 - I x c - 1/2 a^2) (Erfc[(I x - c)/Sqrt[2]] + 
       E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) /. {b -> 5, d -> 0.001/5, 
     c -> 5, a -> 5}], {x, 0, 4, 8}, MaxRecursion -> 22, 
  PrecisionGoal -> 12, AccuracyGoal -> 12]

enter image description here

Using a higher precision goal as suggested above produces much more sampling points:

NIntegrateSamplingPoints@
 NIntegrate[
  Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/
      2 E^(-(1/2) x^2 - I x c - 1/2 a^2) (Erfc[(I x - c)/Sqrt[2]] + 
       E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) /. {b -> 5, d -> 0.001/5, 
     c -> 5, a -> 5}], {x, 0, 4, 8}, MaxRecursion -> 22, 
  PrecisionGoal -> 12, 
  Method -> {"GlobalAdaptive", "SingularityHandler" -> None}]

enter image description here

$\endgroup$
  • $\begingroup$ Wow, such a detailed explanation, thank you very much! $\endgroup$ – Alex Apr 29 '16 at 7:17
  • $\begingroup$ Good, I glad you find the explanation useful. $\endgroup$ – Anton Antonov Apr 29 '16 at 15:20
2
$\begingroup$

I'd suggest a simpler approach that seems to work just fine:

  1. Use arbitrary precision input if you can. In your case, instead of d -> 0.001/5, use d -> 5/1000

    integrand = Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c - 1/2 a^2) 
                 (Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) 
                  /. {b -> 5, d -> 1/5000, c -> 5, a -> 5}];
    
  2. I am not sure why you were trying to calculate 50,000 points for the plot. Letting Plot do its thing seems to work fine here:

    Plot[integrand, {x, 0, 10}, PlotRange -> All]
    

plot

  1. When you set an AccuracyGoal, you should set WorkingPrecision to be at least as large as the AccuracyGoal you request, and generally larger (see docs, under "Details"). You wouldn't have been able to do that in your code, because your value for d was expressed in machine precision.

  2. Set an explicit WorkingPrecision for NIntegrate rather than setting an accuracy goal. I request 20 digits of precision here; you can play around with that value and you will notice that you actually don't require quite as many to receive a correct result; in addition, requiring more digits does not change the value returns (although the calculation takes much longer), which at least makes us hope that we obtained a good stable result.

    NIntegrate[integrand, {x, 0, 8}, WorkingPrecision -> 20]
    (* Out: 9.8856438192273040591*10^-16 *)
    
  3. Notice also that an accuracy goal value of 12 is insufficient and it would give you a worse result:

    NIntegrate[integrand, {x, 0, 8}, WorkingPrecision -> 20, AccuracyGoal -> 12]
    (* Out: 1.8125249203590589373*10^-12 *)
    
$\endgroup$
  • $\begingroup$ I understand, so I was confusing Accuracy with Precision. Thank you so much! $\endgroup$ – Alex Apr 29 '16 at 7:18

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.