Recognitions:
Homework Help

Diffusion-type equations with Crank-Nicolson method

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 occuring when I run the program.

Thanks for any help.
Attached Thumbnails