# Numerical solution of the advection equation with Crank–Nicolson finite difference method

I need to implement a numerical scheme for the solution of the one-dimensional advection equation $$\\\frac{\partial u}{\partial t} + C(x, t) \frac{\partial u}{\partial x} = 0 \\\\$$ $$\\ C(x,t) = \frac{\pi\cos(2 \pi t) + 3.5}{2\pi - \pi\sin(\pi x)} \\\\$$ $$\\ 0 \leq x \leq 1, \quad 0 \leq t \leq 1\\\\$$

with Crank–Nicolson finite difference method $$\\ \frac{u^{j + 1}_{k} - u^{j}_{k}}{\tau} + \frac{C}{4h}(u^{j + 1}_{k + 1} - u^{j + 1}_{k - 1} + u^{j}_{k + 1} - u^{j}_{k - 1}) = 0\\\\$$ and boundary conditions $$\\ \phi(x) = u(x,0), \quad \psi_{0}(t) = u(0,t), \quad \psi_{1}(t) = u(1, t) \\\\$$

it's known that the exact solution has the following form $$\\ U(x,t) = \cos(\pi x) - \frac{1}{2}\sin(2 \pi t) + 2\pi x - 3.5t\\\\$$

The numerical solution on each layer in time is found as the solution of a system of linear equations with a tridiagonal matrix $$\\ Au = b \\\\$$ $$\\ A = \begin{pmatrix} 1 & \frac{C\tau}{4h} & 0 & ...\\ -\frac{C\tau}{4h} & 1 & \frac{C\tau}{4h} & ... \\ 0 & -\frac{C\tau}{4h} & 1 & ... \\ ... & ... & ... & ... \end{pmatrix}\\\\$$

I've write a simple python script, but the result doesn't match the theory:

I'm sure I'm wrong somewhere, but I don't see where. Thanks!

• Speed $C$ is positive at $x=1$. You cannot specify dirichlet conditions at $x=1$. – cfdlab Jan 5 at 13:47

Some remarks about the code, perhaps cleaning these up already solves the question, especially the initialization error:

Floating point errors are a thing, so if you want to get the expected results, give some wiggle room (0.01 is arbitrary, I would not expect defects above 1e-12 in this place)

N = int((x1 - x0) / h + 0.01)
M = int((t1 - t0) / tau + 0.01)

There is no guarantee that $$N·h=x_1-x_0$$ etc., so force this grid condition

h = (x1 - x0) / N
tau = (t1 - t0) / M

There is a honest error in the side boundary initialization, the space component is always first

u[0, i] = psi_0(i * tau)
u[N, i] = psi_1(i * tau)

The implicit trapezoidal method that is the method in time direction of C-N reads as $$\frac{y^{j+1}-y^j}{\tau}=\frac{F(t^j,y^j)+F(t^{j+1},y^{j+1})}2.$$ Accordingly, the coefficients a_h for the banded matrix should be computed at time (j+1)*tau to get the correct order $$O(\tau^2+h^2)$$ of the method.

• Thank you for the remarks, I am little ashamed of such obvious errors. Unfortunately, fixing these doesn't solve the question, perhaps there are more issues. – Winabryr Jan 5 at 19:40