# Crank Nicolson Method with closed boundary conditions

I want to simulate 1D diffusion with a constant diffusion coefficient using the Crank-Nicolson method. $$\frac{\partial u (x,t)}{\partial t} = D \frac{\partial^2 u(x,t)}{\partial x^2}.$$

I take an average of the approximations at the points $$(x_i, t_j)$$ and $$(x_i, t_{j+1}):$$

$$\frac{u_i^{n+1} -u_i^n}{\Delta t} \approx D \frac{(u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1})+ (u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n})}{2 \Delta x^2}$$

Letting $$\alpha = D \frac{\Delta t}{\Delta x^2}$$ yields: $$-\alpha u_{i+1}^{j+1} + 2(1+ \alpha)u_{i}^{j+1} - \alpha u_{i-1}^{j+1}= \alpha u_{i+1}^{j} + 2(1- \alpha)u_{i}^{j} + \alpha u_{i-1}^{j}.$$

In matrix form $$Au^{j+1} = Bu^{j}$$ Where $$A$$ is the tridiagonal matrix $$A = \begin{bmatrix} b_0 & c_0 & & & & & & 0 \\ a_1 & b_1 & c_1 & & & & & \\ & a_2 & b_2 & c_2 & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & a_i & b_i & c_i & & \\ & & & &\ddots & \ddots & \ddots & \\ & & & & & a_{n_x -1} & b_{n_x -1} & c_{n_x -1} \\ 0 & & & & & & a_{n_x} & b_{n_x} \end{bmatrix} = \begin{bmatrix} ? & ? & & & & & & 0 \\ -\alpha & 2(1+\alpha) & -\alpha & & & & & \\ & -\alpha & 2(1+\alpha) & -\alpha & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & -\alpha & 2(1+\alpha) & -\alpha & &\\ & & & &\ddots & \ddots & \ddots & \\ & & & & & -\alpha & 2(1+\alpha) & -\alpha \\ 0 & & & & & & ? & ? \end{bmatrix} =$$ And $$B$$ is the tridiagonal matrix $$B = \begin{bmatrix} ? & ? & & & & & & 0 \\ \alpha & 2(1-\alpha) & \alpha & & & & & \\ & \alpha & 2(1-\alpha) & \alpha & & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & \alpha & 2(1-\alpha) & \alpha & &\\ & & & &\ddots & \ddots & \ddots & \\ & & & & & \alpha & 2(1-\alpha) & \alpha \\ 0 & & & & & & ? & ? \end{bmatrix}$$

I have insulated boundary conditions, so to determine the top rows I first introduce a fictitious points $$u_{-1}$$

$$\frac{u_{1}^{j+1} - u_{-1}^{j+1}}{2 \Delta x} = \text{Flux}(0) =0$$ $$u_{-1}^{j+1} = u_{1}^{j+1}$$

Implying: $$-\alpha u_{1}^{j+1} + 2(1+ \alpha)u_{0}^{j+1} - \alpha u_{-1}^{j+1}= \alpha u_{1}^{j} + 2(1- \alpha)u_{0}^{j} + \alpha u_{-1}^{j}$$ $$-2\alpha u_{1}^{j+1} + 2(1+ \alpha)u_{0}^{j+1} = 2\alpha u_{1}^{j} + 2(1- \alpha)u_{0}^{j}$$

So I set the top row of $$A$$ to be $$\begin{bmatrix} 2(1+\alpha) & -2\alpha & ... \end{bmatrix}$$ and the top row of $$B$$ to be $$\begin{bmatrix} 2(1-\alpha) & 2\alpha & ... \end{bmatrix}$$. A similar derivation applies to the bottom rows.

The problem is that this does not seem to do what I want; it does not conserve flux at the boundary. Here is a sample implementation in Python. I first create a function to multiply by the matrix $$B$$ to save memory:

def mltply_B(u,alpha):
u_i = 2*(1-alpha)*u
u_i=np.insert(u_i,0,0)
u_i=np.append(u_i,0)

u_iplus1=alpha*u
u_iplus1=np.append(u_iplus1,[0,0])

u_iminus1=alpha*u
u_iminus1=np.insert(u_iminus1,0,0)
u_iminus1=np.insert(u_iminus1,0,0)

B_u=u_i+u_iplus1+u_iminus1
B_u=B_u[1:-1]
B_u=(1-alpha)*u
B_u[-1]=(1-alpha)*u[-1]

return B_u


Next I implement the tridiagonal matrix algorithm, or Thomas algorithm, to solve the linear equation:

def TDMAsolver(a, b, c, d):
nf = len(d)  # number of equations
ac, bc, cc, dc = map(np.array, (a, b, c, d))  # copy arrays
for it in range(1, nf):
mc = ac[it-1] / bc[it-1]
bc[it] = bc[it] - mc * cc[it-1]
dc[it] = dc[it] - mc * dc[it-1]

xc = bc
xc[-1] = dc[-1] / bc[-1]

for il in range(nf-2, -1, -1):
xc[il] = (dc[il] - cc[il] * xc[il+1]) / bc[il]

return xc


I create a vector which looks like [4,0,0,...,0,0,4] and observe the output: nx=150 alpha=0.07

u=np.zeros(nx+1)
u=4
u[-1]=4

a=np.full(nx, -alpha)
b=np.full(nx+1, 2*(1+alpha))
b=1+alpha
b[-1]=1+alpha
c=np.full(nx, -alpha)

d=mltply_B(u,alpha)
u1=TDMAsolver(a,b,c,d)


Now the output of

np.sum(u)


Is 8.0, as expected, but the output of

np.sum(u1)


is roughly 7.49, an approximately $$1/16$$ loss of mass. If you run this in a loop the amount of mass gradually dwindles (after about 50 loops you only have half of your mass left).

If I initialize $$u$$ so that the mass is concentrated in the center then there is very little loss, at least initially. Eventually, once a substantial amount of mass reaches the endpoints mass leaves the system. In the eventual application I am interested in, I will be injecting a "source" every time step at one of the boundaries, so it is important for the behavior near the boundaries to be accurate.

Edit: If a right and left difference are used to derive the boundary conditions then we take instead:

$$\frac{u_{0}^{j+1} - u_{-1}^{j+1}}{\Delta x} = \text{Flux}(0) =0$$ $$u_{-1}^{j+1} = u_{0}^{j+1}$$

Implying: $$-\alpha u_{1}^{j+1} + 2(1+ \alpha)u_{0}^{j+1} - \alpha u_{-1}^{j+1}= \alpha u_{1}^{j} + 2(1- \alpha)u_{0}^{j} + \alpha u_{-1}^{j}$$ $$-\alpha u_{1}^{j+1} + (2+ \alpha)u_{0}^{j+1} = \alpha u_{1}^{j} + (2- \alpha)u_{0}^{j}$$

So we set the top row of $$A$$ to be $$\begin{bmatrix} (2+\alpha) & -\alpha & ... \end{bmatrix}$$ and the top row of $$B$$ to be $$\begin{bmatrix} (2-\alpha) & \alpha & ... \end{bmatrix}$$. (And again, a similar derivation applies to the bottom rows.)

This does conserve mass; however, multiple sources imply that that it is more accurate (specifically, accurate to second rather than first order) to use the centered difference at the boundary. So what gives? Is there a way to keep the second order accuracy at the boundary and also conserve mass?

• Can you let me know the boundary condition that you are planning to impose on the problem? If its dirichlet, Impose the BC strongly ( By making the diagonal 1 and other entries zero and substituting the RHS). Jun 14 at 4:24
• @ThivinAnandh I am sorry, I do not understand your comment. The boundary condition is that the flux at the ends is zero, or the derivative at the ends is zero, which is a Neumann boundary condition. Jun 14 at 14:14

I tried to run your code and I guess there might be a small mistake in your derivation,

When you impose $$u_{-1}^{j+1} = u_1^{j+1}$$ into the equations at the ends, you have to equate $$u_{-1}^{j+1} = u_0^{j+1}$$ and not $$u_{-1}^{j+1} = u_1^{j+1}$$

When this modification is applied which basically changes $$2(1 \pm \alpha)$$ in $$A$$ and $$B$$ to be $$(2 \pm \alpha)$$, then there is mass conservation, even for multiple time steps.

EDIT: One tip for mass conservation is to check if the sum of the rows and columns of the resulting matrix gives the same vector of the form $$a [1, 1,..., 1]^\top$$, the resulting matrices with central difference for insulating BC above does not hold this property.