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:
- The evaluation of
I306
is very slow (about 30 minutes on my laptop) 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}]
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.