NALPAL: Not A Livermore Physics Applications Language

Background and Motivation

The problem I am interested in: leverage high-level symbol manipulation capabilities in a computer algebra system (CAS) to generate compile-able code for number crunching. More specifically, I am interested in generating routines for numerical partial differential equation (PDE) solvers. This is not a new idea, significant work was accomplished starting in the early 1980s [1, 2, 3], continued into the late ’80s [4, 5, 6], and early ’90s [7]. However, those capabilities have not been consolidated and formally incorporated into my open-source CAS of choice, Maxima (though there is fairly recent work implementing similar ideas for Mathematica [8, 9]). The early work [1] quotes Hamming’s Turing Award speech [10] for motivation:

We seem not to be able to use the machine, which we all believe is a very powerful tool for manipulating and transforming information, to do our own tasks in this very field. We have compilers, assemblers, monitors, etc. for others, and yet when I examine what the typical software person does, I am often appalled at how little he uses the machine in his own work.

I had been thinking about this problem for a while because of my interest in method of manufactured solutions (MMS) for code verification. I thought things were really going to take off when I came across Cook’s 1990 technical report describing code generation for PDE solvers using A Livermore Physics Applications Language (ALPAL), which was built on top of DOE-MACSYSMA. ALPAL aimed to be a general purpose domain specific language for turning symbolic descriptions of PDEs into code. Here’s a solution to the problem already developed! Of course, it couldn’t be that easy. I asked the Maxima mailing list if anyone knew what became of this effort: ALPAL question to Maxima list, and ended up answering my own question by getting a hold of some of the old ALPAL devs: ALPAL answer to the list. Unfortunately, there is a long history behind the divergence between different versions of MACSYMA (and different groups of developers) that mitigates against this old code ever working with the fork that became Maxima (should an archive of the source ever actually turn up, update: it did, see comment) [11, 12, 13].

As you may have guessed from this post’s title, and the difficulties described in the previous paragraph, I’ve decided to pursue a more toolbox approach rather than (re)implementing a domain specific language and associated decision logic. The working title for my toolbox of utilities will be NALPAL, both as a nod to the valuable historical work in this area, and an indication of the path not taken.

Wirth categorizes three main approaches to systems for numerical PDE solution: black-box systems, subroutine packages, and code generation aids [1]. Black-box systems are intended for the novice user, and tend to be constrained to very specific problem types or incorporate a great deal of decision logic. Subroutine packages provide reusable building blocks, but require more knowledge and coding on the part of the user, and tend to remain problem specific. Code generation aids require the user to be an expert in numerical analysis, and tend not to automate any of the analytical work that precedes the code generation step.

Wirth’s suggested approach is a Unix-like toolbox approach where standard subroutine libraries are used when it makes sense to do so, and utilities for doing custom code generation are written in the CAS environment. The alternative is to create a full language, and incorporate the semantics and decision logic to automate the entire process for the novice user. The work of Wirth, Steinberg and Roache is a good example of the former approach, and Cook’s work on ALPAL is an example of the latter (though his early work [3] is more along the lines of the toolbox approach, going so far as to say “the viewpoint taken is that one should be able to start with integro-differential equations and end up with optimal FORTRAN code segments. The key in achieving this goal was to tackle only the specific algebraic problem-i at hand, and not to attempt to provide tools to cover every conceivable numerical scheme.”). The toolbox approach keeps more of the work for the problem solver in set-up/tear-down/customization, while the language approach loads more of the burden on to the domain specific language developer. Wirth’s second objective sums up well the approach that I think makes the most sense (emphasis original):

Build productivity enhancing tools of broad applicability for the expert user so that efficient, special purpose PDE codes can be built reliably and quickly, rather than attempt to second guess the expert and build general purpose PDE codes (black box systems) of doubtful efficiency and reliability.

There are still useful things to learn from ALPAL even if we have chosen the less ambitious road, since the problem being solved is the same. A basic description of the ALPAL use-case [7]:

1. Take as input a PDE description, along with boundary and initial condition definitions
2. Discretize the PDE
3. Analyze the result (e.g. for stability)
4. Calculate the Jacobian (needed for Newton methods in implicit time-integration or non-linear boundary value problem (BVP)s)
5. Generate code

Even if we don’t have a full-featured language, the user of our set of utilities will still be interested in accomplishing largely these same steps. In fact, Wirth’s early work gives these steps (reproduced here verbatim) [1]:

1. Manipulate the set of partial differential equations to cast them into a form that is amenable to numerical solution. For vector PDEs, this might include vector differential calculus operations and reexpression in scalar (component) form, and the application of a linearization approximation for non-linear PDEs.
2. Discretize the time and space domain, and transform the partial differential operators in the PDEs into finite difference operators. This transforms the partial differential equations into a set of algebraic equations. A multitude of possible transformations for the differential operators are possible and the boundary conditions for the PDEs also must be appropriately handled. The resulting difference equations must be analyzed to see if they form an accurate and numerically stable approximation of the original equation set. For real world problems, this analysis is usually difficult and often intractable.
3. After choosing a solution algorithm from numerical linear algebra, the finite difference equations and boundary conditions are coded in a programing language such as FORTRAN.
4. The numerical algorithm is then integrated with code for file manipulations, operating system interactions, graphics output, etc. forming a complete computer program.
5. The production program is then executed, and its output is analyzed, either in the form of numerical listings or computer-generated graphics.

Wirth goes on to say (emphasis added),

With continuing advances in computer technology, the last step in this process has become easier. For a given class of problems, answers can be calculated more quickly and economically. More importantly, harder problems which require more computational resources can be solved. But the first four steps have not yet benefited from advances in computer performance; in fact, they are aggravated by it.

Also, this additional bit of motivation from Chapter 5,

Taken together with the software described in other chapters, these tools allow the user to quickly generate a FORTRAN code, run numerical experiments, and discard the code without remorse if the numerical results are unsatisfactory.

This is similar in sentiment to the idea of numerical throw away code.

PDE Code Gen Recipes

With the problem background and motivation set, the rest of this post will focus on pulling out useful Maxima recipes demonstrated by folks who have generated FORTRAN from MACSYMA with intent to inflict grievous arithmurgical damage on hapless PDEs.

Much work is done in most of these articles from the 1980s to avoid array references and break up expressions by hand to reduce the computational cost in MACSYMA. Also, MACSYMA’s knowledge of the chain rule is not used for calculating the transformations because of an explosion in the number of terms [5], rather an identity for the derivative of a matrix inverse is used. A lot of this effort seems unnecessary today because speed and memory have improved so much. However, the basic approach and variables used to define the problem is still relevant (compare with this example of Burgers’ equation on a curvilinear grid):

  dep : [f, sigma]; /* the dependent variables */
  curvi : [xi, eta, zeta]; /* the curvilinear coordinates */
  indep : [x, y, z]; /* the independent variables */
  depends(curvi, indep);
  depends(dep, curvi);
  nn : length(indep);
  eqn : sum(diff(sigma * diff(f, indep[i]), indep[i]), i, 1, nn);

The result is rather large, but it illustrates Maxima’s knowledge of the chain rule. Of course, generally it is easiest to compute rather than for an arbitrary grid, so you need to make substitutions based on the inverse of the Jacobian of transformation. In Maxima we might do something like

  J : zeromatrix(3,3);
  for j : 1 thru 3 do (
    for i : 1 thru 3 do (
      J[i,j] : ’diff(indep[j],curvi[i])
    )
  );
  K : zeromatrix(3,3);
  for j : 1 thru 3 do (
    for i : 1 thru 3 do (
      K[i,j] : diff(curvi[j], indep[i])
    )
  );
  grid_trans_subs : matrixmap(”=”, K, invert(J));
  /* making substitutions from a list is easier than from a matrix */
  grid_trans_sublis : flatten(makelist(grid_trans_subs[i],i,1,3));

which gives us a list of nine equations we can use to make substitutions so that all our derivatives are with respect to the computational coordinates.

  trans_eqn : subst(grid_trans_sublis, eqn) $
  /* Evaluation took 0.0510 seconds (0.0553 elapsed) using 265.164 KB. */
  trans_eqn_factor : factor(trans_eqn) $
  /* Evaluation took 2.4486 seconds (2.5040 elapsed) using 48.777 MB. */

Factoring the result starts getting on up there in compute time and memory, but still not untractable, or even that uncomfortable for an interactive session. Of course we still haven’t made any difference expression substitutions, and that will expand the number of terms even further. The assumption that the grid is arbitrary is a decision point that would have to be dealt with in a black-box style implementation. Best to leave it to the (hopefully) knowledgeable user.

As an aside, linearization is another example of an assumption that is probably best left to the user. Wirth assumes that a linearization is required to solve nonlinear PDEs, whereas Cook provides for calculating the system Jacobians needed for Newton methods (of course Wirth does note that this is just one approach that could be used, and the tools are flexible enough to be extended to other methods). Using the Jacobian is a more modern approach made practical by the successful development of Newton-Krylov methods. It is impossible to predict the future development of numerical methods that will change the choices that need to be made in discretizing PDEs (though that doesn’t prevent speculation [14]), so a flexible toolbox approach is again indicated.

Once a PDE is defined substitutions of finite differences for partial derivatives must be made to create the discrete approximation. Wirth uses a function called DISCRETIZE which uses the dependencies list (created by a call to depends) to substitute indexed variables for the independent variables. Then the partial derivatives of the indexed variables are replaced by finite difference operators. The substitutions of difference expressions is controlled by using Maxima’s pattern matching rules. The basic recipe is

1. Define dependencies between independent and dependent (and possibly computational coordinates)
2. Associate a list if indices with the coordinates
3. Define rules that transform differential terms into difference terms with the appropriate index shifts and constant multipliers corresponding to the coordinate which the derivative is with respect to and the selected finite difference expression
4. Apply the rules to the PDE to give a finite difference equation (FDE)
5. Use Maxima’s simplification and factoring capabilities to simplify the FDE
6. Output the FDE in FORTRAN format and wrap with subroutine boilerplate using text processing macros

The boilerplate is a set of standard templates for a particular solution procedure. The example Wirth uses is an alternating direction implicit (ADI) scheme. A more modern example might be a Newton-Krylov scheme. Wirth describes an environment, Fast FORTRAN Programing (FFP), whose goal is to move computational physics out of a “cottage industry” state. He describes it thusly, “The system consists of two major components: a subroutine library and a command language for building, compiling and running codes.” Based on my experience, that sounds a whole lot like Scipy, which is built on top of Python, and packages together many of the classic scientific computing libraries.

Pattern matching in Maxima is relatively straight-forward. For instance, say I have a function that I’d like to use to calculate my derivatives (such as one based on the discrete cosine transform (DCT)), I could declare patterns that replaced derivatives with function calls.

  matchdeclare([fmatch,xmatch,nmatch],all),
  defrule(deriv_subst_1, ’diff(fmatch,xmatch,1), diff_1(fmatch,xmatch)),
  defrule(deriv_subst_2, ’diff(fmatch,xmatch,2), diff_2(fmatch,xmatch)),
  mat_expr : apply1(mat_expr, deriv_subst_1),
  mat_expr : apply1(mat_expr, deriv_subst_2),

This would result in all of the first and second derivatives in mat_expr being replaced by function calls diff_1 and diff_2 respectively.

The initial set-up steps Cook uses in [3] are roughly the same as those used by Wirth, with the addition of a constant list. The focus is more on what to do with the discrete expressions once they are generated. The basic recipe is

1. reduce_cnst()
2. gcfac()
3. optimize()
4. fortran()

The constant list allows reduce_cnst to pull constants out of loops, and optimize pulls out common sub-expressions to speed things up. optimize returns a Maxima block, which is similar to a Fortran subroutine. In fact, turning blocks into subroutines in other languages is a common problem which has been addressed on the Maxima mailing list (e.g. block2c, see also cform), and as Dan Stanger points out, this is the purpose of the GENTRAN package. GENTRAN is not yet successfully ported to work with Maxima [13] (though there is a gentran directory that ships in the Maxima source tarball if you want to take a look).

The only reason Cook went to the Lisp level in [3] was to reduce the cost of the factor routines (which is less of a concern now), and to make optimize work with lists of expressions (which it does now). One of the examples cites a need for an outrageous 500 million words of computer memory for one calculation, but that’s roughly half what’s in my old desktop computer I got used about a year ago (for about $100). The computing power available to the average amateur has come a long way since 1982.

Both Wirth and Cook assume that an arbitrary orthogonal coordinate system (like spherical, polar, cylindrical, Cartesian, etc.) will be defined. A slightly more modern approach is presented by Roache et al. [4, 5, 6]. They assume an arbitrary curvilinear coordinate system, which may not have analytical scale factors and metrics (i.e. they include the numerical grid generation task as part of the problem statement / solution procedure).

The approach demonstrated in [5] focuses on transforming a set of PDEs to arbitrary curvilinear coordinates described by Laplace’s equation. This approach couples the PDE solution and grid generation methods together. Modern approaches in the engineering world generally assume that the grid will be generated by a separate tool (the fully coupled approach can still be useful for adaptive or moving grids), though solving problems on a single, canonical grid seems to still be common in the sciences. The MACSYMA routines presented there are

change()

change variables, convert arbitrary second order differential equation in nn variables to an arbitrary coordinate frame in the variables xi[i]

notate

atomic notation for derivatives

notation(exp,vari)

primitive atomic notation

scheme()

introduce differences of unknowns

difference(u,f,exp)

primitive differences, scheme and difference collect the coefficients of the differences and calculate the stencil of the solver and coordinate transformations

myFORTRAN()

write the FORTRAN code

A lot of this extra effort is avoidable now, because it is tractable to use Maxima’s knowledge of the chain rule, and the built-in indexed variables and pattern matching facilities.

Conclusion

After digging in to the old literature on generating code from MACSYMA, and trying out a few things with current (v5.20.1) Maxima, it seems like all the pieces you need to generate PDE solving code already ship with Maxima (not touched on in this post were the vector calculus and tensor capabilities that Maxima has as well). I kind of already knew this since I'd been generating solver code in an ad hoc sort of way already. Perhaps there is a place for building up a pattern matching rules database, and maybe boilerplate templates. That seems to be the area that the modern efforts are focused on (e.g. Scinapse claims that it’s rules encoding PDE solution knowledge make up half it’s 120kloc and is the fastest growing part of the code-base [9]). The recipes presented here seem more like a documentation, or tutorial item rather than a possible new Maxima package.

References

[1]   Wirth, M. C., On the Automation of Computational Physics, Ph.D. thesis, University of California, Davis, 1980.

[2]   Wirth, M. C., “Automatic generation of finite difference equations and fourier stability analyses,” SYMSAC ’81: Proceedings of the fourth ACM symposium on Symbolic and algebraic computation, ACM, New York, NY, USA, 1981, pp. 73–78.

[3]   Cook, G. O., Development of a Magnetohydrodynamic Code for Axisymmetric, High-β Plasmas with Complex Magnetic Fields, Ph.D. thesis, Brigham Young University, December 1982.

[4]   Florence, M., Steinberg, S., and Roache, P., “Generating subroutine codes with MACSYMA,” Mathematical and Computer Modelling, Vol. 11, 1988, pp. 1107 – 1111.

[5]   Steinberg, S. and Roache, P. J., “Symbolic manipulation and computational fluid dynamics,” Journal of Computational Physics, Vol. 57, No. 2, 1985, pp. 251 – 284.

[6]   Steinber, S. and Roache, P., “Using VAXIMA to Write FORTRAN Code,” Applications of Computer Algebra, edited by R. Pavelle, Kulwer Academic Publishers, August 1984, pp. 74–94.

[7]   Cook, G. O., “Construction of large-scale simulation codes using ALPAL (A Livermore Physics Applications Language),” Technical Report UCRL-102469, Lawrence Livermore National Labs, October 1990.

[8]   Husa, S., Hinder, I., and Lechner, C., “Kranc: a Mathematica package to generate numerical codes for tensorial evolution equations,” Computer Physics Communications, Vol. 174, No. 12, 2006, pp. 983 – 1004.

[9]   Akers, R. L., Kant, E., Randall, C. J., Steinberg, S., and Young, R. L., Enabling Technologies for Computational Science, Vol. 548 of The Springer International Series in Engineering and Computer Science, chap. SCINAPSE: A problem solving environment for partial differential equations, Springer, 2000, pp. 109–122.

[10]   Hamming, R. W., “One Man’s View of Computer Science,” Journal of the Association for Computing Machinery, Vol. 16, No. 1, January 1969, pp. 3–12.

[11]   Cook, G. O., electronic mail communication, August 2010.

[12]   Fateman, R. J., electronic mail communication, August 2010.

[13]   Stanger, D., electronic mail communication, August 2010.

[14]   Houstis, E. N. and Rice, J. R., “Future problem solving environments for computational science,” Mathematics and Computers in Simulation, Vol. 54, No. 4-5, 2000, pp. 243 – 257.

Derivative of Unequally Spaced Points

I posted a while back on solving Burgers' equation on a moving grid. The same coordinate transformation techniques applied there can be used to calculate the derivative of unequally spaced data points. Just as in the moving grid application, the grid transformation takes the points in the physical space to the computational space. Once in computational space, the sky's the limit!

One way to do this would be just a simple, low-order finite difference for the first derivative, maybe using an average of two one-sided differences (that would reduce to a central difference if the points were equally spaced). We can easily do better than that though since there are so many great, optimized transforms available in the FFTW library. Spectral methods give results that are accurate down to the precision of the machine for more than about 20 points (for functions without discontinuities).

The cost of using a grid-transformation is that we have to calculate two derivatives and then multiply them together, and the cost of a spectral method is O(n*log(n)) instead of O(n) as in a finite difference approximation. In one-dimension things are pretty easy, the chain rule shows us the two derivatives we need to calculate:

It is not straight-forward to directly calculate the additional derivative, so we calculate it by inverting the derivative of the physical coordinate with respect to the computational one:

Using the discrete cosine transform from FFTW (Chebyshev psuedospectral method) is straightforward to do by using F2Py and Python:

n = 17 # number of data points
a = 0.0 # left side of the interval
b = 3 * np.pi # right side of the interval
alpha = 0.25 # parameter to scale the perturbation
dx = (b - a) / float(n) # this is the equally spaced delta x
half_interval = (b - a) / 2.0 # so the dct diff is properly normalized
n_perturb = 3 # number of perturbed sets to try

dx_perturb = alpha * dx * (sp.random.random_sample((n,n_perturb)) - 0.5)

x = np.linspace(a, b, n)
dxdxi_perturb = np.zeros((n,n_perturb),dtype=float)
dydxi_perturb = np.zeros((n,n_perturb),dtype=float)
y_perturb = np.zeros((n,n_perturb),dtype=float)
x_perturb = np.zeros((n,n_perturb),dtype=float)

for i in xrange(n_perturb):
x_perturb[:,i] = x + dx_perturb[:,i] # add in the random fluctuations

y = np.cos(x)
y_perturb = np.cos(x_perturb)
dxdxi = cs.dct_diff(x, half_interval)
dy_analytical = - np.sin(x)
dydxi = cs.dct_diff(y, half_interval)
# derivatives with respect to the computational coordinate:
for i in xrange(n_perturb):
dxdxi_perturb[:,i] = cs.dct_diff(x_perturb[:,i], half_interval)
dydxi_perturb[:,i] = cs.dct_diff(y_perturb[:,i], half_interval)

# now actually calculate the derivatives we're after:
dydx = dydxi / dxdxi
dydx_perturb = dydxi_perturb / dxdxi_perturb

The derivative of the physical coordinate with respect to the computational coordinate for the equally spaced case and a couple of cases with some random perturbations added:

The actual derivative that we're after:

This example shows the flexibility of the approach, but it also indicates the importance of smooth changes in grid spacing for accurate derivative approximations. Abrupt changes in grid spacing translates into noise in the derivative.

Burgers' Equation on a Moving Grid: Derivative Forcing

Solving Burgers' equation on a moving grid is a way to handle moving fronts or shocks by adapting the grid to the solution at each time-step. One simple monitor function is just the square of the solution derivative. Since we chose the forced diffusion equation as the governing equation for the grid, we just substitute the central difference approximation for the first derivative squared for the forcing function.

We also include two parameters (s1, s2) so we can adjust the strength of the forcing compared to the smoothing (second derivative term). That gives our semi-discrete equation for the rates as

Here's an example of the algorithm applied to an initial condition made-up of a couple exponential functions.

It shows that the delta-x decreases in areas of high gradients and increases in areas of zero derivative. The wave gets diffused pretty significantly because of the backward Euler time integration scheme.

Another illustrative example of this adaptive grid scheme is shown below. The initial condition is a couple wavelengths of cosine. This example shows that we probably want to force grid points into areas of high curvature (second derivative) as well as high first derivative. Notice the large delta-x around the peaks of each wave, zero derivative in this area, but large second derivative.

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.

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.

Kronecker Product of Sparse Matrices

Gilbert Strang discusses the use of a neat built-in function in Octave, kron(), that's pretty useful for doing multidimensional finite difference methods on PDEs:


kron() takes the Kroenecker product of two matrices. This is especially useful when constructing the large, sparse matrices needed for finite difference approximations.

K is the tridiagonal matrix created for a central difference approximation to the second derivative:

n = 5;
K = spdiag(2*ones(n,1),0) + spdiag(-ones(n-1,1),-1) + spdiag(-ones(n-1,1),1);
spy(K)

K is the differences for a 1D row of points, now we can create the 2-dimensional version for a plane of points with kron():

K2D = kron(speye(n),K) + kron(K,speye(n));
spy(K2D);

The extension to 3 dimensions is exactly as before, now we have a matrix describing the differences between points on the same line in the first sub and superdiagonal, differences across lines in the next bands out from those, and differences across planes in the furthest out bands:

K3D = kron(speye(n),K2D) + kron(K2D,speye(n));
spy(K3D);

Pretty simple! Now you're ready to solve multidimensional partial differential equations.

A Treatise on the Calculus of Finite Differences

Boole, George. "A Treatise on the Calculus of Finite Differences", Macmillan and Co., 1872.

The public domain and Google books is a great combination. You can download the classic by George Boole, A Treatise on the Calculus of Finite Differences (TCFD) and read it at your leisure, for free!

Even though it was written in 1860, this book probably even more relevant today than it was then. Why? Because it provides the mathematical basis for solving higher and higher resolution partial differential equations with faster and faster computers for the complex domains of modern engineering problems. The faster, more accurate and more efficient solution to these types of problems provides the foundation for improving quality of life for everyone.

I like his definition of the calculus of finite differences, it says so much, succinctly:

The Calculus of Finite Differences may be strictly defined as the science which is occupied about the ratios of the simultaneous increments of quantities mutually dependant.

This book has to be understood as a second volume to accompany Boole's first book on Differential Equations. He stresses the importance of the relationship between the two in his preface to TCFD. The discussion of the distinction between a calculus of limits and the calculus of differences highlights this connection.

Like most old math or science books, it is the discussion of important applications that feels dated. The ideas are timeless, but the important applications change with the problems and the technology of the time.

Another timeless idea that bears repeating around all practitioners of the numerical art (from Chapter 8):

... we shall very often need to use the method of Finite Differences for the purpose of shortening numerical calculation, and here the mere knowledge that the series obtained are convergent will not suffice; we must also know the degree of approximation.

To render our results trustworthy and useful we must find the limits of the error produced by taking a given number of terms of the expansion instead of calculating the exact value of the function that gave rise thereto.

Sobel, Scharr and Truncation Error

In image processing the discrete difference operators to calculate gradients and second derivatives are called "masks". They are often referred to as "stencils" in other applications of finite differences for solving PDEs, estimating sensitivities, etc.

Some of the kernels (a mask plus particular values) for calculating gradients are Sobel:

[ 1 2 1]
Gy = [ 0 0 0]
[-1 -2 -1]

[ 1 0 -1]
Gx = [ 2 0 -2]
[ 1 0 -1]

or Sharr:

[ 3 10 3]
Gy = [ 0 0 0]
[-3 -10 -3]

[ 3 0 -3]
Gx = [ 10 0 -10]
[ 3 0 -3]

or the simpler stencil we used in the previous posts, which only used neighbouring pixels in one direction to estimate the gradient. What if we want to design our own masks? A useful tool is the multidimensional Taylor series:

In two dimensions we have:

which includes up to third order terms. We can't actually include the two pure third order terms because we only have three points in either the x or y direction, and we need a minimum of four. We can include the mixed third order terms though, so that gives us 8 constraints on our coefficients. Instead of including all 9 points in the stencil we have to drop one, so in the Maxima code below you'll notice that the central point is commented out.

Lets examine the truncation error for the two masks above. First we define a vector of delta x and delta y corresponding to each point in the stencil:

/* define the points in the stencil or mask: */
steps : matrix( /* dx, dy */
[-1, 1],
[0, 1],
[1, 1],
[-1, 0],
/* [0,0], */
[1, 0],
[-1, -1],
[0,-1],
[1, -1] );

Then vectors describing the coefficients in the two masks (we'll just use the y-component for this example):

/* coefficients in sobel mask */
sobel_y : matrix(
[1],
[2],
[1],
[0],
/*[0],*/
[0],
[-1],
[-2],
[-1]
);
/* coefficients in sharr mask */
sharr_y : matrix(
[3],
[10],
[3],
[0],
/*[0],*/
[0],
[-3],
[-10],
[-3]
);

To evaluate the truncation error we also need to create an operator that is based on the multidimensional Taylor series shown above. This will let us simply multiply those vectors by the operator to see what the masks are actually estimating. Here's the operator:

n : matrix_size( steps );
A : zeromatrix(n[1],n[1]);

for j : 1 thru n[1] do ( A[1,j] : 1 );/* f */
for j : 1 thru n[1] do ( A[2,j] : steps[j,1] );/* df/dx */
for j : 1 thru n[1] do ( A[3,j] : steps[j,2] );/* df/dy */
for j : 1 thru n[1] do ( A[4,j] : steps[j,1]^2/2! );/* d2f/dx2 */
for j : 1 thru n[1] do ( A[5,j] : steps[j,2]^2/2! );/* d2f/dy2 */
for j : 1 thru n[1] do ( A[6,j] : 2*steps[j,1]*steps[j,2]/2! );/* d2f/dxdy */
for j : 1 thru n[1] do ( A[7,j] : 3*steps[j,1]*steps[j,2]^2/3! );/* d3f/dxdy2 */
for j : 1 thru n[1] do ( A[8,j] : 3*steps[j,1]^2*steps[j,2]/3! );/* d3f/dx2dy */

Now the truncation error terms are shown by a couple of matrix vector multiplies:

sobel_te : A.sobel_y;
sharr_te : A.sharr_y;

Rather than just getting a single non-zero entry in the third place of the resulting vector we get a non-zero entry in the last place as well. Looking back at our operator, this is the row for the third order mixed derivative term, (d^3 f)/(dx^2 dy). In other words both the Sobel and Sharr operators introduce an error in their estimation of the gradient that is proportional to this mixed derivative.

How to avoid this error? Well, we can choose a right hand side vector that only has a 1 in the third position and zeros everywhere else, and then invert our operator to solve for what the values in our mask should be.

/* define which derivative we want to estimate, constrain the unwanted
ones to be zero: */
derivs : matrix(
[0], /* f */
[0], /* df/dx */
[1], /* df/dy */
[0], /* d2f/dx2 */
[0], /* d2f/dy2 */
[0], /* d2f/dxdy */
[0], /* d3f/dx2dy */
[0] /* d3f/dxdy2 */
);
inv_A : invert(A);
coeffs : inv_A.derivs;

This gives a result with only non-zero coefficients for the points at [0,1] and [0,-1], in other words, the simple estimate we've already been using! There's really no need to include all those other points for the gradient calculation.

More Edge Detection

I can't really decide on the best way to treat the different channels (red, green, blue) in a color image when doing edge detection. Is it better to just average the magnitude of the gradient in each channel? Treat the channels separately and find red edges, green edges and blue edges? My pragmatic (if ad hoc) choice is to just add them all up, and use a cludge to keep from getting negative gray-scale values.

Here's the Python for gradient-based edges. The magnitude of each channel's gradient is subtracted from 255. Since gray-scale values range from 0 to 255 we use a max() function to set any resulting negative values to 0. The edges detected using this scheme are shown below, on the right.

inpix = img.load()
s = 1.0;
# calculate the gradient in all three channels, sum of squares for mag
for y in xrange(2,h-1):
for x in xrange(2,w-1):
rdx = ( inpix[x+1,y][0] - inpix[x-1,y][0] )/2
gdx = ( inpix[x+1,y][1] - inpix[x-1,y][1] )/2
bdx = ( inpix[x+1,y][2] - inpix[x-1,y][2] )/2
rdy = ( inpix[x,y+1][0] - inpix[x,y-1][0] )/2
gdy = ( inpix[x,y+1][1] - inpix[x,y-1][1] )/2
bdy = ( inpix[x,y+1][2] - inpix[x,y-1][2] )/2
outpix[x,y] = max( 0, 255 - s*( math.sqrt( rdx**2 + rdy**2 ) + math.sqrt( gdx**2 + gdy**2 ) + math.sqrt( bdx**2 + bdy**2 ) ) )

outimg.save(sys.argv[2],"JPEG")

Taking a careful look at the code, you may notice the parameter s floating around in there. This lets you "darken" the edges by choosing something larger than 1. Here's the same image with s=3.0 (on the left):

Like in the previous image processing posts, these examples were made using the Python Imaging Library for image input/output and easy access to the pixels.

Asymmetric Complex Step

To solve Laplace's equation with a complex step finite difference we needed to develop asymmetric difference relations to solve on the +/- imaginary planes.

On the plane in the negative imaginary direction we difference centrally in the real direction and forward two steps in the imaginary direction:

And on the plane in the positive imaginary direction we difference centrally in the real direction and backward two steps in the imaginary direction:

We still want the real part of the answer as we did before, but now we have imaginary coefficients for our function values in the real direction. While we have a compact stencil on the real plane (only nearest neighbors are used), we lose that property on either imaginary plane of the solution. This means our system matrix bandwidth is larger than it would be if we had a purely compact stencil for all parts of the solution.

It's fairly easy to calculate arbitrary derivative approximation formulas using Maxima.

/* calculate the coefficients for a finite difference approximation of
a derivative, can include complex steps */

steps : matrix([1,2,%i,2*%i]); /* multiple of the step size h */

/* sum df d2f d3f d4f */
derivs : matrix([0], [1], [0], [0], [0]); /* the rhs for our system */

n : matrix_size( steps );

A : zeromatrix(n[2]+1,n[2]+1);

for j : 1 thru n[2]+1 do
( A[1,j] : 1,
for i : 2 thru n[2]+1 do
(
if j = 1 then A[i,j] : 0 else A[i,j] : steps[1,j-1]^(i-1)/(i-1)!
)
);

invA : invert(A);

coeffs : expand( invA.derivs );
expand( 20*coeffs );

Where the vector steps tells which points will be included in the stencil, and the vector derivs tells which derivative we want to approximate. This example is calculating a real and complex forward difference approximation for the first derivative.

You can download this Maxima batch file here. This example can be run from the command line:

$ maxima -b fd.mac

or using the batch command from an interactive Maxima session:

(%i1) batch("fd.mac")

Complex Step: Subtractive Cancellation

We showed in an earlier post that approximating first derivatives with complex steps leads to a second order accurate formula that avoids the subtractive cancellation associated with a normal finite difference.

Following the simple convergence study for second derivative approximations we did for cos(). We'll look at the convergence of some first derivative approximation formulas.

The simple complex step approximation is given by

The standard second order finite difference is given by

The 4th order combined complex and real step approximation is

Notice that this formula is not able to avoid subtractive cancellation the way the second order approximation does, but we get a fourth order accurate approximation on a compact stencil.

The convergence of these three formulas for approximating the derivative of cosine is shown in the figure below.

The 4th order scheme runs into subtractive cancellation earlier than the standard 2nd order scheme, but it remains more accurate over the examined range of resolution. The complex step approximation continues to converge as the step size decreases even after the other two formulas diverge.

So for doing sensitivity analysis the second order complex step formula should be used, because only one step needs to be taken, and arbitrarily high precision can be achieved. If you need to approximate a derivative over a region, then the higher order method will pay-off because you can get better accuracy with fewer points.

Planar Magnetron

A magnetron is a device useful for generating lots of microwave radiation, most homes have magnetron's which use roughly 1kW of microwave power for heating food. A magnetron acts like a diode, when the magnetic field between the inner cathode and the outer anode reaches a certain level, electrons starting at the cathode will not be able to reach the anode. This level of magnetic field is known as the Hull cut-off magnetic field of the magnetron.

Suppose we have a simple planar magnetron with constant E and B fields. The path taken by an electron starting from rest at the cathode is governed by the relation of the time rate of change of the particle to the forces on the particle caused by the E and B fields.

Assuming constant particle mass and neglecting relativistic effects gives an equation for particle acceleration:

This second order differential equation can be recast as a system of first order equations:

Where X is the vector of position and velocity components given by

And Y is the vector of velocity and acceleration components given by

The initial conditions are X = 0 and

This system can be integrated numerically with the explicit forward Euler method:

The bisection method is used to find the value of B that causes a maximum of 1 cm
displacement in the y-direction. This 1 cm-high path taken by the electron with Bz = 33.7884
Gauss and Ey = −1/3 stat-volt/cm is shown in the figure below.
Explicit Integration:

Alternatively, an implicit backward Euler method can be used:

The left hand side can be consolidated into a single operator, L, acting on the unkown vector at the future time-step and the additive term due to the electric field:


This operator must be inverted, this can be done once and then each successive time-step
can be calculated with a matrix-vector multiply.


First Order Implicit Integration:


Using a second order accurate approximation for the first derivative (see finite difference for their relationship with derivatives) provides a more accurate solution, and makes the solution for the larger time-step sizes more reasonable:

Collecting all of the n + 1 terms as we did before, gives

Our implicit integration operator then just has a few slight modifications, and the right hand side now includes n and n − 1 terms.




Second Order Implicit Integration:

The second order method achieves something closer the "correct" shape much more quickly than the first-order scheme.

This analysis ignores any losses due to radiation from the accelerated particle (which of course is why magnetrons are useful in the first place).

A pdf and various Octave/Python/Maxima scripts for this problem are on:
https://sites.google.com/site/variousconsequences/Home/magnetron-problem