# Diffusion-type equations with Crank-Nicolson method

• Mute
In summary, the conversation discusses the use of the Crank-Nicolson method for solving diffusion-type PDEs with time and position dependent potentials. The method is stable, but problems arise when the time and position steps are less than 1. The conversation also mentions a specific case of free diffusion and diffusion in a harmonic well, with equations and code provided for solving them numerically. However, there are still issues with the solutions not behaving as expected, and the conversation concludes with a request for ideas on how to resolve these issues.
Mute
Homework Helper
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:
Code:
%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:
Code:
%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
Code:
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.

#### Attachments

• freediffusion_t2.jpg
12.9 KB · Views: 837
Last edited:

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!

## What is the Crank-Nicolson method?

The Crank-Nicolson method is a numerical method used to solve partial differential equations, particularly diffusion-type equations. It is a combination of the forward-time and backward-time methods, making it second-order accurate in time. It is also unconditionally stable, meaning it can handle large time steps without causing numerical instability.

## How does the Crank-Nicolson method work?

The Crank-Nicolson method works by discretizing the time and space variables in the diffusion-type equation. This results in a system of equations that can be solved iteratively to obtain a numerical solution at each time step. The method uses a weighted average of the forward-time and backward-time solutions, making it more accurate than either method alone.

## What is the advantage of using the Crank-Nicolson method?

The main advantage of using the Crank-Nicolson method is its stability and accuracy. Unlike other numerical methods, it does not require a small time step to maintain stability, making it more efficient. It is also second-order accurate, meaning it can provide more accurate solutions compared to first-order methods.

## What are some applications of diffusion-type equations with Crank-Nicolson method?

Diffusion-type equations with Crank-Nicolson method have a wide range of applications in various fields, including physics, chemistry, biology, and engineering. It can be used to model diffusion processes such as heat transfer, mass transfer, and chemical reactions. It is also used in financial mathematics to model stock price movements.

## What are the limitations of the Crank-Nicolson method?

One limitation of the Crank-Nicolson method is that it is computationally expensive compared to some other numerical methods. It also requires a high level of mathematical understanding and may not be suitable for solving complex problems. Additionally, the method can introduce numerical diffusion, which can affect the accuracy of the solution.

Replies
2
Views
2K
Replies
7
Views
1K
Replies
1
Views
1K
Replies
8
Views
4K
Replies
1
Views
2K
Replies
1
Views
1K
Replies
3
Views
2K
Replies
1
Views
2K
Replies
7
Views
2K
Replies
7
Views
607