Friday, March 25, 2011

CFL Condition (Stability in Numerical PDE)

My Traveling to Math Canyon Face
Grab y'all reins and dig in your spurs, were taking a trip down to Math Canyon!

For an lab I had to write a numerical partial differential equation solver for a parabolic equation - basically solve the heat transfer / diffusion equation. Naturally being an inquisitive young man I started to look into the how the error of the solution compare to the step sizes - both in the temporal and spatial domains.  What I observed (Figure __) seemed to be very incorrect - as the number of steps per unit length (on the same interval) increased, the error also increased!

Figure 1 - max error vs x nodes
Let's give that a moment; it should be relished like a fine wine? - the error increased with decreasing step size.  Unless you are encountering round off error from finite machine precision (which would be orders of magnitude lower than what I observed), it is always taught that decreasing the step size increases the accuracy.  This was very puzzling, and I spent a long  time trying to debug my PDE solver - but I found no bug!


Solution with 8 spatial steps


Solution with 64 spatial steps

In looking at the plots of the error between 8 spatial steps and 64 spatial steps I noticed that the magnitude of the error was coming from the boundary conditions, so it could not be due to round off.
  
For  a parabolic PDE signals propagate with infinite speed (actually the speed of light, but diffusion is much
slower than c so the information propagation is essentially infinite). Adding xSteps equates to more computations for a signal to transverse cells as a the conservation law on which the numerical solution is based is driven by the flux (concentration gradient) - line  $$u(j,n) = u(j,n-1) + dt*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1)); $$ in the PDE routine.

Suppose that the x-domain of the PDE is the interval $$[a b]$$, discretized into N intervals of width $$dx = (b-a)/N$$.  Further suppose that a signal (information, i.e. a heat source placed on the boundary) is incident on boundary $$x=a$$.  How many time iterations will it take to cross to the right boundary, $$x=b$$?

From $$u(j,n) = u(j,n-1) + dt*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1)); $$ it is apparent that some amount of the signal will propagate to b - it will be scaled by what is already there, the time step inteval dt, and nodes on the left and right.  For convince, however, let us assume this signal occurs at t=0 and the intial conditions are $$u(x,t=0)=0$$.  Below I calculate the order of $$dt$$ that the term depends on, where n  (rows) is the time vector, and i (columns) is the spatial vector.

Finite Propagation Speed
It is now Flannel Friday, and I need  to catch up on Thirsty Thursday, but why for n=5 are most of the terms proportional to $$dt^7$$?  Does a pattern emerge?

If you really want your mind blown, wear a pair of headphones through a library door and hear the magnets!


1 comment: