Coupled nonlinear PDEs with time dependence on the RHS

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?

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.

• Indeed, the RK I use is explicit (5th order adaptive RK from Numerical Recipes). The BCs are not so crucial, in the sense that, if it makes a big difference in terms of computation time, I can stick with periodic ones. Could you please explain your second paragraph in more detail though? Im quite new to this field :) – Andrei May 30 '14 at 11:22
• I added more detail about operator splitting. Have you had a course in numerical methods? If not, take one. – David Ketcheson May 30 '14 at 18:18
• Thanks a lot :) One last thing, in the splitting approach that you suggest, can one use an adaptive time-step or does $\Delta t$ have to be fixed? – Andrei May 30 '14 at 22:23
• Also, could you please suggest a good book/course material? The information I found on the internet about Strang splitting is quite lacking. – Andrei May 31 '14 at 14:28
• Thanks a lot for all the extra info! Unfortunately I completely forgot I also have a time-dependence on the rhs (see my edit above) so I'm not sure how applicable your answer is right now. – Andrei Jun 4 '14 at 10:54