Ordinary and Partial Differential Equation Routines in C, C++, Fortran, Java®, Maple®, and MATLAB®

Initial value ordinary differential equations (ODEs) and partial differential equations (PDEs) are among the most widely used forms of mathematics in science and engineering. However, insights from ODE/PDE-based models are realized only when solutions to the equations are produced with acceptable accuracy and with reasonable effort.

Most ODE/PDE models are complicated enough (e.g., sets of simultaneous nonlinear equations) to preclude analytical methods of solution; instead, numerical methods must be used, which is the central topic of this book. The calculation of a numerical solution usually requires that wellestablished numerical integration algorithms are implemented in quality library routines. The library routines in turn can be coded (programmed) in a variety of programming languages. Typically, for a scientist or engineer with an ODE/PDE- based mathematical model, finding routines written in a familiar language can be a demanding requirement, and perhaps even impossible (if such routines do not exist).

The purpose of this book, therefore, is to provide a set of ODE/PDE integration routines written in six widely accepted and used languages. Our intention is to facilitate ODE/PDE-based analysis by using the library routines to compute reliable numerical solutions to the ODE/PDE system of interest.

However, the integration of ODE/PDEs is a large subject, and to keep this discussion to reasonable length, we have limited the selection of algorithms and the associated routines. Specifically, we concentrate on explicit (nonstiff) Runge Kutta (RK) embedded pairs. Within this setting, we have provided integrators that are both fixed step and variable step; the latter accept a userspecified error tolerance and attempt to compute a solution to this required accuracy. The discussion of ODE integration includes truncation error monitoring and control, h and p refinement, stability and stiffness, and explicit and implicit algorithms. Extensions to stiff systems are also discussed and illustrated through an ODE application; however, a detailed presentation of stiff (implicit) algorithms and associated software in six languages was judged impractical for a book of reasonable length.

Further, we have illustrated the application of the ODE integration routines to PDEs through the method of lines (MOL). Briefly, the spatial (boundary value) derivatives of the PDEs are approximated algebraically, typically by
finite differences (FDs); the resulting system of initial-value ODEs is then solved numerically by one of the ODE routines.

Thus, we have attempted to provide the reader with a set of computational tools for the convenient solution of ODE/PDE models when programming in any of the six languages. The discussion is introductory with limited mathematical details. Rather, we rely on numerical results to illustrate some basic mathematical properties, and we avoid detailed mathematical analysis (e.g., theorems and proofs), which may not really provide much assistance in the actual calculation of numerical solutions to ODE/PDE problems. Instead, we have attempted to provide useful computational tools in the form of software. The use of the software is illustrated through a small number of ODE/PDE applications; in each case, the complete code is first presented, and then its components are discussed in detail, with particular reference to the concepts of integration, e.g., stability, error monitoring, and control.

Since the algorithms and the associated software have limitations (as do all algorithms and software), we have tried to point out these limitations, and make suggestions for additional methods that could be effective.
Also, we have intentionally avoided using features specific to a particular language, e.g., sparse utilities, object-oriented programming. Rather, we have emphasized the commonality of the programming in the six languages, and thereby illustrate how scientific computation can be done in any of the languages.

Of course, language-specific features can be added to the source code that is provided. We hope this format will allow the reader to understand the basic elements ofODE/PDEintegration, and then proceed expeditiously to a numerical solution of the ODE/PDE system of interest. The applications discussed in detail, two in ODEs and two in PDEs, can be used as a starting point (i.e., as templates) for the development of a spectrum of new applications.

We welcome comments and questions about how we might be of assistance (directed to wes1@lehigh.edu). Information for acquiring (gratis) all the source code in this book is available from http://www.lehigh.edu/˜ wes1/
wes1.html. Additional information about the book and software is available from the CRC PressWeb site, http://www.crcpress.com. Dr. Fred Chapman provided expert assistance with the Maple programming. We note with sadness the passing of Jaeson Lee, father of H. J. Lee, during the completion of H. J. Lee’s graduate studies at Lehigh University.
H. J. Lee
W. E. Schiesser
Bethlehem, PA

Contents

1 Some Basics of ODE Integration
1.1 General Initial Value ODE Problem
1.2 Origin of ODE Integrators in the Taylor Series
1.3 The Runge Kutta Method
1.4 Accuracy of RK Methods
1.5 Embedded RK Algorithms
1.6 Library ODE Integrators
1.7 Stability of RK Methods
2 Solution of a 1×1 ODE System
2.1 Programming in MATLAB
2.2 Programming in C
2.3 Programming in C++
2.4 Programming in Fortran
2.5 Programming in Java
2.6 Programming in Maple
3 Solution of a 2×2 ODE System
3.1 Programming in MATLAB
3.2 Programming in C
3.3 Programming in C++
3.4 Programming in Fortran
3.5 Programming in Java
3.6 Programming in Maple
4 Solution of a Linear PDE
4.1 Programming in MATLAB
4.2 Programming in C
4.3 Programming in C++
4.4 Programming in Fortran
4.5 Programming in Java
4.6 Programming in Maple
5 Solution of a Nonlinear PDE
5.1 Programming in MATLAB
5.2 Programming in C
5.3 Programming in C++
5.4 Programming in Fortran

5.5 Programming in Java
5.6 Programming in Maple
Appendix A Embedded Runge Kutta Pairs
Appendix B Integrals from ODEs
Appendix C Stiff ODE Integration
C.1 The BDF Formulas Applied to the 2×2 ODE System
C.2 MATLAB Program for the Solution of the 2×2 ODE System
C.3 MATLAB Program for the Solution of the 2×2 ODE System Using ode23s and ode15s
Appendix D Alternative Forms of ODEs
Appendix E Spatial p Refinement
Appendix F Testing ODE/PDE Codes