I am trying to solve an equation of the following type $$\partial_zE(z)=-c_0J$$ with $$J=c_1\beta E^3(z)$$ using the boost::odeint-framework and a fixed time stepper, with $c_0$, $c_1$ and $\beta$ constant values.
The right hand side of the function is calculated by the following procedure (depending on the step size $\Delta z$):

  • Multiplication of $E$ with a vector $\exp(i\cdot k_z\cdot\Delta z)$ to move the field from $z=0$ to $z=\Delta z$
  • Calculation of $J$, with the equation given above
  • Multiplication of $J$ with a vector $\exp(-i\cdot k_z\cdot\Delta z)$ to move the result back from $z=\Delta z$ to $z=0$
  • Return the value to the solver

This means that for a common RK4-approach I have to call that function four times, once for $z=0$, twice for $z=0.5\Delta z$ and once for $z=\Delta z$. This approach works fine for $\max(E)\ll E_{lim}$, but for $\max(E)\approx E_{lim}$ I have to reduce my step size $\Delta z$ significantly, else my resulting values for $E$ will explode within a few steps. $E_{lim}$ can be estimated by using $J_{lim}\approx 10^{13}$ (which might depend on the values for $c_0$ and $c_1$). If $J$ stays below that threshold, even large steps are not an issue, but the closer I get to that limit, the smaller I have to choose my steps (which in turn increase my calculation time).
Why does my equation behave like that? When calculating the solution for a similar type equation $$f'(t)=-h\cdot f^3$$ on WolframAlpha I get $$f(t)=\pm\frac{1}{\sqrt{c_0+2ht}}$$

To address the comment from @MaximUmansky:
I am trying to solve a reduced version of the UPPE, written as $$\partial_z\hat{E}(\omega, z)=ik_z\hat{E}(\omega, z)-\frac{\omega}{2\varepsilon_0 c^2k_z}\hat{J}(\omega, z)$$ with $J$ given by the non-linear $K$-photon-absorption ($K=2$): $$\begin{split}J&=\frac{2\varepsilon_0n_0c\sigma_2I^2\hbar\omega_0\varrho_{nt}}{I}E\\ &=\varepsilon_0^2n_0^2c^2\sigma_2\hbar\omega_0E^3\end{split}$$ with $\sigma_2$ the absorption cross section which can be reduced to the equation given above: $$J=c_1\beta E^3(z)$$ For solving the equation I can split it into the linear and the non-linear part, and start by solving the non-linear part using a RK4-solver, and applying the linear part afterwards. Implementation of the right-hand side function is similar as described above, except by the addition of an FFT for transferring $\hat{E}\rightarrow E$ and back, before and after the calculation of $J$.

Still, as soon as $J$ reaches a value $J\geq J_{lim}$ (which occurs for either large values of $E$ or larger values for $\sigma_2$, i.e. a large absorption), the problems described above occur, and the closer $J$ gets to $J_{lim}$, the smaller I have to choose my step size, to prevent $J$ becoming larger than $J_{lim}$, which is a behavior I still do not understand.

  • 1
    $\begingroup$ Is this from electrodynamics in continuous media? It would help to know the physical meaning of those equations, that may help understand the problem you see in the calculations. $\endgroup$ Jul 16 '20 at 15:42
  • $\begingroup$ @MaximUmansky: I extended my question, adding additional information. If those are not sufficient, please let me know. $\endgroup$
    – arc_lupus
    Jul 17 '20 at 9:05
  • $\begingroup$ Have you tried to employ a more stable method than RK4? Most of the stable method are implicit, but there exist explicit alternative. You can have a look at Runge-Kutta-Chebyshev (RKC) methods, or ROCK. $\endgroup$
    – G. Gare
    Jul 17 '20 at 9:37
  • $\begingroup$ @G.Gare: I'd like to avoid implicit methods, usually my matrix for $E$ has $>1e6$ entries, but will take a look at RKC and ROCK methods. $\endgroup$
    – arc_lupus
    Jul 17 '20 at 10:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.