# numerical integration over divergence

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

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.

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

• +1 for solving problem #2 and thanks for pointing out the issue of using I as a variable. – Felix Mar 28 '16 at 16:38
• 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. – Felix Mar 28 '16 at 21:47