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.

Ares/Orion Slip and Classic Mistake

The latest news about the new NASA Constellation program indicates that they are having cost and schedule problems, this is no surprise, everyone has cost and schedule problems. The surprising thing is that they are deleting tests from the planned program development.

CxP attempted to protect the schedule and budgetary pressures by offsetting these additional strains by deleting test items - notably on the Upper Stage. However, this only proved to cause further disconnects throughout the program.

I suppose it should not be too surprising, this is a classic program management mistake. When the schedule gets tight and the costs get too high above estimates, cut testing. Why? Because it is at the end of the program, and program managers have a bit of hubris that leads them to think that the design is sound because they worked so hard on it up-front. Which is obviously the case because they had to spend more time than they thought in that stage of the program, right? Maybe, but also maybe not, that's why the tests are critically important: to validate the design performance. No slide-show, CAD drawing, CFD or FEA simulation can provide that sort of knowledge.

I can't understand why someone would think that under-testing a new product that is probably delayed and expensive because of its new technology and complexity would produce a good result. Rather than removing unneeded technology or complexity, rather than simplifying or reducing the requirements, they are rushing to field a fragile system. It is unfortunate that the burden of this foolishness is borne by the astronaut corps and not the program managers responsible for these short-sited expediencies. The program manager will probably get a reward, move to a new job, and then astronauts will be at undue risk in an already risky business.

Complex Step Jacobian-Free Newton-Krylov

Newton's method are really useful for solving non-linear systems.

Each iteration of the Newton's method requires solving a linear equation:

where the Jacobian is the matrix of partial derivatives

and then applying the update to the solution vector:

For small systems it is tractable to analytically calculate the Jacobian and invert it directly. This rapidly becomes intractable for realistic engineering problems because the number of unknowns can be on the order of 1e6 (if you have a good workstation) to 1e9 (if you have a big cluster). Rather than explicitly creating a 1e6-squared matrix and inverting it at every iteration, an iterative method that only requires the action of the matrix on a vector can be used (such as conjugate gradient or GMRES). Then we can approximate this matrix vector product by

which is a good improvement in storage over the matrix approach, but we can improve it still further by using a complex step

This means we only need one system evaluation to approximate the Jacobian, and much more accurate estimates can be made. As shown in a previous post it is second order accurate, and not as limited by subtractive cancellation as the standard difference approach is.

Here's a simple example of the idea in Maxima

/* trying out complex step finite differences for jacobian free newton
krylov methods */

/* define a simple system of equations */
F : zeromatrix(2,1);
F[1,1] : x[1]*x[2] + 2*x[2];
F[2,1] : 2*x[1] + x[1]*x[2];

/* analytic Jacobian */
J : zeromatrix(2,2);
J[1,1] : diff(F[1,1],x[1]);
J[1,2] : diff(F[1,1],x[2]);
J[2,1] : diff(F[2,1],x[1]);
J[2,2] : diff(F[2,1],x[2]);

u : zeromatrix(2,1);
v : zeromatrix(2,1);
u[1,1] : 1.5;
u[2,1] : 1.0;
v[1,1] : 0.2;
v[2,1] : 0.1;

c_step : 0.0 + 1.0e-200*%i; /* complex step */

/* approximate Jacobian using complex step */
Jv_approx : (
ev(F, x[1] = u[1,1] + c_step*v[1,1], x[2] = u[2,1] + c_step*v[2,1])
) / imagpart(c_step);
imagpart(Jv_approx);

/* evaluate the analytic Jacobian to compare */
Jv_analyt : ev(J.v,x[1]=u[1,1],x[2]=u[2,1]);

A good survey article introducing the JFNK method and describing some applications:
Jacobian-free Newton-Krylov methods: a survey of approaches and applications

Nice short paper about using complex steps in finite differences:
Using Complex Variables to Estimate Derivatives of Real Functions

Three Dimensional Poisson's Equation

This is a follow up to the post I wrote a while back about implementing numerical methods. The goal is to demonstrate a work-flow using tools that come included in the Fedora Linux (and many other) distributions.

1. In Maxima
   
   • Define governing equations
     
   • Substitute finite difference expressions for differential expressions in the governing equations
     
   • Output appropriate difference expressions to Fortran using f90()
     
2. In Fortran: wrap the expression generated by Maxima in appropriate functions, subroutines and modules
   
   
3. In Python
   
   • Compile modules with f2py, which comes packaged along with numpy
     
   • Write the high-level application code, parsing input decks, performing optimization, grid convergence studies, or visualizations that use the compiled Fortran modules
     

Maxima

The governing equation is the three-dimensional Poisson's equation. In Cartesian coordinates The Maxima code to define this equation is straightforward:

depends(u,x); depends(u,y); depends(u,z); /* Poisson's equation: */ p : 'diff(u,x,2) + 'diff(u,y,2) + 'diff(u,z,2) = -f;

Now that we have the continuous differential equation defined we can decide on what sort of discrete approximation we are going to solve. The simple thing to do with second derivatives is to use central differences of second order accuracy, which require only information from closest neighbours in the finite difference grid. Replacing the differential expressions with difference expressions is accomplished in Maxima with the ratsubst() function.

/* substitue difference expressions for differential ones: */
p : ratsubst((u[i-1,j,k] - 2*u[i,j,k] + u[i+1,j,k])/dx**2 , 'diff(u,x,2), p);
p : ratsubst((u[i,j-1,k] - 2*u[i,j,k] + u[i,j+1,k])/dy**2 , 'diff(u,y,2), p);
p : ratsubst((u[i,j,k-1] - 2*u[i,j,k] + u[i,j,k+1])/dz**2 , 'diff(u,z,2), p);
/* substitute the correct array value of f:*/
p : ratsubst(f[i,j,k], f, p);

This assumes that the solution u and the forcing function f are stored in three dimensional arrays, if not then the indexing in the difference expressions becomes more complicated. If we want to apply a Gauss-Seidel method to solve our elliptic problem we just need to solve p for u[i,j,k], and then dump that expression to Fortran format.

gs : solve(p, u[i,j,k]);
del : part(gs[1],2) - u[i,j,k];
/* output for fortran: */
with_stdout("gs.f90", f90(d = expand(del)));

This gives the expression for the update to the solution at each iteration of the solver.

Fortran

The Maxima expressions developed above need control flow added so they get applied properly to our solution array.

do k = 1, nz
do j = 1, ny
do i = 1, nx
d = (dy**2*dz**2*u(i+1,j,k)+dx**2*dz**2*u(i,j+1,k)+dx**2*dy**2* & u(i,j,k+1)+dx**2*dy**2*dz**2*f(i,j,k)+dx**2*dy**2* & u(i,j,k-1)+dx**2*dz**2*u(i,j-1,k)+dy**2*dz**2*u(i-1,j,k))/ & ((2*dy**2+2*dx**2)*dz**2+2*dx**2*dy**2)-u(i,j,k) u(i,j,k) = u(i,j,k) + w*d
end do
end do
end do

The expression for the update was generated from our governing equations in Maxima, one of the things to notice is the w that multiplies our update, this is so we can do successive over relaxation. Also, we'll go ahead and package several of these iterative schemes into a Fortran module

module iter_schemes
implicit none
contains

...bunches of related subroutines and functions ...

end module iter_schemes

Python

Now we want to be able to call our Fortran solution schemes from Python. Using F2py makes this quite simple:

[command-line]$ f2py -m iter_schemes -c iter_schemes.f90

This should result in a python module iter_schemes.so that can be imported just as any other module.

import numpy
from iter_schemes import *

... do some cool stuff with our new module in Python ...

Conclusion

This approach might seem like over-kill for the fairly simple scalar equation we used, but think about the complicated update expressions that can be generated for large vector partial differential equations like the Navier-Stokes or Maxwell's equations. Having these expressions automatically generated and ready to be wrapped in loops can save lots of initial development and subsequent debugging time. It also makes generating unit tests easier, and so encourages good development practice. But does it work? Of course, here's the convergence of the Gauss-Seidel and SOR schemes discussed above along with a Jacobi scheme and a conjugate gradient scheme applied to a 50 by 50 by 50 grid with a manufactured solution. The plot was generated with the Python module matplotlib.