3
$\begingroup$

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
$\endgroup$
9
  • $\begingroup$ @user91411 It could be better to explain your problem in Latex as well. $\endgroup$ Nov 30, 2021 at 2:26
  • $\begingroup$ @MichaelE2 thanks for clearing up the code. I cleaned up the code a bit more, so that it can be executed after copy pasting. $\endgroup$
    – user91411
    Nov 30, 2021 at 12:45
  • $\begingroup$ @AlexTrounev Hi Alex. Thanks for the suggestion. I have added more detail about what I am trying to do. $\endgroup$
    – user91411
    Nov 30, 2021 at 12:45
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Nov 30, 2021 at 14:07
  • $\begingroup$ @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. $\endgroup$
    – user91411
    Nov 30, 2021 at 14:11

2 Answers 2

2
$\begingroup$

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}]

Figure 1

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

Figure 2

{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]"}]}

Figure 3

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"}]

Figure 4

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]

Figure 5 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 Figure 6 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)$ Figure 7

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)$.

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – user91411
    Dec 9, 2021 at 14:46
  • $\begingroup$ @user91411 Ah, sorry, I don't pay attention to second question since it looks like very simple. What is value a in your model? $\endgroup$ Dec 9, 2021 at 16:08
  • $\begingroup$ 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. $\endgroup$
    – user91411
    Dec 9, 2021 at 17:52
  • $\begingroup$ @user91411 Your second problem is not defined well since t is also unknown. Is this also randomly generated number? $\endgroup$ Dec 10, 2021 at 2:24
  • $\begingroup$ Please see the second edit. You can set $t$ to take any value within the domain of integration. $\endgroup$
    – user91411
    Dec 10, 2021 at 8:56
0
$\begingroup$

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
```      
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.