Suppose we have the negative, inhomogeneous advection equation: $$\left(\frac{\partial}{\partial x}-\frac{1}{c}\frac{\partial}{\partial t}\right)v(t,x)=u(t,x)\qquad(t\in\mathbb{R}_{+},x\in\mathbb{R})$$ I want to solve it numerically, backward in time and forward in space.

Now, I have tried a few schemes without avail. For example, the Crank-Nicolson scheme going backwards in time:

In particular, we may discretise as

$$\frac{1}{c}\frac{\mathrm{d}}{\mathrm{d}t}v_\ell(t)=\underbrace{\frac{v_{\ell+1}(t)-v_\ell(t)}{\Delta x}-u_\ell(t)}_{G_\ell(t)}+\mathcal{O}(\Delta x)\qquad\Delta x\to 0$$ where we have taken the forward difference of the spatial derivative.

Then applying the Crank-Nicolson scheme yields

$$\frac{1}{c}\frac{v^{n+1}_\ell-v^n_\ell}{\Delta t}=\frac{1}{2}[G_\ell(t^{n+1})+G_\ell(t^n)]=\frac{1}{2}\left[\frac{v^{n+1}_{\ell+1}-v^{n+1}_\ell+v^n_{\ell+1}-v^n_\ell}{\Delta x}-u^{n+1}_\ell-u^n_\ell\right]$$ Hence $$v^{n+1}_\ell-\frac{c\Delta t}{2\Delta x}(v^{n+1}_{\ell+1}-v^{n+1}_\ell)=v^n_\ell+\frac{c\Delta t}{2}\left[\frac{v^n_{\ell+1}-v^n_\ell}{\Delta x}-u^{n+1}_\ell+u^n_\ell\right]$$ i.e. $$\left(\mathrm{I}-\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)v^{n+1}_\ell=\left(\mathrm{I}+\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)v^n_\ell-\frac{c\Delta t}{2\Delta x}(u^{n+1}_\ell+u^n_\ell)$$ where $$\mathrm{A}=\begin{pmatrix} -1 & 1 & & 0 \\ &\ddots & \ddots \\ & & -1 & 1 \\ 0 & & & -1 \end{pmatrix}$$ Now, since we want to go backwards in time, we replace $\Delta t$ with $-\Delta t$ and we get $$v^{n-1}_\ell=\left(\mathrm{I}+\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)^{-1}\left[\left(\mathrm{I}-\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)v^n_\ell+\frac{c\Delta t}{2\Delta x}(u^{n-1}_\ell+u^n_\ell)\right]$$ Now, I try to implement this on Matlab via

%% Parameters
L = 5; % size of domain 
T = 5; % measurement time
dx = 1e-2; % spatial step
dt = 1e-3; % time step
x0 = 0; % point of measurement
c = 1; % speed of advection
t = 0:dt:T; % time vector
x = [0:dx:L]'; % position vector
nt = length(t); % number of time steps
nx = length(x); % number of position steps
mu = dt/dx;
I = eye(nx,nx+1); % identity matrix
A = spdiags(ones(nx,1)*[-1 1],0:1,nx,nx);
%% Solve Backward Advection Equation
% preallocate the memory
v = zeros(nx,nt);
u = zeros(nx,nt);
% final condition, sinc function centered around x = .5
v(:,nt) = sinc((x-x0)/dx);
for k = nt-1:1
    v(:,k-1) = (I+(c/2)*mu*A)\((I-(c/2)*mu*A)*v(:,k)+(c/2)*mu*(u(:,k-1)+u(:,k)));

My solution simply "blows up" on the boundary of the domain and is zero everywhere else. I have similar results using backwards Euler instead of Crank-Nicolson. Where am I going wrong?

I provide two plots: enter image description here enter image description here

  • $\begingroup$ Since your testcase has $u=0$, try comparing the solution with the exact solution (from the method of characteristics). This could just be a simple programming error: for example, I notice that you aren't using v(:,nt), and you start your recurrence from v(:,nt-1) instead. Also, I think your solution on the boundary doesn't blow up, because you are going backwards in time, so it's equal to the IC at the end, and is zero everywhere else. Incorrect solutions that blow up usually blow up to infinity, not 1. $\endgroup$
    – Kirill
    Jul 4 '17 at 20:41
  • $\begingroup$ @Kirill I start from v(:,nt-1)' instead of v(:,nt) since I "start" from a "final" condition. $\endgroup$
    – Jason Born
    Jul 4 '17 at 21:36
  • $\begingroup$ Your code never uses the IC at all: the rhs only evaluates v(:,nt-1) and down, which are all zero. $\endgroup$
    – Kirill
    Jul 4 '17 at 22:45
  • $\begingroup$ @Kirill I understand, but changing from k=nt-1:1 to k=nt:1 appears to have no effect $\endgroup$
    – Jason Born
    Jul 5 '17 at 9:07
  • $\begingroup$ Okay, so it turns out that I also needed to introduce a step factor, namely k=nt:-1:1, however, I then get the error Subscript indices must either be real positive integers or logicals. $\endgroup$
    – Jason Born
    Jul 5 '17 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.