Ordinary And Partial Differential Equation Routines In C, C++, Fortran, Java, Maple, And Matlab

Ordinary And Partial Differential Equation Routines In C, C++, Fortran, Java, Maple,...

(Parte 1 de 5)

Ordinary and Partial Differential

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

A CRC Press Company Boca Raton London New York Washington, D.C.

Ordinary and Partial Differential

H.J. Lee and W.E. Schiesser

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

Java is a registered trademark of Sun Microsystems, Inc. Maple is a registered trademark of Waterloo Maple, Inc. MATLAB is a registered trademark of The MathWorks, Inc. For product information, please contact:

The MathWorks, Inc. 3 Apple Hill Drive Natick, MA 01760-2098 Tel.: 508-647-7000 Fax: 508-647-7001 e-mail: info@mathworks.com Web: w.mathworks.com http://www.mathworks.com/

This book contains information obtained from authentic and highly regarded sources. Reprinted material is quoted with permission, and sources are indicated. A wide variety of references are listed. Reasonable efforts have been made to publish reliable data and information, but the author and the publisher cannot assume responsibility for the validity of all materials or for the consequences of their use.

Neither this book nor any part may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, microfilming, and recording, or by any information storage or retrieval system, without prior permission in writing from the publisher.

The consent of CRC Press LLC does not extend to copying for general distribution, for promotion, for creating new works, or for resale. Specific permission must be obtained in writing from CRC Press LLC for such copying.

Direct all inquiries to CRC Press LLC, 2000 N.W. Corporate Blvd., Boca Raton, Florida 33431.

Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are used only for identification and explanation, without intent to infringe.

Visit the CRC Press Web site at w.crcpress.com

© 2004 by Chapman & Hall/CRC

No claim to original U.S. Government works

International Standard Book Number 1-58488-423-1

Library of Congress Card Number 2003055809

Printed in the United States of America 1 2 3 4 5 6 7 8 9 0 Printed on acid-free paper

Library of Congress Cataloging-in-Publication Data

Lee, H. J. (Hyun Jin)

Ordinary and partial differential equation routines in C, C++, Fortran, Java, Maple, and

MATLAB / H.J. Lee and W.E. Schiesser. p. cm.

Includes bibliographical references and index. ISBN 1-58488-423-1 (alk. paper) 1. Differential equations—Data processing. 2. Differential equations, Partial—Data processing. I. Schiesser, W. E. I. Title.

QA371.5.D37L44 2003

C231 disclaimer.fm Page 1 Friday, October 17, 2003 9:28 AM

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC


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 anODE/PDE-basedmathematicalmodel,findingroutineswritteninafamiliar 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,wehaveillustratedtheapplicationoftheODEintegrationroutines 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.

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

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 formofsoftware.Theuseofthesoftwareisillustratedthroughasmallnumber 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.

Wehope this format will allow the reader to understand the basic elements ofODE/PDEintegration,andthenproceedexpeditiouslytoanumericalsolution 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 Press Web 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

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC


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 1x1 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 2x2 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

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

5.5 Programming in Java 5.6 Programming in Maple

AppendixA Embedded Runge Kutta Pairs AppendixB Integrals from ODEs

AppendixC Stiff ODE Integration

C.1 The BDF Formulas Applied to the 2x2 ODE System C.2 MATLAB Program for the Solution of the 2x2 ODE System

C.3 MATLAB Program for the Solution of the 2x2 ODE System Using ode23s and ode15s

AppendixD Alternative Forms of ODEs AppendixE Spatial p Refinement AppendixF Testing ODE/PDE Codes

Copyright © 2004 by Chapman & Hall/CRC

(C) 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

1 Some Basics of ODE Integration

The central topic of this book is the programming and use of a set of library routines for the numerical solution (integration) of systems of initial value ordinary differential equations (ODEs). We start by reviewing some of the basic concepts of ODEs, including methods of integration, that are the mathematical foundation for an understanding of the ODE integration routines.

1.1 General Initial Value ODE Problem The general problem for a single initial-value ODE is simply stated as

where y = dependent variable

t0 = initial value of the independent variable y0 = initial value of the dependent variable

Equations 1.1 and 1.2 will be termed a 1x1p roblem (one equation in one unknown). The solution of this 1x1p roblem is the dependent variable as a function of the independent variable, y(t) (this function when substituted into Equations 1.1 and 1.2 satisfies these equations). This solution may be a mathematical function, termed an analytical solution.

Copyright © 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

Toillustrate these ideas, we consider the 1x1p roblem, from Braun1 (which will be discussed subsequently in more detail)

where λ and α are positive constants.

Equation 1.3 is termed a first-order, linear, ordinary differential equation with variable coefficients. These terms are explained below.

Term Explanation

Differential equation Equation 1.3 has a derivative dy/dt = f (y, t) = λe−αty

Ordinary Equation 1.3 has only one independent variable, t,s o that the derivative dy/dt is a total or ordinary derivative

First-order The highest-order derivative is first order (dy/dt is first order)

Linear y and its derivative dy/dt are to the first power; thus,

Equation 1.3 is also termed first degree (do not confuse order and degree)

Variable coefficient The coefficient e−αt is a function of the independent variable, t (if it were a function of the dependent variable, y, Equation 1.3 would be nonlinear or not first degree)

The analytical solution to Equations 1.3 and 1.4 is from Braun:1

where exp(x) = ex. Equation 1.5 is easily verified as the solution to Equations 1.3 and 1.4 by substitution in these equations:

Terms in Substitution of Equation 1.5 Equations 1.3 and 1.4 in Equations 1.3 and 1.4 dy dt y0 exp

thus confirming Equation 1.5 satisfies Equations 1.3 and 1.4.

Copyright © 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

As an example of a nxn problem (n ODEs in n unknowns), we will also subsequently consider in detail the 2x2 system

The solution to Equations 1.6 is again the dependent variables, y1,y 2,a sa function of the independent variable, t. Since Equations 1.6 are linear, constant coefficient ODEs, their solution is easily derived, e.g., by assuming exponential functionsint orbytheLaplacetransform.Ifweassumeexponentialfunctions

where c1,c2, and λ are constants to be determined, substitution of Equations 1.7 in Equations 1.6 gives

Cancellation of eλt gives a system of algebraic equations (this is the reason assumingexponentialsolutionsworksinthecaseoflinear,constantcoefficient ODEs)

Equations 1.8 are the 2x2 case of the linear algebraic eigenvalue problem

Copyright © 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC

n = 2 for Equations 1.8, and we use a bold faced symbol for a matrix or a vector. The preceding matrices and vectors are

A nxn coefficient matrix I nxn identity matrix c nx1 solution vector 0 nx1 zero vector

The reader should confirm that the matrices and vectors in Equation 1.9 have the correct dimensions for all of the indicated operations (e.g., matrix additions, matrix-vector multiples).

Note that Equation 1.9 is a linear, homogeneous algebraic system (homogeneousmeansthattheright-handside(RHS)isthezerovector).Thus,Equation 1.9, or its 2x2 counterpart, Equations 1.8, will have nontrivial solutions (c = 0) if and only if (iff) the determinant of the coefficient matrix is zero, i.e.,

Equation 1.10 is the characteristic equation for Equation 1.9 (note that it is a scalar equation). The values of λ that satisfy Equation 1.10 are the eigenvalues of Equation 1.9. For the 2x2p roblem of Equations 1.8, Equation 1.10 is

Equation 1.1 is the characteristic equation or characteristic polynomial for Equations 1.8; note that since Equations 1.8 are a 2x2 linear homogeneous algebraic system, the characteristic equation (Equation 1.1) is a second-order

Copyright © 2004 by Chapman & Hall/CRC Copyright 2004 by Chapman & Hall/CRC polynomial. Similarly, since Equation 1.9 is a nxn linear homogeneous algebraic system, its characteristic equation is a nth-order polynomial. Equation 1.1 can be factored by the quadratic formula

Thus, as expected, the 2x2 system of Equations 1.8 has two eigenvalues. In general, the nxn algebraic system, Equation 1.9, will have n eigenval-

ues, λ1, λ2,, λn (which may be real or complex conjugates, distinct or


Since Equations 1.6 are linear constant coefficient ODEs, their general solution will be a linear combination of exponential functions, one for each eigenvalue

Equations 1.13 have four constants which occur in pairs, one pair for each eigenvalue. Thus, the pair [c11 c21]T is the eigenvector for eigenvalue λ1 while

[c12 c22]T is the eigenvector for eigenvalue λ2.I n general, the nxn system of Equation 1.9 will have a nx1 eigenvector for each of its n eigenvalues. Note that the naming convention for any constant in an eigenvector, cij ,i st he ith constant for the jth eigenvalue. We can restate the two eigenvectors for

(Parte 1 de 5)