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

528 Pages·2006·1.71 MB·English

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®, and MATLAB®
H.J. Lee and W.E. Schiesser
CHAPMAN & HALL/CRC
A CRC Press Company
Boca Raton London New York Washington, D.C.
C231 disclaimer.fm Page 1 Friday, October 17, 2003 9:28 AM
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.: [phone]
Fax: [phone]
e-mail: [email protected]
Web: www.mathworks.com http://www.mathworks.com/
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. II. Title.
QA371.5.D37L44 2003
515¢.352¢0285—dc22 2003055809
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 www.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
Preface
Initial value ordinary differential equations (ODEs) and partial differential
equations(PDEs)areamongthemostwidelyusedformsofmathematicsin
science and engineering. However, insights from ODE/PDE-based models
arerealizedonlywhensolutionstotheequationsareproducedwithaccept-
ableaccuracyandwithreasonableeffort.
Most ODE/PDE models are complicated enough (e.g., sets of simultane-
ousnonlinearequations)toprecludeanalyticalmethodsofsolution;instead,
numericalmethodsmustbeused,whichisthecentraltopicofthisbook.
The calculation of a numerical solution usually requires that well-
establishednumericalintegrationalgorithmsareimplementedinqualityli-
braryroutines.Thelibraryroutinesinturncanbecoded(programmed)ina
varietyofprogramminglanguages.Typically,forascientistorengineerwith
anODE/PDE-basedmathematicalmodel,findingroutineswritteninafamil-
iarlanguagecanbeademandingrequirement,andperhapsevenimpossible
(ifsuchroutinesdonotexist).
The purpose of this book, therefore, is to provide a set of ODE/PDE in-
tegration routines written in six widely accepted and used languages. Our
intention is to facilitate ODE/PDE-based analysis by using the library rou-
tines to compute reliable numerical solutions to the ODE/PDE system of
interest.
However,theintegrationofODE/PDEsisalargesubject,andtokeepthis
discussiontoreasonablelength,wehavelimitedtheselectionofalgorithms
andtheassociatedroutines.Specifically,weconcentrateonexplicit(nonstiff)
Runge Kutta (RK) embedded pairs. Within this setting, we have provided
integratorsthatarebothfixedstepandvariablestep;thelatteracceptauser-
specifiederrortoleranceandattempttocomputeasolutiontothisrequired
accuracy.ThediscussionofODEintegrationincludestruncationerrormoni-
toringandcontrol,handprefinement,stabilityandstiffness,andexplicitand
implicitalgorithms.Extensionstostiffsystemsarealsodiscussedandillus-
tratedthroughanODEapplication;however,adetailedpresentationofstiff
(implicit) algorithms and associated software in six languages was judged
impracticalforabookofreasonablelength.
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
solvednumericallybyoneoftheODEroutines.
Thus,wehaveattemptedtoprovidethereaderwithasetofcomputational
tools for the convenient solution of ODE/PDE models when programming
inanyofthesixlanguages.Thediscussionisintroductorywithlimitedmath-
ematicaldetails.Rather,werelyonnumericalresultstoillustratesomebasic
mathematicalproperties,andweavoiddetailedmathematicalanalysis(e.g.,
theoremsandproofs),whichmaynotreallyprovidemuchassistanceinthe
actualcalculationofnumericalsolutionstoODE/PDEproblems.
Instead, we have attempted to provide useful computational tools in the
formofsoftware.Theuseofthesoftwareisillustratedthroughasmallnumber
ofODE/PDEapplications;ineachcase,thecompletecodeisfirstpresented,
and then its components are discussed in detail, with particular reference
to the concepts of integration, e.g., stability, error monitoring, and control.
Sincethealgorithmsandtheassociatedsoftwarehavelimitations(asdoall
algorithms and software), we have tried to point out these limitations, and
makesuggestionsforadditionalmethodsthatcouldbeeffective.
Also,wehaveintentionallyavoidedusingfeaturesspecifictoaparticular
language,e.g.,sparseutilities,object-orientedprogramming.Rather,wehave
emphasizedthecommonalityoftheprogramminginthesixlanguages,and
therebyillustratehowscientificcomputationcanbedoneinanyofthelan-
guages.Ofcourse,language-specificfeaturescanbeaddedtothesourcecode
thatisprovided.
Wehopethisformatwillallowthereadertounderstandthebasicelements
ofODE/PDEintegration,andthenproceedexpeditiouslytoanumericalsolu-
tionoftheODE/PDEsystemofinterest.Theapplicationsdiscussedindetail,
two in ODEs and two in PDEs, can be used as a starting point (i.e., as tem-
plates)forthedevelopmentofaspectrumofnewapplications.
We welcome comments and questions about how we might be of assis-
tance([email protected]).Informationforacquiring(gratis)allthe
source code in this book is available from http://www.lehigh.edu/˜wes1/
wes1.html.Additionalinformationaboutthebookandsoftwareisavailable
fromtheCRCPressWebsite,http://www.crcpress.com.
Dr. Fred Chapman provided expert assistance with the Maple program-
ming. We note with sadness the passing of Jaeson Lee, father of H. J. Lee,
duringthecompletionofH.J.Lee’sgraduatestudiesatLehighUniversity.
H. J. Lee
W. E. Schiesser
Bethlehem, PA
Contents
1 SomeBasicsofODEIntegration...................................1
1.1 GeneralInitialValueODEProblem...........................1
1.2 OriginofODEIntegratorsintheTaylorSeries................7
1.3 TheRungeKuttaMethod...................................13
1.4 AccuracyofRKMethods ................................... 18
1.5 EmbeddedRKAlgorithms..................................51
1.6 LibraryODEIntegrators....................................72
1.7 StabilityofRKMethods.....................................95
2 Solutionofa1x1ODESystem...................................107
2.1 ProgramminginMATLAB.................................107
2.2 ProgramminginC.........................................143
2.3 ProgramminginC++......................................174
2.4 ProgramminginFortran...................................206
2.5 ProgramminginJava......................................232
2.6 ProgramminginMaple....................................263
3 Solutionofa2x2ODESystem...................................291
3.1 ProgramminginMATLAB.................................291
3.2 ProgramminginC.........................................298
3.3 ProgramminginC++......................................306
3.4 ProgramminginFortran...................................314
3.5 ProgramminginJava......................................321
3.6 ProgramminginMaple....................................329
4 SolutionofaLinearPDE.........................................339
4.1 ProgramminginMATLAB.................................344
4.2 ProgramminginC.........................................359
4.3 ProgramminginC++......................................366
4.4 ProgramminginFortran...................................374
4.5 ProgramminginJava......................................380
4.6 ProgramminginMaple....................................387
5 SolutionofaNonlinearPDE....................................397
5.1 ProgramminginMATLAB.................................402
5.2 ProgramminginC.........................................411
5.3 ProgramminginC++......................................418
5.4 ProgramminginFortran...................................425
5.5 ProgramminginJava......................................431
5.6 ProgramminginMaple....................................438
AppendixA EmbeddedRungeKuttaPairs........................451
AppendixB IntegralsfromODEs.................................459
AppendixC StiffODEIntegration................................465
C.1 TheBDFFormulasAppliedtothe2x2ODESystem........465
C.2 MATLABProgramfortheSolutionofthe
2x2ODESystem...........................................468
C.3 MATLABProgramfortheSolutionofthe2x2ODESystem
Usingode23s andode15s ..................................477
AppendixD AlternativeFormsofODEs..........................489
AppendixE SpatialpRefinement.................................493
AppendixF TestingODE/PDECodes.............................503
Index......................................................................511
1
Some Basics of ODE Integration
The central topic of this book is the programming and use of a set of li-
brary 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 GeneralInitialValueODEProblem
Thegeneralproblemforasingleinitial-valueODEissimplystatedas
dy
= f(y,t), y(t )= y (1.1)(1.2)
0 0
dt
where
y= dependentvariable
t = independentvariable
f(y,t)= derivativefunction
t = initialvalueoftheindependentvariable
0
y = initialvalueofthedependentvariable
0
Equations1.1and1.2willbetermeda1x1problem(oneequationinoneun-
known).Thesolutionofthis1x1problemisthedependentvariableasafunction
oftheindependentvariable,y(t)(thisfunctionwhensubstitutedintoEquations
1.1 and 1.2 satisfies these equations). This solution may be a mathematical
function,termedananalyticalsolution.
1