# 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. – user14634 Feb 19 at 5:54
• @user14634 You're welcome. Thank you for the accept. – Michael E2 Feb 19 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? – user14634 Feb 17 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? – user14634 Feb 17 at 13:50
• I updated the code by using {n, -50, 50}. – user14634 Feb 18 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. – Michael E2 Feb 18 at 17:48
• Many thanks to Michael E2 for these helpful comments. – user14634 Feb 19 at 6:01