9
$\begingroup$

The toy model is: $$\int_{-2}^{2}\int_{-2}^{2}\frac{1}{1-(x^2+y^2)}\, dx\, dy$$

The integrand have opposite sign across the circle $x^2+y^2=1$, so one would expect that the integral has meaning only under principal value.

( I'm not sure whether I use the terminology "multi-dimensional principal value" mathematically correct, but I think you understand what I mean.)

NIntegrate[1/(1 - x^2), {x, 0, 1, 2}, Method -> "PrincipalValue"]
(*0.549306*)

However, Mathematica can only handle one dimensional problem. The following code doesn't work:

NIntegrate[1/(1 - (x^2 + y^2)), {x, -2, 2}, {y, -2, 2}, 
 Exclusions -> x^2 + y^2 == 1, Method -> "PrincipalValue"]

Mathematica gives:

NIntegrate::pvdim: PrincipalValue can be used for one-dimensional integrals only. 
NIntegrate::pvdim: PrincipalValue can be used for one-dimensional integrals only. 
NIntegrate::pvdim: PrincipalValue can be used for one-dimensional integrals only.  
General::stop: Further output of NIntegrate::pvdim will be suppressed during this calculation. 
(*0.*)

(It would be appreciated if someone tell me why Mathematica returns more than three NIntegrate::pvdim warnings)

In this toy model, via polar coordinates, one can work out the integration (at least this shows the fact that the integration converge): $$\int_{-2}^{2}\int_{-2}^{2}\frac{1}{1-(x^2+y^2)}\, dx\, dy = 2\pi\int_{0}^{2}\frac{r}{1-r^2}\, dr+4\int_{0}^{2}\int_{\sqrt{4-x^2}}^{2}\frac{1}{1-(x^2+y^2)}\, dy\, dx$$

and this time Mathematica can work out the one dimensional principal value integration quickly:

2 Pi NIntegrate[r/(1 - r^2), {r, 0, 1, 2}, Method -> "PrincipalValue"]
4 NIntegrate[1/(1 - (x^2 + y^2)), {x, 0, 2}, {y, Sqrt[4 - x^2], 2}]
%%+%
(*-3.45139,  -0.872413,   -4.32381*)

Seems good. However, the problem is the toy model is so contrived; to deal with realistic problem, the method discribed above is barely useful.

Is there any general method to deal with this kind of problem in Mathematica? I'll be grateful if someone tell me some references related to the topic.

The second question is related to the first:

Integrate[1/(1 - (x^2 + y^2)), {x, -2, 2}, {y, -2, 2}]
N@%

the first line returns:

$$\int_{-2}^2 -\frac{2 \tan ^{-1}\left(\frac{2}{\sqrt{x^2-1}}\right)}{\sqrt{x^2-1}} \, dx$$

during the second evaluation, Mathematica gives some warnings:

NIntegrate::slwcon: Numerical integration converging too slowly...
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections...

and returns:

(*-4.30645 - 9.85225 I*)

Why N@Integrate[...] actually invoke NIntegrate and is it a bug for Mathematica to return a complex number?

Edit

To be concrete, the expression I want to integrate is:

$$\frac{\exp({-\text{p1px}^2-\text{p1py}^2-\text{p1pz}^2-\text{p1x}^2-\text{p1y}^2-\text{p1z}^2})}{\text{p1px}^2+\text{p1py}^2+\text{p1pz}^2+\left(\text{p1x}+\frac{\sqrt{3}}{2}\right)^2+\text{p1y}^2+\left(\text{p1z}+\frac{1}{2}\right)^2-2}$$

The integrand should be integrate over all space (i.e. from $-\infty$to$\infty$)for all the variables. However, because of the Gaussian surpression, one can integrate from -5 to 5.

Naive input such as:

NIntegrate[E^-(p1px^2 + p1py^2 + p1pz^2 + p1x^2 + p1y^2 + p1z^2) 
 1/(-2 + p1px^2 + p1py^2 + p1pz^2 + (Sqrt[3]/2 + p1x)^2 + p1y^2 + (1/2 + p1z)^2),
  {p1x, -5, 5}, {p1y, -5, 5}, {p1z, -5, 5},
  {p1px, -5, 5}, {p1py, -5, 5}, {p1pz, -5, 5}]

will return NIntegrate::slwcon: and NIntegrate::eincr warnings, and the error is much bigger than the result. ( 9.99021 and 156.922 for the integral and error estimates. )

How to improve the output?

Improve this question
$\endgroup$
8
  • $\begingroup$ Who knows, maybe even your actual problem is amenable to a clever substitution. Why not post it? $\endgroup$
    – J. M.'s eventual burnout ♦
    Aug 31, 2015 at 17:33
  • $\begingroup$ @J. M. The actual problem involve lengthy polynomials in the denominator. It's too complicated to post. What I'm looking for is some bulit in method that can handle the multi-dimensional principal value problem or someone can told me which mathematical field does this question belong to. $\endgroup$
    – luyuwuli
    Sep 1, 2015 at 2:57
  • $\begingroup$ @J. M. I append one more realistic example to the post. I hope this would help. Please have a look:) $\endgroup$
    – luyuwuli
    Sep 1, 2015 at 12:12
  • $\begingroup$ Your new example looks like it can be dealt with by a hyperspherical coordinate substitution, suitably shifted. $\endgroup$
    – J. M.'s eventual burnout ♦
    Sep 1, 2015 at 13:28
  • $\begingroup$ You might wish to have a look at my EDIT #1. $\endgroup$
    – Dr. Wolfgang Hintze
    Sep 2, 2015 at 11:41

2 Answers 2

Reset to default
6
$\begingroup$

Using polar coordinates r and f, the region of integration is given by 

{ 0<= r <=2/Cos[f], 0<= f <= 2 \[Pi] }

We can then proceed as follows.

First integration with PrincipalValue

g = 8 Integrate[r/(1 - r^2), {r, 0, 2/Cos[f]}, Assumptions -> 0 < f < \[Pi]/4, PrincipalValue -> True]

(* Out[451]= -4 Log[-1 + 4 Sec[f]^2] *)

Second integration

h = Integrate[g, {f, 0, \[Pi]/4}];

FullSimplify[h]

(* Out[453]= 4 Catalan + 1/2 \[Pi] Log[97 - 56 Sqrt[3]] - 
 2 I (PolyLog[2, 7 I - 4 I Sqrt[3]] - PolyLog[2, I (-7 + 4 Sqrt[3])]) *)

% // N

(* Out[454]= -4.32381 + 0. I *)

This is a real quantity, and it is in agreement with the real part of one of the values mentioned in the OP.

EDIT #1

Your edit provides a mathematically more natural problem since it avoids the hard ("man-made") cutoff.

Your function to be integrated of the whole space is

f = Exp[-(p1px^2 + p1py^2 + p1pz^2 + p1x^2 + p1y^2 + p1z^2) ] 1/(-2 + p1px^2 + p1py^2 + p1pz^2 + (Sqrt[3]/2 + p1x)^2 + p1y^2 + (1/2 + p1z)^2);

Simplifying the name of your variables a bit gives

f1 = f /. {p1px -> x, p1py -> y, p1pz -> z, p1x -> u, p1y -> v, p1z -> w}

(*
Out[6]= E^(-u^2 - v^2 - w^2 - x^2 - y^2 - z^2)/(-2 + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2)
*)

Now let us produce the denominator using

Integrate[Exp[-a h], {a, 0, \[Infinity]}, Assumptions -> h > 0]

(*
Out[97]= 1/h
*)

i.e.

Integrate[Exp[-a h], {a, 0, \[Infinity]}, Assumptions -> h > 0] /. 
 h -> p + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2

(*
Out[99]= 1/(p + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2)
*)

Here we have replaced the number (-2) temporarily by the parameter p.

Now, under the a-intergal we have all terms in the exponent:

f2 = Exp[-u^2 - v^2 - w^2 - x^2 - y^2 - z^2 - 
    a (p + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2)];

and the integral over the whole space is easily done

f3 = Integrate[
  f2, {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, {z, -\
\[Infinity], \[Infinity]}, {u, -\[Infinity], \[Infinity]}, {v, -\[Infinity], \
\[Infinity]}, {w, -\[Infinity], \[Infinity]}, Assumptions -> a > 0]

(*
Out[104]= (E^(-a (1/(1 + a) + p)) \[Pi]^3)/(1 + a)^3
*)

The expression we are looking for is therefore given by

g[p_] = Integrate[f3, {a, 0, \[Infinity]}, Assumptions -> p > 0]

(*
Out[107]= Integrate[(E^(-a (1/(1 + a) + p)) \[Pi]^3)/(1 + a)^3, {a, 0, \[Infinity]}, Assumptions -> p > 0]
*)

In order for the integral to converge we have tried to get an expression with p>0 hoping to then be able to continue it analytically.

But the integral is not done analytically by Mathematica. Hence we try a numeric approach 

gn[p_] := NIntegrate[(
  E^(-a (1/(1 + a) + p)) \[Pi]^3)/(1 + a)^3, {a, 0, \[Infinity]}]

For example

gn[1]

(* Out[112]= 7.54936 *)

Now we would need information about a kind of mapping from this divergent integral if p<0 and the procedure of taking the PrincipalValue.

I'll leave the results of my preliminary studies of this question for a later edit.

Sumarizing

We have reduced the multidimensional integral to a one dimensional integral. It remains to be studied how taking the principal value in the original integral maps to a procedure to give a meaning the divergent one-dimensional integral.

Remark: The reduction to a one dimensional integral described here can of course also be done for the "cut-off" version of the problem.

EDIT #2

The analytic continuation to p < 0 can be done as follows.

After the transformation a -> 1/z-1 the integral g can be written as

g1[p_] = \[Pi]^3 E^(p - 1) Integrate[z E^z Exp[-p/z], {z, 0, 1}];

Expanding E^z this becomes

g2[p_] = \[Pi]^3 E^(p - 1)
   Sum[1/k! Integrate[z^(k + 1) Exp[-p/z], {z, 0, 1}], {k, 0, \[Infinity]}];

But

Integrate[z^(k + 1) Exp[-p/z], {z, 0, 1}, Assumptions -> p > 0] 

(*
Out[212]= ExpIntegralE[3 + k, p]
*)

so that

g3[p_] = \[Pi]^3 E^(p - 1)
   Sum[1/k! ExpIntegralE[3 + k, p], {k, 0, \[Infinity]}];

But now, this expression can be continued analytically to p < 0.

And we obtain numerically to a very good approximation

Table[{k, \[Pi]^3 (E^(-1 + p) ExpIntegralE[3 + k, p])/k! /. p -> -2.}, {k, 0, 
  15}]

(*
Out[213]= {{0, 1.81404 - 9.69943 I}, {1, 5.01155 - 6.46628 I}, {2, 
  2.67871 - 1.61657 I}, {3, 0.73738 - 0.215543 I}, {4, 
  0.140661 - 0.0179619 I}, {5, 0.021617 - 0.00102639 I}, {6, 
  0.00288102 - 0.0000427664 I}, {7, 0.000342928 - 1.35766*10^-6 I}, {8, 
  0.0000368633 - 3.39416*10^-8 I}, {9, 3.6023*10^-6 - 6.85689*10^-10 I}, {10, 
  3.21984*10^-7 - 1.14282*10^-11 I}, {11, 
  2.64847*10^-8 - 1.59834*10^-13 I}, {12, 
  2.01624*10^-9 - 1.90279*10^-15 I}, {13, 
  1.42798*10^-10 - 1.95158*10^-17 I}, {14, 
  9.4526*10^-12 - 1.74248*10^-19 I}, {15, 5.87243*10^-13 - 1.36665*10^-21 I}}
*)

Summing up

Plus @@ %

(*
Out[214]= {120, 10.4072 - 18.0169 I}
*)

Now we make the bold assumptions that taking the principal value is equivalent to taking the real part here.

Then the original integral becomes numerically

gg = 10.407221182257118 ...

Of course we need to justify the bold assumption. I'll try to do this later

EDIT #3

Gaussian decay instead of cut off

Let us first go back to the original toy example, but now with Gaussian decay instead of hard cut off.

The integral is now

f2g = Integrate[Exp[-x^2 - y^2]/(
  1 - x^2 - y^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \
\[Infinity]}]

During evaluation of In[286]:= Integrate::idiv: Integral of E^(-x^2-y^2)/(-1+x^2+y^2) does not converge on {-[Infinity],[Infinity]}. >>

Using PrincipalValue gives a finite result

f2gp = Integrate[Exp[-x^2 - y^2]/(
  1 - x^2 - y^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \
\[Infinity]}, PrincipalValue -> True]

(*
Out[287]= (\[Pi] ExpIntegralEi[1])/E
*)

% // N

(* Out[288]= 2.19024 *)

We can also transform to polar coordinates. The angular part is now trivial and gives a factor 2 [Pi]. The integral is therefore

2 \[Pi] Integrate[(Exp[-r^2] r)/(1 - r^2), {r, 0, \[Infinity]}, 
  PrincipalValue -> True]

(*
Out[289]= (\[Pi] ExpIntegralEi[1])/E
*)

in agreement with the integral in Cartesian coordinates.

Notice that, in contrast to the cut off version, the values is positive.

Now let's do the same integral with replacement of the denominator (as in EDIT #1)

Integrate[Exp[-x^2 - y^2 - 
   a (p - x^2 - 
      y^2)], {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}]

(*
Out[322]= ConditionalExpression[-((E^(-a p) \[Pi])/(-1 + a)), Re[a] < 1]
*)

Integrate[%[[1]], {a, 0, \[Infinity]}, Assumptions -> p > 0, 
 PrincipalValue -> True]

(*
Out[323]= E^-p \[Pi] ExpIntegralEi[p]
*)

% /. p -> 1.

(* Out[324]= 2.19024 *)

Higher dimensions with standard integrand

3 dimensions:

f3 = Integrate[Exp[-x^2 - y^2 - z^2]/(
  p - x^2 - y^2 - 
   z^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, {z, \
-\[Infinity], \[Infinity]}, PrincipalValue -> True, Assumptions -> p > 0]

(*
Out[310]= 2 \[Pi]^(3/2) (-1 + 2 Sqrt[p] DawsonF[Sqrt[p]])
*)

4 dimensions

f4 = Integrate[Exp[-x^2 - y^2 - z^2 - t^2]/(
  p - x^2 - y^2 - z^2 - 
   t^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, {z, \
-\[Infinity], \[Infinity]}, {t, -\[Infinity], \[Infinity]}, 
  PrincipalValue -> True, Assumptions -> p > 0]

(*
Out[336]= \[Pi]^2 (-1 + E^-p p ExpIntegralEi[p])
*)

Arbitrary dimension

ff[n_] := Integrate[Exp[-Sum[x[i]^2, {i, 1, n}]]/(p - Sum[x[i]^2, {i, 1, n}]),
   Sequence @@ Table[{x[i], 0, \[Infinity]}, {i, 1, n}], 
  PrincipalValue -> True, Assumptions -> p > 0]

t = Table[ff[n], {n, 1, 6}]

(*
Out[354]= {
(Sqrt[\[Pi]] DawsonF[Sqrt[p]])/Sqrt[p], 
1/4 E^-p \[Pi] ExpIntegralEi[p], 
 1/4 \[Pi]^(3/2) (-1 + 2 Sqrt[p] DawsonF[Sqrt[p]]), 
 1/16 \[Pi]^2 (-1 + E^-p p ExpIntegralEi[p]), 
 1/48 \[Pi]^(5/2) (-1 - 2 p + 4 p^(3/2) DawsonF[Sqrt[p]]), 
 1/128 \[Pi]^3 (-1 - p + E^-p p^2 ExpIntegralEi[p])}
*)

Plot[{t[[1]], t[[3]], t[[5]]}, {p, -2, 3}]
(* 150902_plot_135.jpg *)

enter image description here

Plot[{t[[2]], t[[4]], t[[6]]}, {p, -2, 3}]
(* 150902_plot_246.jpg *)

enter image description here

EDIT #4

We can get a confirmation of the my bold hypothesis by considering another type of regularization of the Gaussian toy example.

$$\text{toyg}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\frac{e^{-x^2-y^2}}{-x^2-y^2+1}dydx$$

This method considers a general power k of the denominator and applies analytic continution of the result to k->-1

Integrate[Exp[-x^2 - y^2] (1 - x^2 - y^2)^
   k, {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}]

(*
Out[136]= ConditionalExpression[
 E^(-1 - I k \[Pi]) \[Pi] ((-1 + E^(2 I k \[Pi])) Gamma[1 + k] + 
    E Subfactorial[k]), Re[k] > -1]
*)

f = First[%]

(*
Out[137]= E^(-1 - I k \[Pi]) \[Pi] ((-1 + E^(2 I k \[Pi])) Gamma[1 + k] + 
   E Subfactorial[k])
*)

Hence toyg should be 

tyogk = Limit[f, k -> -1]

(*
Out[138]= (\[Pi] (-2 I \[Pi] - Gamma[0, -1]))/E
*)

% // N

(*
Out[139]= 2.19024 - 3.63082 I
*)

We have obtained a complex quantity.

Comparing this with the PV form

toygpv = Integrate[Exp[-x^2 - y^2]/(1 - x^2 - 
    y^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, 
 PrincipalValue -> True]

(*
Out[134]= (\[Pi] ExpIntegralEi[1])/E
*)

% // N

(* Out[135]= 2.19024 *)

shows that the real parts conincide.

Improve this answer
$\endgroup$
12
  • $\begingroup$ Thanks. I didn't realize that I can write the upper limite as a function of other parameters. (I also fixed a stupid mistake in the original post.) $\endgroup$
    – luyuwuli
    Sep 1, 2015 at 2:48
  • $\begingroup$ Do you think that at present all we can do is analysing the symmetry of the integrand and simplify it by hand? Are there any general method to do this? I mean automatically find and handle the "sigular curve" etc? At the first time, I take it for granted that the problem can be solved by some singularity handling algorithms (IMT,DoubleExponential...) and Mathematica should have some built in options to do this. Maybe I've asked too much. $\endgroup$
    – luyuwuli
    Sep 1, 2015 at 3:09
  • $\begingroup$ @ luyuwuli As long as I am concerned, you asked too much. It is probably possible to treat the 3D case. I started it, did the r integral but then I got stuck. $\endgroup$
    – Dr. Wolfgang Hintze
    Sep 1, 2015 at 6:59
  • $\begingroup$ For the 3D case, if you just want the numerical result, I've tried ImplicitRegion, it works (and pretty fast). $\endgroup$
    – luyuwuli
    Sep 1, 2015 at 7:45
  • 1
    $\begingroup$ @luyuwuli, I don't believe those PV-specialized versions have ever been implemented anywhere; I might consider taking a crack at it over the weekend. At least according to my own tests, the current Mathematica implementations of DE and IMT are not equipped for PV use; the convergence rate is markedly poor. $\endgroup$
    – J. M.'s eventual burnout ♦
    Sep 1, 2015 at 13:46
1
$\begingroup$

The following is a numerical solution to the toy model. It's too long to be a comment and not related to the question I really want to solve. Maybe someone can get help from the solution.

$$\int_{-2}^{2}\int_{-2}^{2}\frac{1}{1-(x^2+y^2)}\, dx\, dy = 2\pi\int_{0}^{2}\frac{r}{1-r^2}\, dr+\int_{\mathrm{rest\ region}}\frac{1}{1-(x^2+y^2)}\, dV$$

In Mathematica, by using ImplicitRegion, the formula above can be easily implemented as:

a1=2 Pi NIntegrate[r/(1 - r^2), {r, 0, 1, 2}, Method -> "PrincipalValue"]
a2=NIntegrate[1/( 1 - x^2 - y^2), {x, y} \[Element]
  ImplicitRegion[ x^2 + y^2 >= 4 \[And] -2 <= x <= 2 \[And] -2 <= y <= 2, {x, y}]]
a1+a2
(*-3.45139,-0.872413,-4.32381*)

The advantage of the ImplicitRegion is that the above code can easily extended to 3D or even higher dimension. For 3D case,

$rest=ImplicitRegion[
      x^2 + y^2 + z^2 >= 4 \[And] -2 <= x <= 2 \[And] -2 <= y <= 2 \[And] -2 <= z <= 2, {x, y, z}];  
    a1=4 Pi NIntegrate[r^2/(1 - r^2), {r, 0, 1, 2}, Method -> "PrincipalValue"]
    a2=NIntegrate[1/(1 - x^2 - y^2 - z^2), {x, y, z} \[Element] $rest]
a1+a2
(*-18.23,-6.88868,-25.1186*)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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