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

$$\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?

• This is just the spatial average of the forward differences, and should still be first order accurate. What paper are you reading?
– Jeff
Commented Aug 16, 2016 at 17:01
• @Jeff The paper I was reading is: Cook, P. G., and Solomon, D. K., Transport of atmospheric trace gases to the water table: Implications for groundwater dating with chlorofluorocarbons and krypton 85, Water Resources Research, Vol. 31, No. 2, 263-270, Feb 1995. Commented Aug 16, 2016 at 22:14
• @Jeff Do you know why the average is weighted the way it is 1/6, 4/6, 1/6? Commented Aug 16, 2016 at 22:15
• Does $D_e$ stand here for constant? Commented Aug 17, 2016 at 16:02
• It's not particularly important exactly what the weights are. I'll add a more detailed answer below.
– Jeff
Commented Aug 17, 2016 at 16:10

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.