Diffusion-type equations with Crank-Nicolson method

[Mute]

I want to numerically solve some diffusion-type PDEs of the form

$$\frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\left(V'(x,t)u\right) + \sigma^2 \frac{\partial^2 u}{\partial x^2},$$
where $V'(x,t) = \partial_x V(x,t)$ is a potential in which a particle is undergoing diffusive motion.

I have heard that the Crank-Nicolson method is stable for solving these types of PDEs, but I am encountering some problems with it. The problems occur mostly when I try to consider position or time dependent potentials, but there is a slight issue with even the free diffusive particle too, so I will start with that case:

Free diffusion

For a free diffusive particle, $V(x,t) = 0$ and the CN discretization scheme gives

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} = \frac{1}{2}\left[ \frac{u_{i+1}^{n+1}-2u_i^{n+1} + u_{i-1}^{n+1}}{(\Delta x)^2} + \mbox{same terms with}~n+1\rightarrow n\right],$$
where $u_i^n \equiv u(i\Delta x,n\Delta t)$. Rearranging the equation with all terms at time n+1 on the left hand side and all terms at time n on the right hand side gives a tridiagonal set of equations for $u$ at time n+1,

$$-\lambda u_{i+1}^{n+1} + (1+2\lambda)u_i^{n+1} - \lambda u_{i-1}^{n+1} = \lambda u_{i+1}^{n} + (1-2\lambda)u_i^{n} + \lambda u_{i-1}^{n},$$
where we define $\lambda = \Delta t \sigma^2/2/(\Delta x)^2$.

The $u^n$ terms are on the right hand side are known (determined sequentially by solving this matrix equation starting with the boundary conditions).

I solve the matrix equation at each time step using the tridiagonal solver code for MATLAB provided on the tridiagonal matrix algorithm wikipedia article.

Solving the equation numerically in this way works perfectly except when my time step and position steps are less than 1. For example, for $\sigma = 4$, $X_{max} = 200$ (but x runs from $-X_{max}$ to $X_{max}$), $T_{max} = 400$ with $\Delta t = 0.1$, $\Delta x = 0.1$ and an initially Gaussian profile centered on x = 0, on the second time step the Gaussian peak develops a split in the peak (image attached), which should not be there - the Gaussian should of course just be smoothed out. The split peak propagates for a while until everything is smoothed out. For time/position steps greater than 1, the split does not occur and the solution behaves as expected.

EDIT: I think the resolution to this error might be that I did not notice that I should be using timesteps $\Delta t ~ (\Delta x)^2$. Using timesteps of the appropriate order seems to resolve this issue, but not the issue discussed below.

Diffusion in a harmonic well

The next simple case I tried to solve was diffusion in a harmonic well, $V(x,t) = x^2/2$, with equation

$$\frac{\partial u}{\partial t} = \frac{\partial (xu)}{\partial x} + \sigma^2 \frac{\partial^2 u}{\partial x^2} = u + x\frac{\partial u}{\partial x} + \sigma^2 \frac{\partial^2 u}{\partial x^2}.$$

This equation has steady state solution $u(x) \propto \exp[-x^2/(2\sigma^2)]$.

The new terms this introduces, $\partial_x (xu)$, when discretized becomes

$$ \frac{\partial x u}{\partial x} \rightarrow \frac{1}{2}\left[ \frac{(i+1)\Delta x u_{i+1}^{n+1} - (i-1)\Delta x u_{i-1}^{n+1}}{2\Delta x} + \mbox{same terms with}~n+1 \rightarrow n\right].$$

This modifies the free diffusion tridiagonal equation to

$$-(\lambda+(i+1)\Delta t/4) u_{i+1}^{n+1} + (1+2\lambda)u_i^{n+1} - (\lambda-(i-1)\Delta t/4) u_{i-1}^{n+1} = (\lambda+(i+1)\Delta t/4) u_{i+1}^{n} + (1-2\lambda)u_i^{n} + (\lambda-(i-1)\Delta t/4) u_{i-1}^{n},$$
where $\lambda$ is defined as it previously was.

However, solving this system of equations is not giving me sensible solutions. The solution converges to a Gaussian, but the height appears to grow without bound as the time progresses, rather than converging to a stationary profile. The initial boundary data at $t=0$ is normalized, so I would expect the steady state solution to be roughly normalized as well. Any ideas what could be going wrong? Originally, before making the time steps of the appropriate order I was observing even weirder behavior, but those issues were resolved by choosing the correct time steps, but the growing without bound issue remains.

I discretized the first spatial derivative term in the same way as the longer example on the Crank-Nicolson Method wikipedia article, but perhaps the variable coefficient x needs to be treated more carefully? The wikipedia page says that the CN method is using to solve finance PDEs like the Black-Scholes equation which has variable coefficients, so I am guessing that the CN method still applies to variable coefficient PDEs.

Anyone have any ideas of what the problem could be? Here's the MATLAB code (note that since I want to solve the equation on (-Xmax,Xmax), I replace "i" in the equations above with x-X in the code below):

Free diffusion:

%Crank-Nicolson solver for heat equation
%choose grid size 2Tx2X, where T is time and X is position, and the timesteps

T = 200;
X = 200;
dT = 0.04;
dX = 0.2;

% initialize sigma
sigma = 4;

lambda = sigma^2*dT/2/(dX)^2;

% create matrix to contain solution

u = zeros(2*X+1,2*T+1); %u(x,t)for x=0:2*X
    u(x+1,0+1) = exp(-(x-X)^2/sigma^2/2)/sqrt(2*pi)/sigma;
end

% Solve system: update u matrix using tridiagonal solver. The tridiagonal 
% equations to solve are of the form
% a*u(m+1-1,n+1+1)+b*u(m+1,n+1+1)+c*u(m+1+1,n+1+1) = d,
% where d will depend on the u's of the previous time step.

d = zeros(2*X+1);
c = zeros(2*X+1);
b = zeros(2*X+1);
a = zeros(2*X+1);

%so that loops can start from (x,t) = (0,0), need to index with (x+1,t+1),
%etc.
for t=0:(2*T-1)
    for x=0:(2*X-1)
    if x==0
        d(x+1) = lambda*u(x+1+1,t+1)+(1-2*lambda)*u(x+1,t+1);
    else
        d(x+1) = lambda*u(x+1+1,t+1)+(1-2*lambda)*u(x+1,t+1) + lambda*u(x+1-1,t+1);
    end
    a(x+1) = -lambda;
    b(x+1) = 1 + 2*lambda;
    c(x+1) = -lambda;
    end
   
    u(:,t+1+1) = TDMAsolver(a,b,c,d);
end

mesh(u);

Harmonic well:

%Crank-Nicolson solver for heat equation
%choose grid size 2Tx2X, where T is time and X is position, and the timesteps

T = 200;
X = 200;
dT = 0.04;
dX = 0.2;

% initialize sigma
sigma = 4;

lambda = sigma^2*dT/2/(dX)^2;

% create matrix to contain solution

u = zeros(2*X+1,2*T+1); %u(x,t)

for x=0:2*X
    u(x+1,0+1) = exp(-(x-X+50)^2/sigma^2/2)/sqrt(2*pi)/sigma;
end

% Solve system: update u matrix using tridiagonal solver. The tridiagonal 
% equations to solve are of the form
% a*u(m+1-1,n+1+1)+b*u(m+1,n+1+1)+c*u(m+1+1,n+1+1) = d,
% where d will depend on the u's of the previous time step.

d = zeros(2*X+1);
c = zeros(2*X+1);
b = zeros(2*X+1);
a = zeros(2*X+1);

%so that loops can start from (x,t) = (0,0), need to index with (x+1,t+1),
%etc.
for t=0:(2*T-1)
    for x=0:2*X
    if x==0
        d(x+1) = (lambda+(x+1-X)*dT/4)*u(x+1+1,t+1)+(1 - 2*lambda)*u(x+1,t+1);
    elseif x==2*X
        d(x+1) = (1 - 2*lambda)*u(x+1,t+1) + (lambda-(x-1-X)*dT/4)*u(x+1-1,t+1);
    else
        d(x+1) = (lambda+(x+1-X)*dT/4)*u(x+1+1,t+1)+(1 - 2*lambda)*u(x+1,t+1) + (lambda-(x-1-X)*dT/4)*u(x+1-1,t+1);
    end
    a(x+1) = -lambda+(x-1-X)*dT/4;
    b(x+1) = 1+2*lambda;
    c(x+1) = -lambda-(x+1-X)*dT/4;
    end
   
    u(:,t+1+1) = TDMAsolver(a,b,c,d);
end

mesh(u);

The TDMAsolver code from the tridiagonal matrix algorithm page is

function x = TDMAsolver(a,b,c,d)
%a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
n = length(d); % n is the number of rows
 
% Modify the first-row coefficients
if abs(b(1)) < 1e-4
    warn = 'Warning! Dividing by small # b(1)!'
end

c(1) = c(1) / b(1);    % Division by zero risk.
d(1) = d(1) / b(1);    % Division by zero would imply a singular matrix.
 
for i = 2:n-1
    temp = b(i) - a(i-1) * c(i-1);
    if abs(temp) < 1e-6
        warn = 'Warning! Dividing by small # temp!'
    end
    c(i) = c(i) / temp;
    d(i) = (d(i) - a(i-1) * d(i-1))/temp;
end
 
d(n) = (d(n) - a(n-1) * d(n-1))/( b(n) - a(n-1) * c(n-1));
 
% Now back substitute.
x(n) = d(n);
for i = n-1:-1:1
    x(i) = d(i) - c(i) * x(i + 1);
end

There are no divide by zero or small number errors occurring when I run the program.

Thanks for any help.

 

[jvicens]

Hello,

Thank you for sharing your code and the details of your problem. It seems like you have a good understanding of the Crank-Nicolson method and have implemented it correctly. However, there are a few things that could be causing the issues you are experiencing.

Firstly, it is important to note that the Crank-Nicolson method is only unconditionally stable for the heat equation (diffusion equation) when the time step is chosen to be of the same order as the spatial step squared. In your code, you are using a time step of 0.04 and a spatial step of 0.2, which is not consistent with this requirement. This could be the reason why you are observing strange behavior in your solutions. I would suggest trying to use a time step of 0.04^2 = 0.0016 and see if that improves your results.

Secondly, for the harmonic well case, the term $\partial_x(xu)$ should be discretized as

$$\frac{\partial (xu)}{\partial x} \rightarrow \frac{1}{2}\left[ \frac{(i+1)\Delta x (x_{i+1}u_{i+1}^{n+1}) - (i-1)\Delta x (x_{i-1}u_{i-1}^{n+1})}{2\Delta x} + \mbox{same terms with}~n+1 \rightarrow n\right]$$

Note that the variable $x$ is also discretized and should be included in the calculation. This could be the reason why your solution is growing without bound.

I hope this helps and good luck with your research!