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
NIntegrate[AiryAi[y],...]with the result ofIntegrate[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 theWorkingPrecisionoption ofNIntegrate,NDSolve, and other numerical solvers. $\endgroup$