Here we extend our discussion and implementation of the Crank- Nicolson (CN) method to convection-diffusion systems.

To clarify nomenclature, there is a physically important difference between
convection and advection.
Since we are interested in the transport of a protein suspended in a
(semi)liquid medium, we will use
the term *advection* (as opposed to *convection*) in the following discussion.

We study the following advection-diffusion equation:

with concentration $u$, diffusion coefficient $D$, advection velocity $a$ (velocity of the medium), reaction term $f$, and domain length $L$. Our Neumann boundary conditions make certain that no amount of $u$ enters or leaves our domain.

We use the same grid as before. Also see our earlier discussion of finite difference stencils.

As described, for instance, in Trefethen Table 3.2.1 the Crank-Nicolson stencil for the advection term takes the following form (assuming that $a$ is constant)

Extending our previous CN stencil with this expression and applying the resulting, extended stencil to our advection-diffusion equation in grid point $(j,n)$ we obtain:

As before, we define $\sigma = \frac{D \Delta t}{2 \Delta x^2}$ and $\rho = \frac{a \Delta t}{4 \Delta x}$ and rearrange the above approximation of our advection-diffusion equation:

This equation holds for all grid points with spatial indices $j=1,\ldots,J-2$. On the two boundaries of our spatial grid we apply the above Neumann boundary conditions that dictate

resulting in the following amended expressions for $j=0$ and $j=J-1$:

We use the same vector notation, $\mathbf{U}^n$ and $\mathbf{f}^n$, and matrix names $A$ and $B$ as before and can now write the above equations as a linear system in matrix notation:

where the tridiagonal matrix $A$ has the following vector of length $J$ on its diagonal

the following vector of length $J-1$ on its first superdiagonal

and the following vector of length $J-1$ on its first subdiagonal:

The tridiagonal matrix $B$ has the following vector of length $J$ on its diagonal

the following vector of length $J-1$ on its first superdiagonal

the following vector of length $J-1$ on its first subdiagonal

Let us implement the above numerical scheme for the same example we discussed previously.

Construction of the linear system we solve to integrate our advection-diffusion system numerically over one time step is very similar to our previous linear system.

We therefore expect, that we can reuse our previous code
and amend merely matrices `A_u`

, `B_u`

, `A_v`

, and `B_v`

as described above.

Here is just a copy-and-paste of our previous code. Refer to that post for explanations.

All that changes in this code segment compared with before
are our definition of `rho`

and construction of matrices `A_u`

, `B_u`

, `A_v`

,
`B_v`

.
(We actually also change `J`

and `N`

- see comment below)

Note that `rho`

is the same for both protein species `U`

and `V`

.

As advection velocity `a`

we choose a negative value so that the fluid carrying
proteins `U`

and `V`

flows from right to left.

```
import numpy
from matplotlib import pyplot
numpy.set_printoptions(precision=3)
L = 1.
J = 500
dx = float(L)/float(J-1)
x_grid = numpy.array([j*dx for j in range(J)])
T = 200
N = 2000
dt = float(T)/float(N-1)
t_grid = numpy.array([n*dt for n in range(N)])
D_v = float(10.)/float(100.)
D_u = 0.01 * D_v
k0 = 0.067
f_vec = lambda U, V: numpy.multiply(dt, numpy.subtract(numpy.multiply(V,
numpy.add(k0, numpy.divide(numpy.multiply(U,U), numpy.add(1., numpy.multiply(U,U))))), U))
sigma_u = float(D_u*dt)/float(2.*dx*dx)
sigma_v = float(D_v*dt)/float(2.*dx*dx)
a = -0.0003
rho = float(a*dt)/float(4.*dx)
total_protein = 2.26
no_high = 10
U = numpy.array([0.1 for i in range(no_high,J)] + [2. for i in range(0,no_high)])
V = numpy.array([float(total_protein-dx*sum(U))/float(J*dx) for i in range(0,J)])
A_u = numpy.diagflat([-sigma_u+rho for i in range(J-1)], -1) +\
numpy.diagflat([1.+sigma_u+rho]+[1.+2.*sigma_u for i in range(J-2)]+[1.+sigma_u-rho]) +\
numpy.diagflat([-(sigma_u+rho) for i in range(J-1)], 1)
B_u = numpy.diagflat([sigma_u-rho for i in range(J-1)], -1) +\
numpy.diagflat([1.-sigma_u-rho]+[1.-2.*sigma_u for i in range(J-2)]+[1.-sigma_u+rho]) +\
numpy.diagflat([sigma_u+rho for i in range(J-1)], 1)
A_v = numpy.diagflat([-sigma_v+rho for i in range(J-1)], -1) +\
numpy.diagflat([1.+sigma_v+rho]+[1.+2.*sigma_v for i in range(J-2)]+[1.+sigma_v-rho]) +\
numpy.diagflat([-(sigma_v+rho) for i in range(J-1)], 1)
B_v = numpy.diagflat([sigma_v-rho for i in range(J-1)], -1) +\
numpy.diagflat([1.-sigma_v-rho]+[1.-2.*sigma_v for i in range(J-2)]+[1.-sigma_v+rho]) +\
numpy.diagflat([sigma_v+rho for i in range(J-1)], 1)
U_record = []
V_record = []
U_record.append(U)
V_record.append(V)
for ti in range(1,N):
U_new = numpy.linalg.solve(A_u, B_u.dot(U) + f_vec(U,V))
V_new = numpy.linalg.solve(A_v, B_v.dot(V) - f_vec(U,V))
U = U_new
V = V_new
U_record.append(U)
V_record.append(V)
```

We notice that the numerical stability of this scheme is greatly dependent on
the magnitude of `a`

.

As was discussed by Walther *et
al.*,
preservation of total protein mass in a certain regime near ```
total_protein =
2.26
```

is important to
preserve the pattern-forming pattern of the
original system (excluding the advection term).

Notice that even though we increased `J = 500`

and `N = 2000`

for a much finer
numerical grid
compared with
before
(`J = 100`

, `N = 1000`

), we still lose some protein mass to numerical errors.

Initial protein mass:

```
sum(numpy.multiply(dx,U_record[0]) + numpy.multiply(dx,V_record[0]))
2.259999999999998
```

Protein mass at the end of the simulation:

```
sum(numpy.multiply(dx,U_record[-1]) + numpy.multiply(dx,V_record[-1]))
2.2064888602225103
```

We can reduce the loss of protein mass by reducing the magnitude of `a`

or
refining our grid even further
by increasing both `J`

and `N`

.

An important question to answer would be how much numerical accuracy we gain by
increasing `J`

and `N`

.

And here is a kymograph showing the time and space behaviour of `U`

.

```
U_record = numpy.array(U_record)
V_record = numpy.array(V_record)
fig, ax = pyplot.subplots()
pyplot.xlabel('x'); pyplot.ylabel('t')
heatmap = ax.pcolor(x_grid, t_grid, U_record, vmin=0., vmax=1.2)
colorbar = pyplot.colorbar(heatmap)
colorbar.set_label('concentration U')
```