2
$\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 '21 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 '21 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 '21 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 '21 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 '21 at 14:11
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 '21 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 '21 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 '21 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 '21 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 '21 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, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.