Crank-Nicholson method for spherical diffusion

[hunt_mat]

TL;DR Summary

Code for Crank-Nicholson method is highly dependent upon discretisation

The code I have for solving the diffusion equation on the spherical domain. The solution seems to differ drastically depending on the refinement of the mesh which obviously shouldn't be the case. I have checked and double checked my derivation and code and I can't seem to find an error. One thing is for sure, the fineness of the points I use shouldn't make the solution change wildly, and yet it does. I'm at a complete loss as to what is going on. The main function I use to code up is:

function u=spherical_fd(u_0,q,D,r,t)
global X;
n=length(r);m=length(t);  %The dimensions of u are defined by the length os r and t.
u=zeros(m,n); %Initialise u
dr=r(2); %as both t and r start at zero, dr and dt can be defined in this way
dt=t(2);
alpha=D*dt/dr^2; %compute alpha and beta
beta=D*dt/(2*dr);
A_main=-(1+alpha)*ones(n,1); %Begin constructing the vectors to put in the diagonals
A_u=0.5*alpha+beta./r(1:n-1);
A_l=0.5*alpha-beta./r(2:n);
A=diag(A_u,1)+diag(A_l,-1)+diag(A_main,0); %Construct A
B_main=(1-alpha)*ones(n,1); %Only the main diagonal for differes by anyting more than a sign
B=diag(A_u,1)+diag(B_main,0)+diag(A_l,-1); %Construct B
A(1,1)=1-3*alpha; %Insert the inner boundary values for A and B
A(1,2)=3*alpha;
B(1,1)=1+3*alpha;
B(1,2)=3*alpha;
A(n,n)=-(1+alpha); %Insert the outer boundary values for A and B
A(n,n-1)=alpha;
B(n,n)=-alpha;
B(n,n-1)=alpha-1;
p=zeros(n,1); %Used in the consutruction of c
p(n)=1;
B_1=A\B;
w=A\p;
gamma=-2*(dr/D)*(0.5*alpha+beta/r(n)); %For computational efficency.
u(1,:)=u_0;
for i=2:m
c=gamma*w*(q(i-1)+q(i));
sol=-B_1*u(i-1,:)'+c;
u(i,:)=sol';
X(i)=sum((r'.^2).*(A*(u(i,:)'))+(r'.^2).*(B*(u(i-1,:)')));
end

and I use a wrapper function:

%This is a wrapper to test the main spherical diffusion.
global X
clear;clc;
N=800;
M=400;
r=linspace(0,1,N);
t=linspace(0,20,M);
X=zeros(M,1);
u_av=zeros(M,1);
u_0=0.95*ones(length(r),1);
u_0_av=trapz(r,(r.^2).*u_0');
D=0.0001;
q=0.15*ones(M,1);
M_1=floor(M/5);
q(M_1:end)=0;
y=u_0_av+0.3*t;
y(M_1+1:end)=y(M_1);
u=spherical_fd(u_0,q,D,r,t);
for i=1:M
    v=(r.^2).*u(i,:);
    u_av(i)=trapz(r,v);
end
u_bdy=u(:,end);
plot(t,u_av);
hold on
plot(t,y);
xlabel('Time');
ylabel('Average Value')
axis([0 18 0 2])

Am I making some trivial error?

 

[pasmith]

I'm not sure you are implementing the boundary conditions correctly.

At $r = 0$ you want $\frac{\partial u}{\partial r} = 0.$ You are using a central difference, so that means that $u_{i,-1} = u_{i,1}$. So at $r = 0$ you should have (a single subscript on $u$ is a spatial index) $$\frac{d u_0}{d t} = 3D\frac{2u_1 - 2u_0}{\delta r^2} = \frac{6\alpha}{\delta t}(u_1 - u_0)$$ so that $$u_{i+1,0} - u_{i,0} = 3\alpha(u_{i+1,1} - u_{i+1,0} + u_{i,1} - u_{i,0})$$ and $$(1 + 3\alpha)u_{i+1,0} - 3\alpha u_{i+1,1} = (1 - 3\alpha)u_{i,0} + 3\alpha u_{i,1}.$$ This not consistent with equation (10) of your paper; I haven't checked whether it's consistent with your code.

At $r = r_N = R$, you have $$\frac{\partial u}{\partial r} = \frac{q(t)}{\lambda}.$$ This can be subsituted directly into the second term on the right hand side of governing PDE; in the approximation of the second derivative you must eliminate $u_{N+1}$ using $$u_{N+1} = u_{N-1} + 2\delta r \frac{q(t)}{\lambda}.$$
Putting these together, I get $$\begin{split} \frac{d u_N}{dt} &= D\frac{u_{N+1} + u_{N-1} - 2u_N}{\delta r ^2} + \frac{2Dq(t)}{R\lambda} \\ &= D\frac{2u_{N-1} - 2u_N}{\delta r ^2} + 2D\frac{q(t)}{\lambda \delta r}\left(1 + \frac{\delta r}{R}\right) \\ &= \frac{2\alpha}{\delta t}(u_{N-1} - u_N) + \frac{4\beta q(t)}{\delta t}\left(1 + \frac{1}{N}\right). \end{split}$$ I'm not sure this is consistent with your equation (13).

(I think your paper would be easier to follow if you first discretise in space to obtain $$\frac{du_j}{dt} = \dots$$ and then apply Crank-Nicholson to obtain $$u_{i+1,j} - \frac{\delta t}{2} \left.\frac{du_{j}}{dt}\right|_{t=t_{i+1}} = u_{i,j} + \frac{\delta t}2 \left.\frac{du_{j}}{dt}\right|_{t=t_i}$$ as I have done above.)

 

[hunt_mat]

I agree with your correction for the inner boundary condition but I have no idea what you're doing with the outer boundary condition. Can you explain what you're doing a little more carefully? I'm the opposite of you, as what your doing is very confusing. I'm not doing the method of lines or anything like that, just simple finite differences, nothing fancy.