Coupled nonlinear PDEs with time dependence on the RHS

Question

I would like to numerically solve the following system of 2 coupled partial differential equations for the unknown functions $\psi_X(x,y,t)$ and $\psi_C(x,y,t)$:

$\partial_t \psi_X = -i\psi_C - i\vert\psi_X\vert^2 \psi_X$

$\partial_t \psi_C = -i F(x,y)\exp{[i(k^p_x x+k^p_y y - \omega_p t)]} -i\psi_X -i [\delta - \Delta -i\kappa_C]\psi_C$

with $F(x,y) > 0$ (typically a gaussian function), $\kappa_C >0$ and $\delta \in \mathbb{R}$ constants.

The domain of integration is a square box of length 2L, with PBC.

So far I have tried a 5th order Runge Kutta approach, with adaptive time step, and calculating the laplacian term in fourier space.

Could anyone suggest a more robust/faster alternative?

Answer 1

How are you choosing your time step? If your Runge-Kutta method is explicit, you will need to take $\Delta t = {\mathcal O}((\Delta x)^2)$, which is the real reason that this method is slow.

In order to march faster in time, I would recommend using an operator-splitting approach, in which you alternate between integrating the ODE terms and integrating the spatial derivatives:

Step 1: Integrate

$$ \partial_t \psi_X = -i\psi_C - i\vert\psi_X\vert^2 \psi_X \\ \partial_t \psi_C = -i F(x,y) -i\psi_X -i [\delta -i\kappa_C]\psi_C $$ over a time step $\Delta t$ using an explicit Runge-Kutta method.

Step 2: Integrate $$\partial_t \psi_C = i\Delta \psi_C$$ exactly (over a time step $\Delta t$) by taking the FFT, multiplying by the appropriate phase factor, and taking the inverse FFT.

Since the time integration in step 2 is exact, the time step size will be restricted only based on step 1 (and that part doesn't look stiff). The overall accuracy will be first order. You can apply higher order splitting methods if necessary (look up Strang splitting for second order).

In order to impose the correct boundary conditions, you should use an implicit time integration method for the spatial part. The implicit trapezoidal or midpoint method should work well, and will again allow you to take arbitrarily large time steps. But you'll need to use an efficient linear solver.

A detailed example (including Python code) of how to implement such methods for a different problem can be found in this IPython notebook (written by me).

A very basic but useful introduction to low-order splitting methods and other alternatives can be found, for instance, in Randall LeVeque's finite difference book.