Burgers' Equation on a Moving Grid

The familiar inviscid Burgers' equation on a non-moving, physical-space coordinate system:

Now we want to transform the equation, which is easily done in Maxima:

depends([u, x], [t]);
depends([xi], [x]);
depends([u], [xi]);

burgers : 'diff(u, t) + u * 'diff(u, x);

ev(burgers, nouns);

This gives us Burgers' equation transformed to computational space and a moving grid:

Now we've got some extra terms. The grid velocity:

And the grid transformation (to go from the computational space to the physical space):

The grid velocity is the troublesome term, because while we have an equation governing the time evolution of the unknown, we don't have a governing equation for the grid. This is where moving grids get fun, we get to pick the governing equation. How about a forced, unsteady, linear diffusion equation?

The method of lines is a popular and practical way to formulate PDEs for numerical solution. In Maxima, it is very straightforward to substitute in some finite difference expressions for the spatial derivatives, to get our semi-discrete equation for the time-derivatives.

/* substitute central differences for the spatial derivatives: */
burgers_discrete_space : ratsubst(u[i+1] - u[i-1], 'diff(u, xi), burgers_trans);
burgers_discrete_space : ratsubst(1/(x[i+1] - x[i-1]), 'diff(xi, x), burgers_discrete_space);

/* governing equation for the grid (forced diffusion equation): */
grid : 'diff(x,t) + 'diff(x, xi, 2) = f;
/* substitute central differences for the spatial derivative: */
grid_discrete_space : ratsubst(x[i+1] - 2*x[i] + x[i-1], 'diff(x,xi,2), grid);

The reason for transforming to computational coordinates is because the nodes are (usually) equally spaced in the computational domain. That way we can use the same difference formula for every point no matter what the spacing is in the physical space.

Now we can use the augcoefmatrix() function to get an operator and right-hand-side vector for our rate equation.

Rate_Operator : augcoefmatrix([burgers_discrete_space, grid_discrete_space], ['diff(u,t), 'diff(x,t)]);
Rate_RHS : -submatrix(Rate_Operator, 1, 2);
Rate_RHS : ratsubst(u[i], u, Rate_RHS);
Rate_Operator : submatrix(Rate_Operator, 3);
Rate_Operator_inv : invert(Rate_Operator);
Rate_Vector : transpose(['diff(u,t), 'diff(x,t)]);
Rate_expression : fullratsimp(Rate_Operator_inv.Rate_RHS);

Here's the resulting semi-discrete equations:

Because our choice of governing equation for the grid doesn't depend on the time-derivative of the solution, the operator is easily inverted:

and then applied, to give a vector of nonlinear expressions for our rates:

Note that in the absence of any forcing, the steady state grid is equally distributed grid points. The forcing function is what allows us to cluster grid points where higher resolution is needed. This forcing function should be based on some measure of the local error in the numerical solution. It is convenient just to use magnitude of the solution derivative since we've already calculated it, and areas with high gradients tend to need more points to resolve well.

Equal spacing in computational space is standard for using finite differences, but what about spectral methods where we usually choose nodes according to the Chebyshev polynomial roots? In this case it is convenient to choose our computational spacing to be the roots spacing rather than equal spacing. Here's the derivative of the coordinate with respect to computational coordinate using the DCT derivative routine.

Which shows that the roots nodes give a constant grid transform while the equally spaced nodes have a non-constant transform (opposite the normal situation of equally spaced computational coordinates).

PowerPoint Engineering Considered Harmful

I have always been inspired by Edward Tufte's work. As an engineer in the aerospace industry, successfully working in large organizations (read bureaucracies) is a necessary skill. Unfortunately, bureaucracy has its own particular pathology. The trick is to not catch the diseases yourself while at the same time thriving in close proximity. In that vein, I have lately been considering the Cognitive Style of PowerPoint, which is a popular tool in the industry.

Everybody likes to beat up on NASA, and there's plenty to find fault with (as there is in any bureaucracy). Honest past criticism from a much admired physicist:

It appears that there are enormous differences of opinion as to the probability of a failure with loss of vehicle and of human life. The estimates range from roughly 1 in 100 to 1 in 100,000. The higher figures come from the working engineers, and the very low figures from management. What are the causes and consequences of this lack of agreement?
...
Let us make recommendations to ensure that NASA officials deal in a world of reality in understanding technological weaknesses and imperfections well enough to be actively trying to eliminate them. They must live in reality in comparing the costs and utility of the Shuttle to other methods of entering space. And they must be realistic in making contracts, in estimating costs, and the difficulty of the projects. Only realistic flight schedules should be proposed, schedules that have a reasonable chance of being met. If in this way the government would not support them, then so be it. NASA owes it to the citizens from whom it asks support to be frank, honest, and informative, so that these citizens can make the wisest decisions for the use of their limited resources.

For a successful technology, reality must take precedence over public relations, for nature cannot be fooled.
-- Rogers Commission, Appendix F, by R. P. Feynman

The same "debt of honesty" that Feynman claims NASA owes the public can be reasonably attributed to any organization towards their stakeholders. Pretending the problems experienced by this particular organization are isolated would be foolish. NASA does provide an accessible case study of an organization being firmly infected by the pitch culture though. Others might put it as the triumph of the personality-ethic over the character-ethic, style over substance.

From a more recent report precipitated by disaster:

The Mission Management Team Chairʼs position in the hierarchy governed what information she would or would not receive. Information was lost as it traveled up the hierarchy. A demoralized Debris Assessment Team did not include a slide about the need for better imagery in their presentation to the Mission Evaluation Room. Their presentation included the Crater analysis, which they reported as incomplete and uncertain. However, the Mission Evaluation Room manager perceived the Boeing analysis as rigorous and quantitative. The choice of headings, arrangement of information, and size of bullets on the key chart served to highlight what manage- ment already believed. The uncertainties and assumptions that signaled danger dropped out of the information chain when the Mission Evaluation Room manager condensed the Debris Assessment Teamʼs formal presentation to an informal verbal brief at the Mission Management Team meeting.
-- CAIB Report, Vol. 1, pp 201

Any technical organization will have trouble clearly and accurately communicating technical information to decision makers, who most likely are not experts in the particular domain being addressed, and often are not even technically competent in any speciality. Even when they are experts, using an inappropriate tool never helps.

From Thomas Ricks' book Fiasco:

[Army Lt. General David] McKiernan had another, smaller but nagging issue: He couldn't get Franks to issue clear orders that stated explicitly what he wanted done, how he wanted to do it, and why. Rather, Franks passed along PowerPoint briefing slides that he had shown to Rumsfeld: "It's quite frustrating the way this works, but the way we do things nowadays is combatant commanders brief their products in PowerPoint up in Washington to OSD and Secretary of Defense...In lieu of an order, or a frag [fragmentary order], or plan, you get a bunch of PowerPoint slides...[T]hat is frustrating, because nobody wants to plan against PowerPoint slides." That reliance on slides rather than formal written orders seemed to some military professionals to capture the essence of Rumsfeld's amateurish approach to war planning. "Here may be the clearest manifestation of OSD's contempt for the accumulated wisdom of the military profession and of the assumption among forward thinkers that technology--above all information technology--has rendered obsolete the conventions traditionall governing the preparation and conduct of war," commented retired Army Col. Andrew Bacevich, a former commander of an armored cavalry regiment. "To imagine that PowerPoint slides can substitute for such means is really the height of recklessness." It was like telling an automobile mechanic to use a manufacturer's glossy sales brochure to figure out how to repair an engine.

This raises some of the same issues discussed in the report by members of the NASA Return to Flight Task Force in the contribution at the top of this thread.
-- Edward Tufte, August 11, 2006

Writing a good technical report or white-paper takes thought and hard work. Clearly communicating intent with effective mission orders takes initiative and competence. Throwing bullets on a slide is the path of less resistance (I say this as someone who has regrettably succumb to this temptation on more than one occasion). It is also the refuge of the bureaucrat "of an entirely different character" planning to profit from ambiguity and disingenuous politicking.

Knowing a thing is bad is nice in a theoretical way, but what about actually doing well because of that knowledge. How to effectively communicate technical results to hurried managers? I think the Tufte-inspired handout format is a nice start. Doing a good job of tightly incorporating relevant graphics into succinct narrative is probably the best an engineer can do. Improving the report to the point that it isn't a long-winded, dry, snooze-inducer means that management may actually read it, and their decisions may benefit from it. Of course, the hand-out can reference more in-depth technical documents for the competent, curious and critical reader.

Chebyshev Spectral Method Shock Filtering

Now that I've got a nice handy spectral derivative routine, time to test it out on some practical problems (ie non-smooth).

A useful scalar test function is a square wave with a little Gaussian blip added to it:

This is good test function because it contains discontinuities as well as small smooth features that we'd like to be able to resolve.

There are tons of options to filter out high frequency "ringing" around discontinuities (Gibbs phenomenon). One nice option is the Butterworth filter, the gain for which is:

Where the cut-off frequency and the order of the filter are free parameters you can "fiddle" with. For this set of derivatives I'm going to choose a cut-off frequency of half the highest frequency in my grid and see what the effect of the order parameter is. Here's a plot of the gain for a couple different order filters:

Notice they all have the same gain at the cut-off frequency, they just pass through it differently.

The un-smoothed spectral derivative compared to the derivative after different order Butterworth filters are applied to the Chebyshev coefficients before the inverse transform are shown below.

The presence of the discontinuities causes "ringing" throughout the entire grid, but the filters do a pretty good job of containing it to the area immediately around the shocks. Here's a zoomed in view of the smooth feature area:

Without filtering the derivatives of the smooth feature are totally swamped by ringing, but you actually get a nice result out after some filtering.

FFTW Discrete Cosine Transform Derivative

I've done several posts lately about finite difference methods. Spectral methods are the ultimate extension of these same basic ideas used from Boole on up to the present. For an interesting read on this progression (or lack thereof) check out the preface of Boyd's spectral methods book.

Instead of using a series of overlapping local derivative expressions, the derivative at each point is calculated with an interpolating polynomial that spans the entire grid. This leads to exponential (rather than the usual first, second, etc) order of convergence. For folks used to finite differences, Fornberg has written a concise guide that focuses especially on the connection between finite difference approximations and pseudospectral methods.

There are a couple of things that make spectral methods practical for a wide variety of engineering problems (not just periodic boundaries). The most important being Chebyshev Polynomials, which through a little change of variables can be calculated with Fourier transforms, or discrete cosine transform (DCT) for a factor of four speed-up over the FFT. Luckily the excellent FFTW library has specialized transforms for real, even data which we can use to calculate our DCTs with nice, optimized O(n*log(n)) algorithms for arbitrary length data. Here is a thread discussing the usage of FFTW's DCT for calculating derivatives.

Chebyshev polynomial in terms of cos():

From the FFTW docs the forward and inverse discrete cosine transforms (DCTs) are


Notice the similarity between the Chebyshev polynomial and the formula for the coefficients of the DCT.

Here's a Fortran routine that uses the DCT Type II and Type III from FFTW to calculate the derivative of a non-periodic function. An important thing to do when using spectral methods on non-periodic functions is to choose grid points that are at the roots of the Chebyshev polynomials. This makes them equally spaced in the angular coordinate, so we can use the standard transform libraries (which expect equally spaced data). Choosing nodes this way also helps with Runge's Phenomenon.

subroutine dct_diff(dydx, y, half_interval, n)
! use the discrete cosine transform from fftw to calculate the derivative
integer, intent(in) :: n ! number of data points
double precision, intent(in), dimension(n) :: y ! function values
double precision, intent(out), dimension(n) :: dydx ! derivative values
double precision, intent(in) :: half_interval ! half the interval length
double precision, dimension(n) :: beta ! derivative coefficients
double precision, dimension(n) :: alpha ! function coefficients
integer :: plan1, plan2, i
integer :: n_logical ! the logical size of the transform, size of
! the corresponding real symmetric DFT

! the logical size depends on the type of transform, check the docs:
! http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
n_logical = 2*(n)

! forward DCT:
call dfftw_plan_r2r_1d(plan1, n, y, alpha, FFTW_REDFT10, FFTW_ESTIMATE)
call dfftw_execute_r2r(plan1, y, alpha)
call dfftw_destroy_plan(plan1)
alpha = alpha / half_interval

! recurrence for the derivative coefficients:
beta(n) = 0D0
beta(n-1) = 2D0 * dble(n-1) * alpha(n)
do i = n-1, 2, -1
beta(i-1) = beta(i+1) + 2D0 * dble(i-1) * alpha(i)
end do
beta = - beta ! this makes it work, but not sure why!

! inverse DCT:
call dfftw_plan_r2r_1d(plan2, n, beta, dydx, FFTW_REDFT01, FFTW_ESTIMATE)
call dfftw_execute_r2r(plan2, beta, dydx)
call dfftw_destroy_plan(plan2)

! FFTW computes the un-normalized transforms, normalize by logical size
dydx = dydx / dble(n_logical)

end subroutine dct_diff

I compile this with F2Py so I can use all the nice facilities of Numpy and Scipy. Here's a Python convenience function for giving nodes based on the Chebyshev roots:

def arb_int_opt_nodes(a,b,n): # function to give Chebyshev points on an interval
y = numpy.zeros(n)
for i in xrange(n):
y[n-i-1] = numpy.cos( (2*i + 1)*numpy.pi/(2*n) )
return(((b - a) * y + b + a) / 2.0)

Try it out on a test function:

with a known analytical derivative:

The test function and its derivative using six grid points:

Convergence graph:

Pretty awesome, down to the round-off plateau with less than 20 points on the interval. This recent book on Chebyshev polynomials has a good tabulated example of integrating and differentiating a Chebyshev series on pp 33 - 35 that was very useful in debugging.

Perhaps a more practically useful test function is one-wavelength of cos(). This will give us an idea of how many nodes per wavelength we need for high accuracy solutions.


Which looks to give a rule of thumb of about twenty nodes per wavelength to get down to the round-off plateau (32 bit double precision).

Scroll back up and take a look at the Fortran again, see that sign change after the recurrence? Know why it's there? Leave a comment and let me know if you think you do.

Maxima Functions to Generate Finite Difference Expressions

I've written a couple posts demonstrating the use of Maxima to calculate finite difference expressions of various varieties (complex step derivatives, compact differences, plain old differences). The batch files associated with those aren't very user friendly, if you want to calculate different expressions you have to go in and edit them in several places. So I decided to write a couple functions that would take two arguments and return the finite difference expression you asked for.

Maxima file for explicit finite difference expressions.

Maxima file for implicit or compact finite difference expressions.

Usage is just:

load("explicit_fde.mac");
foo : explicit_fde([-1,0,1],1);

which would put a central difference expression for the first derivative into foo. The first argument is a list representing the stencil to use, the second argument is the derivative you want to estimate.

The implicit expressions are similar:

load("implicit_fde.mac");
foo : implicit_fde([-1,0,1],1,1);

The extra argument is so that you don't get the trivial solution back (see this post for details). Also, the return value for the implicit expression is the augmented coefficient matrix for the first derivative. That means the last column is the explicit side of the expression (terms involving the function value), and the first n columns are the coefficients for the first derivative at the stencil points.

Thanks to the helpful folks on the Maxima List for fixing some of my mistakes and giving me some good debugging tips.

Open Math Tools Custom Search Engine

I found a really simple how-to for creating custom Google search engines and used it to make my own little Open Source Math Tools search engine. I find it really annoying when I'm searching for documentation or examples in Maxima that the results are full of Nissan Maxima pages. Things are even worse when trying to search for stuff relevant to R.

Using the custom search engine I don't have to specify that I'm searching for Maxima, Octave, R or SciPy stuff, because I already told Google to favor results from those pages. Right now the custom search engine (CSE) only has about twenty sites included here's the big ones (list updated 10 Nov 2009):

• ntrs.nasa.gov
  
• mit.edu
  
• umich.edu
  
• projects.scipy.org/pipermail/*
  
• nabble.com
  
• nist.gov
  
• http://www.csulb.edu/~woollett/*
  
• http://beshenov.ru/maxima/*
  
• axiom-developer.org
  
• clusterbuilder.org
  
• linuxhpc.org
  
• clustermonkey.net
  
• matplotlib.sourceforge.net
  
• f2py.org
  
• http://cens.ioc.ee/projects/f2py2e/*
  
• scipy.org
  
• http://maxima.sourceforge.net/email-archives/*
  
• http://www.math.utexas.edu/pipermail/maxima/*
  
• maxima.sourceforge.net
  
• www.gnu.org/software/octave
  
• octave.org
  
• r-project.org
  

And of course the mailing list archives for those projects as well. The idea of this CSE is to help the user find articles, tutorials and documentation about using open source math tools to do specific tasks.

The results are pretty good. I choose the option of including broader web search results, so high ranking results from Wikipedia or MathWorld still show up, but the first page is full of relevant links to the favoured sites. For instance the results for 'bayes model averaging' have a mix of about half and half papers from academic websites and links to various R project pages, or the results for the Maxima function augcoefmatrix() which links to the Maxima manual as well as several web sites with usage examples.

It really improves the signal to noise ratio. Leave a comment if there's a site you think I'm missing.

Deriving Compact Differences with Maxima

I wrote a while back about using Maxima to derive asymmetric complex step finite difference expressions for derivatives. The basic approach was sound, but it didn't really use much of the power of the built-in symbolic computation capability. This post will explore deriving compact difference relations with Maxima and avoid all the for-loops by using Maxima's built-in knowledge of Taylor series and Linear Algebra.

First we'll use the taylor() function to define truncated Taylor series expressions for the function and it's first derivative. Then we'll us the ev() function to evaluate the expressions with different values of x: x0, x0+h, and x0-h. Where h is the step size.

f_taylor : taylor(f(x),x,x0,5);
fp_taylor : taylor('diff(f(x),x),x,x0,4);

rhs_list : [at(f(x),x=x0),
'at('diff(f(x),x), x=x0),
'at('diff(f(x),x,2), x=x0),
'at('diff(f(x),x,3), x=x0),
'at('diff(f(x),x,4), x=x0),
'at('diff(f(x),x,5), x=x0)
];

f_coefs : submatrix(factor(augcoefmatrix([f_taylor], rhs_list)), 7);
fp_coefs : submatrix(factor(augcoefmatrix([fp_taylor], rhs_list)), 7);

Using the augcoefmatrix() function is an easy way to extract the coefficients of the various derivatives in the Taylor series. Coefficients of the function expansion:

Coefficients in the first derivative expansion:

Now we need to pack the coefficients for each point into a matrix.

/* add the columns for the points with function values: */
A : addcol(transpose(ev(f_coefs, x=x0-h)));
A : addcol(A, transpose(ev(f_coefs, x=x0)));
A : addcol(A, transpose(ev(f_coefs, x=x0+h)));

/* add the columns for the points with function derivative values: */
A : addcol(A, transpose(ev(fp_coefs, x=x0-h)));
A : addcol(A, transpose(ev(fp_coefs, x=x0)));
A : addcol(A, transpose(ev(fp_coefs, x=x0+h)));

The matrix A represents the linear system we need to solve to find the coefficients of the function value and its derivative to use in our finite difference approximation.

Since we might want to solve for more than one expression, it's a good idea to factor A rather than just inverting it. The matrix should look familiar if you've used a Taylor Table before.

/* we need to invert A to get finite difference formulas so factor it: */
M : lu_factor(A);
[P,L,U] : get_lu_factors(M);

Now it is a straight forward matter to generate a right-hand-side vector for our linear system and solve for the unknown coefficients. The numerical rhs-vector represents the derivatives of the function.

Putting a one in the appropriate position of the vector gives the corresponding finite difference expression for that derivative.

/* right hand side vector for getting the first derivative: */
b_fp : matrix([0],[1],[0],[0],[0],[1]);
/* coefficients for the compact difference formula: */
a_fp : fullratsimp(lu_backsub(M, b_fp));
/* values vector (function and derivative at the points): */
vv : matrix([at(f(x),x=x0-h),
at(f(x),x=x0),
at(f(x),x=x0+h),
'at('diff(f(x),x),x=x0-h),
'at('diff(f(x),x),x=x0),
'at('diff(f(x),x),x=x0+h)
]);
/* compact difference formula for the first derivative: */
cd_fp : vv.a_fp;
cd_fp_reduced : expand((h**4) * (cd_fp) / 30);

You might notice the extra 1 in the 6th place of the right-hand-side vector, it's there so we don't get the trivial solution of using the value for the derivative at x=x0 for our estimate of the value of the derivative (we obviously don't know it yet). This results in the following difference expression.

Which can be multiplied by an arbitrary constant to achieve a more computation ready form.

The h**4 term can be neglected, resulting in the well-known fourth order compact (implicit) expression for the first derivative.

Here's the Maxima batch file used for this.

Maxima Latex Listings

The listings package is a really nice way to include typeset source-code listings in Latex documents. I wanted to include an appendix in one of my documents to show the Maxima calculations I used to arrive at some results shown in the paper. Unfortunately Maxima is not one of the (many, many) built-in languages recognized by the package. I tried using C, because the comments in Maxima are like C comments, but it doesn't handle single quotes (used to get noun forms like 'diff()) well. So, how to include listings of Maxima code in Latex?

It turns out to be really easy to define your own language in the preamble of the document.

\lstdefinelanguage{Maxima}{
keywords={addrow,addcol,zeromatrix,ident,augcoefmatrix,ratsubst,diff,ev,tex,%
with_stdout,nouns,express,depends,load,submatrix,div,grad,curl,%
rootscontract,solve,part,assume,sqrt,integrate,abs,inf,exp},
sensitive=true,
comment=[n][\itshape]{/*}{*/}
}

All this definition does is italicize comments and bold face the keywords, which is sufficient to make the code more readable and pleasant looking. Here's an example screenshot of the results:

Performance optimization (allocation inside a for loop)

An interesting discussion about an old topic just popped up on the Octave lists. It's not interesting because someone noticed that failing to pre-allocate a vector before entering a loop is slow (doesn't anyone RTFM?), but because of the follow-up discussion on indexing and range objects.

Here's the interesting bit:

octave:1> tic(); n=1e5; retval=1:n; toc()
Elapsed time is 0.000528962 seconds.
octave:2> tic(); n=1e5; retval = (1:n)(1:n); toc
Elapsed time is 0.00593709 seconds.
octave:3> tic();n=1e5;retval=[1:n]; toc
Elapsed time is 0.010952 seconds.

Why the significant difference in performance? According to jwe:

In Octave, an expression like 1:n creates a range object, which contains only the base, limit, and increment as double precision values, so no matter how many elements are in the range, it only takes a few bytes of storage (24 for the data plus some overhead for the internal octave_value object itself).

If you write [1:n], you force a Matrix object with N elements to be created. It will require 8*N bytes of storage, plus the overhead for the internal octave_value object itself.

Similar to the difference between range() and xrange() in Python.