hunt_mat
Homework Helper
 1,686
 12
 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 coordinates 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 CrankNicholson 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_r3,1);b=zeros(N_r2,1);c=zeros(N_r3,1);
for i=1:N_r3
a(i)=(s_1+s_2/r(i+1));
b(i)=(1+2*s_1);
c(i)=(s_1s_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,end1)=(2*s_14*s_2/r(end1))/3;A(end,end)=(3+2*s_14*s_2/r(end1))/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_r2,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,end1)u(i,end2))/3;
%Construct the RHS for the solution
for j=2:N_r1
A_1=(s_1+s_2/r(j));
A_2=(2*s_11);
A_3=(s_1s_2/r(j));
v(j1)=A_1*u(i1,j+1)+A_2*u(i1,j)+A_3*u(i1,j1);
end
v(end)=(s_1+s_2/r(end1))*u(i1,end)+(2*s_11)*u(i1,end1)(s_1s_2/r(end1))*u(i1,end2)2*dr*(s_1+s_2/r(end1))*J(i)/3;
%Solve for the new time step
w=A\v;
u(i,2:N_r1)=w';
end
figure
for i=1:N_t
plot(r,u(i,:));
pause(0.1);
end
Attachments

109.9 KB Views: 32
Last edited by a moderator: