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[0]=(1-alpha)*u[0]
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[0]=4
u[-1]=4
a=np.full(nx, -alpha)
b=np.full(nx+1, 2*(1+alpha))
b[0]=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?
