## Differential approaches

*Ray Girvan* on two alternative packages reflecting a common quest for PDE solutions.

If programs that solve differential equations have been a recurring theme over the last few issues of *Scientific Computing World*, it reflects their central importance to scientific computing. Most physical phenomena are described by differential equations, from simple first order linear ordinary differential equations (ODEs), to the classically hard-to-solve coupled systems of non-linear partial differential equations (PDEs) arising in fluid flow problems.

This theme links a pair of products that illustrate two approaches to the problem of solving "difficult" differential equations. Diffpack, from Numerical Objects AS of Oslo, shows the power of the numerical method. Maple 7, the latest release of Waterloo Maple's well-established mathematics package, shows the major inroads made by research into the problem of analytical solution.

Diffpack is not a program, but a development environment for writing finite element PDE solvers in C++. Diffpack's main feature is to give the programmer new C++ classes of concise high-level constructs - object-oriented programming - for the matrix, grid, and calculus operations needed for this application. Although Diffpack is also available for Linux and Unix, the 32-bit Windows version reviewed here needs Microsoft Visual C++ 6.0, and needs around 200Mb on top of this to install the main Kernel.

It needs to be said straight away that Diffpack is aimed at users who have acquired, or are prepared for hard work to acquire, an extensive understanding of C++ and the mathematics of the finite element method. It is possible, however, to approach Diffpack with less knowledge. For Windows users, the demo CD-Rom provides a very clear entry point, with about 50 code samples ("Jump Start Applications") covering Diffpack’s repertoire: from a simple program to print "Hello World!", via heat flow, Poisson and various wave equations, to difficult problems such as coupled heat conduction and flow in a non-Newtonian fluid. Of these, 17 are precompiled and can be run with only a small System Kit installed.

Using Diffpack from Visual C++ needs little preamble beyond adding paths to Diffpack’s Library and Include files (so that #include can call Diffpack tools). Programs need an **initDiffpack** line and there are some dialect differences - for instance, output **cout** in C++ becomes **s_o** in Diffpack - but otherwise coding and compilation are unchanged. Programs compiled with Diffpack read an input file defining the simulation parameters and output a **SimRes** file that can be piped through "filters" to plotting systems such as MATLAB or Iris Explorer. For prototyping, however, the Windows version has its own graphical user interface, a friendly control panel that can be incorporated in compiled FEM applications. From this you can view and edit the input file, launch and log the simulation, build simple grids, plot graphs, browse the fields of the results file, and display the output. The display, based on the public domain Vtk (Visualisation Toolkit), produces nicely rendered output with detailed options for angle, zoom lighting, contouring, animation and so on.

According to Numerical Objects, a major strength of Diffpack is its flexibility: it can be applied to non-standard PDE problems, including coupled multiphysics PDEs, outside mainstream engineering. Its philosophy is that Diffpack should take over most of the solution process, so that the application programmer need only be concerned with necessary points of data input or decision. Nor does it blindly apply object-oriented techniques, but reserves these constructs "for higher level administrative tasks" while most of the computation happens in fast low-level code. The result runs fast - the makers quote little speed difference from optimised Fortran - while remaining flexible for prototyping.

Facilities in the Kernel, now at v3.5, include vectors, matrices, report generation, pre-processors and postprocessors, batch runs, solvers for sparse linear and non-linear systems, finite element and finite difference meshes, scalar and vector fields, finite elements and FEM algorithms, numerical integration, random number generators, real and complex arithmetic, and ODE solution. Beyond this, Diffpack can be extended with optional Toolboxes: Adaptivity for adaptive meshing; Datafilter, which can import CAD meshes in the major Abaqus, Ansys and Nastran formats; Multilevel, for developing multigrid/multilevel linear solvers; and Parallel, for implementing models on multiprocessor architectures.

Diffpack is well documented, but in forms that need prior knowledge. The main reference is a 680-page hardback book, Hans Petter Langtangen’s *Computational Partial Differential Equations*. Thematically linked to the example programs, it is the text for a University of Oslo course using Diffpack to teach numerical PDE methods, and needs a grasp of C++ along with linear algebra and multivariable calculus. The other documentation, an HTML guide to the C++ classes and the internal documentation of the C++ header files, will mainly be useful to practiced programmers.

However, the companion book does say that Diffpack is not entirely out of bounds to users like myself, with a little C++ and some knowledge of mathematics and general programming, and I had no trouble exploring and risking modifications to the simpler models. Obviously, though, I cannot claim to have explored Diffpack to its limits, and professional numerical programmers are its real intended market. For instance, the Numerical Objects customer list includes corporates such as DaimlerChrysler, Intel and Nasa, and universities such as Bristol, Cambridge, Cornell and Stanford.

The demo CD-Rom features Windows Media Player animations of Diffpack output in many fields: deformation of a metal bracket, electrical activity in the heart, and groundwater flow. Waves are speciality: a flypast of waves in Hasvik harbour, 3D non-linear waves around oil rigs legs, and a wave power plant. These impressive images, I think, are the best indication of what Diffpack can do in the right hands.

**Maple 7** In contrast to Diffpack, Maple is a standalone program. One of the major mathematics packages, it is often characterised as emphasising pure maths, but this picture has begun to change over recent releases. The expansion of Waterloo Maple's online Application Center into engineering, general sciences, statistics and finance, and its establishment of IQTraders.com (recently sold to Engineering.com) to develop Maple-derived Web technology for engineering industry projects, show a solid foothold in applied mathematics.

Maple 7, launched last June, has no real surprises in the interface. Briefly, for readers unfamiliar with Maple, it uses the worksheet format that is ubiquitous for general-purpose maths packages. Essentially this is a command line system, but with more flexibility over what each input-output region - the "execution group" - can hold: arithmetic, algebra, rotatable 3D plots, spreadsheets, non-executable text, and so on. Maple comprises a small kernel that calls code from a larger library, divided into a main section, always available, and packages that need explicitly loading. It also includes a programming language for writing stored "procedures".

However, Waterloo Maple has announced significant additions under three central themes: units, connectivity, and differential equation solution. The first, continuing Maple's shift toward applied mathematics, is a new Units package containing more than 500 unit names, usable either in "standard" or "natural" forms: e.g. **3*unit(feet)** or **3*feet**. All unit calculations output in the default unit set, such as SI, but are easily changed back with a **convert** function. With plenty of synonyms (e.g. **foot = feet = ft**), units are well implemented, bringing Maple into line with Mathematica and Mathcad.

Connectivity is an active area for maths packages, now that MathML, the W3C (World Wide Web Consortium) Mathematical Markup Language, enables online display that is machine-readable for further processing. Maple has followed Mathcad 2001 and Mathematica 4 in supporting MathML, with both import and export of the latest MathML 2.0 standard. A related Maple facility is its new Sockets package for downloading Web data in real time using simple, **OPEN("address", port)** then **GET page.htm\n** commands. Though the Maple documentation is rather abstract, a demo worksheet downloaded from Waterloo Maple's applications library shows how to pull raw HTML from the Lycos.com weather page, then strip out city temperature data using string functions. This worked straight away over a dialup connection, and I easily modified it to read from BBC News Online and the *Scientific Computing World* site. Doing more would need familiarity with TCP/IP protocols, and Maple warns against using Sockets where security is crucial.

There are other package additions such as String, List, CurveFitting and XMLTools, and internal algorithm speedups, particularly numerical integration routines by NAG (the Numerical Algorithms Group). But the most radical overhaul is to the DEtools and PDEtools packages of differential equation solvers. On the ODE front, there’s a new boundary value problem (BVP) solver, and the IVP numerical routines in the main solver **dsolve** have been updated both for stiff and non-stiff equations. Maple 7 now has a powerful ability to attempt symmetry transformations and classify equations against known forms, and this especially improves **dsolve's** handling of non-linear ODEs. The worked examples draw heavily on the ODE lists in Kamke's Differentialgleichungen, and Waterloo Maple’s press releases claim that "Maple 7 solves virtually 100 per cent of standard benchmark test suites, by far the highest achievement in the market".

The most striking advance is in the **pdsolve** for PDEs and PDE systems. Normally trying for an analytical PDE solution from a maths package is at the level of "a slight hope it'll work if the PDE is very easy". With Maple 7 I found myself working through textbooks trying to catch it out.

Even if a general solution isn't found, Maple often returns a partial result separated into soluble ODEs for the individual variables. A **build** command converts this to an explicit form that can be checked with **pdetest**, but it may not be the general solution - and that's where the learning curve steepens. If you do understand, say, the theory of Lie groups, Maple 7 provides the manual options to use this knowledge. If not, it doesn't matter; the default settings are friendly and robust, and I think most scientists will find something useful in this remarkable package of solvers.

**Analysis: Pushing the boundaries** Despite the success of numerical methods, there are strong advantages to finding analytical solutions to differential equations. One reason is that it characterises the solution exactly; another is the potential unreliability of numerics. Even with the most straightforward type of ODE, the initial value problem (IVP), "stiffness" can make numerical solvers go unstable. Users need to watch out for this if the ODE has solution time-scales or coefficients that differ by orders of magnitude, and so choosing a stiff or non-stiff solver is an important first step in ODE numerical solution.

PDEs are even more prone to peculiarities. A characteristically difficult case is the "singularly perturbed" PDE, such as the Navier-Stokes fluid flow equations, where a small factor multiplying the highest order derivative makes numerical methods founder. But Lothar Collatz, in *The Numerical Treatment of Differential Equations*, warns that even for "quite simple examples", approximate methods can converge to false solutions "in a disarmingly innocuous manner". This makes analytical solutions helpful for testing numerical PDE solvers. Even though object geometry may not allow it in general, an exact PDE solution may be derivable for, say, heat flow in a disk to demonstrate that the FEM algorithm is reliable.

The difficulty is that the more exotic the differential equation, the less the chance of an analytical solution. First order ODEs and linearity bode well; non-linearity is a nuisance for derivation because superposition doesn’t work as a means to creating general solutions from special forms. Second order ODEs are already starting to involve special functions such as Airy and Bessel. Stephen Wolfram’s *The Mathematica Book* points out that many such functions were historically designed to represent these ODE solutions. Third order often leads to hypergeometric functions, and for non-linear and higher order ODEs, analytical solutions become the exception.

This scarcity of exact solutions applies even more to PDEs. However, a good many second order linear PDEs with symmetrical forms are soluble, and fortunately these include some classic formulae describing physical phenomena, such as the Laplace, Poisson, diffusion, and wave equations.

The traditional approach to the search for an analytical solution to an ODE or PDE is to attempt to find - or find a transformation that leads to - one of the known differential equation forms tabulated in books such as Erich Kamke’s *Differentialgleichungen: Losungsmethoden und Losungen* or Daniel Zwillinger’s *Handbook of Differential Equations*.

This approach can be automated by computer algebra systems, but there are promising alternatives such as Lie group methods. This dates from the late 19th century when Sophus Lie developed the theory for solving ODEs by identifying symmetries. Impractical in Lie's time, it's being successfully implemented now that symbolic programs can do the laborious calculations.

Zwillinger's book is, incidentally, highly worth reading for its categorised analysis of solutions and their methods, both exact and numeric, right across the spectrum of ODE and PDE solution.

**You can use the online Reader Enquiry service at Scientific Computing World to make contact with organisations referred in this article, or to visit relevant websites.**