# Bug in the "PrincipalValue" method of NIntegrate?

## Update 2

The technical support suggests setting AccuracyGoal. As I quote here

It does appear that the program is returning an incorrect result along with a warning message when evaluating the following command:

NIntegrate[f[x], {x, 0, 1, 3}, Method -> "PrincipalValue"]

As stated in the warning message, we should specify a finite value for the AccuracyGoal option. As such the following command will work smoothly and gives correct output:

NIntegrate[f[x], {x, 0, 1, 3}, Method -> "PrincipalValue", AccuracyGoal -> 10]

In the meantime, I've forwarded your example to our developers who will take a closer look and consider how to fix any issues they find.

By setting AccuracyGoal, all warnings are gone (even the warnings about integral and error estimates are 0 for 1/(x-1) over the region $[1-\epsilon, 1+\epsilon]$ are gone, which is quite strange). One may even set AccuracyGoal->0 and still get the correct result.

Despite of the fix, the mechanism of PrincipalValue still remains a puzzle to me.

## Update 1

As pointed by Emilio Pisanty, one should pay attention to EvaluationMonitor to monitor principal value integration. I update the codes, and find some more weird output.

Block[{lst1, i, f, ans1, ans2},
lst1 = {};
i = 0;
f[x_?NumericQ] := (AppendTo[lst1, {x, ++i}]; 1/(x - 1));
ans1 = NIntegrate[f[x], {x, 0, 1, 3}, Method -> "PrincipalValue"];
ans2 = NIntegrate[1/(x - 1), {x, 0, 1, 3}, Method -> "PrincipalValue"];
Show[{Plot[1/(x - 1), {x, 0, 3}, PlotRange -> {-10, 600}],
ListPlot[lst1, PlotRange -> All, PlotStyle -> Red]},
PlotLabel -> {ans1, ans2}]
]

>NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {9.46688*10^-16}. NIntegrate obtained -0.928109 and 1.0151159883865095 for the integral and error estimates. >>

>NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>


The answer is $\ln2 \approx 0.693$. However, one can clearly see that by monitoring the evaluation in a more rigorously way, one will get a totally wrong answer, even though that the sampling points seems correctly distribute around the singularity.

If one doesn't monitor the evaluation, only defining the function as f[x_?NumericQ] := 1/(x - 1), one will still get wrong answer.

It seems that shutting down the symbolic processing will change the result. This is really disturbing. Since this technique is used by the official. And one has to shut down the symbolic processing to evaluate the multi-dimensional principal integration.

Any ideas?

## OP

This question is related to this one, I encountered some weird outputs when trying to solve it.

Here is the code:

Block[{i, ansWrong, ansCorrect},
i = 1;
ansWrong =
Reap[NIntegrate[1/(10000 x - 2025), {x, 0, 81/400, 1/2},
Method -> "PrincipalValue", EvaluationMonitor :> Sow[{x, i++}]]];
i = 1;
ansCorrect =
Reap[NIntegrate[1/10000 1/(x - 81/400), {x, 0, 81/400, 1/2},
Method -> "PrincipalValue", EvaluationMonitor :> Sow[{x, i++}]]];
GraphicsGrid[{{
ListPlot[ansWrong[[2, 1]], PlotLabel -> First@ansWrong],
ListPlot[ansCorrect[[2, 1]], PlotLabel -> First@ansCorrect]}}
]
]


NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {2.04383*10^-16}. NIntegrate obtained 0.000028159986301613934and 0.00013893352072413546 for the integral and error estimates. >>

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>

The two integrands are mathematically identical, however, from the left plot, one can clearly see that the region is wrongly sampled and it is only one of the regions to be calculated. (The empty part of the sampling region is $[x_0-\epsilon, x_0+\epsilon]$ which includes the singularity.) What's worse, the program is stuck at origin, giving totally wrong answer.

So the question is how to understand these behaviors? (v10.4.0, Windows10)

PS: the above question can be bypassed by making the integrand unsymmetrical with respect to the singularity.

NIntegrate[(1 - Exp[-x])/(10000 x - 2025), {x, 0, 81/400, 1/2},
Method -> "PrincipalValue"] +
NIntegrate[Exp[-x]/(10000 x - 2025), {x, 0, 81/400, 1/2},
Method -> "PrincipalValue"]

(*==>0.0000384674*)


And the solution to that post is

ClearAll[f]
f[y_?NumericQ] :=
NIntegrate[(1 - Exp[-x])/(
100^2 x - 90^2 y (1 - y)), {x, 0, 81/100 y (1 - y), 1 - y},
Method -> "PrincipalValue"] +
NIntegrate[Exp[-x]/(
100^2 x - 90^2 y (1 - y)), {x, 0, 81/100 y (1 - y), 1 - y},
Method -> "PrincipalValue"]
NIntegrate[f[y], {y, 0, 1}]

(*==>0.0000600275*)


which agrees with the analytical solution Log[(100*10^(38/81))/(81*19^(19/81))]/10000.

• Note that the samplings are not as different as they look in your plots, which can be brought out by setting PlotRange→Full in your plots. The sampling near the singularity appears to match pretty much exactly, regardless of the value of SingularPointIntegrationRadius. The difference is that the first integral does a lot of sampling around $x=0$ (not that that is any less mysterious, of course), and this persists even if the lower limit of integration is negative. Jun 30, 2016 at 17:28
• ... and, also, it's important to note that EvaluationMonitor is not that reliable when Method→"PrincipalValue" is used, and it differs from the explicit evaluation tracking presented in the CPV section of the Integration Strategies page. Jun 30, 2016 at 17:38
• @EmilioPisanty I'm so misguided by the ListPlot. Unlike Plot, I think for the case it didn't show all the data points, it should have given us some warnings. You are right about the EvaluationMonitor for Principalvalue`, although the werid behavior persists. Jul 1, 2016 at 2:07
• @luyuwuli what did you mean by "the above question can be bypassed by making the integrand unsymmetrical with respect to the singularity"? The original integrand is already unsymmetric w.r.t. $x_0=81/400=0.2025$ in the domain $[0, 1/2]$.
– lxy
Jun 27, 2019 at 6:43
• @jsxs What I mean is that the integrand is symmetric around the singularity, i.e. on $(x_0-\delta,x_0+\delta)$ where $x_0$ is the singular point and $\delta$ doesn't have to be small. By making the integrand slightly distorted, the problem can be bypassed. Jun 27, 2019 at 15:48