Crank-Nicolson Finite Difference: Why use a forward time difference weighted across multiple positions?

Question

I was reading through a paper about the following PDE

$$ \frac{\partial u}{\partial t} = D_e \frac{\partial^2 u}{\partial z^2} - V_e \frac{\partial u}{\partial z}. $$

This PDE models the transport of gas tracers through unsaturated zones in soils.

Along with the paper I had a numerical solver for this PDE written by one of the paper's authors. In an attempt to understand the solver I wrote my own using the Crank-Nicolson method.

When I wrote my solver I approximated $ \frac{\partial u}{\partial t} $ using the forward difference approximation,

$$ \frac{\partial u}{\partial t} \approx \frac{u_{m + 1,\ j} - u_{m,\ j}}{\Delta t}, $$

where $u_{m,\ j}$ is the value of $u$ at the $m$th step in time and the $j$th step in position on the solution grid, and $\Delta t$ is the length of the time step which is constant.

This method obtained a solution close to the author's solver but after reading through the author's solver I noticed that he has used a different forward difference approximation for $ \frac{\partial u}{\partial t} $,

$$ \frac{\partial u}{\partial t} \approx \frac{u_{m + 1,\ j + 1} - u_{m,\ j + 1}}{6\Delta t} + 2\frac{u_{m + 1,\ j} - u_{m,\ j}}{3\Delta t} + \frac{u_{m + 1,\ j - 1} - u_{m,\ j - 1}}{6\Delta t}. $$

Is this approximation more accurate than the one I used? Should it always be used instead of the one I used?

It is not mentioned in the text I used as a reference for the Crank-Nicolson method. What is this weighted approximation called so I can look for a reference?

Why are the different forward difference approximations weighted as they are: one sixth for the positions before and after and two thirds for the current position?

Answer 1

To check the order of accuracy, you can simply plug a Taylor expansion into the finite difference scheme. If $u(x,t)$ is smooth and $u_{m,j} = u(m\Delta t,j\Delta x)$ then

$$u_{m+1,j+1}-u_{m,j+1} = u_t\Delta t + \frac{1}{2}u_{tt}\Delta t^2 + u_{xt} \Delta x\Delta t + O((\Delta t + \Delta x)^3),$$ $$u_{m+1,j} - u_{m,j} = u_t\Delta t + \frac{1}{2}u_{tt}\Delta t + O(\Delta t^3),$$ $$u_{m+1,j-1} - u_{m,j-1} = u_t\Delta t + \frac{1}{2}u_{tt}\Delta t^2 - u_{xt} \Delta x\Delta t + O((\Delta t + \Delta x)^3).$$

Suppose we average with weights $a,b,c$, so we form

$$A:= a\frac{u_{m+1,j+1}-u_{m,j+1}}{\Delta t} + b\frac{u_{m+1,j}-u_{m,j}}{\Delta t} + c\frac{u_{m+1,j-1}-u_{m,j-1}}{\Delta t}.$$

Then

$$A = (a+b+c)u_t + \frac{1}{2}(a+b+c)u_{tt}\Delta t + (a-c)u_{xt}\Delta x + O((\Delta t + \Delta x)^3/\Delta t).$$

The scheme is consistent with $u_t$ if $a+b+c=1$. It is first order accurate in time for any choices of weights. If $a\neq c$ it is first order in $x$ as well. If $a=c$, then the $\Delta x$ term drops out and you get $O(\Delta x^2 + \Delta x^3/\Delta t)$ accuracy in $x$.

So all that is important is $a=c$ and $a+b+c=1$. Your choice was $a=c=0$ and $b=1$, while theirs is $a=c=1/6$ and $b=2/3$. Both have the same accuracy.

The reason for choosing nonzero $a$ and $c$ could be related to stability. If you perform a stability analysis (i.e., compute the eigenvalues of the right hand side), you might find that averaging increases the values of $\Delta t$ and $\Delta x$ that yield a stable scheme. This is just a guess though.