Numerical Solution of the Advection Dispersion equation

I am facing a simple (at first glance) problem. I need to implement a numerical scheme for the solution of the first order wave propagation equation with chromatic dispersion included. My original problem is (for a forward propagating wave):

$$\frac{1}{c} \frac{\partial u(x,t)}{\partial t} = -\frac{ \partial u}{ \partial x} - \frac{i \beta_2}{2} \frac{ \partial^2 u}{ \partial t^2},$$ where $c$ is the velocity of light, $u$ is the (complex) envelope of the field, $\beta_2$ is the 2nd order dispersion coefficient. Assume also that the wave is propagating inside a ring cavity of length , say, L where I take periodic boundary conditions: $u(x+L,t) = u(x,t)$ and also that at $t=0$ we know $u(x,0)$ and $u_t(x,0)$.

I am trying to implement a time-stepping numerical scheme and in the process I tried the following:

1) MOL approach, where I do semidiscretization along $x$, reduce the set of equations to a system of first order ODEs (by setting $v = \dot{u}$) and I establish a system: $$\begin{bmatrix} \dot{v} \\ \dot{u} \end{bmatrix} = A \begin{bmatrix} v \\ u \end{bmatrix} .$$

When I solve the corresponding ODEs via 4th Runge-Kutta, Crank-Nicholson , or simply precomputing the matrix exponential, unfortunatelly, all my solutions eventually blow up to Inf. I implemented the periodic boundary conditions by modifying the matrix $A$ as $A \leftarrow PA$ , where $P$ is the identity matrix with the first row identical copy of the last row.

I also tried a simple finite differences approach where the spatial derivative is approximated via an upwind FD (first order) but to no avail.

Lastly I tried a strang splitting approach based on the two equations: \begin{align} \frac{1}{c}\dot{u} &= -\frac{\delta u }{\delta x} \\ \frac{1}{c}\dot{u} &= -\frac{i\beta_2}{2}\frac{\delta^2 u }{\delta t^2} , \end{align}

where this time the solution does not blow up but it looks somehow unphysical.

Does someone here know a stable and possibly higher than first order time-stepping scheme for this equation? Please, note that a solution based on a Fourier transform in $x$ is also not a good option for me because I would like to have the flexibility to implement different non-periodic boundary conditions. I would also dislike substituting the second order derivative in time with one in space due to the fact that this complicates the implementation of boundary conditions.

Thanks.

• Is it supposed to be $\partial^2/\partial t^2$ or $\partial^2/\partial x^2$? – Kirill Jul 25 '15 at 21:58
• The second derivative is in time. – kenny Jul 26 '15 at 7:46
• Are you familiar with von Neumann stability analysis? That would give you some insight under which conditions your discretization is stable – nluigi Dec 9 '15 at 11:19
• this question has been also asked on physics SE before any answers were posted here. – Anton Menshov Jun 5 '19 at 5:05

If $\beta_2$ is piece-wise constant the equation is solvable by Fesnel integrals so that you can develop approximte solution based on the boundary conditions. You have to introduce a co-moving frame and solve from there.