Finite differencing I

From Interactive System for Ice sheet Simulation
Revision as of 12:00, 4 August 2009 by Hoffman (Talk | contribs)

Jump to: navigation, search



In one dimension, the general form of the convection diffusion equation is

\frac{\partial u(x,t) }{\partial t} - \frac{\partial}{\partial x} D(x) \frac{\partial}{\partial x} u(x,t)  - C(x)\frac{\partial}{\partial x}  u(x,t)  = S(x,t),

u is a general variable, D is a spatially-varying diffusivity, C is a spatially-varying convection rate, and S is a source term. The second term on the left represents diffusion of a solute or other material property, the third term represent convection.

This equation can be used to model a wide range of phenomena, including the distribution of temperatures (or energy conservation) in an ice sheet. It also bears similarity to the equations expressing conservation of momentum, and analysis of the numerical solutions to this equation are representative of the analysis of numerous numerical treatments in computational fluid dynamics. For these reasons, we will make convection diffusion our first "model problem", or problem to solve in order to strengthen intuition.

We will take a step-wise approach to solving this equation, first solving for the diffusion, or parabolic equation, then solving for the convection portion. Finally, we will solve the complete equation. Through this process, we will be looking at the stability of the numerical methods used to solve the equation.

You should think of this as a starting point for both your learning to program, as well as your learning to solve PDEs with programs.

Diffusion and explicit solution

First, we will solve a simplified version of the equation explicitly. Explicit here refers to the way (or what) the differentiation operators are applied to. In this situation they are directly applied to the solution at the present time step in order to determine the next time.

The Stencil for the most common explicit method for the parabolic equation.

To better understand, apply the idea to what is called the parabolic, diffusion, or sometimes heat equation. In terms of convection diffusion this is D(x,t) = 1, C(x,t) = 0 and S(x,t) = 0,

 \frac{\partial u(x,t) }{\partial t} = \frac{\partial ^2 u(x,t)}{\partial x^2},

The finite difference approximation of the equation is

 \frac{u(x,t+\Delta t) - u(x,t)}{\Delta t} = \frac{u(x+\Delta x,t) - 2u(x,t) + u(x-\Delta x,t)}{\Delta x^2}.

Where both derivative approximations are known from the previous lesson. One is called the 'forward Euler' approximation of the time derivative, and the other is the second order accurate, centered second derivative.

The equation is then algebraically solved for  u(x,t)

u(x,t + \Delta t) = u(x,t) + \Delta t \frac{u(x+\Delta x,t) - 2u(x,t) +
u(x-\Delta x,t)}{\Delta x^2}

There you have it, a way to compute the future, using the present. Your first task will be to change this into an algorithm.

One final note, the stencil at the left makes great sense, and understanding it will help make other algorithms clear. However, to understand it we must modify the notation a little. Make super-scripts n refer to time and subscripts j refer to space. The previous equation becomes

u_j^{n+1} = u_j^n + \frac{\Delta t}{\Delta x ^2} \left( u_{j-1}^{n} - 2 u_j^{n} + u_{j+1}^{n} \right).

Make sure you recognize how this corresponds with the diagram of the explicit stencil.

Numerical solution

One can see from the above equation that one way to numerically solve the parabolic equation is to use stencils or operators for computing second derivatives. Once the derivative is computed, finding the solution corresponding to the next time step, u(x,t+\Delta t) is just a matter of multiplying the derivative by \Delta t and adding  u(x,t).

In psuedocode, the solution looks something like

Initialize variables
loop t over time:
  loop i over space:
     u(t,i) = u(t-1,i) +  delta_t * (u_old(t-1,i-1) - 2*u_old(t-1,i) + u_old(t-1,i+1) ) / delta_x**2
     store solution as needed
  end loop over space
end loop over time

Note that in psuedo-code the final result is an T by N array (T is time steps, N space points). This is typical, the data structures used to store output is often as complex as the algorithm. Sometimes more.

You'll need to do this in fortran 90, and you'll find Gethin's Pragmatic Programming very helpful. Plotting is also a key to understanding simulation output. Again, Gethin's got what you need, although you'll also find the matplotlib documentation helpful (should be familiar if you've used Matlab).

Plotting results

The above line

      store solution as needed

is perhaps, frustratingly vague. Gethin's Fortran1: Fortran for beginners discusses input/output, IO, so the writing should not be a problem. Even

  write(*,*) u

will suffice, if you can re-direct program output to a file. Assuming your program is called cd_prg, this is done with

>>.\cd_prg >> data

where the >> operator will replace whatever was previously in the file data.

As for reading it into Python and plotting it, the main challenges are reading in the data and animating the time series data. Here is some code to do that

!/usr/bin/env python
# Import only what is needed
from numpy import loadtxt,shape,linspace
from pylab import plot,show,clf,show,ion
# Import data file, called 'data'
#Determine how much data came in
dims = shape(d)
clf()    # Clears the screen
ion()    # Interactive plot mode, critical for animation
# x data, note that this must correspond to program's domain
x = linspace(0,1,dims[1])  
# Initial plot, very Matlab(ish), note return of plot handle that allows plot to
# be altered elsewhere in code.
ph,=plot(x,d[0,:],'k')           # matplot lib requires show to be called
# Loop to plot each time step
for i in range(1,dims[0]):
    ph.set_ydata(d[i,:])   # Only update y data (faster than replot)


  1. Using the algorithm for the 'explicit' method, find a numerical solution to this heat conduction problem:
 \frac{\partial u(x,t) }{\partial t} = \frac{\partial ^2 u(x,t)}{\partial x^2},
 u(x,0) = sin(\pi x)
 u(0,t) = u(1,t) = 0

Use \Delta x = 0.1, \Delta t = 0.005125, and t_{end} = 1.025. Compare the computed solution to the exact solution  u(x,t) = \exp(-\pi^2 t) \sin(\pi x). Repeat the experiment with \Delta t = 0.006 and t_{end} = 1.026.

Group one, parabolic, explicit

Group six, parabolic, explicit

Convection and numerical stability

For this unit consider the first-order hyperbolic PDE

 \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}.

Mathematically, this statement is saying that a quantity u(x,t) exists on some grid, and is being carried along by a wind with velocity v. Before applying finite difference operators, clean up the notation so that super-scripts (u^n) refer to time and subscripts (u_j) refer to space.

Continuing to work with explicit schemes, the machinery of descritization allows us to quickly move to the form

 \frac{u_j^{n+1} - u_j^n}{\Delta t} = -v \left( \frac{u^n_{j+1} - u_{j-1}^n}{2\Delta x} \right )

and solve to give a recurrence relation

 u_j^{n+1} =u_j^n  - \frac{v \Delta t}{2\Delta x} \left(u^n_{j+1} - u_{j-1}^n \right )

von Nuemann Stability Analysis

Before implementing this, consider the stability of the solutions by assuming a very generic form of solution

 u_j^n = A(k)^n e^{ijk}.

Maybe you can recall a course in differential equations where you spent the better part of a semester making similar substitutions into equations to find solutions? This complex exponential is the Swiss Army knife of functions, and satisfies many equations.

In our assumed solutions the amplitude is A^n(k) (exponentiated to higher powers with time) and the wave number is k = \frac{2 \pi}{\lambda}. Said in words, we assume that the solution will be oscillatory (recall  e^{ikx} = cos(kx) + i sin(kx)) and that the solution's amplitude will depend on the frequency, or k. In our discrete case j serves as a proxy for space, x.

 A^{n+1} e^{ijk} = A^{n} e^{ijk}  - \frac{v \Delta t}{2\Delta x} \left( A^{n} e^{i(j+1)k} - A^{n} e^{i(j-1)k} \right )

divide through by  A^n e^{ijk}

 A = 1 - \frac{v \Delta t}{2\Delta x} (e^{ik} - e^{-ik}) = 1 - i\frac{v \Delta t}{\Delta x} sin k

What does that mean? If |A|^2>1, the solution grows without bound in time, because each time step applies an higher exponent to  A(k)^n. So, this solution scheme is unstable for all time steps, and space steps. Bummer.

Getting stability

Now, let's try that again, and when discretizing do what will become a favorite trick, average or smear the values of the function. A new discretization will be

 u_j^{n+1} = \frac{1}{2} (u_{j+1}^n + u_{j-1}^n) - \frac{v \Delta t}{2\Delta x} \left(u^n_{j+1} - u_{j-1}^n \right ),

this is called the Lax method. Now consider stability in the same way. Omitting some algebra

 A = cos k - i \frac{v \Delta t}{2\Delta x} sin k,

is the amplitude. Requiring

 |A|^2 = cos^2 k +  \left(\frac{v \Delta t}{2\Delta x}\right)^2 sin^2 k \leq 1,

to avoid unbound growth, yields

 \frac{|v|\Delta t}{\Delta x} \leq 1.

This is called the Courant-Friedrichs-Levy stability criterion. It states that the information on a grid has a velocity of  \frac{\Delta x}{\Delta t} and that the velocity in the system can not be exceeded by it (causing the ratio to exceed one). For such a thing to happen would be completely unphysical. Consider what happens when an object exceeds the velocity of waves in the media that carries it, a sonic boom. This is a "numerical boom".


  1. Implement the Lax method for a linear system. Is this method explicit or implicit? Use a 10 unit domain and begin with a height 1.0 square wave between 4.5 \leq x \leq 5.5. Fix the ends at 0. let v be 1.0. Also track the sum of the solution before and after the simulation. End the simulation after 4 seconds. Report the behavior with and without the CFL being satisfied. If the CFL is very small, do things improve. What about when it's just under 1.0. What's going on here. Try subtracting u^j_n from both sides of the discretization and inspect for differences between the original discretization and Lax. See an extra term?
  2. Try the leapfrog method for descretization
u^{n+1}_j = u^{n-1}_j - \frac{v\Delta t}{\Delta x} (u^n_{j+1} - u^n_{j-1}).
Does this improve the numerical diffusion?

Final program and exercise

Bring the descretization schemes for diffusion and convection to solve the convection-diffusion equation explicitly, with finite differences. Consider a non-dimensional form of the equation.

\frac{\partial \phi}{\partial t} + u \phi -  \frac{\partial }{\partial x}\mathrm{Pe}^{-1} \frac{\partial }{\partial x} \phi = q

The \mathrm{Pe}^{-1} , the inverse of the Peclet number, the ratio of the velocity scale U times the length scale L to the diffusivity D,

 \mathrm{Pe} = \frac{UL}{D}.
  • On a unit domain, specify \phi_(t,0)= .2, \phi_(t,1)=1, and Pe = 10.
  • Compare the solution to c_a(x), the analytic solution of the equation, is
 c_a(x) = a + (b-a)*\frac{\exp((x-1)\mathrm{Pe}) - \exp(\mathrm{-Pe})}{1-\exp(\mathrm{-Pe})}
  • Experiment with the Peclet number and mesh resolution to determine how stable your numerical scheme is.