# Numerical Solution of the Propagation-Dispersion equation

I have asked this question on Computational Science and also on Mathoverflow, but no satisfactory answers so far. I thought maybe the physics community could shed some insight on the issue.

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):

\begin{equation} \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}, \end{equation} 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{equation} \begin{bmatrix} \dot{v} \\ \dot{u} \end{bmatrix} = A \begin{bmatrix} v \\ u \end{bmatrix} . \end{equation} 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.

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

3. 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 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 am starting to think that the problem is in general not well posed, since more or less all my numerical solutions so far agree in the way they diverge to infinity. Is there something wrong with the boundary conditions?

Such an equation is usually found in the context of ligth propagation inside a transmission line (optical fiber), where people usually solve the non-linear Schroedinger equation via a split-step Fourier scheme.

**********EDIT********

As Kyle suggested the above equation somehow resembles the Dampened Harmonic oscillator equation. Rewrite in the following manner: \begin{equation} M\ddot{U} + C\dot{U} + DU = 0, \end{equation} where $M = \frac{i \beta_2}{2}$, $C = \frac{1}{c}$ and $D_{ij} = (\delta_{i,j}-\delta_{i,j-1})/\Delta x$ - for upwind differentiation. Applying the Velocity-Verlet method (setting $A(t) = \ddot{U}(t)$ and $V(t) = \dot{U}(t)$ ) yields: \begin{align} U(t+\Delta t) &= U(t) + \Delta{t}V(t) +\frac{{\Delta t }^2}{2}A(t) \\ V(t+\Delta t) &= V(t) + \frac{\Delta t}{2}\left ( A(t) + A(t+\Delta t) \right) \\ &= V(t) + \frac{\Delta t}{2} A(t) - \frac{\Delta t}{2} M^{-1} \left( CV(t+\Delta t) + DU(t+\Delta t) \right ) \end{align}

This yields an implicit equation for $V(t+\Delta t)$. Solving it yields: \begin{align} V(t+\Delta t) &= \left [1+\frac{\Delta t}{2}M^{-1}C \right ]^{-1} \left( V(t) + \frac{\Delta t}{2} A(t) - \frac{\Delta t}{2} M^{-1} DU(t+\Delta t) \right) . \end{align} Which looks reasonable. In order to finalize the scheme, a boundary condition is needed. Assuming (for simplicity's sake might be another kind of BC, say Dirichlet) periodic boundaries, i.e. $U_1 = U_N$, this imposes the following modification of the differentiation matrix $D_{ij}$ (when the wave is propagating in positive x direction): \begin{align} D_{1,j} &= (\delta_{1,j} - \delta_{N-1,j})/\Delta(x) \\ \end{align}

The solutions I obtained look somewhat more reliable but still unstable. Especially at the boundaries. Take a look at the attached matlab code:

clear;clc; close all;
%Velocity, domain boundaries, dispersion coeff, grid size, grid, grid spacing and time step
c = +2; a = -10; b = 10;
beta = +3; N = 250;  x = linspace(a,b,N)';
dx = x(2)-x(1); dt = 0.01*dx/abs(c);

% %differentiation matrix.
Dx = zeros(N,N); Dx(1,1) = 1/dx; Dx(1,N-1) = -1/dx;
for i = 2:N Dx(i,i) = 1/dx; Dx(i,i-1) = -1/dx; end;

% % Different initial conditions
z0 = (b+a)/20;

% % 1) A sech pulse.
tau =1; %FWHM duration of the pulse
k = 5; %arbitrary phase shift...
aE = @(z,t) sech((t-(z-z0)/c)/tau).*exp(1i*k*z);
aE_t =  @(z,t) -1/tau.*sech((t-(z-z0)/c)/tau).*tanh((t-(z-z0)/c)/tau).*exp(1i*k*z);

% % 2) A gaussian
% tau =1; %sec
% k = 5; %arbitrary phase shift...
% aE = @(z,t)exp(-((t-(z-z0)/c)/tau).^2).*exp(1i*k*z);
% aE_t =  @(z,t) 2*(t-(z-z0)/c)/tau.*exp(-((t-(z-z0)/c)/tau).^2).*exp(1i*k*z);

% % 3) A travelling sine wave!
% m = 5; %order of the oscillations
% kz = 2*pi/(b-a)*m % wave number!
% w = kz*c;
% aE = @(z,t)sin(w*t-kz*z);
% aE_t =  @(z,t) w*cos(w*t-kz*z);

U = aE(x,0);  V = aE_t(x,0);
U(1) = U(end); V(1) = V(end);

%mass;  damping; "elasticity" matrix!
M = 1i*beta/2; C = (1/c); D = Dx;
ctr = 1; t = dt;

tend = 2*(2*(b-a)/abs(c));
Unew = U; Vnew = V;
while t<tend

if(mod(ctr-1,100)==0)
plot(x,real(U),x,real(V),x,real(aE(x,t)),x,real(aE_t(x,t)));
%         plot(x,real(aE(x,t)),x,real(aE_t(x,t))); %<- plot the analytical expressions!
legend('U^n','V^n','U(t)','V(t)');
title(['t = ' num2str(t)]);
xlim([a-1 b+1])
ylim([-10 10])
getframe;
end

%acceleration, U and velocity !
A = -(C*V+D*U)/M;
Unew = U + V*dt+0.5*(dt^2)*A;
Vnew = (V+dt/2*A-dt/2*D*Unew/M)/(1+dt/2*C/M);
U = Unew; V = Vnew;  ctr =ctr+1;     t=t+dt;
end

• Are you always evaluating "left to right"? With periodic boundary conditions I can see that could give problems. What would happen if you stepped in alternate (space) directions? One iteration clockwise, the next anti clockwise. It might make the periodic BC better behaved. – Floris Jul 28 '15 at 11:55
• I am always evaluating "left to right", yes. Sorry but I am not sure I understand what you mean. The equation above is supposed to describe a wave packet propagating in a single direction (depending on the sign of c). What do you mean to make one iteration counter clockwise? – kenny Jul 28 '15 at 12:42
• Have you tried Verlet integration? – Kyle Kanos Jul 28 '15 at 12:47
• @Floris I am not evaluating the time evolution iteratively as you suggested, but rather simultaneously (that is I store everything in a mesh function and do the relevant operations on that mesh function.) So forward or backward iteration do not play a role here. The only place this is cruicial is the approximation of the difference operator d/dx . This one I tried with an upwind scheme or with a global Chebyshev differentiation matrix, but again (as mentioned already) without success. – kenny Jul 28 '15 at 15:11
• This is not the equation of a wave. The spatial derivative must be of second order. If your problem if one-dimensional, you probably should consider the telegrapher's equations. – Tom-Tom Sep 21 '15 at 13:09