Using mathematica version 10.3, I observe that an integral over a piecewise function evaluates a factor 10 slower than the same integral without the piecewise distinction. Here is the code:
Originally, I have defined the function g1
and the integral over g1
, fs
, as
Kaold[k_] := EllipticK[k^2/(-1 + k^2)]/Sqrt[1 - k^2]
kold[p_, y_, z_] := Sqrt[4 y z/(p^2 + (y + z)^2)]
g1old[x_, y_] :=
1/Pi Sqrt[1/(
x y)] (kold[0, x, y] Kaold[kold[0, x, y]] -
kold[1, x, y] Kaold[kold[1, x, y]])
fsold[r_, r0_, d_] := NIntegrate[x g1old[x d, r d], {x, 0, 20}]
However, g1
in this definition has issues with floating point precision when x-y
becomes small. Therefore, I developed an alternative representation
Ka[k_] := EllipticK[k^2/(-1 + k^2)]/Sqrt[1 - k^2]
Kaomks[x_] := EllipticK[(x - 1)/x]/Sqrt[x]
k0[z_] := Sqrt[(4 z)/(1 + z)^2]
k1[y_, z_] := Sqrt[(4 y z)/(1 + (y + z)^2)]
omks[p_, y_, z_] := (p^2 + (y - z)^2)/(p^2 + (y + z)^2)
g1[x_, y_] := \[Piecewise] {
{1/Pi Sqrt[1/(
x y)] (k0[x/y] Kaomks[omks[0, x, y]] -
k1[x, y] Kaomks[omks[1, x, y]]), k0[x/y] > 0.5},
{1/Pi Sqrt[1/(x y)] (k0[x/y] Ka[k0[x/y]] - k1[x, y] Ka[k1[x, y]]),
True}
}
fs[r_, r0_, d_] := NIntegrate[x g1[x d, r d], {x, 0, 20}]
Both functions yield the same result, if evaluated with sufficient precision. The second version is much more robust against finite numerical precision. However, unfortunately, the integral fs
of the second approach evaluated much slower than in the initial representation. On my PC,
fs[25., 25., 0.01] // Timing
takes 0.14 seconds and
fsold[25., 25., 0.01] // Timing
takes 0.01 seconds. I believe that this is due to the piecewise distinction in the function g1
because if I choose only one of the piecewise cases
g2[x_, y_] :=
1/Pi Sqrt[1/(
x y)] (k0[x/y] Kaomks[omks[0, x, y]] -
k1[x, y] Kaomks[omks[1, x, y]])
fs2[r_, r0_, d_] := NIntegrate[x g2[x d, r d], {x, 0, 20}]
g3[x_, y_] :=
1/Pi Sqrt[1/(x y)] (k0[x/y] Ka[k0[x/y]] - k1[x, y] Ka[k1[x, y]])
fs3[r_, r0_, d_] := NIntegrate[x g3[x d, r d], {x, 0, 20}]
then both fs2
and fs3
evaluate within 0.01 seconds (and yield the same result).
How can I improve the performance of the piecewise function? I would like to avoid splitting fs
in two separate integrals, where each integral contains only one of the two representations of g1
, because I am trying to evaluate fs
with a given precision goal and that is much easier if it is only one term.
Method->{Automatic,“SymbolicProcessing“->0}
to yourNIntegrate
speed things up? $\endgroup$ – Lukas Apr 17 '16 at 6:38