- #1

- 1,798

- 33

- 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:

and I use a wrapper function:

Am I making some trivial error?

Code:

```
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:

Code:

```
%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?