I am implementing the advection equation $u_x+(1/c)u_t=0$ following a Crank-Nicholson finite difference scheme. The equation for this is \begin{eqnarray*} -\frac{\gamma}{4} w_{n-3 j+1} + w_{n-2 j+1} + \frac{\gamma}{4} w_{n-1 j+1} = \frac{\gamma}{4} w_{n-3 j} + w_{n-2 j} - \frac{\gamma}{4} w_{n-1 j}. \end{eqnarray*}
In matrix form it can be written as:
\begin{eqnarray*} A \bf{w}_{j+1} = \it B {\bf{w}_j } + \bf{b} \end{eqnarray*} where
\begin{eqnarray*} A = \begin{pmatrix} 1 & \frac{\gamma}{4} & 0 & 0 & \cdots & 0 \\ -\frac{\gamma}{4} & 1 & \frac{\gamma}{4} & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & -\frac{\gamma}{4} & 1 & \frac{\gamma}{4} \\ 0 & 0 & \cdots & \cdots & -\frac{\gamma}{4} & 1 \end{pmatrix} \; , \; \bf{w}_{j+1} = \begin{pmatrix} w_{1 j+1} \\ w_{2 j+1} \\ \vdots \\ \vdots \\ \vdots \\ w_{n-2 j+1} \\ \end{pmatrix} \end{eqnarray*}
\begin{eqnarray*} B = \begin{pmatrix} 1 & -\frac{\gamma}{4} & 0 & 0 & \cdots & 0 \\ \frac{\gamma}{4} & 1 & -\frac{\gamma}{4} & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{\gamma}{4} & 1 & -\frac{\gamma}{4} \\ 0 & 0 & \cdots & \cdots & \frac{\gamma}{4} & 1 \end{pmatrix} \; , \; \bf{w}_{j} = \begin{pmatrix} w_{1 j} \\ w_{2 j} \\ \vdots \\ \vdots \\ \vdots \\ w_{n-2 j} \\ \end{pmatrix} \; , \; \bf{b} = \begin{pmatrix} \frac{\gamma}{4} w_{0 j+1} + \frac{\gamma}{4} w_{0j} \\ 0 \\ \vdots \\ 0 \\ -\frac{\gamma}{4} w_{n-1 j+1} - \frac{\gamma}{4} w_{n-1 j} \end{pmatrix} \end{eqnarray*}
As an initial condition I have a Gaussian pulse $u(x,0)=\mathrm{e}^{-400(x-1/2)^2}$. The spacial domain is $x \in [0,1]$. Constant velocity $c=1$, $dt=0.001$, $dx=0.002$ and $\gamma=c dt/dx$.
My implementation is in Python and I can provide it upon request.
I attach a few screen captures of the "wavefield".
Can someone please help me understand what is going on after the pulse hits the left boundary? I am using Dirichlet boundary conditions with 0 on the boundary.
Thanks.
*****Update*****
I found someone else asking the same question and some interesting answers which I am still trying to digest. Still I would like to have more discussion on this issue.
Here is the link on this
