# Runge Kutta solution blows up for a first order ODE with very large coefficients

I am solving a first-order ODE:

$\frac{\partial \rho }{\partial t} = -a \rho^2 + b |A(t)|^2 \rho +c|A(t)|^{2m}$

This is the evolution of the plasma density in the presence of a laser pulse (complex amplitude A, which is normalized). The coefficients a, b and c depend on the material and the initial laser intensity. Typically, these coefficients are very large (for example a = 1e13, b=2e16 c=1e7).

Now I solve this equation (in Python, and I also played with C) with the conventional Runge-Kutta 4th order method (also tried 6th order); here is a little bit of code (quite standard, nothing fancy):

for l in range(Nr):
mag= np.abs(Uo[:,l] )**2

for k in range(1,Nt):
k1 = dT*(b*rho[k-1,l]*mag[k-1] + c*mag[k-1]**m - a*rho[k-1,l **2.0 )

tmp1=0.5  * ( mag[k]+mag[k-1])
tmp2=rho[k-1,l] + 0.5*k1
k2 = dT * (b*tmp2 * tmp1 + c*tmp1**m- a* tmp2**2.0  )

tmp2=rho[k-1,l] + 0.5*k2
k3 = dT * (b*tmp2 * tmp1 + c*tmp1**m- a* tmp2**2.0 )

tmp2=rho[k-1,l] + k3
k4 = dT * (b*tmp2 * mag[k] + c*mag[k]**m- a* tmp2**2.0 )

rho[k,l]=   rho[k-1,l] + (k1/6. + k2/3. + k3/3. + k4/6.)


I also loop over space (Nr) because I have a time and space dependent plasma density. Since I discretized my time axis (with a total of Nt points, dT is usually 1e-16s), I do not use adaptive step size control (besides, I want this code to be fast).

Now, if I do this integration with low power, everything works fine, and I get a nice result (image showing the temporal evolution of plasma density $\rho$ for low power, at one space point):

But if I increase the pulse power, my Runge Kutta solution blows up, and I get NaNs. Reason is, that for high powers (say 1e9 Watts), the coefficients become huge (for example a = 1e13, b=2e52, c=1e13).

My question is: Do you have an idea how I can still get a result, without using step size control or using a finer grid? I wouldn't mind using approximations. I am really thankful for any help:)

• Though the following will not necessarily solve your problem, it is highly recommended to implement the ode dimensionless such that all relevant coefficients are around 1 to stabilize solvers. The result is then scaled back accordingly.
– Bort
Oct 26, 2016 at 9:41
• Hi Bort, thanks. Yes, you are right, of course. I also implemented it to be dimensionless. It helped in a way, and I can go up to 1e6W until it blows up. But still not the high power levels I need. Oct 26, 2016 at 9:48
• The reason for the blow up might be some form of stability limit. Did you try using an implicit method? Since your ODE is scalar, it should be relatively straightforward to implement a Newton iteration to solve the arising nonlinear problem. Jun 18, 2017 at 12:27

Though the following will not necessarily solve your problem, it is highly recommended to implement the ode dimensionless such that all relevant quantities are around 1 to stabilize solvers. The result is to be scaled back accordingly. This will also give you a hint what goes wrong in terms of the size of the coefficients or step size in the discretization.

So in your case, try: $\rho = \hat{\rho} \rho_0, t = \hat{t} t_0$

This gives you:

$\frac{\partial \hat{\rho}}{\partial \hat{t}}=-\hat{a}\hat{\rho}^2+\hat{b}|A(t_0\hat{t})|^2 \hat{\rho}+\hat{c}|A(t_0\hat{t})|^{2m}$

with $\hat{a}=at_0\rho_0,\hat{b}=b t_0, \hat{c}=c \frac{t_0}{\rho_0}$

This way, when your expected plasma density increases because of larger laser powers, you pick larger $\rho_0$ and scale down $\hat{a}, \hat{b}, \hat{c}$.

E.g. You can always pick $\hat{b}=1\Rightarrow t_0=\frac{1}{b}, \hat{a}=1 \Rightarrow \rho_0=\frac{b}{a}, \hat{c}=\frac{c a}{b^2}$.

So for your given values (a = 1e13, b=2e52, c=1e13), you will probably need very small time steps, but all $\hat{a}=\hat{b}=1,\hat{c}=25\cdot 10^{-80}$, in the scaled version. So this should not blow up and the c term can be neglected.

In case of the limit $c\rightarrow 0$, you can even find a analytical solution for verification of your solver.

$\hat{\rho}(\hat{t})=\frac{e^{\int_1^\hat{t} |A(t_0 \xi)|d\xi}}{C_1 -\int_1^\hat{t}-e^{\int_1^\zeta |A(t_0 \xi)|d\xi}d\zeta}$

Here $C_1$ is the integration constant.