Message Boards Message Boards

Solve convergence problem in two-dimentional NIntegrate?

Posted 7 years ago

I need to calculate the principal value of following integral

 l=1/10;
 k=1/10;
 NIntegrate[Sqrt[w^2 E^(-w/l)]Sqrt[v^2 E^(-v/k)](1-E^(1000I(w+v)))/((w+v)(v-0.01)),{w,0,?},{v,0,?}]

But it gives the slwcon error:

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

I searched all over the internet. I tried different numerical integration methods (GlobalAdaptive, LocalAdaptive, ...) with different options (AccuracyGoal,PrecisionGoal, ...) but none of them get rid of the error. The integral for very small values of k and l converges and gives the result without error. This means that the integral is convergent. But I need the result for the larger values of k and l as specified in the code above. How should I resolve the error?

POSTED BY: F TQ
3 Replies

MMA can find symbolic solution for this double integral:

ClearAll["Global`*"];
Remove["Global`*"];

f[w_, v_] :=  Sqrt[w^2 E^(-w/l)] Sqrt[v^2 E^(-v/k)] (1 - E^(1000 I (w + v))/(w + v) (v - 1/100))

sol = Integrate[f[w, v], {w, 0, Infinity}, {v, 0, Infinity}, Assumptions -> {l, k} > 0, PrincipalValue -> True]

(* (2 k^2 l^2 (800000000 k^6 (I + 2000 l) - 
     l^3 (I + 200 I l + 400000 l^2) - 
     800000 k^5 (1 + 2000 I l + 8000000 l^2) + 
     k l^2 (I - (6000 - 1200 I) l + 800000 l^2 + 1600000000 I l^3) + 
     2000 k^4 (-1 + (800 + 4000 I) l + 1600000 I l^2 + 
        4800000000 l^3) + 
     k^2 l (I + (10000 - 600 I) l + (4400000 + 8000000 I) l^2 - 
        5600000000 I l^3 + 1600000000000 l^4) - 
     k^3 (I + (2000 + 400 I) l + (5600000 + 16000000 I) l^2 - 
        1600000000 I l^3 + 6400000000000 l^4) - 
     2 k (I + 2000 k)^2 l (I + 2000 l) (k - l + 600 k l) Log[(
       k (I + 2000 l))/((I + 2000 k) l)]))/(25 (I + 2000 k)^2 (k - 
     l)^4 (I + 2000 l))*)


     sol1 = Limit[(sol /. k -> 1/10), l -> 1/10]
     (* 153615360523216720031/96009600360006000037500 - (3190240234 I)/960096003600060000375*)
      N[sol1, 39]
      (* 0.00159999999945012915551108225908344037372 - 3.32283461449437974587676677293*10^-12 I  *)

Maple solution:

enter image description here

Numeric solution:

ClearAll["Global`*"];
Remove["Global`*"];
l = 1/10; 
(* Dont use 0.1. Better is using  exact number or alternative
 SetPrecision[0.1, 100] ,or 0.1`100 you may increase the value 100 for better result *)
k = 1/10;
f[w_, v_] := Sqrt[w^2 E^(-w/l)] Sqrt[v^2 E^(-v/k)] (1 - E^(1000 I (w + v))/(w + v) (v - 1/100))
NIntegrate[f[w, v], {w, 0, Infinity}, {v, 0, Infinity}, Method -> "GlobalAdaptive", WorkingPrecision -> 19]
(*0.001599999999406708480 - 3.323499234519002782*10^-12 I*)

For better result increase WorkingPrecision.

POSTED BY: Mariusz Iwaniuk
Posted 7 years ago

Sorry. By looking at the Maple code I realized that the input code was incorrect. I corrected the code. Also, I posted the code in the informal Mathematica forum, here. Please analyze this code.

POSTED BY: F TQ

MMA can find symbolic solution to this function. :)

ClearAll["Global`*"];
Remove["Global`*"];

f[x_, y_] := Sqrt[x Exp[-x/a]] Sqrt[y Exp[-y/a]] (1 - Exp[c*I (x + y)])/((x + y) (y - b))

solving for x:

sol = Integrate[f[x, y], {x, 0, Infinity}, Assumptions -> {{a, b, c, y} > 0, {a, b, c, y} \[Element] Reals}, PrincipalValue -> True]
(*(1/(Sqrt[1 - 2 I a c] (-b + y)))Sqrt[
 E^(-(y/a)) y] (Sqrt[a] Sqrt[1 - 2 I a c] Sqrt[2 \[Pi]] - 
   Sqrt[a] E^(I c y) Sqrt[2 \[Pi]] + 
   E^(y/(2 a)) \[Pi] Sqrt[y - 2 I a c y] - 
   E^(y/(2 a)) \[Pi] Sqrt[y - 2 I a c y]
     Erf[(Sqrt[1/a - 2 I c] Sqrt[y])/Sqrt[2]] - 
   Sqrt[1 - 2 I a c] E^(y/(2 a)) \[Pi] Sqrt[y]
     Erfc[Sqrt[y/a]/Sqrt[2]])*)
sol1 = FullSimplify[sol, Assumptions -> {a, b, y, c} > 0]
(*(1/(Sqrt[1 - 2 I a c] (-b + y)))Sqrt[
 E^(-(y/a)) y] (Sqrt[a] (Sqrt[1 - 2 I a c] - E^(I c y)) Sqrt[
    2 \[Pi]] + 
   E^(y/(2 a)) \[Pi] (Sqrt[y - 2 I a c y]
        Erfc[(Sqrt[1/a - 2 I c] Sqrt[y])/Sqrt[2]] - 
      Sqrt[1 - 2 I a c] Sqrt[y] Erfc[Sqrt[y/a]/Sqrt[2]]))*)

solving for y:

  sol2 = Integrate[sol1, {y, 0, Infinity}, Assumptions -> {{a, b, c} > 0, {a, b, c} \[Element] Reals}, PrincipalValue -> True]
 (* 2 a \[Pi] - 
  2 Sqrt[2] Sqrt[a] Sqrt[b] \[Pi] DawsonF[Sqrt[b/a]/Sqrt[2]] + (
  Sqrt[a] Sqrt[
   2 \[Pi]] (-(Sqrt[2 \[Pi]]/Sqrt[1/a - 2 I c]) + 
     Sqrt[b] E^(-(b/(2 a)) + I b c) \[Pi] Erfi[(
       Sqrt[b] Sqrt[1/a - 2 I c])/Sqrt[2]]))/Sqrt[
  1 - 2 I a c] - (\[Pi] (b^2 HypergeometricPFQ[{1, 1}, {3/2, 2}, -(b/(
        2 a))] + 
     a (a - b EulerGamma - b Log[2] + b Log[a] - b Log[b])))/a + (1/(
  a (1 - 2 I a c)))\[Pi] (-b^2 (I + 2 a c)^2 HypergeometricPFQ[{1, 
        1}, {3/2, 2}, -(b/(2 a)) + I b c] + 
     a (a - b EulerGamma + 2 I a b c EulerGamma - b Log[2] + 
        2 I a b c Log[2] + (b - 2 I a b c) Log[a] + 
        b (-1 + 2 I a c) Log[b] - b Log[1 - 2 I a c] + 
        2 I a b c Log[1 - 2 I a c]))*)
sol3 = FullSimplify[sol2, Assumptions -> {{a, b, c} > 0, {a, b, c} \[Element] Reals}]
(*\[Pi] ((2 a (I + a c - (I Sqrt[a] Sqrt[1/a - 2 I c])/Sqrt[
      1 - 2 I a c]))/(I + 2 a c) - 
   2 Sqrt[2] a Sqrt[b/a] DawsonF[Sqrt[b]/(Sqrt[2] Sqrt[a])] + (
   2 Sqrt[2] Sqrt[a] Sqrt[b]
     DawsonF[(Sqrt[b] Sqrt[1/a - 2 I c])/Sqrt[2]])/Sqrt[
   1 - 2 I a c] - (
   b^2 HypergeometricPFQ[{1, 1}, {3/2, 2}, -(b/(2 a))])/a + (
   b^2 (1 - 2 I a c) HypergeometricPFQ[{1, 1}, {3/2, 2}, -(b/(2 a)) + 
      I b c])/a - b Log[1 - 2 I a c])*)

enter image description here

sol3 /. a -> 100 /. b -> 1/1000 /. c -> 1000 // N
(* 314.121 + 0.00134748 I *)
POSTED BY: Mariusz Iwaniuk
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract