Symbolic Generation of Finite Difference Formulae

Robert M. Corless
Department of Applied Mathematics
and the Ontario Research Centre for Computer Algebra
University of Western Ontario

This work is joint work with Jichao Zhao, a Ph.D. student working under my direction, and is based in part on work done with Jacek Rokicki (T.U.Warsaw).  The applications in mathematical finance (not talked about here) are joint work with Matt Davison (Applied Math, Western) and Jichao Zhao.

Introduction

The finite difference calculus is very old, going back at least to Newton; even as a topic in modern numerical analysis it is very old.  It would be surprising if there were something really new to say...today we will have a tour around some old results, and do some bookkeeping & programming, and make a neat, though minor, refinement on a classical but recently more popular finite difference technique, now called "compact finite differences".  We will end (time permitting) with a discussion of some stability aspects.

From "A Treatise on the Calculus of Finite Differences" by George Boole, 1st ed. 1860, 2nd ed. 1872 edited by J. F. Moulton:

"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 dependent."

From "The Calculus of Finite Differences" by L. M. Milne-Thompson, 1933:

"Let f(x) be a given function of the variable x.  The Differential Calculus is concerned with the properties of which is still a function of the single variable

which is a function of the two variables
and ."

Familiar finite differences include the forward difference above, which approximates the derivative to order h:

 > Delta[h](f)(x) = (f(x+h)-f(x))/h;

 > lhs() = series(rhs(),h,3);

The central difference operator which approximates the derivative to order h^2:

 > delta[h](f)(x) = (f(x+h/2)-f(x-h/2))/h;

 > lhs() = series(rhs(),h,4);

The old-fashioned usages of finite difference formulae were different from today's needs (indeed the first use of the Method of Lines was to derive a partial differential equation, not solve it! (see Hairer, Norsett & Wanner, Solving ODE I, section I.6, for a discussion of the contributions of Brook Taylor 1715, Johann Bernoulli 1727, and d'Alembert 1747), and by and large their operational methods for generation of formulae have fallen by the wayside: we know pretty much all the formulae we need.  Or do we?

 >

Symbolic Generation, the easy way

Generation of finite difference formulae is, nowadays, pretty easy: just look them up, many of them are known already.  But, how can we find them ourselves?

Operator Formalism

The older texts make the connection between operator expansions and series:

 > restart;

 > f(x+h) = series( f(x+h),h );

The following is the "shift" operator; it takes a function, and returns a function with shifted argument:

 > E := f -> (t->f(t+h));

 > E(sin);

 > E(sin)(x);

Composition of operators is represented by the symbol @.  Repeated composition is by @@, analogous to the old Fortran * for times and ** for power.

 > (E@@3)(sin)(x);

"It is obvious that" Taylor expansion means that the operators E and D may be related by the following series, which looks like the exponential series:

 > E = I + h*D + h^2*(D@@2)/2! + h^3*(D@@3)/3! + O(h^4);

In other words, for a slightly modified meaning of the exponential function,

 > E = exp(h*D);

Of course, that statement is a shorthand notation for an operator relation:

 > Delta = E - I;

 > Delta = (exp(h*D) - 1)/h;

 > Dln := solve(%,D);

Wait a minute, that only LOOKED like an exponential---now this one LOOKS LIKE a logarithm...

 > series( %, h );

 > D = %;

Let's check this, that seems somewhat strange:

 > D(f)(x) - (f(x+h)-f(x))/h + (f(x+2*h)-2*f(x+h)+f(x))/2/h ;

 > series(%,h,4);

So this technique appears to work!  Can we soup it up?

 > (f(x+h/2) - f(x-h/2))/h;

 > series( %, h, 8 );

Of course everyone immediately recognizes that series as being a hyperbolic sine...

 > series( 2*sinh(h*D/2)/h, h, 8 );

 > (exp(h*D/2) - exp(-h*D/2))/h = convert( (exp(h*D/2) - exp(-h*D/2))/h, trigh );

 > delta = 2*sinh( h*D/2 )/h;

 > Ddelta := solve(%,D);

 > D = series(%,h,8);

 > delta := f -> (t->(f(t+h/2)-f(t-h/2))/h);

The following is a "Maple trick" to ensure that h is treated as a constant in the following applications of operators.

 > h := t -> h;

 > D6 := delta - h^2*(delta@@3)/24 + 3/640*h^4*(delta@@5);

 > expand(D6(f)(t));

 > series(%,h,8);

Hey, rather neat!  This operator--series correspondence is fun stuff!

One last trick: using Pade approximants instead of series.

 > series(Ddelta,h,6);

 > convert(%,ratpoly);

 > small := expand( (80*(t->t)+9*h^2*(delta@@2))(D(f))(t) - 240/3*delta(f)(t) - 17/3*h^2*(delta@@3)(f)(t) );

 > series( small, h, 8 ):

 > simplify(%);

We have just seen our first "compact" finite difference formula, also known as a Pade method.  Let's try another:

 > series(Dln,h,6);

 > convert(%,ratpoly, 2,2 );

 > Delta := f -> (t->(f(t+h)-f(t))/h);

 > small := expand((10*(t->t) + 12*h*Delta + 3*h^2*(Delta@@2))(D(f))(t) - (10*Delta + 7*h*(Delta@@2) + h^2*(Delta@@3)/3)(f)(t));

 > series( small, h, 7 );

Good, so we have a fifth-order accurate formula that uses four sample points...maybe we can do better??  And, wait a minute, that formula gives us a relation between the derivative of f, at three points, and the values of f, at four points.  We'll have to do some work to sort that out!

 >

Brute Force

Maxim #31, from http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/maxims.html

Computational mathematics is mainly based on two ideas: Taylor series, and linear algebra.

We don't have to be clever about this kind of thing, any more, though: we can simply use brute force.

 > stencil := [[0,0],[h,0],[-h,0],[0,h],[0,-h]];

The stencil is the description of where (which grid points) we know information at.

 > N := nops( stencil );

 > anszatz := add( b[i]*f(stencil[i][1],stencil[i][2]),i=1..N ) + b[N+1]*h^2*(D[1,1])(f)(-h,0) + b[N+2]*h^2*(D[1,1])(f)(h,0) ;

The anszatz is a description of what we know at the stencil, and reflects the desired form of the finite difference approximation.  In the way I have written it above, the anszatz has "zero weight".

 > difexp := h*(D[1])(f)(0,0);

The differential expression is what we wish to approximate by the anszatz.  Now we take Taylor series of the difference between the anszatz and the differential expression, and set the coefficients to zero.

 > zero := series( difexp - anszatz, h, N+3 );

 >

 >

 >

 >

 >

 > sol := solve( eqs, vars );

 >

 > formula := map( factor, subs( sol, anszatz ) );

 > err := map(factor, series( difexp - formula, h, 7 ) );

This worksheet procedure can be modified for whatever purpose you desire...in some sense this is all one needs.  It's a bit clumsy, sort of like re-writing your Fortran every time you want to run the program again, and to make it work one has to ensure that all weights are correct, and that enough series terms are taken, but the final step (generation of the error coefficients) ensures that the results are correct.

 >

 >

Programmatic Generation

The brute force approach discussed above has been turned into a reasonably easy to use program.

 > fd,err := FINDIF( f, [h], [[0],[h],[2*h]], [[[0]],[[0]],[[0]]], h*D[1](f)(h), 3, d );

 > stencil := [[0,0,0], [Delta[x],0,0], [-Delta[x],0,0], [0,Delta[y],0], [0,-Delta[y],0], [0,0,Delta[z]], [0,0,-Delta[z]]];

 > mask := [seq([[0,0,0],[1,0,0],[0,1,0],[0,0,1]],i=1..nops(stencil))];

 > intrp := FINDIF( f, [Delta[x],Delta[y],Delta[z]], stencil, mask, Delta[x]*D[1](f)(a*Delta[x],b*Delta[y],c*Delta[z]), 3, d );

 >

 >

 >

 >

 >

Minimizing error coefficients

 >

 >

 >

 >

 >

 >

 >

 >

 >

 >

Compact Finite Difference Formulae

Maxim #18: Symbolic computing is mainly useful when you want a symbolic answer.

We have used this approach to generate a number of new formulae; new chiefly in the issue of closure at the boundary.  I give a brief tutorial on the use of compact finite difference formulae here, and then discuss the refinements.  First, let us look at a compact finite difference formula for the first derivative.

 > restart;

 > Digits := trunc(evalhf(Digits));

 > f := x -> sin(Pi*x)+cos(2*x);

 > df := (D)(f);

FDMP = First Derivative Matrix Pencil

We use this program to have a look at a particular family of compact finite differences, generated using FINDIF in order for fast solution.  Each row of the matrix pencil represents one invocation of a compact finite difference formula; the first row and the last row are special because they cannot involve points to the left (for the first row) or points to the right (for the last row), and special formulae were generated for those equations.

 >

 >

 > M, C := FDMP( N, h, 0, 0, float[8] );

 >

Why such a peculiar M matrix?  We have that M D(f) = C f, and so in order to find our derivatives we must solve the system.  We see above that M is tridiagonal, which makes it easy...but in fact the first entry of M was chosen to make it even easier!

 >

 > L,U :=LinearAlgebra:-LUDecomposition( M, output=['L','U'] );

 >

 > 2.0+sqrt(3.0);  1/(2.0+sqrt(3.0));

The LU factoring of M can be done exactly, analytically, and consists of factors with constant diagonals...therefore the solution of M f[x] = C f can be carried out in a matrix-free, division-free manner!  See the file FirstBox.mpl.  This is very close to as fast as evaluation of an explicit finite difference formula.  The matrix L is well-conditioned if we choose c = 2 + sqrt(3) (but not if we choose c=2 - sqrt(3), which also formally allows exact factoring).

 >

 > X := Vector(N,i->evalf(i/N));

 >

 >

The following can be done very quickly (almost as fast as central differences) if the analytical factoring of M is used (in fact, of course, neither M, nor C, nor even the explicit factors, are formed explicitly: instead the differences are done in a black-box manner; for this talk we show the matrices explicitly for comprehension, not execution).

 >

 >

 >

In the following, we just check that our procedure really does generate fourth order accuracy:

 >

 > #Digits := 30;

 > Ns := [seq( combinat[fibonacci](k), k=8..14 ) ];

 > err := table();

 > for i to nops(Ns) do  N := Ns[i];  h := 1.0/N;  M,C := FDMP(N,h,-10,-10,float[8]);  Vals := Vector( N, datatype=float[8], i -> evalf(f(i*h)) );  C0   := LinearAlgebra:-LinearSolve(M,C.Vals);  Z0 := C0 - Vector( N, datatype=float[8], i-> evalf( df(i*h)) );  err[i] := LinearAlgebra:-Norm(Z0,2); end do:

 > p1 := plots[loglogplot]( [seq([Ns[i],err[i]], i=1..nops(Ns))], style=POINT, symbol=DIAMOND, symbolsize=30 ):

 > p2 := plots[loglogplot]( [seq([Ns[i], err[3]*(Ns[3]/Ns[i])^4], i=1..nops(Ns))] ): plots[display]( {p1,p2} );

 >

 >

 >

 >

Stability Issues

A good reference here is Carpenter, Gottlieb and Abarbanel, J. Comp. Phys, 108, 1993.  Studying semidiscretized PDE by the method of lines where the spatial derivatives are done by these compact methods leads to studying the eigenvalues (generalized eigenvalues) of the matrix pencils presented here.

Eigenvalues of First Derivative Matrix Pencil

with(LinearAlgebra):

 > with(LinearAlgebra):

 > N := 377; h := 1.0/N;

 > M, C := FDMP( N, h, 0, 0, float[8] );

 >

 > st := time(): E,V:= Eigenvectors( A, output=['values','vectors'] ); time()-st;

 > for i to N do if Re(E[i])>0 then break end if end do;

 > i;

 > E[i];

 > Norm( V[1..-1,i] );

 > plot( [seq([Re(V[k,i]),Im(V[k,i])], k=1..N)], style=POINT );

 > #

 > plot([seq([Re(E[i]),Im(E[i])],i=1..N)],style=POINT,view=[0..7,-50..50]);

 > E;

Eigenvalues of Bound Second Derivative Matrix Pencil

In solving (say) with one wants to discretize and replace the second spatial derivative with a compact finite difference operator so that we may instead solve the vector differential equation  The eigenvalues of then become of interest; they should all be negative, to preserve the stability of this equation for positive It turns out that the use of the boundary conditions is important here, and if we do use them, then the eigenvalues are indeed all real and negative.  The pseudospectra are also quite nice, though is not a normal matrix.

 > Digits := trunc(evalhf(Digits));

 > Ns := [seq( combinat[fibonacci](k), k=8..14 ) ];

 > err := table();

 > f := x -> sin(Pi*x);

 > d2f := (D@@2)(f);

 > for i to nops(Ns) do  N := Ns[i];  h := 1.0/(N+1);  M,C := SDMP(N,h,float[8]);  Vals := Vector( N, datatype=float[8], i -> evalf(f(i*h)) );  C0   := LinearAlgebra:-LinearSolve(M,C.Vals);  Z0 := C0 - Vector( N, datatype=float[8], i-> evalf( d2f(i*h)) );  err[i] := LinearAlgebra:-Norm(Z0,2); end do:

 > p1 := plots[loglogplot]( [seq([Ns[i],err[i]], i=1..nops(Ns))], style=POINT, symbol=BOX ):

 > p2 := plots[loglogplot]( [seq([Ns[i], err[1]*(Ns[1]/Ns[i])^4], i=1..nops(Ns))] ): plots[display]( {p1,p2} );

 > N := Ns[2];  h := 1.0/(N+1);  M,C := SDMP(N,h,float[8]);

 > with(LinearAlgebra):

 > E,V := Eigenvectors( C, M, output=['values','vectors']);

 > bigE := max( seq(Re(E[i]),i=1..N) );

 > plot([seq(h^2/12*[Re(E[i]),Im(E[i])],i=1..N)],style=POINT);

From EigTool (see the Pseudospectra Gateway) we can see that the largest eigenvalue is bounded away from zero, on the left, as are pseudospectra for perturbations that are of reasonably small size.

Future work

Of course, it is the pseudospectra that probably really matters here.  We wish to investigate (more fully) the pseudospectra of these analytically factorable matrices (In the bound 2nd derivative case, the pseudospectra are very well behaved, even though the matrix is quite far from normal).   A detailed description of just how much the computation really costs, and how much we really save by pre-factoring the compact finite difference matrix, would be useful.

 >