Crank Nicholson method for cylindrical co-ordinates

I am trying to solve the heat equation in cylindrical co-ordinates using the Crank-Nicholson method, the basic equation along with boundary/initial conditions are:

The numerical algorithm is contained in the document. The code I wrote for it is:

%This solves the heat equation with source term in cylindrical
%co-ordinates.

%T_t=(k/r)*(rT_r)_r+Q(t,r)
N_r=300; %Radial steps
N_t=200; %Time steps

r=linspace(0,1,N_r); %Radial co-ordinate
t=linspace(0,2,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments

%Solution matrix
T=zeros(N_t,N_r); 

%Physical characteristics
k=0.01; %Thermal conductivity
h=0.01; %Thermal exchange constant
theta=0.1;
A=k*dt/(2*dr*2); B=k*dt/(2*dr); %Useful constants
Q=-0.1*t.^2.*exp(-t);
%Construct the matrices

%(i+1)th matrix
a_1=-(theta+2*A)*ones(N_r,1);
a_1(1)=-(4*A-theta);
a_1(N_r)=-(2*h*dr*(A+B/r(N_r))+2*A+theta);
a_2=A+B./r(1:N_r-1);
a_2(1)=4*A;
a_3=A-B./r(2:N_r);
a_3(N_r-1)=2*A;
alpha=diag(a_1)+diag(a_2,1)+diag(a_3,-1);

%ith matrix
b_1=(2*A-theta)*ones(N_r,1);
b_1(1)=4*A-theta;
b_1(N_r)=2*A-theta+2*h*dr*(A+B/r(N_r));
b_2=-a_2;
b_3=-a_3;
beta=diag(b_1)+diag(b_2,1)+diag(b_3,-1);

%Constructing the solution:
b=zeros(1,N_r);

for i=2:N_t
    b=-0.5*(Q(i,:)'+Q(i-1,:)')*dt;
    b(N_r)=2*dr*(Q(i-1)+Q(i))*(A+B/r(N_r));
    u=T(i-1,:)';
    C=beta*T(i-1,:)'+b';
    x=alpha\C;
    T(i,:)=x';
end

for i=1:N_t
    plot(r,T(i,:));
    pause(0.1);
end

There is a great deal of oscillations in the region of r=0 and I'm not sure how to get rid of it. Any suggestions?

[J. Alex Lee]

are you oscillations in time or space or both? What happens if you refine your r-steps toward r=0?

[Bill Greene]

I believe your BC at r=0 is wrong. dT/dr should equal zero there. Handling this BC is somewhat tricky. Equation 2.61 in Morton and Mayers shows how to handle this.

As an aside, is there a reason you don't simply use the pdepe function to solve this equation?

[Matthew Hunt]

Bill, If you look at the attached pdf, you'll see that is indeed the BC I have at t=0. To handle it, look at equation 11 in the attached pdf. I don't use pdepe because I want to make this a function to use in and optimiser.

J. Alex Lee, this happens in space. When I solved the heat equation in spherical co-ordinates, I needed no refinement of r-steps.

[J. Alex Lee]

Maybe there is something funky about your BCs...your final condition seems to imply that you want the temperature difference with the outside where outside has 0 temperature, but your intial temperature is 0. Don't you want to either start with a non-zero temperature field, or modify your outside boundary condition as

or something like that, where T_0 is a non-zero constant?

It doesn't solve your problem, but maybe helpful to understand the steady state solution, and see if your FD scheme will suffer oscillations with the steady state solution.

[Matthew Hunt]

There should be a heat source term in the equation.