# NDSolve very slow

I am trying to run this code to get the output from NDSolve for 4 dynamical equations. However, I have two issues:

(1) It is taking a very long time (can reach several hours for a given instance) although the problem is not that difficult and shouldn't take that long (!). Also it seems to act slower/faster if I vary the value of Scalar1, but it is not clear if there is any linear trend in this.

Any advice would be appreciated on how to speed it up.

In the code below, please keep k=0 in function input variables, and import (example, using InterpEr=Import["InterpEr.wdx"]) the three uploaded interpolation functions that will be called in the code below (download from here).

For example try Scalar1 values between 1 and 1000, and observe the long delays in NDSolve.

c = 3*10^8;
λ = c/(1.1424*10^10);
Mu0 = 4 *Pi* 10^-7;
m0 = 9.1093837015*10^-31;
e = 1.602176634*10^-19;
n = 10;
m = 100;
omegaHz = 2*Pi (1.1424*10^10);
Δϕ = 2 Pi/m;

InterpEr = Import["InterpEr.wdx"]
InterpEz = Import["InterpEz.wdx"]
InterpHphi = Import["InterpHphi.wdx"]

Er[r_?NumberQ, t_?NumberQ, z_?NumberQ, k_?NumberQ] := Module[{er},
er = InterpEr[r, z];
er Cos[omegaHz t + k Δϕ]
]
Ez[r_?NumberQ, t_?NumberQ, z_?NumberQ, k_?NumberQ] := Module[{ez},
ez = InterpEz[r, z];
ez Cos[omegaHz t + k Δϕ]
]
Hphi[r_?NumberQ, t_?NumberQ, z_?NumberQ, k_?NumberQ] := Module[{hphi},
hphi = InterpHphi[r, z];
hphi Cos[omegaHz t + k Δϕ]
]

Scalar1 = 1;
ErScaled[r_?NumberQ, t_?NumberQ, z_?NumberQ, k_?NumberQ] :=
Scalar1 Er[r, t, z, k]/((m0*c^2)/(e λ));
EzScaled[r_?NumberQ, t_?NumberQ, z_?NumberQ, k_?NumberQ] :=
Scalar1 Ez[r, t, z, k]/((m0*c^2)/(e λ));
BphiScaled[r_?NumberQ, t_?NumberQ, z_?NumberQ, k_?NumberQ] :=
Scalar1 (Mu0*Hphi[r, t, z, k])/((m0*c)/(e λ));

ss = NDSolve[{D[tt[zz], zz] == γ[zz]/
Sqrt[γ[zz]^2 - pr[zz]^2 - 1],
D[rr[zz], zz] == pr[zz]/Sqrt[γ[zz]^2 - pr[zz]^2 - 1],
D[γ[zz], zz] ==
D[rr[zz], zz] ErScaled[rr[zz], tt[zz], zz, 0] +
EzScaled[rr[zz], tt[zz], zz, 0],
D[pr[zz], zz] ==
D[tt[zz], zz] ErScaled[rr[zz], tt[zz], zz, 0] -
BphiScaled[rr[zz], tt[zz], zz, 0], rr[0.0] == 0,
tt[0.0] == 0, γ[0.0] == 1.1174, pr[0.0] == 0}, {rr, tt,
pr, γ}, {zz, 0.0, 1.76}, MaxSteps -> ∞]


Please note that this is my attempt to use NDSolve to solve the following set of equations:

$$\dfrac{\mathrm{d}t(z)}{\mathrm{d}z} = \dfrac{\gamma(r, t; z)}{\sqrt{\gamma(r, t; z)^{2}-p_{r}(r, t; z)^{2}-1}}$$

$$\dfrac{\mathrm{d}r(z)}{\mathrm{d}z} = \dfrac{p_{r}}{\sqrt{\gamma(r, t; z)^{2}-p_{r}(r, t; z)^{2}-1}}$$

$$\dfrac{\mathrm{d}\gamma(r, t; z)}{\mathrm{d}z} = \dfrac{\mathrm{d}r(z)}{\mathrm{d}z} E_{r}(r, t; z) + E_{z}(r, t; z)$$

$$\dfrac{\mathrm{d}p_{r}(r, t; z)}{\mathrm{d}z} = \dfrac{\mathrm{d}t(z)}{\mathrm{d}z} E_{r}(r, t; z) - B_{\phi}(r, t; z)$$

which have $$E_{r}, E_{z}, B_{\phi}$$ as the given (uploaded) functions and the following 4 dependent variables ($$z$$ is the independent variable):

$$t = t(z)$$

$$r = r(z)$$

$$\gamma = \gamma(r, t; z)$$

$$p_{r} = p_{r}(r, t; z)$$

(2) Was my choice to solve this set using NDSolve a good choice? Are there any other/better methods (other than NDSolve) to solve this set of equations efficiently?

So far, the output from NDSolve seems to me like numerical noise. For example, after NDSolve, if we plot the function $$\gamma$$, which had initial value 1.1174 at $$z$$=0.0, we see that its value is fluctuating around the same value (it should be increasing or decreasing with $$z$$):

Plot[γ[zz] /. ss[[1]], {zz, 0.01, 1.75},
AxesLabel -> {"z", "γ[z]"}, PlotLabel -> "Scalar" Scalar1]


I am trying to increase Scalar1 in the hope (convincing myself?) to find out that it may be a matter of scaling, but I am starting to wonder whether something else could be wrong in this code/method?

Thanks for any suggestions

• It takes about 15 seconds on my machine. Also you don't need to use Module - for example write: er Cos[omegaHz t + k Δϕ] as InterpEz[r, z] Cos[omegaHz t + k Δϕ] – flinty May 25 '20 at 15:48
• @flinty Just 15 seconds! That's quite different from what I get. Did you try to do it a few times with different Scalar1 values between 1 and 1000? And did you mean that Module is a potential reason for slow down? I used it just to make sure the variables inside it are local and not mixed with any other global variables used before or later. – user135626 May 25 '20 at 23:22
• Which version of Mathematica are you using? I'm using 12.1. Module will be a little slower - you don't need to localize variables as you can inline them into the expressions. – flinty May 26 '20 at 10:55