Numerical Solution to system of ODEs coupled to a PDE

In a nutshell: I want to numerically (and maybe efficiently) solve a system of ODEs that are coupled to a PDE.

I want to do a simulation of a pulse of light entering a collection of atoms which can be coupled by two sets of differential equations, the Maxwell-Schrodinger equations with the optical Bloch equations.

The Optical-Bloch equations are a system of complex differential equations (describing how the energy levels of the atoms change as they are hit by laser light). The Maxwell-Schrodinger equations are simply the wave-equation with the approximation that the envelope varies slowly.

Here is what one such system of ODEs looks like: \begin{align*} \dot p_{11} &= E_p p_{13} - E_p p_{31} \\ \dot p_{12} &= E_p p_{13} - E_p p_{32} - p_{12} (\Delta_{31}-\Delta_{32}) \\ \dot p_{13} &= E_p p_{11} + E_c p_{12} - E_p p_{33} + p_{13} \Delta_{31} \\ \dot p_{21} &= E_p p_{23} - E_c p_{31} + p_{21}(-\Delta_{31} + \Delta_{32}) \\ \dot p_{22} &= E_c p_{23} - E_c p_{32} \\ \dot p_{23} &= E_p p_{21} + E_c p_{22} - E_c p_{33} + p_{23} \Delta_{31} + p_{23} (-\Delta_{31} + \Delta_{32}) \\ \dot p_{31} &= - E_p p_{11} - E_c p_{21} + E_p p_{33} - p_{31} \Delta_{31} \\ \dot p_{32} &= - E_p p_{21} - E_c p_{22} - E_c p_{33} + p_{23} \Delta_{31} + p_{23} (-\Delta_{31} + \Delta_{32}) \\ \dot p_{33} &= E_p p_{13} - E_c p_{23} + E_p p_{31} + E_c p_{32} \end{align*}

where $p_{i j}(t, z)$ represents the populations of atoms at different energy levels. The initial condition assumes $p_{11}(0, z) = 1$ and all other $p_{ij}(0, z)=0$.

$\Delta_{ij}$ (for all ij) and $E_c$ are constants. (i is the imaginary number, c is the speed of light, and k is a constant)

$E_p(t, z)$ is a function of t and z, and evolves according to the PDE called the Maxwell-Schrodinger equation (it represents the laser light that is cycling the population of the atom's energy levels): $$ \frac{\partial E_p}{\partial z } + \frac{1}{c} \frac{\partial E_p}{\partial t } = i k p_{13}(t, z) $$

We apply the boundary condition $E(0, t) = e^{-(t-t_0)^2} $, which implies that our wave starts as an unobstructed gaussian pulse that then enters the medium over time.

Typically these solutions are solved by analytically solving them at steadystate (taking their time derivatives equal to zero and simply plugging them in).

I'm interested in understanding what different methods can be used to solve such a dynamical system. (and when such methods break down).

The standard method for solving the Optical Bloch Equations (the system of linear equations) is simply putting these equations in matrix form $\dot P = MP$ and then recognizing that the solution simply has the form $P_0 e^{M*t}$.

With the inclusion of the slowly-varying Maxwell equation (the PDE), it's not obvious to me that there is a simple solution like before. In fact, things seem to be much more complicated, because if one tries to make the PDE into a set of ODEs, then we identify that the previous ODEs now are nonlinear ODEs since there are terms proportional to $E_p p_{13}$.

Are there any obvious tricks that can be used here for finding a numerical solution? Is it necessary to rewrite the PDE as an ODE and try to numerically solve the nonlinear system of ODEs?

Solution 1:

The usual trick for handling such coupled PDE-ODE systems with different time scales is splitting. If we consider the vector of ten unknowns $\mathbf{U} = (E_p,p_{11}\dots p_{33})^\top$, then the system of differential equations can be written $$ \frac{\partial}{\partial t} \mathbf{U} + \mathbf{A} \frac{\partial}{\partial z} \mathbf{U} = \mathbf{R}(\mathbf{U}) \, , $$ where $\mathbf{A}$ is a $10\times 10$ matrix with only nonzero coefficient $\mathbf{A}_{11} = c$, and \begin{aligned} \mathbf{R}(\mathbf{U}) &= \left(i k c p_{13},\dot{p}_{11}\dots \dot{p}_{33}\right)^\top\\ &= \left(i k c p_{13},E_p(p_{13} - p_{31})\dots E_p(p_{13} + p_{31}) + E_c (p_{32}-p_{23})\right)^\top . \end{aligned} In the present case, the two involved physical time scales relate to the wave speed $c$, and to the relaxation spectrum deduced from $\mathbf{R}(\mathbf{U})$. Splitting consists in solving successively $$ \frac{\partial}{\partial t} \mathbf{U} + \mathbf{A} \frac{\partial}{\partial z} \mathbf{U} = \mathbf{0}\, , \qquad (\mathbf{H}_a) $$ and $$ \frac{\partial}{\partial t} \mathbf{U} = \mathbf{R}(\mathbf{U})\, , \qquad\qquad (\mathbf{H}_b) $$ with a dedicated method for each operator. For instance, the Strang splitting scheme writes $$ \mathbf{U}_i^{n+1} = \mathbf{H}_b(\Delta t/2)\circ \mathbf{H}_a(\Delta t) \circ \mathbf{H}_b(\Delta t/2)\, \mathbf{U}_i^n \, , $$ where $\mathbf{U}_i^n$ approximates $\mathbf{U}(z_i,t_n)$ and $\Delta t$ is the time step. This scheme is second-order accurate if both suboperators are stable and second-order accurate. Possible dedicated methods for each operator are

• a second-order accurate solver for the linear transport equations $\mathbf{H}_a$, e.g. a Lax-Wendroff or ADER scheme. Such a scheme is stable under the CFL condition.

• a second-order accurate solver for the nonlinear ODE $\mathbf{H}_b$, e.g. a Runge-Kutta method. If $\mathbf{H}_b$ is stiff, then the stability condition for such a solver can be penalizing.