- #1

- 1,796

- 33

- TL;DR Summary
- I am solving the spherical diffusion equation in 1D using CN and for some reason it diverges.

I'm trying to solve the diffusion equation in spherical co-ordinates with spherical symmetry. I have included the PDE in question and the scheme I'm using and although it works, it diverges which I don't understand as Crank-Nicholson should be unconditionally stable for the diffusion. The code I'm using is:

Any suggestions are welcome,

Matlab:

```
%This program solves the spherical diffusion equation with constant
%diffusion coefficient, it has boundary conditions c_r(t,0)=0 and
%c_r(t,R)=J.
%The equation is T_t=(D/r^2)(r^2T_r)_r
clear;
N_r=50;N_t=1000; %N_r is the number of points in the electrodes, N_t is the
D=0.01; %This is the diffusion coefficient.
r=linspace(0,1,N_r);
t=linspace(0,5,N_t);% set time and distance scales
dr=r(2)-r(1);dt=t(2)-t(1); %Increments
s_1=D*dt/(2*dr^2);%Useful for the solution.
s_2=D*dt/(4*dr);
%The main problem is the matrix problem Au=v, the matrix A is constant, so
%we do this before the main loop.
u=zeros(N_t,N_r); %set up solution matrix
a=zeros(N_r-3,1);b=zeros(N_r-2,1);c=zeros(N_r-3,1);
for i=1:N_r-3
a(i)=(s_1+s_2/r(i+1));
b(i)=-(1+2*s_1);
c(i)=(s_1-s_2/r(i+1));
end
b(end)=-(1+2*s_1);
%The beginning and end rows are different to the others due to the boundary
%conditions
A=diag(a,1)+diag(c,-1)+diag(b);
%Add in the corrections for the boundary conditions
A(1,1)=(2*s_1+4*s_2/r(2))/3;A(1,2)=-(3+2*s_1+4*s_2/r(2))/3; %Boundary conditions ar r=0;
A(end,end-1)=(2*s_1-4*s_2/r(end-1))/3;A(end,end)=-(3+2*s_1-4*s_2/r(end-1))/3; %Boundary conditions at r=R
%Add in intial condition:
u(1,:)=0;
%This is the RHS for the linear system of equations.
v=zeros(N_r-2,1);
%This is the boundary condition on r=R
J=0.02*exp(0.5*t);
%Solve the system
for i=2:N_t
%Add in boundary conditions
u(i,1)=(4*u(i,2)-u(i,3))/3;
u(i,end)=(J(i)+4*u(i,end-1)-u(i,end-2))/3;
%Construct the RHS for the solution
for j=2:N_r-1
A_1=-(s_1+s_2/r(j));
A_2=(2*s_1-1);
A_3=(s_1-s_2/r(j));
v(j-1)=A_1*u(i-1,j+1)+A_2*u(i-1,j)+A_3*u(i-1,j-1);
end
v(end)=-(s_1+s_2/r(end-1))*u(i-1,end)+(2*s_1-1)*u(i-1,end-1)-(s_1-s_2/r(end-1))*u(i-1,end-2)-2*dr*(s_1+s_2/r(end-1))*J(i)/3;
%Solve for the new time step
w=A\v;
u(i,2:N_r-1)=w';
end
figure
for i=1:N_t
plot(r,u(i,:));
pause(0.1);
end
```

Any suggestions are welcome,

#### Attachments

Last edited by a moderator: