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

Question

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:)

Answer 1

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.