4
$\begingroup$

I have the following test function which I want to evaluate at various points from $x=0$ to $5$:

$f(x) = \int^{20}_{1}\frac{1}{x-y} dy$

There's clearly a singularity here in the integral at $x=y$, which is why in the code sample below I specified the principal value method & exclusion.

TestList = Range[0, 5, 0.5]
TestFunc[x_] :=  NIntegrate[1/(x - y), {y, 1, 20}, Method -> "PrincipalValue", Exclusions -> {x == y}]

Mathematica is able to evaluate TestFunc at various points, e.g. if I input

TestFunc[2]

I get $-2.89037$, which is correct ($\int^{20}_{1}\frac{1}{2-y} dy=-\mathrm{ln}(18)= -2.89037$). However if I input the list

TestFunc[TestList]

Mathematica refuses to evaluate, returning the error 

NIntegrate::pvexeq: PrincipalValue cannot work with the specified exclusions.

Why? Strictly, when I evaluate TestFunc[2], I get this warning message 

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."

However the result is still accurate, and increasing MinRecursion to 80 as well as MaxRecursion to 90 doesn't change anything.

Improve this question
$\endgroup$
1
  • 4
    $\begingroup$ You are feeding the entire list to the Exclusions option in a single Integrate command. Either use Map (/@): TestFunc /@ TestList or define your function as Listable. $\endgroup$
    – SPPearce
    Commented Dec 13, 2018 at 5:42

2 Answers 2

Reset to default
5
$\begingroup$

Interpreting the OP's example as an MWE and not the actual problem, there are three issues: (1) evaluating the integral for a list of test values; (2) the integral and error being zero; and (3) the example integral being divergent even in the principal value sense if the singularity is at an endpoint of the interval of integration.

Listability

The easiest way to make testFunc operate on a list of input values is to give it the attribute Listable.

testList = Range[0, 5, 0.5];
ClearAll[testFunc];
SetAttributes[testFunc, Listable];
testFunc[x_] := 
 NIntegrate[1/(x - y), {y, 1, 20}, Method -> "PrincipalValue", 
  Exclusions -> {x == y}, AccuracyGoal -> 100];

The following gives some errors because x == 1 is in testList and the principal value does not exist for it (3rd entry).

testFunc[testList]
(*
{-2.99573, -3.66356, -39.0857, -3.61092, -2.89037, -2.45674, 
 -2.14007, -1.88707, -1.67398, -1.48808, -1.32176}
*)

Zero integrals

Integrals that are zero cause an error in the error norm NIntegrate uses by default and the NIntegrate::izero warning message to be issued.

NIntegrate[0., {x, 0, 1}]

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

(*  0.  *)

For an integral that is zero, the solution is to set a finite value for AccuracyGoal as in testFunc[] above. The default setting is AccuracyGoal -> Infinity, and in that case NIntegrate uses only relative error in determining sufficient convergence. However, the relative error when the estimate and absolute error estimate are both zero is undefined. Setting a finite value for AccuracyGoal indicates how close to zero the estimate has to be for a result of zero to be considered acceptable. A result of zero for both the integral and error might occur because underflow occurred at all sampling points, or the sampling points accidentally all landed on roots of the integrand, and perhaps other rare possibilities.

Here is an example of a positive integral and integrand that evaluates to zero:

NIntegrate[Exp[-1000 (Exp[1000 (x - Sqrt[2.])^2] - 1)], {x, 0, 2}]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option....

(*  0.  *)

NIntegrate[Exp[-1000 (Exp[1000 (x - Sqrt[2])^2] - 1)], {x, 0, Sqrt[2], 2}]
(*  0.000885895  *)

Instead of using MinRecursion as suggested, I manually subdivided the interval at a point x == Sqrt[2] where the function does not underflow to zero. (Setting a finite value for AccuracyGoal in the first integral, even as high as AccuracyGoal -> 100 makes the warning go away but still returns the inaccurate value of 0.) Further, the second, positive result is still inaccurate. Increasing MinRecursion helps. Increasing WorkingPrecision can also help with underflow, but in this case, one also has to raise MaxRecursion:

NIntegrate[[-1000 (Exp[1000 (x - Sqrt[2])^2] - 1)], {x, 0, Sqrt[2], 2}, 
 MinRecursion -> 2]
(*  0.00177179  *)

NIntegrate[Exp[-1000 (Exp[1000 (x - Sqrt[2])^2] - 1)], {x, 0, 2}, 
 MaxRecursion -> 20, WorkingPrecision -> 100]
(*  0.00177179007867237284978305329244440725705339189727922053943991469...  *)

Why does the OP's integral evaluate to zero? It does not even come close to underflowing!

This has to do with how the "PrincipalValue" method is implemented. It takes a symmetric subinterval around each singularity, divides it at the singularity, reflects the integrand around the singularity and adds. If the integral converges in the principal value sense, the singularities will cancel and there will be a finite residual integral. In the OP's case, it is zero. For $x=2$ in the OP's case, the transformation is as follows: $$ \int_{3/2}^{5/2} {dy \over 2-y} = \int_{3/2}^{2} {dy \over 2-y} + \int_{2}^{5/2} {dy \over 2-y} = \int_0^{1/2} {du \over u} + \int_0^{1/2} -{du \over u} = \int_0^{1/2} 0\; dy \,, $$ where the first substitution is $2 - y = u$ and the second substitution is $2 - y = -u$.

We can examine the setup NIntegrate uses with IntegrationMonitor:

Quiet@Reap[       (* ignore error messages *)
   Block[{x = 2},
    NIntegrate[1/(x - y), {y, 1, 20}, Method -> "PrincipalValue", 
     Exclusions -> {x == y}, MaxRecursion -> 0,
     IntegrationMonitor -> (Sow[
         Map[{#1["Integrand"]["FunctionExpression"], #1@
             "Boundaries", #1@"Integral", #1@"Error"} &, #1]] &)]
    ]
   ][[-1, -1]]

(* integrand     interval      integral est.  error est.
  {{{1/(2 - y),  {{1, 1.5}},    0.693147,     2.27069*10^-8}},
   {{1/(2 - y),  {{2.5, 20}},  -3.58393,      0.132922}},
   {{0,          {{0, 0.5}},    0.,           0.}}}
*)

The third entry is the principal value sub-integral.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for answer, very helpful. Question: what does MWE stand for? $\endgroup$
    – Allure
    Commented Dec 14, 2018 at 2:52
  • $\begingroup$ @Allure MWE = minimal working example. Users are encouraged to post the simplest code that illustrates the problem instead of a large complicated mess of code, most of which has nothing to do with the issue. $\endgroup$
    – Michael E2
    Commented Dec 14, 2018 at 3:45
3
$\begingroup$

This integral can be calculated exactly

TestList = Range[0, 5, 1/2];
    TestFunc[x_List] := Integrate[1/(x - y), {y, 1, 20}, PrincipalValue -> True]
    TestFunc[TestList]
    {-Log[20], -Log[39], -\[Infinity], -Log[37], -Log[18], -Log[35/
       3], -Log[17/2], -Log[33/5], -Log[16/3], -Log[31/7], -Log[15/4]}
Improve this answer
$\endgroup$

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