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]]), 
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.

  • $\begingroup$ Unable to test right now, but does adding Method->{Automatic,“SymbolicProcessing“->0} to your NIntegrate speed things up? $\endgroup$ – Lukas Apr 17 '16 at 6:38
  • $\begingroup$ Yes, it does. Thank you! If you post this as a reply, I will accept it. $\endgroup$ – Felix Apr 17 '16 at 15:08

It appears that by far most time is spent on symbolic precomputations during NIntegrate. You can avoid this by calling the NIntegratecommand with the option Method -> {Automatic, "SymbolicProcessing" -> 0}.

fs[25., 25., 0.01] // AbsoluteTiming
fsold[25., 25., 0.01] // AbsoluteTiming
(* {0., 695.762} *)
(* {0., 695.762} *)

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.