Implementing Crank Nicholson in spherical setting

I am trying to solve the spherical diffusion equation using the Crank-Nicholson method. I think that I understand the method, I have included a derivation of the method in the file FD_diffusion.pdf and I have included my files for the code, spherical_fd.m is the main code used to solve the system and test_wrap.m is simply a wrapper function.

When you run the code to see what it does it seems to work relatively well, there isn't anything which seems to be wrong. However there seems to be two major issues which are apparent:

1. The solution changes dramatically depending on the fineness of the discretisation taken.
2. The conservation property doesn't seem to be obeyed

I have gone over my derivations time and time again, and I've looked at my code and there doesn't seem to by anything wrong but there obviously is but for the life of me, I can't figure it out.

[J. Alex Lee]

If you look at your u(r) solutions in time, you will see that your boundary condition at r=0 is not met based on just a visual check (a spike vs zero slope), and your boundary condition at r=R also does not appear to be met

% pick some random time
time_index = 10

% get the spatial solution
u_check = u(time_index,:)

% visual check
plot(r,u_check)

dr = r(2)
lam = D

% slope at r = R according to BC
q(time_index)/lam

% slope at r = R using FD
(u_check(end)-u_check(end-1))/dr

[Matthew Hunt]

I'm pretty sure that it is. I encoded it into the scheme.

[Torsten]

Yes, but the results show that the implementation must somehow be wrong.

But since this was your question, it doesn't help, I know.

Your u_average is wrong, by the way, as you can easily see if you integrate u=1 over r.

The average u is

u_avg = 3/R^3 * integral(r^2*u(r)  dr)

(comes from

1/(4/3*pi*R^3) * integral(4*pi*r^2*u(r) dr) 

).

[J. Alex Lee]

By the way, is this just an example for you to study the CN method and numerical solution schemes, or is it to solve a real problem? If latter, I would recommend looking into pdepe.

[Matthew Hunt]

I know how to use pdepe but it's too slow for my requirements whereas the direct CN sceme is fast.

[J. Alex Lee]

on my computer, running your algorithm spherical_fd takes ~0.85 seconds, whereas a pdepe implementation of your problem takes ~1.3 seconds...is your real application significantly different in scale that this difference is important?

[Matthew Hunt]

I've since got this code to work. I'm now working on the cylindrical case.

Answer 1

[Matthew Hunt]

Question resolved.

Answer 2

[J. Alex Lee]

I believe this should be the right formulation for CN (taking a few liberties interpreting your problem).

If you want speed, probably better to use sparse matrices and do the linear solve every iteration rather than pre-solve a full matrix and doing only the forward multiplication in time loop.

However, keep in mind that CN has a fixed accuracy in time whereas pdepe seems to use ode15s to advance so that it can have a higher order accuracy (and also will also adapt time step size), so probably compared to CN on the basis of same accuracy, it will probably be a lot faster.

D = 1e-3;
q0 = 0.15;
tc = 4;

xmesh = linspace(0,1,800).';
tspan = linspace(0,20,50);

N = numel(xmesh);

r = xmesh;

dr = xmesh(2)-xmesh(1);
dt = tspan(2)-tspan(1);

% helper
g = D*dt/2/dr/dr;

% vectorized sparse construction of matrices
tic
offdiags = ...
	+ sparse(2:(N-1),1:(N-2),+(dr./r(2:(end-1))-1)*g , N,N) ...
	+ sparse(2:(N-1),3:(N  ),-(dr./r(2:(end-1))+1)*g , N,N);

ML = sparse(2:(N-1),2:(N-1), ones(N-2,1) + 2*g , N,N) + offdiags;
MR = sparse(2:(N-1),2:(N-1), ones(N-2,1) - 2*g , N,N) - offdiags;

% boundary conditions
ML(1,  1:3  ) = [-3/2 ,  2 , -1/2]/dr;
ML(N,(N-2):N) = [ 1/2 , -2 ,  3/2]/dr;


u = zeros(N,1);

figure; hold on
for k = 1:(numel(tspan)-1)
	up = u;
	
	d = MR*up;
	d(N) = qfun(tspan(k+1),q0,tc)/D;
	
	u = ML\d;
	plot(r,u)	
end

%% transient flux function

function q = qfun(t,q0,tc)
q = q0;
if t >= tc
	q = 0;
end
end

[Matthew Hunt]

I have no idea what you're doing here. It might be what I need or it might not but unless you get rid of the clever code, I can't say for sure.

[Torsten]

Alex invested quite a lot of time to understand your code.

I think now it's only fair that you invest a little time to understand his improvements.

Answer 3

[J. Alex Lee]

For one thing, your boundary conditions (as developed on your pdf) don't look right.

At r=0, you can just directly use Eq 8 LHS as your equation. In the scheme

A*u_(i+1) = B*u_(i) + c

this means your first row in B will be zeroed out, and the first row in A will just be your FD coefficients. You have j=-1 and j=+1, but that is confusing since there is no -1th node...you can either use 0 and 1 (1st and 2nd columns), or use a 3 point forward FD, with coefficients -3/2, 2, -1/2.

[Matthew Hunt]

As I said, you take the limit as r->0 to get the equation necessary.

[J. Alex Lee]

take it or leave it, man.

[Matthew Hunt]

It's the usual way that people get rid of the problem at r=0. I think your method will not implement the actual inner boundary condition. This is, I believe well known.