- 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,
%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
109.9 KB Views: 32
Last edited by a moderator: