Crank-Nicholson for diffusion-advection vs diffusion equation

Question

Let's consider the following 1D diffusion equation:

$\frac{\partial u}{\partial t} = xk \frac{\partial}{\partial x}(\frac{1}{x}\frac{\partial u}{\partial x})$

where we assume that the diffusion coefficient $k$ is constant.

In order to discretize this equation, let us define the following:

$x_i = 1 + i\Delta x, \forall i=0,1,...,I \\ t_n = n\Delta t, \forall n=0,1,...,N \\ r = \frac{k\Delta t}{2 \Delta x^2} \\ s = \frac{k\Delta t}{4 \Delta x} $

where $\Delta t$ and $\Delta x$ are the time step and the spacing, respectively.

Let's discretize the diffusion equation by means of the Crank-Nicholson method:

$u_i^{n+1}-u_i^n = r x_i [ \frac{u_{i+1}^n - u_i^n}{x_{i+\frac{1}{2}}} - \frac{u_i^n - u_{i-1}^n}{x_{i-\frac{1}{2}}}] + r x_i [ \frac{u_{i+1}^{n+1} - u_i^{n+1}}{x_{i+\frac{1}{2}}} - \frac{u_i^{n+1} - u_{i-1}^{n+1}}{x_{i-\frac{1}{2}}}]$

where upper indexes refer to time, lower indexes refer to space.

So far, so good. Now, using the product rule, the initial diffusion equation can be rewritten as:

$\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2} - \frac{k}{x}\frac{\partial u}{\partial x}$

The latter equation is now a diffusion-advection equation, but should be equivalent to the original diffusion equation.

Let us discretize the latter equation by means of the Crank-Nicholson method:

$u_i^{n+1}-u_i^n = r [ u_{i+1}^n - 2 u_i^n + u_{i-1}^n ] - \frac{s}{x_i} [ u_{i+1}^n - u_{i-1}^n] + r [ u_{i+1}^{n+1} - 2 u_i^{n+1} + u_{i-1}^{n+1}] - \frac{s}{x_i} [u_{i+1}^{n+1} - u_{i-1}^{n+1}]$

Those 2 discretizations are different. So, here is my question:

Which discretization should I choose?

Answer 1

Take the first expression and start to reduce the $x$ values to $x_i$, \begin{align} \frac{u_{i+1}^n - u_i^n}{x_{i+\frac{1}{2}}} - \frac{u_i^n - u_{i-1}^n}{x_{i-\frac{1}{2}}} &= \frac{(x_i-\frac12Δx)(u_{i+1}^n - u_i^n) - (x_i+\frac12Δx)(u_i^n - u_{i-1}^n)}{x_{i-\frac{1}{2}}x_{i+\frac{1}{2}}} \\ &=\frac{x_i}{x_i^2-\frac14Δx^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n) \\&\qquad - \frac{Δx}{2(x_i^2-\frac14Δx^2)}(u_{i+1}^n - u_{i-1}^n) \end{align} After multiplying with $rx_i$ you find that the difference to the second formula are only in terms of size $O(Δx^2)$, that is, within the error order of the method.