Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to NIntegrate this expression. It is very close to zero, and so gives all possible types of errors.

More precisely, I am actually solely interested in the Sign of the outcome. However, this sign changes for different options:

NIntegrate[(1/(108 π))E^(-((2 Pq^2)/9)-1/18 (53-16 Pq+2 Pq^2)^2) (12+
    E^(1/18 (53-16 Pq+2 Pq^2)^2) Sqrt[2 π] (9-4 Pq^2)-E^(1/18 (53-
     16 Pq+2 Pq^2)^2) Sqrt[2 π] (-9+4 Pq^2) Erf[(-53+16 Pq-2 Pq^2)/(3 Sqrt[2])]),
 {Pq,-∞,∞},
 WorkingPrecision -> 20, AccuracyGoal -> 10]
(* -1.7126730471588308508*10^-16 *)

NIntegrate[(1/(108 π))E^(-((2 Pq^2)/9)-1/18 (53-16 Pq+2 Pq^2)^2) (12+
    E^(1/18 (53-16 Pq+2 Pq^2)^2) Sqrt[2 π] (9-4 Pq^2)-E^(1/18 (53-
     16 Pq+2 Pq^2)^2) Sqrt[2 π] (-9+4 Pq^2) Erf[(-53+16 Pq-2 Pq^2)/(3 Sqrt[2])]),
 {Pq,-∞,∞},
 WorkingPrecision -> 30, AccuracyGoal -> 10]
(* 8.61998254998083821843412792489*10^-16 *)

Any ideas how to make this NIntegrate work without errors, while at the same time being confident in the Sign of the outcome?

Thanks! Laurens

share|improve this question

2 Answers 2

up vote 1 down vote accepted

No errors or warnings:

NIntegrate[(1/(108 π)) E^(-((2 Pq^2)/9) - 
     1/18 (53 - 16 Pq + 2 Pq^2)^2) (12 + E^(1/18 (53 - 16 Pq +
           2 Pq^2)^2) Sqrt[2 π] (9 - 4 Pq^2) - 
    E^(1/18 (53 - 16 Pq + 2 Pq^2)^2) Sqrt[
      2 π] (-9 + 4 Pq^2) Erf[(-53 + 16 Pq -  2 Pq^2)/(3 Sqrt[2])]), {Pq, -∞, ∞}, 
 MinRecursion -> 6, MaxRecursion -> 20, Method -> "LocalAdaptive", 
 WorkingPrecision -> 20, AccuracyGoal -> 10]

(* -3.2972035487866904746*10^-15 *)
share|improve this answer
    
That's great Rasher, thanks, but the answer you get is different from both the previous two answers I got. How can I be confident about this result? –  LBogaardt Apr 24 at 11:05

Changing the parameters in NIntegrate as the OP does can show there are difficulties with the numerics.

Checking the numerics

First let int1 be the OP's integrand and int3 be a simplified version to be discussed later.

int1 = (1/(108 π)) E^(-((2 Pq^2)/9) - 1/18 (53 - 16 Pq + 2 Pq^2)^2) (12 + 
     E^(1/18 (53 - 16 Pq + 2 Pq^2)^2) Sqrt[2 π] (9 - 4 Pq^2) - 
     E^(1/18 (53 - 16 Pq + 2 Pq^2)^2) Sqrt[2 π] *
       (-9 + 4 Pq^2) Erf[(-53 + 16 Pq - 2 Pq^2)/(3 Sqrt[2])]);
int3 = 1/108 E^(-((2 Pq^2)/9)) ((12 E^(-(1/18) (53 - 16 Pq + 2 Pq^2)^2))/π + 
     Sqrt[2/π] (9 - 4 Pq^2) Erfc[-((-53 + 16 Pq - 2 Pq^2)/(3 Sqrt[2]))]);

The OP's integrand evaluates accurately with 20 digits of precision only over a fairly small domain. For instance, it takes 40 digits of precision for Pq -> 1 to get one digits of precision in the output:

int1 /. Pq -> 1.`20
(* 0.*10^-20 *)

int1 /. Pq -> 1.`40
Precision[%]
(*
  6.*10^-39
  1.06324
*)

The problem is, as often is the case, subtraction of like-sized terms. The integrand can easily be simplified to

$$\frac{e^{-\frac{2 \text{Pq}^2}{9}-\frac{1}{18} \left(2 \text{Pq}^2-16 \text{Pq}+53\right)^2} \left(\sqrt{2 \pi } e^{\frac{1}{18} \left(2 \text{Pq}^2-16 \text{Pq}+53\right)^2} \left(9-4 \text{Pq}^2\right) \left(\text{erf}\left(\frac{-2 \text{Pq}^2+16 \text{Pq}-53}{3 \sqrt{2}}\right)+1\right)+12\right)}{108 \pi }\,.$$ Here the problem can be seen: The error function is nearly equal to $-1$ for all Pq. The best thing to do is to use Erfc. This leads to int3, which is set equal to $$\frac{1}{108} e^{-\frac{2 \text{Pq}^2}{9}} \left(\sqrt{\frac{2}{\pi }} \left(9-4 \text{Pq}^2\right) \text{erfc}\left(-\frac{-2 \text{Pq}^2+16 \text{Pq}-53}{3 \sqrt{2}}\right)+\frac{12 e^{-\frac{1}{18} \left(2 \text{Pq}^2-16 \text{Pq}+53\right)^2}}{\pi }\right)\,.$$

int1 - int3 // FullSimplify
(* 0 *)

Below are comparisons of the numerics for 20 and 100 digits of precision (code given for 20 digits only). One can see the great difficulty that the formula int1 presents for Mathematica. On the other hand, Mathematica handles int3 fairly easily.

With[{prec = 20},
 Labeled[
  TableForm[Table[{
     Round@pq,
     Precision[int1 /. Pq -> (N[pq, prec] /. 0 -> SetAccuracy[0, prec])],
     Precision[int3 /. Pq -> (N[pq, prec] /. 0 -> SetAccuracy[0, prec])]},
    {pq, 4 - 6, 4 + 6}], 
   TableHeadings -> {None, {"Pq", "int1", "int3", "N@int1", "N@int3"}}],
  Row[{"Precision ", prec}, BaseStyle -> {Italic}]
  ]
 ]

The tables show the digits of precision in the output value:

Mathematica graphics Mathematica graphics Mathematica graphics

Integrating

The obvious thing to try is

NIntegrate[int3, {Pq, -Infinity, Infinity}]
(* -3.29719*10^-15 *)

It evaluates quickly and without messages. But having gotten into a muddle with the original attempts, we might wonder whether this integral is the accurate one.

Verification

Let's start by observing [ref] that $$\text{erfc}(x)<\frac{e^{-x^2}}{\sqrt{\pi }}, \quad\hbox{for}\ x > 1\,.$$ We can use this to calculate (exactly) bounds for the tails of our integrals. If we restrict the interval of integration to $2\le {\rm Pq} \le 6$, we are concerned with bounding the integral outside this interval. Outside this interval, the argument to the error function is greater than 6 and

Abs[D[-((53 - 16 Pq + 2 Pq^2)^2/18) - (2 Pq^2)/9, Pq]] > 1

Thus the error introduced by integrating from 2 to 6 is bounded by the sum of the following two integrals:

intupper = 
  1/108 E^(-((2 Pq^2)/9)) ((12 E^(-(1/18) (53 - 16 Pq + 2 Pq^2)^2))/π + 
     Sqrt[2] E^(-(1/18) (53 - 16 Pq + 2 Pq^2)^2)/Pi);
du = D[(53 - 16 Pq + 2 Pq^2)^2/18 + (2 Pq^2)/9, Pq];

err1 = Integrate[intupper*du, {Pq, 6, Infinity}]
err2 = Integrate[intupper*du, {Pq, -Infinity, 2}]
N@Abs[err2]
N@Abs[err1]
(*
  (12 + Sqrt[2])/(108 E^(985/18) π)
 -((12 + Sqrt[2])/(108 E^(857/18) π))

  8.31283*10^-23
  6.78317*10^-26
*)

The integral from 2 to 6 is

NIntegrate[int3, {Pq, 2, 6}]
(* -3.29719*10^-15 *)

NIntegrate[int3, {Pq, 2, 6}, WorkingPrecision -> 100, 
 MaxRecursion -> 20, AccuracyGoal -> Infinity, PrecisionGoal -> 20]
(* -3.2971877967...*10^-15 *)

This differs from the whole integral by less than 10^-23:

NIntegrate[int3, {Pq, -Infinity, Infinity}, WorkingPrecision -> 100, 
 MaxRecursion -> 20, AccuracyGoal -> Infinity, PrecisionGoal -> 20]
(* -3.2971877942...*10^-15 *)

It's hard (impossible?) to be completely certain about a numerical approximation.

Another verification

For a partition of the interval from 2 to 6, we can approximate the upper and lower sums, which bound the integral. The integrand has two extrema, one a local maximum and one a local minimum:

NSolve[D[int3, Pq] == 0 && 2 < Pq < 6, Pq]
D[int3, {Pq, 2}] /. %

(* {{Pq -> 3.39448}, {Pq -> 4.05311}} *)
(* {-5.46427*10^-14, 1.02747*10^-13} *)

When a subinterval of the partition does not contain these numbers, the maximum and minimum of the integrand will be at the endpoints. When it does, we can extremum if appropriate. Below are the two functions that compute the lower and upper sums.

Clear[max, min, lower, upper];

max[prec_] :=
  ({max[prec], min[prec]} =
     Pq /. NSolve[D[int3, Pq] == 0 && 2 < Pq < 6, Pq, WorkingPrecision -> prec];
     max[prec]);
min[prec_] := (max[prec]; min[prec]);

Options[lower] = Options[upper] = {WorkingPrecision -> 30};
lower[n_Integer, OptionsPattern[]] := 
  Module[{vals = 
     int3 /. Pq -> Range[SetPrecision[2, OptionValue[WorkingPrecision]], 6, 4/n], 
    pq1 = 2 + Floor[min[OptionValue[WorkingPrecision]] - 2, 4/n]},
     (4/n) (Total[
       Min /@ Transpose[{Most[vals], Rest[vals]}]] + (Min[#] - Min@Most[#] &[int1 /.
           Pq -> {pq1, pq1 + 4/n, min[OptionValue[WorkingPrecision]]}]))
   ];
upper[n_Integer, OptionsPattern[]] := 
  Module[{vals = 
     int3 /. Pq -> Range[SetPrecision[2, OptionValue[WorkingPrecision]], 6, 4/n], 
    pq1 = 2 + Floor[max[OptionValue[WorkingPrecision]] - 2, 4/n]},
   (4/n) (Total[
       Max /@ Transpose[{Most[vals], Rest[vals]}]] +
         (Max[#] - Max@Most[#] &[int1 /.
           Pq -> {pq1, pq1 + 4/n, max[OptionValue[WorkingPrecision]]}]))
   ];

This is sufficient to prove the integral is negative, since the magnitude is greater than the error from truncating the interval of integration:

upper[16]
(* -1.20840289638770*10^-15 *)

Here are bounds on the truncated integral:

lower[1000000, WorkingPrecision -> MachinePrecision]
upper[1000000, WorkingPrecision -> MachinePrecision]

(* -3.29722*10^-15 *)
(* -3.29715*10^-15 *)

Conclusion

The integral is accurate:

NIntegrate[int3, {Pq, -Infinity, Infinity}]
(* -3.29719*10^-15 *)

It is well-behaved under increasing the precision:

NIntegrate[int3, {Pq, -Infinity, Infinity}, WorkingPrecision -> #, 
   AccuracyGoal -> Infinity, PrecisionGoal -> Min[20, # - 1]] & /@ {10, 20, 30, 40}
(*
  {-3.297187795*10^-15,
   -3.2971877942145651745*10^-15, 
   -3.29718779421456517136335352150*10^-15, 
   -3.297187794214565171363353521467832469704*10^-15}
*)

One of the take-aways here is that sometimes it is helpful to apply mathematics to Mathematica.

share|improve this answer
    
+1 for taking the time, very nice. –  rasher Apr 27 at 20:54
    
@rasher Thanks. It's kind of a slow day and doing mathematics is fun, even if it's fairly elementary. –  Michael E2 Apr 27 at 21:02
    
Cool! Thanks a lot of all the effort :) –  LBogaardt Apr 29 at 22:38
    
@LBogaardt You're welcome. –  Michael E2 Apr 30 at 0:56

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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