numerical integration over divergence

Question

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}]

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

Answer 1

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:

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