Finite Difference Heat Equation using NumPy

The problem we are solving is the heat equation

\frac{ \partial u}{\partial t} = \frac{ \partial^2 u}{\partial x^2}

with Dirichlet Boundary Conditions ( u(0,t) = u(1,t) = 0 ) over the domain 0 \le x \le 1 with the initial conditions

u(x,0) = 10~sin( \pi x )

You can think of the problem as solving for the temperature in a one-dimensional metal rod when the ends of the rod is kept at 0 degrees. Intuitively, you know that the temperature is going to go to zero as time goes to infinite.

To solve this problem using a finite difference method, we need to discretize in space first. I will be using a second-order centered difference to approximate \frac{ \partial^2 u}{\partial x^2} . This will give the following semi-discrete problem:

\frac{ \partial U}{\partial t} = A~\vec{U}

where

A = \frac{1}{\Delta x^2} \left( \begin{array}{cccccccc}-2 & 1 & 0 & 0 & \ldots & 0 & 0 & 0 \\1 & -2 & 1 & 0 & \ldots & 0 & 0 & 0 \\0 & 1 & -2 & 1 & \ldots & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & \ldots & 1 & -2 & 1 \\0 & 0 & 0 & 0 & \ldots & 0 & 1 & -2 \end{array} \right)

and

\vec{U}= \left( \begin{array}{c} u(0,t) \\ u(\Delta x,t) \\ u( 2 \Delta x, t) \\ \vdots \\ u(1-\Delta x,t) \\ u(1,t) \end{array} \right)

The next step is to discretize in time. We can do this by using the Crank-Nicolson method which is

\frac{\partial y}{\partial t} \approx \frac{ u_{n+1}-u_n }{\Delta t} = \frac{1}{2} ( f_{n+1}(x) + f_n(x) )

where

\frac{\partial y}{\partial t} = f(x)

If we apply this method to the semi-discrete problem, we will get

\frac{ \vec{U}_{n+1}-\vec{U}_n }{\Delta t} = \frac{1}{2} ( A~\vec{U}_{n+1} + A~\vec{U}_n)

If we want to solve for \vec{U}_{n+1}, we get the following system of equations

(I -\frac{\Delta t}{2}~A)~\vec{U}_{n+1} = ( I + \frac{\Delta t}{2}~A)~\vec{U}_n

We can implement this method using the following python code. This code will then generate the following movie.

Please note that the movie was sped up by a factor of ten.