Crank-Nicholson solution of 1D heat equation

  • #1
hunt_mat
Homework Helper
1,741
25
I wish to numerically compute solutions of the 1D heat equation using the Crank-Nicholson scheme:

The equation is:
[tex]
\partial_{t}u=\partial^{2}_{x}u
[/tex]

I use the discretisation:
[tex]
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})
[/tex]

Where [itex]s=\delta t/2\delta x[/itex]. This can be rearranged into the following:
[tex]
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}
[/tex]

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 [itex]\partial_{x}u(t,0)=\partial_{x}u(t,L)=0[/itex] how would I go about doing it? I've heard about people going from an [itex]n\times n[/itex] matrix to an [itex](n-2)\times (n-2)[/itex], 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:

Code:
%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
 

Answers and Replies

  • #2
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.
 
  • #3
hunt_mat
Homework Helper
1,741
25
It shouldn't matter with the transient problem as we're just dealing with boundary conditions. so you're suggesting I use something like:

[tex]
\frac{\partial u}{\partial x}\approx\frac{u_{i,j+1}-u_{i,j-1}}{2\delta x}
[/tex]

That's the accuracy you're referring to right?
 
  • #4
hunt_mat
Homework Helper
1,741
25
The code is now working and I have also solved the following extension of the simple heat equation:

[tex]
\rho c\frac{\partial T}{\partial t}=\frac{\partial}{\partial x}\left(k(x)\frac{\partial T}{\partial x}\right)+\dot{Q}
[/tex]
 

Related Threads on Crank-Nicholson solution of 1D heat equation

Replies
1
Views
2K
Replies
0
Views
4K
Replies
1
Views
2K
  • Last Post
Replies
5
Views
3K
  • Last Post
Replies
5
Views
10K
Replies
0
Views
7K
Replies
1
Views
1K
Replies
0
Views
11K
Top