Simulation of ODE/PDE Models with MATLAB�,OCTAVE and SCILAB

Alain Vande Wouwer • Philippe SaucezCarlos Vilas

Simulation of ODE/PDEModels with MATLAB�,OCTAVE and SCILAB

Scientific and Engineering Applications

123

Alain Vande WouwerService d’AutomatiqueUniversité de MonsMonsBelgium

Philippe SaucezService de Mathématique et Recherche

OpérationnelleUniversité de MonsMonsBelgium

Carlos Vilas(Bio)Process Engineering GroupInstituto de Investigaciones Marinas

(CSIC)VigoSpain

ISBN 978-3-319-06789-6 ISBN 978-3-319-06790-2 (eBook)DOI 10.1007/978-3-319-06790-2Springer Cham Heidelberg New York Dordrecht London

Library of Congress Control Number: 2014939391

� Springer International Publishing Switzerland 2014This work is subject to copyright. All rights are reserved by the Publisher, whether the whole or part ofthe material is concerned, specifically the rights of translation, reprinting, reuse of illustrations,recitation, broadcasting, reproduction on microfilms or in any other physical way, and transmission orinformation storage and retrieval, electronic adaptation, computer software, or by similar or dissimilarmethodology now known or hereafter developed. Exempted from this legal reservation are briefexcerpts in connection with reviews or scholarly analysis or material supplied specifically for thepurpose of being entered and executed on a computer system, for exclusive use by the purchaser of thework. Duplication of this publication or parts thereof is permitted only under the provisions ofthe Copyright Law of the Publisher’s location, in its current version, and permission for use mustalways be obtained from Springer. Permissions for use may be obtained through RightsLink at theCopyright Clearance Center. Violations are liable to prosecution under the respective Copyright Law.The use of general descriptive names, registered names, trademarks, service marks, etc. in thispublication does not imply, even in the absence of a specific statement, that such names are exemptfrom the relevant protective laws and regulations and therefore free for general use.While the advice and information in this book are believed to be true and accurate at the date ofpublication, neither the authors nor the editors nor the publisher can accept any legal responsibility forany errors or omissions that may be made. The publisher makes no warranty, express or implied, withrespect to the material contained herein.

Printed on acid-free paper

Springer is part of Springer Science+Business Media (www.springer.com)

Additional material to this book can be downloaded from http://extras.springer.com

To V. L. E.

A.

To Daniel Druart

Philippe Saucez

Foreword

With the availability of computers of increasing power and lower cost, computer-based modeling is now a widespread approach to the analysis of complex scientificand engineering systems. First-principles models and numerical simulation can beused to investigate system dynamics, to perform sensitivity analysis, to estimateunknown parameters or state variables, and to design model-based controlschemes. However, due to the usual gap between research and common practice,scientists and engineers still often resort to conventional tools (e.g. low-orderapproximate solutions), and do not make use of the full array of readily availablenumerical methods.

Many systems from science and engineering are distributed parameter systems,i.e., systems characterized by state variables (or dependent variables) in two or morecoordinates (or independent variables). Time and space are the most frequentcombination of independent variables, as is the case of the following (time-varying,transient, or unsteady state) examples:

• temperature profiles in a heat exchanger• concentration profiles in a sorptive packed column• temperature and concentration profiles in a tubular reactor• car density along a highway• deflection profile of a beam subject to external forces• shape and velocity of a water wave• distribution of a disease in a population (spread of epidemics)

but other combinations of independent variables are possible as well. Forinstance, time and individual size (or another characteristic such as age) occur inpopulation models used in ecology, or to describe some important industrialprocesses such as polymerization, crystallization, or material grinding. In thesemodels, space can also be required to represent the distribution of individuals (ofvarious sizes) in a spatial region or in a nonhomogeneous reactor medium (due tononideal mixing conditions in a batch reactor, or to continuous operation in atubular reactor).

The preceding examples show that there exists a great variety of distributedparameter systems, arising from different areas of science and engineering, whichare characterized by time-varying distributions of dependent variables. In view of

vii

the system complexity, a mathematical model, i.e., a mathematical description ofthe physical (chemical, biological, mechanical, electrical, etc.) phenomena takingplace in the system, is often a prerequisite to system analysis and control. Such amodel consists of partial differential equations (PDEs), boundary conditions (BCs),and initial conditions (ICs) describing the evolution of the state variables. Inaddition, distributed parameter systems can interact with lumped parameter sys-tems, whose state variables are described by ordinary differential equations(ODEs), and supplementary algebraic equations (AEs) can be used to expressphenomena such as thermodynamic equilibria, heat and mass transfer, and reactionkinetics (combinations of AEs and ODEs are also frequently termed differential-algebraic equations, or DAEs). Hence, a distributed parameter model is usuallydescribed by a mixed set of nonlinear AE/ODE/PDEs or PDAEs. Most PDAEsmodels are derived from first principles, i.e., conservation of mass, energy, andmomentum, and are given in a state space representation which is the basis forsystem analysis.

This book is dedicated to numerical simulation of distributed parameter systemsdescribed by mixed systems of PDAEs. Special attention is paid to the numericalmethod of lines (MOL), a popular approach to the solution of time-dependentPDEs, which proceeds in two basic steps. First, spatial derivatives are approxi-mated using finite difference, element, or volume approximations. Second, theresulting system of semi-discrete (discrete in space continuous in time) equationsis integrated in time using an available solver. Besides conventional finite dif-ference, element, and volume techniques, which are of high practical value, moreadvanced spatial approximation techniques are examined in some detail, includingfinite element and finite volume approaches.

Although the MOL has attracted considerable attention and several general-purpose libraries or specific software packages have been developed, there is still aneed for basic, introductory, yet efficient, tools for the simulation of distributedparameter systems, i.e., software tools that can be easily used by practicing sci-entists and engineers, and that provide up-to-date numerical algorithms.

Consequently, a MOL toolbox has been developed within MATLAB/OCTAVE/SCILAB. These environments conveniently demonstrate the usefulnessand effectiveness of the above-mentioned techniques and provide high-qualitymathematical libraries, e.g., ODE solvers that can be used advantageously incombination with the proposed toolbox. In addition to a set of spatial approxi-mations and time integrators, this toolbox includes a library of applicationexamples, in specific areas, which can serve as templates for developing newprograms. The idea here is that a simple code template is often more compre-hensible and flexible than a software environment with specific user interfaces.This way, various problems including coupled systems of AEs, ODEs, and PDEsin one or more spatial dimensions can easily be developed, modified, and tested.

This text, which provides an introduction to some advanced computationaltechniques for dynamic system simulation, is suitable as a final year undergraduatecourse or at the graduate level. It can also be used for self-study by practicingscientists and engineers.

viii Foreword

Preface

Our initial objective in developing this book was to report on our experience innumerical techniques for solving partial differential equation problems, usingsimple programming environments such as MATLAB, OCTAVE, or SCILAB.Computational tools and numerical simulation are particularly important forengineers, but the specialized literature on numerical analysis is sometimes toodense or too difficult to explore due to a gap in the mathematical background. Thisbook is intended to provide an accessible introduction to the field of dynamicsimulation, with emphasis on practical methods, yet including a few advancedtopics that find an increasing number of engineering applications. At the origin ofthis book project, some years ago, we were teaming up with Bill Schiesser (LehighUniversity) with whom we had completed a collective book on Adaptive Methodof Lines. Unfortunately, this previous work had taken too much of our energy, andthe project faded away, at least for the time being.

Time passed, and the book idea got a revival at the time of the post-doctoralstay of Carlos in the Control Group of the University of Mons. Carlos had justachieved a doctoral work at the University of Vigo, involving partial differentialequation models, finite element techniques, and the proper orthogonal decompo-sition, ingredients, which all were excellent complements to our backgroundmaterial.

The three of us then decided to join our forces to develop a manuscript with anemphasis on practical implementation of numerical methods for ordinary andpartial differential equation problems, mixing introductory material to numericalmethods, a variety of illustrative examples from science and engineering, and acollection of codes that can be reused for the fast prototyping of new simulationcodes.

All in one, the book material is based on past research activities, literaturereview, as well as courses taught at the University of Mons, especially introductorynumerical analysis courses for engineering students. As a complement to the text, awebsite (www.matmol.org) has been set up to provide a convenient platform fordownloading codes and method tutorials.

Writing a book is definitely a delicate exercise, and we would like to seize thisopportunity to thank Bill for his support in the initial phase of this project. Many ofhis insightful suggestions are still present in the current manuscript, which hasdefinitely benefited from our discussions and nice collaboration.

ix

Of course, we also would like to express our gratitude to our colleagues atUMONS and at IMM-CSIC (Vigo), and particularly Marcel, Christine, Antonio,Julio, Eva, and Míriam, and all the former and current research teams, for the niceworking environment and for the research work achieved together, which was asource of inspiration in developing this material. We are also grateful to a numberof colleagues in other universities for the nice collaboration, fruitful exchanges atseveral conferences, or insightful comments on some of our developments:Michael Zeitz, Achim Kienle, Paul Zegeling, Gerd Steinebach, Keith Miller, SkipThompson, Larry Shampine, Ken Anselmo, Filip Logist, to just name a few.

In addition, we acknowledge the support of the Belgian Science Policy Office(BELSPO), which through the Interuniversity Attraction Program DynamicalSystems, Control and Optimization (DYSCO) supported part of this research workand made possible several mutual visits and research stays at both institutions(UMONS and IIM-CSIC) over the past several years.

Finally, we would like to stress the excellent collaboration with Oliver Jackson,Editor in Engineering at Springer, with whom we had the initial contact for thepublication of this manuscript and who guided us in the review process andselection of a suitable book series. In the same way, we would like to thankCharlotte Cross, Senior editorial assistant at Springer, for the timely publicationprocess, and for her help and patience in the difficult manuscript completion phase.

Mons, March 2014 Alain Vande WouwerVigo Philippe Saucez

Carlos Vilas

x Preface

Contents

1 An Introductory Tour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Some ODE Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 An ODE/DAE Application . . . . . . . . . . . . . . . . . . . . . . . . . . 321.3 A PDE Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 371.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

2 More on ODE Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.1 A Basic Fixed Step ODE Integrator . . . . . . . . . . . . . . . . . . . . 452.2 A Basic Variable-Step Nonstiff ODE Integrator. . . . . . . . . . . . 512.3 A Basic Variable Step Implicit ODE Integrator . . . . . . . . . . . . 672.4 MATLAB ODE Suite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 852.5 Some Additional ODE Applications . . . . . . . . . . . . . . . . . . . . 86

2.5.1 Spruce Budworm Dynamics . . . . . . . . . . . . . . . . . . . 862.5.2 Liming to Remediate Acid Rain . . . . . . . . . . . . . . . . 93

2.6 On the Use of SCILAB and OCTAVE . . . . . . . . . . . . . . . . . . 1092.7 How to Use Your Favorite Solvers in MATLAB? . . . . . . . . . . 115

2.7.1 A Simple Example: Matrix Multiplication. . . . . . . . . . 1172.7.2 MEX-Files for ODE Solvers . . . . . . . . . . . . . . . . . . . 122

2.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

3 Finite Differences and the Method of Lines . . . . . . . . . . . . . . . . . . 1253.1 Basic Finite Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1263.2 Basic MOL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1273.3 Numerical Stability: Von Neumann and the Matrix Methods. . . 1293.4 Numerical Study of the Advection Equation . . . . . . . . . . . . . . 1363.5 Numerical Study of the Advection-Diffusion Equation . . . . . . . 1423.6 Numerical Study of the Advection-Diffusion-Reaction

Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1503.7 Is it Possible to Enhance Stability? . . . . . . . . . . . . . . . . . . . . 1513.8 Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1533.9 Accuracy and the Concept of Differentiation Matrices . . . . . . . 157

xi

3.10 Various Ways of Translating the Boundary Conditions. . . . . . . 1673.10.1 Elimination of Unknown Variables . . . . . . . . . . . . . . 1703.10.2 Fictitious Nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1743.10.3 Solving Algebraic Equations . . . . . . . . . . . . . . . . . . . 1763.10.4 Tricks Inspired by the Previous Methods . . . . . . . . . . 1793.10.5 An Illustrative Example (with Several

Boundary Conditions). . . . . . . . . . . . . . . . . . . . . . . . 1813.11 Computing the Jacobian Matrix of the ODE System . . . . . . . . 1903.12 Solving PDEs Using SCILAB and OCTAVE . . . . . . . . . . . . . 1973.13 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201

4 Finite Elements and Spectral Methods . . . . . . . . . . . . . . . . . . . . . 2034.1 The Methods of Weighted Residuals . . . . . . . . . . . . . . . . . . . 211

4.1.1 Interior Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2134.1.2 Boundary Method . . . . . . . . . . . . . . . . . . . . . . . . . . 2134.1.3 Mixed Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2134.1.4 Galerkin Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 2144.1.5 Collocation Method . . . . . . . . . . . . . . . . . . . . . . . . . 2144.1.6 Orthogonal Collocation Method . . . . . . . . . . . . . . . . . 214

4.2 The Basics of the Finite Element Method . . . . . . . . . . . . . . . . 2154.3 Galerkin Method Over Linear Lagrangian Elements . . . . . . . . . 216

4.3.1 LHS of the Weighted Residual Solution . . . . . . . . . . . 2184.3.2 First Term in the RHS of the Weighted

Residual Solution. . . . . . . . . . . . . . . . . . . . . . . . . . . 2194.3.3 Second Term in the RHS of the Weighted

Residual Solution. . . . . . . . . . . . . . . . . . . . . . . . . . . 2204.3.4 Third Term in the RHS of the Weighted

Residual Solution. . . . . . . . . . . . . . . . . . . . . . . . . . . 2214.3.5 Fourth Term in the RHS of the Weighted

Residual Solution. . . . . . . . . . . . . . . . . . . . . . . . . . . 2244.4 Galerkin Method Over Linear Lagrangian Elements:

Contribution of the Boundary Conditions . . . . . . . . . . . . . . . . 2274.4.1 Dirichlet Boundary Conditions. . . . . . . . . . . . . . . . . . 2284.4.2 Neumann Boundary Conditions . . . . . . . . . . . . . . . . . 228

4.5 The Finite Element Method in Action . . . . . . . . . . . . . . . . . . 2294.6 The Finite Element Method Applied to Systems of PDEs . . . . . 2354.7 Galerkin Method Over Hermitian Elements. . . . . . . . . . . . . . . 237

4.7.1 LHS Term of the Weighted Residual Solution. . . . . . . 2374.7.2 First and Second Terms of the RHS Term

of the Weighted Residual Solution . . . . . . . . . . . . . . . 2394.7.3 Third Term of the RHS Term

of the Weighted Residual Solution . . . . . . . . . . . . . . . 241

xii Contents

4.7.4 Fourth Term of the RHS Termof the Weighted Residual Solution . . . . . . . . . . . . . . . 242

4.7.5 Galerkin Method Over Hermitian Elements:Contribution of the Boundary Conditions . . . . . . . . . . 243

4.8 An Illustrative Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2454.9 The Orthogonal Collocation Method. . . . . . . . . . . . . . . . . . . . 249

4.9.1 LHS Term of the Collocation Residual Equation . . . . . 2514.9.2 First Three Terms of Collocation

Residual Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 2524.9.3 Fourth Term of the RHS of the Collocation

Residual Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 2544.9.4 Contribution of the Boundary Conditions . . . . . . . . . . 2554.9.5 A Brief Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . 256

4.10 Chebyshev Collocation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2564.11 The Proper Orthogonal Decomposition . . . . . . . . . . . . . . . . . . 262

4.11.1 The Method of Snapshots . . . . . . . . . . . . . . . . . . . . . 2654.11.2 Example: The Heat Equation. . . . . . . . . . . . . . . . . . . 2684.11.3 Example: The Brusselator . . . . . . . . . . . . . . . . . . . . . 273

4.12 On the Use of SCILAB and OCTAVE . . . . . . . . . . . . . . . . . . 2774.13 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282

5 How to Handle Steep Moving Fronts? . . . . . . . . . . . . . . . . . . . . . 2855.1 Conservation Laws. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2865.2 The Methods of Characteristics and of Vanishing Viscosity . . . 2885.3 Transformation-Based Methods . . . . . . . . . . . . . . . . . . . . . . . 2935.4 Upwind Finite Difference and Finite Volume Schemes. . . . . . . 2955.5 A Divide and Conquer Approach . . . . . . . . . . . . . . . . . . . . . . 2985.6 Finite Volume Methods and Slope Limiters . . . . . . . . . . . . . . 3035.7 Grid Refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3215.8 An Additional PDE Application. . . . . . . . . . . . . . . . . . . . . . . 3315.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335

6 Two Dimensional and Time Varying Spatial Domains . . . . . . . . . . 3396.1 Solution of Partial Differential Equations in More

than 1D Using Finite Differences. . . . . . . . . . . . . . . . . . . . . . 3396.1.1 The Heat Equation on a Rectangle . . . . . . . . . . . . . . . 3406.1.2 Graetz Problem with Constant Wall Temperature . . . . 3446.1.3 Tubular Chemical Reactor. . . . . . . . . . . . . . . . . . . . . 3486.1.4 Heat Equation on a Convex Quadrilateral . . . . . . . . . . 3526.1.5 A Convection-Diffusion Equation on a Square . . . . . . 3576.1.6 Burgers Equation on a Square . . . . . . . . . . . . . . . . . . 360

Contents xiii

6.2 Solution of 2D PDEs Using Finite Element Techniques . . . . . . 3616.2.1 FitzHugh-Nagumo’s Model . . . . . . . . . . . . . . . . . . . . 3646.2.2 Reduced-Order Model for FitzHugh-Nagumo

Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3826.3 Solution of PDEs on Time-Varying Domains . . . . . . . . . . . . . 388

6.3.1 The Freeze-Drying Model . . . . . . . . . . . . . . . . . . . . . 3896.3.2 The Landau Transform . . . . . . . . . . . . . . . . . . . . . . . 3916.3.3 The Finite Element Representation. . . . . . . . . . . . . . . 394

6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 399References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 400

Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403

xiv Contents

Acronyms

AE Algebraic EquationALE Algebraic Lagrangian EulerianBC Boundary ConditionsBDF Backward Differentiation FormulaDAE Differential Algebraic EquationFD Finite DifferencesFEM Finite Element MethodIBVP Initial Boundary Value ProblemIC Initial ConditionsIVP Initial Value ProblemLHS Left Hand SideMOL Method of LinesMC Monotonized CentralODE Ordinary Differential EquationPDAE Partial Differential Algebraic EquationPDE Partial Differential EquationPOD Proper Orthogonal DecompositionRHS Right-Hand SideRK Runge–KuttaWRM Weighted Residuals Methods

xv