# Optimization for numerical integration of Airy function and using NIntegrate inside NDSolve

Sorry for the long post. I need to integrate a function which uses the trajectory data given by the Lorentz force equation. We have the electromagnetic field:

$$\vec{E}=E_0 \cos{\omega(t+\bar{z})}\text{exp}\left[-\frac{(t+\bar{z})^2}{\tau^2}\right]\hat{y}\\ \vec{B}=\frac{E_0}{c} \cos{\omega(t+\bar{z})}\text{exp}\left[-\frac{(t+\bar{z})^2}{\tau^2}\right]\hat{x}$$ with $$\bar{z}=z/c$$ and the force equation:

$$\dot{\bar{p}}_y =a_0\omega\cos{\omega (t+\bar{z})}\text{exp}\left[-\frac{(t+\bar{z})^2}{\tau^2}\right] \left(1+ \frac{\bar{p}_z}{\sqrt{1+\bar{p}^2}}\right)\\ \dot{\bar{p}}_z =-a_0\omega\cos{\omega (t+\bar{z})}\text{exp}\left[-\frac{(t+\bar{z})^2}{\tau^2}\right]\frac{\bar{p}_y}{\sqrt{1+\bar{p}^2}}\\ \dot{\bar{z}}=\frac{\bar{p}_z}{\sqrt{1+\bar{p}^2}}$$ The normalized momentum is denoted by $$\bar{p}=p/mc$$ and $$a_0$$ is the field strength parameter. I need to extract 3 parameters as a function of time and $$\epsilon_{\gamma}$$ from the numeric solution of above equations: $$\chi_1=\frac{\alpha}{\hbar \, \bar{p}^2},\quad \bar{p}^2=\bar{p}_y^2+\bar{p}_z^2\\ \chi_2=\left(\frac{E_s}{|\vec{E}|} \frac{\epsilon_{\gamma}}{\bar{p}-\epsilon_{\gamma}}\frac{1}{\bar{p}+\bar{p}_z}\right)^{2/3}\\ \chi_3=\frac{|\vec{E}|}{E_s}\left(1+\frac{\bar{p}_z}{\bar{p}}\right)\epsilon_{\gamma}$$ The numerical solution given by NDSolve

χ[ϵγ_?NumericQ, t_?NumericQ] := χ[ϵγ, t] =
With[{α := 1/137 , h := 105 10^(-36), Es := 13 10^(17),
E0 := 27 10^(13), a0 := 70,
ω := 22 10^(14), τ := 22 10^(-15), py0 := 0,
pz0 := 1200,
z0 := -100 10^(-15), t0 := 0}, {-α /(h (py[t]^2 + pz[t]^2)), ((
Es Exp[(t + z[t])^2/τ^2] ϵγ)/(
E0 Abs[Cos[(t + z[t]) ω]] (-ϵγ + Sqrt[
py[t]^2 + pz[t]^2]) (pz[t] + Sqrt[py[t]^2 + pz[t]^2])))^(
2/3), (ϵγ E0/Es Exp[-(t + z[t])^2/τ^2] Abs[
Cos[ω (t + z[t])]]  (1 +
pz[t]/Sqrt[py[t]^2 + pz[t]^2])) } /. NDSolve[{
py'[s] == -a0 ω  Cos[ω (s +
z[s])] Exp[- (s + z[s])^2/τ^2] (1 +
pz[s]/Sqrt[1 + py[s]^2 + pz[s]^2]),
pz'[s] ==
a0  ω  Cos[ω (s +
z[s])] Exp[-(s + z[s])^2/τ^2] py[s]/
Sqrt[1 + py[s]^2 + pz[s]^2],
z'[s] == pz[s]/Sqrt[1 + py[s]^2 + pz[s]^2], py[t0] == py0,
pz[t0] == pz0, z[t0] == z0}, {py, pz, z}, {s, t0,  t}][[1]]]



These $$\chi$$ functions are needed for the integration of Airy function:

$$dW(\epsilon_{\gamma},\, t)=\chi_1 \left(\int^{\infty}_{\chi_{2}}\text{Ai}(y)\, dy + \left(\frac{2}{\chi_2} + \chi_{3}\sqrt{\chi_2}\right) Ai'(\chi_2) \right)$$ The mathematica code for this is:

dWγ[ϵγ_?NumericQ, t_?NumericQ] :=
χ[ϵγ, t][[1]] ( NIntegrate[ AiryAi[s], {s, χ[ϵγ, t][[2]], Infinity}, AccuracyGoal -> 10] +
(2/χ[ϵγ, t][[2]] + χ[ϵγ, t][[3]] (χ[ϵγ, t][[2]])^(1/2)) AiryAiPrime[χ[ϵγ, t][[2]]])


The time plot for $$dW(\epsilon_{\gamma}, t)$$ is given by:

ListPlot[ParallelTable[
dWγ[100, t], {t, 0, 100 10^(-15), 1 10^(-16)}],
Joined -> True] // AbsoluteTiming


with time increment $$\delta t= 10^{-16}$$, it takes about 10 seconds to sample 1000 data points with the error message:

NIntegrate::izero :  Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option.

General::munfl :  Exp[-808.38] is too small to represent as a normalized machine number; precision may be lost.



I think this is to be expected, integral of Ai(x) evaluates to ridiculously small numbers, beyond what machine precision is capable of handling

Question 1: I need to numerically integrate $$dW(\epsilon_{\gamma}, t)$$ over the energy $$\epsilon_{\gamma}$$: $$W(\epsilon_m,\,\epsilon, \, t)=\int_{\epsilon_{m}}^{\epsilon} dW(\epsilon_{\gamma},\, t)\, d \epsilon_{\gamma}$$ Here $$\epsilon_m$$ is fixed and $$\epsilon$$ is variable. For this I used Nintegrate:

Wγ[ϵm_, ϵ_, t_] := 819*10^(-16) NIntegrate[dWγ[ϵγ, t], {ϵγ, ϵm, ϵ}]


$$819 \times 10^{-16}$$ is just a scaling factor. Now if plot the values of $$W(1,1200,t)$$ over the time domain:

ListPlot[ParallelTable[
Wγ[1, 1200, t], {t, 0, 100 10^(-15), 1 10^(-16)}],
Joined -> True] // AbsoluteTiming


with the same time increment $$\delta t= 10^{-16}$$. It takes about 1000 seconds to sample 1000 data points. Is there a way to improve this computation time ? Any suggestion for the improving the performance is greatly appreciated!

Question 2: I need to solve the equation $$W(\epsilon_m,\,\epsilon,\, t)= a$$ for $$\epsilon$$, where $$a$$ is constant and the values of $$t$$ and $$\epsilon_{m}$$ are fixed.

Edit 1: The scaling of $$a$$ is basically given as: $$a= W_{\gamma}(1,1200,t) r$$ where $$0 \leq r \leq 1$$ and it is a randomly generated number.

Edit 2: Value of $$t$$ lies within the integration region $$0 \leq t \leq 1\times 10^{-13}$$. It can take any value within these limits. I thank Alex for bringing up these points.

I thought I could use whenevent command in NDSolve if I can express $$W(\epsilon_{m},\epsilon,t)$$ in terms NDSolve . So I tried:

Wγ2[ϵm_?NumericQ, ϵ_?NumericQ, t_?NumericQ] :=
W[ϵ] /.
NDSolve[{W'[ϵγ] ==
819 10^(-16) dWγ[eγ, t],
W[ϵm] ==
0}, {W}, {ϵγ, ϵm , ϵ}][[1]]



When evaluated, I get the following errors

In[21]:= Wγ2[1, 100, 50 10^(-15)]

During evaluation of In[21]:= NDSolve::ndnum: Encountered non-numerical value for a derivative at ϵγ == 1..

During evaluation of In[21]:= ReplaceAll::reps: {W'[ϵγ]==(819 dWγ[eγ,1/20000000000000])/10000000000000000,W[1]==0} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

Out[21]= W[100] /. {W'[ϵγ] == (819
dWγ[eγ, 1/20000000000000])/10000000000000000,
W[1] == 0}


Sorry for the noob question but why am I getting no evaluation ? Is this because $$dW(\epsilon_{\gamma},t)$$ is failing to get numeric values when placed inside NDSolve ?, although it can be still evaluated outside:

In[7]:= dWγ[100, 50 10^(-15)]
Out[7]= 6.2098*10^25

• @user91411 It could be better to explain your problem in Latex as well. Nov 30, 2021 at 2:26
• @MichaelE2 thanks for clearing up the code. I cleaned up the code a bit more, so that it can be executed after copy pasting. Nov 30, 2021 at 12:45
• @AlexTrounev Hi Alex. Thanks for the suggestion. I have added more detail about what I am trying to do. Nov 30, 2021 at 12:45
• You could replace NIntegrate[AiryAi[y],...] with the result of Integrate[AiryAi[y], {y, z, Infinity}, Assumptions -> z > 0]. -- You will need to use arbitrary precision numbers to avoid underflow (but it won't speed things up, probably). See reference.wolfram.com/language/tutorial/Numbers.html and the WorkingPrecision option of NIntegrate, NDSolve, and other numerical solvers. Nov 30, 2021 at 14:07
• @MichaelE2 Hi Michael. I have tried your suggestion. It returns hypergeometric function with the argument $\chi_2^3$. Its evaluation takes much longer time than Nintegrate so I've decided to stick with Nintegrate option. Nov 30, 2021 at 14:11

To solve this problem we rescale variables and functions as follows $$t\rightarrow \omega t, z\rightarrow \omega z$$, then we have

\[Alpha] = 1/137; h = 105 10^(-36); Es = 13 10^(17); E0 =
27 10^(13); a0 = 70; \[Omega] = 22 10^(14); \[Tau] =
22 10^(-15); k = \[Omega] \[Tau]; py0 = 0; pz0 = 1200; z0 = -100 \
10^(-15) \[Omega]; t0 = 0; tmax = \[Omega] 10^(-13); sol =
NDSolve[{py'[
s] == -a0  Cos[(s + z[s])] Exp[-(s + z[s])^2/k^2] (1 +
pz[s]/Sqrt[1 + py[s]^2 + pz[s]^2]),
pz'[s] ==
a0  Cos[ (s + z[s])] Exp[-(s + z[s])^2/k^2] py[s]/
Sqrt[1 + py[s]^2 + pz[s]^2],
z'[s] == pz[s]/Sqrt[1 + py[s]^2 + pz[s]^2], py[t0] == py0,
pz[t0] == pz0, z[t0] == z0}, {py, pz, z}, {s, t0, tmax}]


Visualization of numerical solution

Table[Plot[
Evaluate[{py[t], pz[t], z[t]} /. sol[[1]]][[i]], {t, 0, tmax},
PlotRange -> All], {i, 3}]


Now we define functions $$\chi ,dW\gamma$$ using approximation for $$\int_t^\infty Ai(t) dt$$ in the form $$f(t)$$

f[t_] := If[t <= 8,
1/18 (6 - (
3 3^(5/6)
t Gamma[1/3] HypergeometricPFQ[{1/3}, {2/3, 4/3}, t^3/
9])/\[Pi] + (
3^(2/3) t^2 HypergeometricPFQ[{2/3}, {4/3, 5/3}, t^3/9])/
Gamma[4/3]), -((5 E^(-((2 t^(3/2))/3)))/(
72 Sqrt[\[Pi]] t^(3/4))) + (41 Erfc[Sqrt[2/3] t^(3/4)])/(
36 Sqrt[6])]; \[Chi] = {-\[Alpha]/(h (py[t]^2 +
pz[t]^2)), ((Es Exp[(t + z[t])^2/
k^2] \[Epsilon]\[Gamma])/(E0 Abs[
Cos[(t + z[t]) ]] (-\[Epsilon]\[Gamma] +
Sqrt[py[t]^2 + pz[t]^2]) (pz[t] +
Sqrt[py[t]^2 + pz[t]^2])))^(2/3), (\[Epsilon]\[Gamma] E0/
Es Exp[-(t + z[t])^2/k^2] Abs[
Cos[ (t + z[t])]] (1 + pz[t]/Sqrt[py[t]^2 + pz[t]^2]))};

dW\[Gamma] = \[Chi][[
1]] (f[\[Chi][[
2]]] + (2/\[Chi][[2]] + \[Chi][[
3]] (\[Chi][[2]])^(1/2)) AiryAiPrime[\[Chi][[2]]]);


Using these definitions we can compute in one second

lst = Table[{ts/\[Omega],
dW\[Gamma] /. {\[Epsilon]\[Gamma] -> 100, t -> ts} /.
sol[[1]]}, {ts, 0, 100 10^(-15) \[Omega],
1 10^(-16) \[Omega]}] // Quiet; // AbsoluteTiming


Also we can plot $$dW\gamma$$ on different intervals as

{Plot[dW\[Gamma] /. {\[Epsilon]\[Gamma] -> 100, t -> ts} /.
sol[[1]], {ts, 0, 100 10^(-15) \[Omega]}, PlotRange -> All],
Plot[dW\[Gamma] /. {\[Epsilon]\[Gamma] -> 100, t -> ts} /.
sol[[1]], {ts, 50, 150}],
Plot[dW\[Gamma] /. {\[Epsilon]\[Gamma] -> 100, t -> ts} /.
sol[[1]], {ts, 0, 30}],
Plot[dW\[Gamma] /. {\[Epsilon]\[Gamma] -> 100, t -> ts} /.
sol[[1]], {ts, 160, 220}]} // Quiet

Table[Plot[
dW\[Gamma] /. {\[Epsilon]\[Gamma] -> eps, t -> ts} /.
sol[[1]], {ts, 0, 100 10^(-15) \[Omega]}, PlotRange -> All,
PlotPoints -> 100], {eps, {1, 500, 1000}}] // Quiet


{Plot[dW\[Gamma] /. {\[Epsilon]\[Gamma] -> eps, t -> 110} /.
sol[[1]], {eps, 1, 1000}, PlotRange -> All,
AxesLabel -> {"\[Epsilon]\[Gamma]", "dW\[Gamma]"}],
LogLogPlot[
dW\[Gamma] /. {\[Epsilon]\[Gamma] -> eps, t -> 110} /.
sol[[1]], {eps, 1, 1000}, PlotRange -> All,
AxesLabel -> {"\[Epsilon]\[Gamma]", "dW\[Gamma]"}]}


Finally we plot  Wγ[1, 1000, t] (it takes about 179.6 s)

lst1 = Table[{ts/\[Omega],
819*10^(-16)*
NIntegrate[
dW\[Gamma] /. {\[Epsilon]\[Gamma] -> eps, t -> ts} /.
sol[[1]], {eps, 1, 1000}] // Quiet}, {ts, 0, tmax,
tmax/1000.}] // AbsoluteTiming
ListPlot[lst1[[2]], PlotRange -> All, AxesLabel -> {"t", "W"}]


To solve problem 2, $$W(\epsilon_m,\epsilon, t)=a$$, with time dependent scaling $$a=W[1,1200,t]RandomReal[]$$ or with scaling a=1.5*10^15 RandomReal[] we use same algorithm as follows

lst = Table[{{ts/\[Omega], epsm},
819*10^(-16)*
NIntegrate[
dW\[Gamma] /. {\[Epsilon]\[Gamma] -> eps, t -> ts} /.
sol[[1]], {eps, 1, epsm}] // Quiet}, {epsm, 1., 1201,
10}, {ts, 0, tmax, tmax/10}];
e = Interpolation[Flatten[lst, 1], InterpolationOrder -> 3];

Plot3D[e[x, y], {x, 0, 10^-13}, {y, 1, 1.09*10^3},
ColorFunction -> Hue, Mesh -> None, PlotPoints -> 50]


To compute solution $$\epsilon(t)$$ with a given r=RandomReal[] we use in the first case

Table[plot =
ContourPlot[
e[x/\[Omega], y] - e[x/\[Omega], 1200] RandomReal[] == 0, {x, 0,
tmax}, {y, 1, 1.09*10^3}, PlotRange -> All], {8}]


This solution looks like In a case of second scaling we have

rd = RandomReal[{0, 1}, 10]

Table[pl[i] =
ContourPlot[
e[x/\[Omega], y] - 1.5 10^15 rd[[i]] == 0, {x, 0, tmax}, {y, 1,
1.2*10^3}, PlotRange -> All, PlotLabel -> rd[[i]]], {i, 10}]


In this case function $$\epsilon(t)$$ looks like a smooth function, and we can recover function $$\epsilon(t,r)$$

We can retrieve numerical data from pl as follows

Do[et[i] = First@Cases[Normal@pl[i], Line[data_] :> data, -1];,{i,Length[rd]}];


Finally we use rd and et to interpolate function $$\epsilon (t,r)$$.

• Hi Alex. Thank you for the detailed answer. The problem is about photon emission by an electron placed in an external field. I am trying to form the backbone of the code that would simulate such reaction. The optimization of Airy integral is neat. As I understand it uses the asymptotic expansion of hypergeometric function as the argument grows large. I will take a closer look at your post. I would also appreciate if you could comment on the second question. Dec 9, 2021 at 14:46
• @user91411 Ah, sorry, I don't pay attention to second question since it looks like very simple. What is value a in your model? Dec 9, 2021 at 16:08
• I apologize, there was a mistake in the previous comment. The number 'a' is essentially a random number but its scaling is different. Please see the edit underneath the question 2. Dec 9, 2021 at 17:52
• @user91411 Your second problem is not defined well since t is also unknown. Is this also randomly generated number? Dec 10, 2021 at 2:24
• Please see the second edit. You can set $t$ to take any value within the domain of integration. Dec 10, 2021 at 8:56

This is a partial answer for now as an addendum to Alex's answer.

Integration can be further optimized if we use the single parameter, $$\chi_2$$ and express all the relevant quantities in terms of it. For this we solve $$\epsilon_{\gamma}$$ using the definition of $$\chi_2$$: $$\epsilon_{\gamma}=\frac{\bar{p} \chi^{3/2}_2\chi_{e}}{(1+\chi^{3/2} _{2}\chi_e)}, \,\, d\epsilon_{\gamma}=\frac{3\bar{p} \chi^{1/2}_2\chi_{e}}{2(1+\chi^{3/2} _{2}\chi_e)^2},\,\, \chi_e=\frac{|\vec{E}|}{E_s}\left(\bar{p}+\bar{p}_z\right)$$ Using the above relations once can also express $$\chi_3$$ as: $$\chi_3=\frac{\chi^{3/2}_2\chi^2_e}{1+\chi_2^{3/2}\chi_{e}}$$

Rewriting $$dW(\epsilon_{\gamma},\, t)$$ in terms of new variables: $$dW(\chi_2,\, t)= -\frac{\alpha}{\hbar\bar{p}^2}\left(\int^{\infty}_{\chi_{2}}\text{Ai}(y)\, dy + \left(\frac{2}{\chi_2} + \frac{\chi^{3/2}_2\chi^2_e}{1+\chi^{3/2}_{2}\chi_{e}}\sqrt{\chi_2}\right) Ai'(\chi_2) \right)$$ with:


dWγ[pe_?NumericQ, χ2_?NumericQ, χe_?NumericQ] :=
-6951 10^(28)/
pe^2 (NIntegrate[
AiryAi[y], {y, χ2,
Infinity}] + (2/χ2 + (χ2^(3/2) χe^2)/(1 + (χ2)^(3/2) χe)
χ2^(1/2)) AiryAiPrime[χ2])


The code for the trajectory is now modified to extract the values of $$\bar{p}, \chi_2$$ and $$\chi_{e}$$:

χ[ϵγ_?NumericQ, t_?NumericQ] := χ[ϵγ, t] =
With[{ Es := 13 10^(17),
E0 := 27 10^(13), a0 := 70,
ω := 22 10^(14), τ := 22 10^(-15), py0 := 0,
pz0 := 1200,
z0 := -100 10^(-15), t0 := 0}, {Sqrt[py[t]^2 + pz[t]^2],((
Es Exp[(t + z[t])^2/τ^2] ϵγ)/(
E0 Abs[Cos[(t + z[t]) ω]] (-ϵγ + Sqrt[
py[t]^2 + pz[t]^2]) (pz[t] + Sqrt[py[t]^2 + pz[t]^2])))^(
2/3), ( E0/Es Exp[-(t + z[t])^2/τ^2] Abs[
Cos[ω (t + z[t])]]  (pz[t] + Sqrt[py[t]^2 + pz[t]^2])) } /. NDSolve[{
py'[s] == -a0 ω  Cos[ω (s +
z[s])] Exp[- (s + z[s])^2/τ^2] (1 +
pz[s]/Sqrt[1 + py[s]^2 + pz[s]^2]),
pz'[s] ==
a0  ω  Cos[ω (s +
z[s])] Exp[-(s + z[s])^2/τ^2] py[s]/
Sqrt[1 + py[s]^2 + pz[s]^2],
z'[s] == pz[s]/Sqrt[1 + py[s]^2 + pz[s]^2], py[t0] == py0,
pz[t0] == pz0, z[t0] == z0}, {py, pz, z}, {s, t0,  t}][[1]]]


Rewriting the integral over energy $$\epsilon_{\gamma}$$ as: $$W(\epsilon_m,\,\epsilon, \, t)=\int_{\epsilon_{m}}^{\epsilon} dW(\epsilon_{\gamma},\, t)d \epsilon_{\gamma} \\ =W(\chi^{m}_{2},\,\chi_{2}, \, t) =\int_{\chi^{m}_{2}}^{{\chi_{2}}} \frac{3\bar{p} (\chi^{'}_{2})^{1/2}\chi_{e}}{2(1+ (\chi^{'}_{2})^{3/2}\chi_e)^2} dW(\chi'_2,\, t) d\chi_2'$$ with

Wγ[
pe_?NumericQ, χ2i_?NumericQ, χ2f_?NumericQ, χe_?
NumericQ] :=
819*10^(-16) (NIntegrate[(3 pe Sqrt[χ2] χe)/(
2 (1 + χ2^(3/2) χe)^2)
dWγ[
pe, χ2, χe], {χ2, χ2i, χ2f}])


Doing so, samples 1000 points around 130 seconds

ParallelTable[{N[t 10^(15)],
Wγ[χ[1, t][[1]], χ[1, t][[2]], χ[1200 - 1/100, t][[2]], χ[1, t][[3]]]}, {t, 0,
100 10^(-15), 1 10^(-16)}] // Quiet // AbsoluteTiming


If I include Alex's definition for the Airy integral:

f[t_] := If[t <= 8,
1/18 (6 - (
3 3^(5/6)
t Gamma[1/3] HypergeometricPFQ[{1/3}, {2/3, 4/3}, t^3/
9])/\[Pi] + (
3^(2/3) t^2 HypergeometricPFQ[{2/3}, {4/3, 5/3}, t^3/9])/
Gamma[4/3]), -((5 E^(-((2 t^(3/2))/3)))/(
72 Sqrt[\[Pi]] t^(3/4))) + (41 Erfc[Sqrt[2/3] t^(3/4)])/(
36 Sqrt[6])];

dWγ2[pe_?NumericQ, χ2_?NumericQ, χe_?NumericQ] :=
-6951 10^(28)/
pe^2 (f[χ2] + (2/χ2 + (χ2^(3/2) χe^2)/(1 + (χ2)^(3/2) χe))AiryAiPrime[χ2])

Wγ2[
pe_?NumericQ, χ2i_?NumericQ, χ2f_?NumericQ, χe_?
NumericQ] :=
819*10^(-16) (NIntegrate[(3 pe Sqrt[χ2] χe)/(
2 (1 + χ2^(3/2) χe)^2)dWγ2[pe, χ2, χe], {χ2, χ2i, χ2f}])


1000 points takes about 17 seconds:

ParallelTable[{N[t 10^(15)],
Wγ2[χ[1, t][[1]], χ[1, t][[2]], χ[1200 - 1/100, t][[2]], χ[1, t][[3]]]}, {t, 0,
100 10^(-15), 1 10^(-16)}] // Quiet // AbsoluteTiming
$$$$