Finite Difference Heat Equation using NumPy

January 2nd, 2010

To expand my experience with NumPy during the Christmas holiday, I decided to solve the heat equation using a finite difference technique. I know this is a trivial problem to solve. However, like I said in my previous post, the goal is to learn NumPy and not to solve complicated mathematical problems. Likewise, the purpose of this post is not to explain how finite difference methods work, but how to implement it using Python and 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.

Comments(0) - Numerical Calculations

Creating Movies with PyPlot/PyLab

January 2nd, 2010

Over the Christmas holidays, I started to dive more into Python. I have been trying to learn more about NumPy, SciPy, and PyLab. So far, I love these libraries. However, there doesn’t seem to be a good way to create movies/animations.

The best way to create a movie using Python is to generate plots using PyLab or PyPlot and save them as pictures. You can then use another program to stitch them together (ie. mencoder or ffmpeg). This process takes substantial amount of time and seems a little bit hacky at best. Please leave comments if you know of a better way.

It took some time to figure out how to get everything configured on MacOS X. My first problem was trying to get a copy of mencoder or ffmpeg. We need one of these programs in order to stitch all the pictures together. I found an interesting post which allowed me to easily install ffmpeg onto my mac. You can download my code and a copy of ffmpeg here.

Once you downloaded the zip file, you will need to

Copy ffmpeg to /usr/local/bin
Execute the following commands from your terminal:
sudo chown root:wheel /usr/local/bin/ffmpeg
sudo chmod 755 /usr/local/bin/ffmpeg

Once you have installed ffmpeg, we can look at how to create a movie. I created a function called CreateMovie which will easily create a movie

def CreateMovie(plotter, numberOfFrames, fps=10):
	import os, sys
	import matplotlib.pyplot as plt
 
	for i in range(numberOfFrames):
		plotter(i)
		fname = '_tmp%05d.png'%i
 
		plt.savefig(fname)
		plt.clf()
 
	os.system("rm movie.mp4")
	os.system("ffmpeg -r "+str(fps)+" -b 1800 -i _tmp%05d.png movie.mp4")
	os.system("rm _tmp*.png")

This function creates a movie called movie.mp4 in the current working directory which has a default frame rate of 10 frames per second. Make sure that you don’t have any files called movie.mp4 and _tmp*.png in the current working directory because they will be deleted.

The function takes the following parameters:

plotter: This parameter is a function with the following form:

def plotter(frame_number)

where frame_number is the current frame that needs to be plotted using the matplotlib.pyplot library.

numberOfFrames: The total number of frames in the movie.

fps: The frames per second. The default is 10.

In order to illustrate the use of this function, check out this code for a moving gaussian.

import CreateMovie as movie
import matplotlib.pyplot as plt
import numpy as np
 
x = np.linspace(-.5,4,500)
 
# Plots a given frame
def plotFunction( frame ):
	plt.plot(x, np.exp( -10*(x - frame/10.0)**2) )
	plt.axis((-.5,4,0,1.1))
 
movie.CreateMovie(plotFunction, 50)

This code will generate the following movie.

Comments(1) - Programming

Poster Presentation on Traffic Flow

December 17th, 2009

Lighthill-Whitham Model

The continuity equation for traffic flow is defined as:

$\frac{\partial \rho}{\partial t} + \frac{\partial q}{\partial x} = 0$

where $\rho$ is the car density, is the car speed, and is the traffic flow ( $q = \rho u$ ). This equation corresponds to a one dimensional fluid-dynamic model where no vehicles will enter or leave the traffic flow except at the boundary.

Let us assume that the car’s velocity depends only on the density ( $u=u(\rho)$ ). The simplest model
would be to consider the following linear relationship:

$u(\rho) = u_{max} ( 1-\frac{\rho}{\rho_{max}})$

where $\rho_{max}$ is the traffic density in which you have bumper-to-bumper traffic and $u_{max}$ is the maximum speed that a vehicle can travel. This usually corresponds to the speed limit.

You will notice that when there are no vehicles on the road an individual vehicle would travel at $u_{max}$ . Likewise, when you have bumper-to-bumper traffic, the vehicles on the road won’t be moving.

The relationship between density and traffic flow can be shown by the Fundamental Diagram of Road Traffic which has a maximum flow at $\frac{\rho_{max}}{2}$ .

If we put our velocity model into the continuity equation, we get the following:

$\frac{\partial \rho}{\partial t} + u_{max} ( 1-2 \frac{\rho}{\rho_{max}} ) \frac{\partial \rho}{\partial x} = 0$

Shock Waves

The Lighthill-Whitham Model is a non-linear wave equation. This equation can be solved using the method of characteristic lines. One of the defining features of this equation is the creation of expansion and shock waves. Shock waves are formed when characteristic lines intersect which causes $\rho (x)$ to become discontinuous. It can be shown that shock waves will have speed

$S(\rho_+,\rho_-) = \frac{q(\rho_+)-q(\rho_-)}{\rho_+-\rho_-}$

where $\rho_+$ and $\rho_-$ is the density when approaching from the right and left of the discontinuity.

Dealing with shock waves numerically has serious difficulties. To avoid shock waves, we can add a small diffusive term to the LW Model. If we define the traffic flow as

$q = u_{max} ( 1-\frac{\rho}{\rho_{max}} ) \rho-\epsilon \frac{\partial \rho}{\partial x}$

we can get a variation of burgers equation

$\frac{\partial \rho}{\partial t} + u_{max} ( 1-2 \frac{\rho}{\rho_{max}} ) \frac{\partial \rho}{\partial x} = \epsilon \frac{\partial^2 \rho}{\partial x^2}$

Fourier-Galerkin Spectral Method

The diffusive LW model can be solved using a Fourier-Galerkin Spectral Method. If we assume the problem has periodic boundary conditions and $x \in [0, 2\pi]$ , we can expand $\rho$ using $e^{inx}$ as a basis. This will lead us to the following system of equations

$\frac{d a_n(t)}{dt} + i~n~u_{max}~a_n(t) + \epsilon~n^2~a_n(t)$ $= \frac{2~i~u_{max}}{\rho_{max}} \displaystyle\sum_{j=-N/2}^{N/2} j~a_j(t)~a_{n-j}(t)$

$a_n(0) = \begin{cases}\frac{1}{2 \pi} \displaystyle\int_{0}^{2\pi} \rho(x,0)~e^{inx}~dx & \text{for $n \in [-\frac{N}{2}, \frac{N}{2}]$} \\ 0 & \text{otherwise} \end{cases}$

where

$\rho(x,t) = \displaystyle\sum_{j=-N/2}^{N/2} a_j(t)~e^{inx}$

After a Traffic Light Turns Green

Let’s analyze the behavior of traffic after a light turns green. We can do this by using the following initial conditions:

$\rho(x,0) = \begin{cases}\rho_{max} & \text{for $0 \leq x \leq \pi$} \\0 & \text{for $\pi < x \leq 2\pi$}\end{cases}$

If we use the Fourier-Galerkin Spectral Method, we get the following results

when $\rho_{max} = 3$ , $\epsilon = .01$ , $u_{max} = .5$ , and N=100 . You will notice the Gibbs phenomenon in our initial conditions. This is because a step-function can’t be exactly represented by our Fourier basis functions. However, the diffusive term will quickly damp any highly oscillating term.

If we look at the exact solution to this problem for the Non-Diffusive LW Model, we will get an expansion wave. The length of expansion can be shown to be $2 u_{max} t$ .

If we look at the Diffusive LW Model, we can see that the expansion length is dependant on $\epsilon$ .

After a Traffic Light Turns Red

Let’s analyze the behavior of traffic after a light turns red. We can do this by using the following initial conditions:

$\rho(x,0) = \begin{cases} \rho_0 & \text{for $0 \leq x \leq \pi$} \\ \rho_{max} & \text{for $\pi < x \leq 2\pi$}\end{cases}$

If we use the Fourier-Galerkin Spectral Method, we get the following results

when $\rho_{max} = 3$ , $\rho_0 = 1$ , $\epsilon = .01$ , $u_{max} = .5$ , and N=100 . After the Gibbs phenomenon dies down from our diffusive term, you will notice a shock that moves at speed -0.1685 .

If we look at the exact solution to this problem for the Non-Diffusive LW Model, we will get a shock that moves at speed

$u = – \frac{\rho_0 ~ u_{max}}{\rho_{max}-\rho_0} ( 1-\frac{\rho_0}{\rho_{max}} )$

If we use the same parameters as the above numerics, we get a shock speed of -0.166667 .

If we look at the Diffusive LW Model, we can see that the shock speed doesn’t vary much for different values of $\epsilon$ .

References:

[1] Kopriva, D. (2009). Implementing Spectral Methods for Partial Differential Equations. Springer Science
[2] Gottlieb, D. & Orszag, S.A. (1977). Numerical analysis of spectral methods: Theory and Applications. Montpelier,Vermont: Capital City Press
[3] Haberman, R. (1977). Mathematical Models: Mechanical Vibrations, Population Dynamics, and Traffic Flow. Englewood Cliffs, NJ: Prentice-Hall
[4] Helbing, D. (2001). Traffic and related self-driven many-particle systems. Rev. Mod. Phys. 73 (4), 1067–1141.

Comments(1) - Uncategorized

Stochastic Traffic Cellular Automaton

November 17th, 2009

Lately, I have started to study Traffic Flow Models. This post will hopefully be a start of a long series of posts regarding the field. Typically, there are two types of approaches when modeling traffic flow. The first approach is to treat traffic as a fluid. In this model, the behavior of individual cars are not directly modeled. Instead of modeling cars, this approach models the behavior of traffic density using Constitutive Laws. The second approach attempts to model traffic flow using a Cellular Automaton. The basic idea is that the driver’s behavior only depends on the car(s) immediately in-front and behind him. In this post, I am going to construct a simple model to describe this later approach.

First, let us divide the road into equal segments of length $\Delta x$ and divide time likewise into equal time-steps $\Delta t$ . The road would look something like this:

At any given timestep, a car can only occupy one grid space. This space is represented by the integer . After each time-step, we will move the car according to it’s current velocity. We will define the integer as the number of steps that the car will take at each time-step. Hence, the new car position will follow the iterative formula $x_{new} = x_{old} + m_{new}$ . In order to understand Stochastic Traffic Cellular Automaton, we will also have to define two more variables. The first is which the number of steps between the car and the car in-front of it. The second variable is which corresponds to the max distance that any car can move on a time-step.

The Stochastic Traffic Cellular Automaton model will then follow the following steps for each car on the road for each time-step:

Speed-up: If $m_{old} \neq M$ , then $m_{new} = m_{old} + 1$

Do not Overrun: If $m_{new} > g_{old}$ , then $m_{new} = g_{old}$

Randomization: If $m_{new} \neq 0$ , then with probability

take $m_{new} = m_{new}-1$ .

Move Car: $x_{new} = x_{old} +m_{new}$ .

Without any traffic, a driver would speed up to a velocity of $\frac{M}{\Delta t}$ with an acceleration of $\frac{1}{\Delta t^2}$ . If there is a car in front, it makes sure that it doesn’t crash into the car (which is usually a good idea). Lastly, there is a random factor placed into the equation. This is to make the model non-deterministic. For example, a driver could over-brake or do other non-deterministic behaviors. The randomization is created to encapsulate all these possible behaviors. As you can see, the model is very basic. However, we can add further rules if we want a more accurate description of a driver.

If you want to try coding this up, I would try simulating traffic behind a red light. What happens when the light turns green?

Holmes, Mark, Introduction to the Foundations of Applied Mathematics, Chapter 5
K. Nagel and M. Schreckenberg, J. Physique I, 2, 2221 (1992), A cellular automaton model for freeway traffic

Comments(1) - Discrete Math

The God Paradox

September 13th, 2009

If God can do anything, can He make a mountain which is too heavy for Him to lift?

If God could make a mountain which he couldn’t lift, there is something that he cannot do. Likewise, if God cannot make a mountain which he couldn’t lift, then there is something which he cannot do. What does this say about the limits of an omnipotent being? Is even God subjected to the rules of logic? Or perhaps nonsense statements remain nonsense statements, even when you append ‘God’ to it.

Comments(2) - Philosophy

Interesting Books

Interesting Blogs

Interesting Articles

Finite Difference Heat Equation using NumPy

Creating Movies with PyPlot/PyLab

Poster Presentation on Traffic Flow

Lighthill-Whitham Model

Shock Waves

Fourier-Galerkin Spectral Method

After a Traffic Light Turns Green

After a Traffic Light Turns Red

Stochastic Traffic Cellular Automaton

The God Paradox