# How to calculate the principal value for a two-dimensional integral numerically? [closed]

I want to calculate the principal value of following two-dimensional integral numerically $$\int_{0}^{\infty}dx\int_{0}^{\infty}dy\sqrt{xe^{-100x}}\sqrt{ye^{-100y}}\frac{1-e^{1000\imath(x+y)}}{(x+y)(y-0.001)}$$ The mathematica code is

NIntegrate[Sqrt[x E^{-100x}]Sqrt[y E^{-100x}]((1-E^{1000I(x+y)})/(x+y))(1/(y-0.001)),{x,0,∞},{y,0,∞}]

The mathematica gives the following error

Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small

NIntegrate failed to converge to prescribed accuracy after 18 recursive bisections in y near {y,x} ={0.000999996,0.00121927}. NIntegrate obtained 0.0509438 +0.00910859I and 0.05892024847005836` for the integral and error estimates.

How can I calculate the principal value of integral without error with NIntegrate? Actually I am encountering with this types of integrals and errors in my line of work frequently. I need a solid method to deal with such integrals.

P.S.: I examined this integral in Matlab and it gave an answer without error. What is the difference between Mathematica and Matlab in this context?

## closed as off-topic by Daniel Lichtblau, MarcoB, m_goldberg, Yves Klett, FeyreNov 26 '16 at 15:47

• The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

• Did you get any answer at all? I got the error but also got the answer {1.79492*10^10} – Jose Enrique Calderon Nov 25 '16 at 7:28
• @JoseECalderon Yes. I Do. But I am not sure that is the correct answer. – Farhad Nov 25 '16 at 7:34
• The integrand in your code is different from that in $\LaTeX$, which one is correct? (Just press Ctrl+Shift+N and check your code, the (y-0.001) term is in the numerator. ) – xzczd Nov 25 '16 at 12:17
• @xzczd I corrected the mathematica code. – Farhad Nov 25 '16 at 12:44
• I'm voting to close this question as off-topic because the behavior indicated is correct: the integral diverges, and for a reason that is fairly clear. – Daniel Lichtblau Nov 25 '16 at 17:39

The integral under consideration diverges in view of the singular line $y=1/1000$ of the integrand. The Integrate command detects it:

Integrate[ Sqrt[x Exp[-100 x]] Sqrt[ y Exp[-100 x]] (1 - Exp[1000 I (x + y)])/((x + y) (y - 0.001)), {x, 0, \[Infinity]}, {y, 0, \[Infinity]}]

Integrate::idiv: Integral of (0.177245 Sqrt[y])/(-0.001+y)-((0.0414551 +0.0375163 I) E^(1000 I y) Sqrt[y])/(-0.001+y)+(3.14159 E^(100 y) y Erf[10 Sqrt[y]])/(-0.001+y)-(3.14159 E^(100 y) y Erf[10 Sqrt[(1-10 I) y]])/(-0.001+y) does not converge on {0,[Infinity]}.

Few syntax errors (e.g. E^{-100x}) are corrected in the above.

Even with PrincipalValue -> True

Integrate[ Sqrt[x Exp[-100 x]] Sqrt[y Exp[-100 x]] (1 - Exp[1000 I (x + y)])/((x + y) (y - 0.001)), {y,0, \[Infinity]}, {x, 0, \[Infinity]}, PrincipalValue -> True]

Integrate[ 1/(0.001 + 1. x) E^(-100 x) ((0. - 0.0993459 I) Sqrt[x] + 3.14159 x + E^(1000 I x) ((0.106478 - 0.0123898 I) Sqrt[x] - 3.14159 x Cos[1000. x] + (0. + 3.14159 I) x Sin[1000. x] + FresnelC[ 25.2313 Sqrt[ x]] ((3.14159 - 3.14159 I) x Cos[ 1000. x] - (3.14159 + 3.14159 I) x Sin[1000. x]) + FresnelS[ 25.2313 Sqrt[ x]] ((3.14159 + 3.14159 I) x Cos[ 1000. x] + (3.14159 - 3.14159 I) x Sin[1000. x]))), {x, 0, [Infinity]}, PrincipalValue -> True] Then

NIntegrate[1/(0.001 + 1. x) E^(-100 x) ((0. - 0.0993459 I) Sqrt[x] + 3.14159 x +E^(1000 I x) ((0.106478 - 0.0123898 I) Sqrt[x] -
3.14159 x Cos[1000. x] + (0. + 3.14159 I) x Sin[1000. x] +
FresnelC[25.2313 Sqrt[x]] ((3.14159 - 3.14159 I) x Cos[1000. x] - (3.14159 + 3.14159 I) x Sin[1000. x]) + FresnelS[ 25.2313 Sqrt[
x]] ((3.14159 + 3.14159 I) x Cos[1000. x] + (3.14159 - 3.14159 I) x Sin[1000. x]))), {x,  0, \[Infinity]}]

0.0257157 - 0.00922848 I

• Verify.. There is a mistake in the equation. Left the y-0.001 in the numerator. – Jose Enrique Calderon Nov 25 '16 at 9:13
• @Jose E Calderon: Thank you. The incorrect addition was deleted by me. – user64494 Nov 25 '16 at 9:23
• Which version are you in? In v9.0.1 and v8.0.4 the 2nd integral returns unevaluated with a few warning. – xzczd Nov 25 '16 at 14:33
• @xczd: I use v. 11. You are right. Indeed, the true output is Integrate[ 1/(0.001 + 1. x) E^(-100 x) ((0. - 0.0993459 I) Sqrt[x] + 3.14159 x + E^(1000 I x) ((0.106478 - 0.0123898 I) Sqrt[x] - 3.14159 x Cos[1000. x] + (0. + 3.14159 I) x Sin[1000. x] + FresnelC[ 25.2313 Sqrt[ x]] ((3.14159 - 3.14159 I) x Cos[1000. x] - (3.14159 + 3.14159 I) x Sin[1000. x]) + FresnelS[25.2313 Sqrt[ x]] ((3.14159 + 3.14159 I) x Cos[ 1000. x] + (3.14159 - 3.14159 I) x Sin[1000. x]))), {x, 0, [Infinity]}, PrincipalValue -> True] Thank you. I edited my answer. – user64494 Nov 25 '16 at 14:56
• Thank you. It seems that the "PrincipalValue" option works for two-dimensional Integrate not NIntegrate. This would do. – Farhad Nov 25 '16 at 18:21

The integrand has a singularity in $y=0.001$, but also have terms than can be separated and integrated firstly wrt x. The integrad is:

Sqrt[x Exp[-100 x]] Sqrt[y Exp[-100 y]] (1 - Exp[1000 I (x + y)])/((x + y) (y - 0.001))

Now, and with $y\geq0$, first integrate the real part:

Assuming[y >= 0, Integrate[(1 - Cos[1000 (x + y)]) Sqrt[x] Exp[-50. x], {x, 0, \[Infinity]}]]

(*0.00250663 + 0.0000182432 Cos[1000. y] + 0.0000212049 Sin[1000. y]*)

with respect to $y$ and taking into account the singularity at 0.001:

NIntegrate[1/(-0.001+ y) Sqrt[y] Exp[-50 y] (0.0025066282746310014 +
0.000018243204055544903 Cos[1000. y] +
0.000021204895736310885 Sin[1000. y]), {y, 0, .001, \[Infinity]},
Method-> "DoubleExponential", MaxRecursion -> 200]

(*0.000567425*)

If we do the same for the imaginary part:

Assuming[y >=0, -I Integrate[( Sin[1000 (x + y)]) Sqrt[ x] Exp[-50. (x)], {x, 0, \[Infinity]}]]

(*-I (0.0000212049 Cos[1000 y] - 0.0000182432 Sin[1000 y])*)

NIntegrate[1/(-0.001 + y) Sqrt[y]
Exp[-50 y] (-I (0.000021204895736310888 Cos[1000 y] -
0.00001824320405554491 Sin[1000 y])), {y, 0, .001, \[Infinity]},
Method -> "DoubleExponential", MaxRecursion -> 200]

(*2.21313*10^-6 I*)

The value seems to be 0.000567425+2.21313$10^{-6}i$

• @ José Antonio Díaz Nava : This is another task than the stated one because you integrate with respect to $y$ from $0.001$ to $\infty$. – user64494 Nov 25 '16 at 16:09
• @user64494 You are right. I have just edited to take into account the singularity in $y$ – José Antonio Díaz Navas Nov 25 '16 at 16:21
• @ José Antonio Díaz Nava : Can you kindly comment your code, especially the following places {y, 0, 0.001, [Infinity]} and "SymbolicProcessing" -> 5? A good code is a commented code. I am not an expert in Mathematica. TIA. – user64494 Nov 25 '16 at 16:46
• Basically you can include the singularities in the integration domain. The number in SymbolicProcessing option tells Mathematica how many seconds for preprocessing the integrand symbolically. Please write my name correctly. – José Antonio Díaz Navas Nov 25 '16 at 16:51
• @ José Antonio Díaz Navas : I am very sorry for the incorrect spelling of your name. Thank you. I did find it (i.e. "Basically you can include the singularities in the integration domain" ) in reference.wolfram.com/language/tutorial/… . However, the integral under consideration diverges. Is a certain generalized value of it calculated in your answer? – user64494 Nov 25 '16 at 17:05