How to speed up this code with Bessel functions

Question

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.

Answer 1

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.

Answer 2

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