# 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?

• do the $p$'s depend on $z$ as well, or only on $t$? – themaker Jun 13 '17 at 17:54
• thanks for pointing that out. yes the p's are a function of t! I'll edit the post to clarify. (Also, @themaker I see this question hasn't gotten much attention; maybe there are some other things I can more clear so that it can get more attention?) – Steven Sagona Jun 16 '17 at 0:12
• but the question is about $z$. I ask because, if $E_p$ depends on $z$, and you have equations like $\dot p_{33} = E_p p_{13} + ...$ then the $p$'s need also depend on $z$. Also would be good to know what kind of boundary conditions are you considering for the PDE. – themaker Jun 16 '17 at 0:20
• sorry that was a typo (I've been working pretty late!) I meant to say that the p's are a function of z – Steven Sagona Jun 16 '17 at 2:50
• I've updated my question to try to make it more clear. @themaker please let me know if this answers your questions – Steven Sagona Jun 16 '17 at 3:44

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.
• @StevenSagona The CFL condition imposes $\Delta t < \Delta x/c$, and the stability of a RK scheme imposes something like $\Delta t < \alpha/\varrho$, where $\varrho$ is the spectral radius of the relaxation function's Jacobian matrix $\partial \mathbf{R}/\partial \mathbf{U}$. – Harry49 Jun 21 '17 at 4:33
• So these conditions you mention are for guaranteeing the individual solutions for each individual operator solutions for $H_a$ and $H_b$? But is there a condition necessary for Strang splitting? That is, what is necessary about this PDE-ODE system for it to be able to be split up into a form like $H_b$ $H_a$ $H_b$ in which $H_a$ and $H_b$ can be individually solved? – Steven Sagona Jun 23 '17 at 21:01