ACCURATE PARALLEL INTEGRATION OF LARGE SPARSE SYSTEMS OF DIFFERENTIAL EQUATIONS

Abstract

We describe a MIMD parallel code to solve a general class of ordinary differential equations, with particular emphasis on the large, sparse systems arising from space discretization of systems of parabolic partial differential equations. The main goals of this work are sharp bounds on the accuracy of the computed solution and flexibility of the software. We discuss the sources of error in solving differential equations, and the resulting constraints on time steps. We also discuss the theory of a posteriori error analysis for the Galerkin finite element methods, and its implementation in error control and estimation. The software is designed in a matrix-free fashion, which enables the solver to effectively tackle large sparse systems with minimal memory consumption and an easy and natural transition to MIMD (distributed memory) parallelism. In addition, there is no need for the choice of a particular representation of a sparse matrix. All memory is dynamically allocated, with a new expandable array object used for archiving results. The implicit solution of the discrete equations is carried out by replaceable modules: the nonlinear solver module may be a full Newton scheme or a quasi-Newton; these in turn are implemented with a linear solver, for which we have used both a direct solver and QMR, an iterative (Krylov space) method. Three computations are presented: the Lorenz system, which has dimension three and the discretized versions of the (partial-differential) bistable equation in one and two dimensions. The Lorenz system demonstrates the quality of the error estimation. The discretized bistable examples provide large sparse systems, and our precise error estimation shows, contrary to standard error estimates, that reliable computation is possible for large times.