Mathematical problems described by partial differential equations (PDEs) are ubiquitous in science and engineering. Examples range from the simple (but very common) diffusion equation, through the wave and Laplace equations, to the nonlinear equations of fluid mechanics, elasticity, and chaos theory. However, few PDEs have closed-form analytical solutions, making numerical methods necessary. The numerical task is made difficult by the dimensionality and geometry of the independent variables, the nonlinearities in the system, sensitivity to boundary conditions, a lack of formal understanding of the kind of solution method to be employed for a particular problem, and so on.
A number of methods have been developed to deal with the numerical solution of PDEs. These fall into two broad categories: the finite-difference methods and the finite-element methods. Roughly speaking, both transform a PDE problem to the problem of solving a system of coupled algebraic equations. In finite-difference methods, the domain of the independent variables is approximated by a discrete set of points called a grid, and the dependent variables are defined only at these points. Any derivative then automatically acquires the meaning of a certain kind of difference between dependent variable values at the grid points. This identification helps us transform the PDE problem to a system of coupled algebraic equations. In the finite-element method, the dependent variable is approximated by an interpolation polynomial. The domain is subdivided into a set of non-overlapping regions that completely cover it; each such region is called a finite element. The interpolation polynomial is determined by a set of coefficients in each element; the small size of the elements generally ensures that a low-degree polynomial is sufficient. A weighted residual method is then used to set up a system of algebraic equations for these coefficients. The numerical problem in both methods therefore is one of solving a system of algebraic equations.
In this article, we present the Mathematica package MathPDE that implements the finite-difference method for time-dependent PDE problems. We give examples of the way in which symbolic programming lends a degree of generality to MathPDE. Before doing so, we would like to put our work in perspective. We first review, very briefly, the extant PDE packages that we are familiar with that have been developed over the years.
Mathematica has a built-in function NDSolve that can numerically solve a variety of PDEs and offers the user a very simple and quick route to solving PDEs numerically. So a natural question that arises is why another solver is needed. In this regard, we emphasize two points. First, MathPDE is able to handle irregular spatial domains, not just rectangular or cubic domains. Second, MathPDE can produce a standalone executable that runs independently of Mathematica. providing portability.
The following systems use either a compiled language directly, or have a high-level interface language that preprocesses the input and employs subroutines in a compiled language. None of them uses symbolic programming.
Diffpack (see [1 ]) presents an object-oriented problem-solving environment for the numerical solution of PDEs. It implements finite-difference as well as finite-element methods and provides C++ modules with a wide selection of interchangeable and application-independent components. ELLPACK (see [2 ]) presents a high-level interface language for formulating elliptic PDE problems and presents more than 50 problem-solving modules for handling complex elliptic boundary value problems. It is implemented as a FORTRAN preprocessor and can handle a variety of system geometries in two dimensions (both finite differences and finite elements) and rectangular geometries in three dimensions (finite differences only). Cogito and COMPOSE were developed at Uppsala University, Sweden (see [3. 4 ]). Both implement finite-difference schemes for time-dependent problems and exploit object-oriented technology to develop a new kind of software library with parts that can be flexibly combined, enhancing easy construction and modification of programs. Cogito is a high-performance solver that comprises the
three layers Parallel, Grid, and Solver, the lower two layers being Fortran 90 parallel code, while the Solver is written in C++. COMPOSE is a C++ object-oriented system that exploits Overture [5 ] for grid generation.
The work on numerical PDE libraries is far too extensive for us to be able to present an exhaustive review here (see [6. 7. 8. 9 ] for overviews). We would like to emphasize that while MathPDE has the generality to handle a wide range of PDE problems, it is not useful, at least in its present form, for certain problems where very specialized finite-difference decoupling schemes are known to work much better than the straightforward finite differencing that our approach uses. The Navier-Stokes equation is one such example (see [4 ] for an approach to this problem). In general, the work discussed in this paper is not adequate for problems where specialized numerical libraries have been known to perform effectively. Instead, it is suitable for PDE problems that are often encountered in engineering modeling situations. Even for such problems, several static subroutine libraries have been developed that perform efficiently, as discussed above. On the other hand, MathPDE can treat a wide range of spatial domains (see Section 3). There is also the great potential to combine the enormous benefits of the interactive Mathematica environment with its large library of built-in mathematical functions, which enables the user to experiment with a large class of PDE problems. Finally, MathPDE generates C++ source code and a standalone executable, addressing the important issue of portability.
Our preliminary results were reported in the references [10. 11. and 12 ]. We now give a brief summary of the present paper. MathPDE accepts the PDE and the initial and boundary conditions, along with the geometry of the system, in Mathematica ’s list format. It supports geometries of fairly general shapes and dimensions. The symbolic engine applies explicit and implicit difference methods and discretizes the geometry to a grid. MathPDE then generates a program for solving the numerical problem, which is essentially an algebraic equation system on the grid; if the system is nonlinear, then a multidimensional Newton-Raphson method [7 ] is used to solve it. This program makes use of the external library SuperLU to efficiently handle sparse linear algebraic systems [13 ]. MathPDE internally calls MathCode [14 ] to translate the numerical program into C++ and then to generate an executable. All these operations are done by invoking a single function (SetupMathPDE ) provided by MathPDE. The executable can then be run either from within Mathematica or externally to obtain the numerical solution of the PDE problem.
Here is an outline of this article. In Section 2, we provide a quick introduction to MathPDE with the example of a one-dimensional diffusion equation.
We discuss the ideas behind the package in Section 3. Section 3.1 is about the numerical algorithms used, such as the discretization of derivatives, the solution of an algebraic system, and Fourier stability. In Section 3.2, we discuss how the spatial domain is discretized. Section 3.3 discusses code generation and its translation into C++ using MathCode .
In Section 4, we discuss application examples that illustrate the strengths and limitations of MathPDE. The examples chosen are divided into problems in one, two, and three dimensions, and there is a section on two nonlinear problems.
In Section 5, we give a summary of our work and make some concluding remarks.
2. A Quick Tour of MathPDE
In this section we provide a brief tour of MathPDE. taking the example of a one-dimensional diffusion equation. We illustrate how to define the PDE problem, set up the numerical problem using MathPDE. and generate C++ code using MathCode. The package MathPDE is loaded in the usual way.
2.1 Problem Definition