# How to speed up this code with Bessel functions

I need to perform a ContourPlot on a function with BesselJ included. But my code is very slow, especially when I set the PlotPoints to be 100. I would like to know how to speed up this code. Any help or suggestion will be highly appreciated! The code is shown below：

λ = 0.500;
k = (2 π)/λ;
ρ[x_, y_] := Sqrt[x^2 + y^2];
ϕ[x_, y_] := ArcTan[x, y];
Ε2[n_, x_, y_, z_] := Exp[I k z + I n ϕ[x, y] + (I k ρ[x, y]^2)/(4 z)]* Sqrt[π/2]*(-I)^(Abs[n]/2)*Sqrt[(k*ρ[x, y]^2)/(4 z)]*(BesselJ[(Abs[n] - 1)/2, (k ρ[x, y]^2)/(4 z)] - I*BesselJ[(Abs[n] + 1)/2, (k ρ[x, y]^2)/(4 z)]);
U[α_, x_, y_, z_] := (Exp[I π α]*Sin[π α])/π*\!$$\*UnderoverscriptBox[\(∑$$, $$n = \(-50$$\), $$50$$]\*FractionBox[$$Ε2[n, x, y, z]$$, $$α - n$$]\);
ContourPlot[Arg[U[4.5, x, y, 30]], {x, -15, 15}, {y, -15, 15},   Frame -> True, PlotPoints -> 10, PlotRange -> All,   ContourLines -> False, ColorFunction -> "Rainbow",    PlotLegends -> Automatic, Contours -> 20,    AspectRatio -> Automatic] // Timing


The time cost is about 69 seconds when the PlotPoints is set as 10 and n set as {-50, 50} , as shown in the figure below. However, the time cost is over 10 minitues, if the PlotPoints is 100.

Here's a refactoring of the OP's code & Rolf's compile idea. Basically I tried to implement the idea in my comment:

Also, BesselJ is not compilable, so it slows down the compiled function with call-backs to the main kernel. You could pass a list of Bessel function values as an argument to Compile. This would reduce the number of Bessel function calls by almost half, too.

It actually reduces the number of Bessel function calls by about three quarters.

I also computed the Bessel function values from the backward recursion, https://dlmf.nist.gov/10.6#E1. (The forward recursion is numerically unstable.) This is four to five times faster for the OP's n = 50. If you would rather use BesselJ instead of the recursively computed values, set $useBesselJ equal to True (see the definitions of bj[]). (* Compute J[n,x] via backward recursion dlmf.nist.gov/10.6#E1 *) padJC = Compile[{{jarray, _Real, 1}, {x, _Real}, {nmax, _Real}}, Join[ Reverse@Module[{n = nmax - Length@jarray + 1}, Rest[ NestList[(n--; #) &[{{2. n /x, -1.}.#, First@#}] &, jarray[[{1, 2}]], (* initial values for recursion *) Ceiling@n (* round up if half integer *) ] ][[All, 1]] ], jarray ], (*CompilationTarget -> "C",*) (* minor speed improvement *) RuntimeOptions -> "Speed" ]; ClearAll[bx, bj]; (* OP's Bessel function argument *) bx[x_?NumericQ, y_?NumericQ,z_?NumericQ] := (k ρ[x, y]^2)/(4*z); (* Computes list of Bessel function values for E2[] *) bj[x_?NumericQ, y_?NumericQ, z_?NumericQ] /; TrueQ[$useBesselJ] :=
BesselJ[Range[-1, nn + 1]/2, (k ρ[x, y]^2)/(4*z)];
(* Computes list of Bessel function values for E2[] *)
bj[x_?NumericQ, y_?NumericQ, z_?NumericQ] := bj@bx[x, y, z];
bj[xx_] /; EvenQ[nn] := Riffle[
padJC[BesselJ[Range[nn - 1, nn + 1, 2]/2, xx], (* half-integer orders *)
xx, (nn + 1)/2],
padJC[BesselJ[Range[nn - 2, nn, 2]/2, xx], xx, nn/2] (* integer orders *)
];
bj[xx_] /; OddQ[nn] := Riffle[
padJC[BesselJ[Range[nn - 2, nn, 2]/2, xx],     (* half-integer orders *)
xx, nn/2],
padJC[BesselJ[Range[nn - 1, nn + 1, 2]/2, xx], xx, (nn + 1)/2]
];


Check the Bessel function values. The max relative error is a little over $$10^{-14}$$, which is good enough for plotting.

xx = bx[1., 2., 30]; (* test value *)
nn = 50;             (* range for n *)
bj[1., 2., 30] // AbsoluteTiming;
BesselJ[Range[-1, nn + 1]/2, xx] // AbsoluteTiming;
(% - %%)/%% // Last // Abs // Max
First /@ {%%%, %%}
(*
1.93447*10^-14       <-- Max relative error
{0.00016, 0.000702}  <-- Timings {bj, BesselJ}
*)


The maximum relative error illustrated with the $useBesselJ flag: (bj[1., 2., 30] - #)/# &@ Block[{$useBesselJ = True}, bj[1., 2., 30]] // Abs // Max
(*  1.93447*10^-14  *)


I had to refactor E2[] a bit. I separated the real/complex and the scalar/vector operations.

U = Compile[
{{α, _Real}, {x, _Real}, {y, _Real}, {z, _Real}, {jn, _Real, 1}},
Block[{n = Range[-(Length@jn - 3), Length@jn - 3], e1 = 0.,
e2 = 0. I, jm = {0.}, jp = {0.}},
With[{a = Drop[jn, -2]},   (* order (Abs[n]-1)/2 *)
jm = Reverse@a ~Join~ Rest@a];
With[{a = Drop[jn, 2]},    (* order (Abs[n]+1)/2 *)
jp = Reverse@a ~Join~ Rest@a];
e1 = (Sin[Pi*α])*Sqrt[π/2]*
Sqrt[(k*(x^2 + y^2))/(4 z)]/Pi;
e2 = Total[
Exp[
I k z + I n ArcTan[x, y] + (I k (x^2 + y^2))/(4 z)]*(-I)^(Abs[
n]/2.)*(jm - I*jp)/(α - n)];
e1*Exp[I*Pi*α]*e2
],
(*CompilationTarget -> "C",*)
RuntimeOptions -> "Speed",
CompilationOptions -> {"InlineExternalDefinitions" -> True}];


The plot for n equal to 50, MaxRecursion -> 3:

nn = 50;
ContourPlot[
Arg[U[4.5, x, y, 30, bj[x, y, 30]]],
{x, -15, 15}, {y, -15, 15}, Frame -> True, PlotPoints -> 15,
MaxRecursion -> 3, PlotRange -> All, ContourLines -> False,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
Contours -> 20, AspectRatio -> Automatic,
ImageSize -> 350] // AbsoluteTiming The OP's 70-sec. plot takes 2 sec. with this approach.

• Many thanks to Michael E2. The new code is really very useful. Feb 19, 2020 at 5:54
• @user14634 You're welcome. Thank you for the accept. Feb 19, 2020 at 13:23

Using Compile speeds things up:

λ = 0.500;
k = (2 π)/λ;
ρ[x_, y_] := Sqrt[x^2 + y^2];
ϕ[x_, y_] := ArcTan[x, y];
Ε2[n_, x_, y_, z_] := Exp[I k z + I n ϕ[x, y] + (I k ρ[x, y]^2)/(4 z)]* Sqrt[π/2]*(-I)^(Abs[n]/2)*Sqrt[(k*ρ[x, y]^2)/(4 z)]*(BesselJ[(Abs[n] - 1)/2, (k ρ[x, y]^2)/(4 z)] - I*BesselJ[(Abs[n] + 1)/2, (k ρ[x, y]^2)/(4 z)]);

U = With[{expr = ((Exp[I*Pi*\[Alpha]]*Sin[Pi*\[Alpha]])*
Sum[\[CapitalEpsilon]2[n, x, y, z]/(\[Alpha] - n), {n, -10, 10}])/Pi
},
Compile[{
{\[Alpha], _Real}, {x, _Real}, {y, _Real}, {z, _Real}},
expr]
];
ContourPlot[Arg[U[4.5, x, y, 30]], {x, -15, 15}, {y, -15, 15},   Frame -> True, PlotPoints -> 10, PlotRange -> All,   ContourLines -> False, ColorFunction -> "Rainbow",    PlotLegends -> Automatic, Contours -> 20,    AspectRatio -> Automatic] // Timing

• Many thanks to Rolf Mertig. Compile really speeds things up. However, when I change the PlotPoints to 100, the code still runs very slow. Is it possile to increase the speed when PlotPoints -> 100? Feb 17, 2020 at 12:20
• A severe problem with compile code: If the number of n is increased from {n, -10, 10} to {n, -50, 50}, the ContourPlot figure is obviously different from the figure obtained using the code without complie. Is it possiple to overcome this problem? Feb 17, 2020 at 13:50
• I updated the code by using {n, -50, 50}. Feb 18, 2020 at 14:22
• @user14634 It may be a problem with Sum, which tries symbolic methods; try replacing Sum with Total@Table. Also, BesselJ is not compilable, so it slows down the compiled function with call-backs to the main kernel. You could pass a list of Bessel function values as an argument to Compile. This would reduce the number of Bessel function calls by almost half, too. Feb 18, 2020 at 17:48
• Many thanks to Michael E2 for these helpful comments. Feb 19, 2020 at 6:01