1
$\begingroup$

Consider the following function g1:

Ka[k_] := EllipticK[k^2/(-1 + k^2)]/Sqrt[1 - k^2]
k[p_, y_, z_] := Sqrt[4 y z/(p^2 + (y + z)^2)]
g1[x_, y_] := 
 1/Pi Sqrt[1/(
   x y)] (k[0, x, y] Ka[k[0, x, y]] - k[1, x, y] Ka[k[1, x, y]])

which has a divergence at x=y. The divergence is weak, i.e., the integral over the divergence is finite.

Now I would like to calculate the two-fold integral over g1 multiplied by the following function mz:

t[x_, y_] := 2 ArcTan[Exp[(x - y)]] + 2 ArcTan[Exp[(x + y)]]
mz[x_, y_] := Cos[t[x, y]]

The integrals are defined as f1 (first integral) and I306 (second integral):

f1[x_?NumericQ, y_?NumericQ, d_?NumericQ] := 
 d^2 (NIntegrate[z g1[z d, x d] (1 - mz[z, y]), {z, 0, x}, 
     PrecisionGoal -> 6, 
     Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100000, 
       Method -> {"GaussKronrodRule", "Points" -> 3}}, 
     MaxRecursion -> 20, WorkingPrecision -> 30] + 
    NIntegrate[z g1[z d, x d] (1 - mz[z, y]), {z, x, y + 15}, 
     PrecisionGoal -> 6, 
     Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100000, 
       Method -> {"GaussKronrodRule", "Points" -> 3}}, 
     MaxRecursion -> 20, WorkingPrecision -> 30])

h[x_, y_, d_] := x (Re[f1[x, y, d]] + mz[x, y] - 1) mz[x, y]
I306[y_, d_] := 
 NIntegrate[h[x, y, d], {x, 0, y + 15}, PrecisionGoal -> 3, 
  Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100000, 
    Method -> {"GaussKronrodRule", "Points" -> 3}}, 
  MaxRecursion -> 20, WorkingPrecision -> 30]

My problems are, illustrated for y=10 and d=100:

  1. The evaluation of I306 is very slow (about 30 minutes on my laptop)
  2. The evaluation of f1 yields very noisy data, despite the requested precision of 6 digits. This can be see from the following plot:

    Plot[h[x, 10, 100], {x, 0, 25}]

h

  1. f1 has sometimes a finite imaginary part, despite the purely numerical integration of a real function.

My attempts to address the problem are reflected in my definitions of f1 and I306. Specifically, I have noticed an increased calculation speed by splitting the f1 integral into two parts such that the divergence is at the integral boundaries. Also, GaussKronrodRule with 3 points seemed to perform best of all methods.

My specific question is: Are there better strategies to solve this integral?

Edit: Renamed I to I306 to avoid conflicts with mathematica's internal variables.

$\endgroup$
3
$\begingroup$

The problem with your Plot is actually unrelated to the later integration. Plot supplies a machine precision value, so your integral in f1 can not be done at the requested precision. If you do this (takes a while)

ListPlot[{#, h[#, 10, 100]} & /@ (Range[1, 20, 5/200]), 
 PlotRange -> All, Joined -> True]

you get a nice plot:

enter image description here

As another observation, for some reason you need to specify a higher working precision in the outer NIntegrate than in the inner ones. WorkingPrecision->32 seems to help. Its still mighty slow though.

Just using default method:

int[y_, d_] := 
 NIntegrate[h[x, y, d], {x, 0, y + 15}, PrecisionGoal -> 3, 
  WorkingPrecision -> 32]

takes about 5 minutes. If you really only want a low precision approximation you can set maxrecursion to 1 or 2 and cut the time in half.

It also helps a little if you first simplify the integrand:

integrand[z_, d_, x_, y_] = Simplify[z g1[z d, x d] (1 - mz[z, y])];

This must be before (not inside) f1 and is not a SetDelayed := .

Also using I is asking for trouble..

$\endgroup$
  • $\begingroup$ +1 for solving problem #2 and thanks for pointing out the issue of using I as a variable. $\endgroup$ – Felix Mar 28 '16 at 16:38
  • $\begingroup$ Actually, with the increased working precision, the integral I306 can be evaluated using an automatic choice of method: I306[y_, d_] := NIntegrate[h[x, y, d], {x, 0, y + 15}, PrecisionGoal -> 3, WorkingPrecision -> 32], which evaluates in 3 minutes instead of 30. Hence, I consider problem #1 solved. $\endgroup$ – Felix Mar 28 '16 at 21:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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