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