Crank-Nicholson solution of 1D heat equation

[hunt_mat]

I wish to numerically compute solutions of the 1D heat equation using the Crank-Nicholson scheme:

The equation is:
<br /> \partial_{t}u=\partial^{2}_{x}u<br />

I use the discretisation:
<br /> u_{i+1,j}-u_{i,j}=s(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1})+s(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1})<br />

Where s=\delta t/2\delta x. This can be rearranged into the following:
<br /> su_{i+1,j+1}-(2s+1)u_{i+1,j}+su_{i+1,j-1}=-su_{i,j+1}+(2s-1)u_{i,j}-su_{i,j-1}<br />

This is a matrix equation which is tridiagonal and can be solved very easily using matlab's many inbuilt functions.

So suppose I wish to use the following boundary conditions \partial_{x}u(t,0)=\partial_{x}u(t,L)=0 how would I go about doing it? I've heard about people going from an n\times n matrix to an (n-2)\times (n-2), solving that and just adding in the boundary conditions after that? I tried that method but my solution blew up annoyingly.

Any suggestions on what I am doing wrong? My code is here:

%This program is meant to test out the Crank-Nicolson scheme using a simple
%nonlinear diffusion scheme.
n=2000;m=30;
t=linspace(0,20,n);% set time and distance scales
x=linspace(0,1,m);
dx=x(2)-x(1);dt=t(2)-t(1); %Increments
s=dt/(2*dx^2);%Useful for the solution.
u=zeros(n,m); %set up solution matrix
p=s*ones(1,m-1); q=-(1+2*s)*ones(1,m);
Z=diag(p,1)+diag(p,-1)+diag(q);
A=Z(2:m-1,2:m-1);
%Add in intial condition:
u(1,:)=exp(-5*(x-0.5).^2);
v=zeros(m-2,1);
%Solve the system
for i=2:n-1
    %Construct the RHS for the solution
    for j=2:m-1
        v(j-1)=s*u(i-1,j+1)-(2*s+1)*u(i-1,j)-s*u(i-1,j-1);
    end
    %Solve for the new time step
    w=A\v;
    u(i,2:m-1)=w;
    u(i,1)=u(i,2);
    u(i,end)=u(i,end-1);
end

 

[Aladin Choura]

A common fall is to use first order finite difference for the derivative boundary condition. Use a 3 points formula instead. I don't know how that can be embedded in the solution for the transient problem though.

 

[hunt_mat]

It shouldn't matter with the transient problem as we're just dealing with boundary conditions. so you're suggesting I use something like:

<br /> \frac{\partial u}{\partial x}\approx\frac{u_{i,j+1}-u_{i,j-1}}{2\delta x}<br />

That's the accuracy you're referring to right?

 

[hunt_mat]

The code is now working and I have also solved the following extension of the simple heat equation:

<br /> \rho c\frac{\partial T}{\partial t}=\frac{\partial}{\partial x}\left(k(x)\frac{\partial T}{\partial x}\right)+\dot{Q}<br />