Numerical methods with matlab implementations and applications pdf

Applied Numerical Methods with MATLAB® for Engineers and Scientists

Fourth Edition

Steven C. Chapra Berger Chair in Computing and Engineering

Tufts University

cha97962_fm_i-xvi.indd 1 07/11/16 8:25 pm

APPLIED NUMERICAL METHODS WITH MATLAB® FOR ENGINEERS AND SCIENTISTS, FOURTH EDITION

Published by McGraw-Hill Education, 2 Penn Plaza, New York, NY 10121. Copyright © 2018 by McGraw-Hill Education. All rights reserved. Printed in the United States of America. Previous editions © 2012, 2008, and 2005. No part of this publication may be reproduced or distributed in any form or by any means, or stored in a database or retrieval system, without the prior written consent of McGraw-Hill Education, including, but not limited to, in any network or other electronic storage or transmission, or broadcast for distance learning.

Some ancillaries, including electronic and print components, may not be available to customers outside the United States.

This book is printed on acid-free paper.

1 2 3 4 5 6 7 8 9 LCR 21 20 19 18 17

ISBN 978-0-07-339796-2 MHID 0-07-339796-2

Chief Product Officer, SVP Products & Markets: G. Scott Virkler Vice President, General Manager, Products & Markets: Marty Lange Vice President, Content Design & Delivery: Betsy Whalen Managing Director: Thomas Timp Brand Manager: Thomas M. Scaife, Ph.D. Director, Product Development: Rose Koos Product Developer: Jolynn Kilburg Marketing Manager: Nick McFadden Market Development Manager: Shannon O’Donnell Digital Product Developer: Joan Weber Director, Content Design & Delivery: Linda Avenarius Program Manager: Lora Neyens Content Project Managers: Jeni McAtee; Sandy Schnee Buyer: Susan K. Culbertson Content Licensing Specialists: Lorraine Buczek Cover Image: Vitalii Nesterchuk/123RF© Compositor: MPS Limited Printer: LSC Communications

All credits appearing on page or at the end of the book are considered to be an extension of the copyright page.

Library of Congress Cataloging-in-Publication Data Chapra, Steven C., author. Applied numerical methods with MATLAB for engineers and scientists / Steven C. Chapra, Berger Chair in Computing and Engineering, Tufts University. Fourth edition. | New York, NY : McGraw-Hill Education, [2018] | Includes bibliographical references and index. LCCN 2016038044 | ISBN 9780073397962 (alk. paper) | ISBN 0073397962 (alk. paper) LCSH: Numerical analysis—Data processing—Textbooks. | Engineering mathematics—Textbooks. | MATLAB—Textbooks. LCC QA297 .C4185 2018 | DDC 518—dc23 LC record available at https://lccn.loc.gov/2016038044 The Internet addresses listed in the text were accurate at the time of publication. The inclusion of a website does not indicate an endorsement by the authors or McGraw-Hill Education, and McGraw-Hill Education does not guarantee the accuracy of the information presented at these sites.

mheducation.com/highered

cha97962_fm_i-xvi.indd 2 07/11/16 8:25 pm

To

My brothers,

John and Bob Chapra

and

Fred Berger (1947–2015)

who I miss as a good friend, a good man.

and a comrade in bringing the light of engineering

to some of world’s darker corners.

cha97962_fm_i-xvi.indd 3 07/11/16 8:25 pm

iv

ABOUT THE AUTHOR

Steve Chapra teaches in the Civil and Environmental Engineering Department at Tufts University, where he holds the Louis Berger Chair in Computing and Engineering. His other books include Numerical Methods for Engineers and Surface Water-Quality Modeling.

Steve received engineering degrees from Manhattan College and the University of Michigan. Before joining the faculty at Tufts, he worked for the Environmental Protection Agency and the National Oceanic and Atmospheric Administration, and taught at Texas A&M University and the University of Colorado. His general research interests focus on surface water-quality modeling and advanced computer applications in environmental engineering.

He has received a number of awards for his scholarly contributions, including the Rudolph Hering Medal, the Meriam/Wiley Distinguished Author Award, and the Chandler- Misener Award. He has also been recognized as the outstanding teacher at Texas A&M University (1986 Tenneco Award), the University of Colorado (1992 Hutchinson Award), and Tufts University (2011 Professor of the Year Award).

Steve was originally drawn to environmental engineering and science because of his love of the outdoors. He is an avid fly fisherman and hiker. An unapologetic nerd, his love affair with computing began when he was first introduced to Fortran programming as an undergraduate in 1966. Today, he feels truly blessed to be able to meld his love of math- ematics, science, and computing with his passion for the natural environment. In addition, he gets the bonus of sharing it with others through his teaching and writing!

Beyond his professional interests, he enjoys art, music (especially classical music, jazz, and bluegrass), and reading history. Despite unfounded rumors to the contrary, he never has, and never will, voluntarily bungee jump or sky dive.

If you would like to contact Steve, or learn more about him, visit his home page at http://engineering.tufts.edu/cee/people/chapra/ or e-mail him at steven.chapra@tufts.edu.

cha97962_fm_i-xvi.indd 4 07/11/16 8:25 pm

v

CONTENTS

About the Author iv

Preface xiv

Part One Modeling, Computers, and Error Analysis 1

1.1 Motivation 1 1.2 Part Organization 2

CHAPTER 1

Mathematical Modeling, Numerical Methods, and Problem Solving 4

1.1 A Simple Mathematical Model 5 1.2 Conservation Laws in Engineering and Science 12 1.3 Numerical Methods Covered in This Book 13 1.4 Case Study: It’s a Real Drag 17 Problems 20

CHAPTER 2

MATLAB Fundamentals 27

2.1 The MATLAB Environment 28 2.2 Assignment 29 2.3 Mathematical Operations 36 2.4 Use of Built-In Functions 39 2.5 Graphics 42 2.6 Other Resources 46 2.7 Case Study: Exploratory Data Analysis 46 Problems 49

CHAPTER 3

Programming with MATLAB 53

3.1 M-Files 54 3.2 Input-Output 61

cha97962_fm_i-xvi.indd 5 07/11/16 8:25 pm

vi CONTENTS

3.3 Structured Programming 65 3.4 Nesting and Indentation 79 3.5 Passing Functions to M-Files 81 3.6 Case Study: Bungee Jumper Velocity 87 Problems 91

CHAPTER 4

Roundoff and Truncation Errors 99

4.1 Errors 100 4.2 Roundoff Errors 106 4.3 Truncation Errors 114 4.4 Total Numerical Error 125 4.5 Blunders, Model Errors, and Data Uncertainty 130 Problems 131

Part Two Roots and Optimization 135

2.1 Overview 135 2.2 Part Organization 136

CHAPTER 5

Roots: Bracketing Methods 138

5.1 Roots in Engineering and Science 139 5.2 Graphical Methods 140 5.3 Bracketing Methods and Initial Guesses 141 5.4 Bisection 146 5.5 False Position 152 5.6 Case Study: Greenhouse Gases and Rainwater 156 Problems 159

CHAPTER 6

Roots: Open Methods 164

6.1 Simple Fixed-Point Iteration 165 6.2 Newton-Raphson 169 6.3 Secant Methods 174 6.4 Brent’s Method 176 6.5 MATLAB Function: fzero 181 6.6 Polynomials 183 6.7 Case Study: Pipe Friction 186 Problems 191

cha97962_fm_i-xvi.indd 6 07/11/16 8:25 pm

CONTENTS vii

CHAPTER 7

Optimization 198

7.1 Introduction and Background 199 7.2 One-Dimensional Optimization 202 7.3 Multidimensional Optimization 211 7.4 Case Study: Equilibrium and Minimum Potential Energy 213 Problems 215

Part Three Linear Systems 223

3.1 Overview 223 3.2 Part Organization 225

CHAPTER 8

Linear Algebraic Equations and Matrices 227

8.1 Matrix Algebra Overview 229 8.2 Solving Linear Algebraic Equations with MATLAB 238 8.3 Case Study: Currents and Voltages in Circuits 240 Problems 244

CHAPTER 9

Gauss Elimination 248

9.1 Solving Small Numbers of Equations 249 9.2 Naive Gauss Elimination 254 9.3 Pivoting 261 9.4 Tridiagonal Systems 264 9.5 Case Study: Model of a Heated Rod 266 Problems 270

CHAPTER 10

LU Factorization 274

10.1 Overview of LU Factorization 275 10.2 Gauss Elimination as LU Factorization 276 10.3 Cholesky Factorization 283 10.4 MATLAB Left Division 286 Problems 287

cha97962_fm_i-xvi.indd 7 07/11/16 8:25 pm

viii CONTENTS

CHAPTER 11

Matrix Inverse and Condition 288

11.1 The Matrix Inverse 288 11.2 Error Analysis and System Condition 292 11.3 Case Study: Indoor Air Pollution 297 Problems 300

CHAPTER 12

Iterative Methods 305

12.1 Linear Systems: Gauss-Seidel 305 12.2 Nonlinear Systems 312 12.3 Case Study: Chemical Reactions 320 Problems 323

CHAPTER 13

Eigenvalues 326

13.1 Mathematical Background 328 13.2 Physical Background 331 13.3 The Power Method 333 13.4 MATLAB Function: eig 336 13.5 Case Study: Eigenvalues and Earthquakes 337 Problems 340

Part Four Curve Fitting 343

4.1 Overview 343 4.2 Part Organization 345

CHAPTER 14

Linear Regression 346

14.1 Statistics Review 348 14.2 Random Numbers and Simulation 353 14.3 Linear Least-Squares Regression 358 14.4 Linearization of Nonlinear Relationships 366 14.5 Computer Applications 370 14.6 Case Study: Enzyme Kinetics 373 Problems 378

cha97962_fm_i-xvi.indd 8 07/11/16 8:25 pm

CONTENTS ix

CHAPTER 15

General Linear Least-Squares and Nonlinear Regression 385

15.1 Polynomial Regression 385 15.2 Multiple Linear Regression 389 15.3 General Linear Least Squares 391 15.4 QR Factorization and the Backslash Operator 394 15.5 Nonlinear Regression 395 15.6 Case Study: Fitting Experimental Data 397 Problems 399

CHAPTER 16

Fourier Analysis 404

16.1 Curve Fitting with Sinusoidal Functions 405 16.2 Continuous Fourier Series 411 16.3 Frequency and Time Domains 414 16.4 Fourier Integral and Transform 415 16.5 Discrete Fourier Transform (DFT) 418 16.6 The Power Spectrum 423 16.7 Case Study: Sunspots 425 Problems 426

CHAPTER 17

Polynomial Interpolation 429

17.1 Introduction to Interpolation 430 17.2 Newton Interpolating Polynomial 433 17.3 Lagrange Interpolating Polynomial 441 17.4 Inverse Interpolation 444 17.5 Extrapolation and Oscillations 445 Problems 449

CHAPTER 18

Splines and Piecewise Interpolation 453

18.1 Introduction to Splines 453 18.2 Linear Splines 455 18.3 Quadratic Splines 459 18.4 Cubic Splines 462 18.5 Piecewise Interpolation in MATLAB 468 18.6 Multidimensional Interpolation 473 18.7 Case Study: Heat Transfer 476 Problems 480

cha97962_fm_i-xvi.indd 9 07/11/16 8:25 pm

x CONTENTS

Part Five Integration and Differentiation 485

5.1 Overview 485 5.2 Part Organization 486

CHAPTER 19

Numerical Integration Formulas 488

19.1 Introduction and Background 489 19.2 Newton-Cotes Formulas 492 19.3 The Trapezoidal Rule 494 19.4 Simpson’s Rules 501 19.5 Higher-Order Newton-Cotes Formulas 507 19.6 Integration with Unequal Segments 508 19.7 Open Methods 512 19.8 Multiple Integrals 512 19.9 Case Study: Computing Work with Numerical Integration 515 Problems 518

CHAPTER 20

Numerical Integration of Functions 524

20.1 Introduction 524 20.2 Romberg Integration 525 20.3 Gauss Quadrature 530 20.4 Adaptive Quadrature 537 20.5 Case Study: Root-Mean-Square Current 540 Problems 544

CHAPTER 21

Numerical Differentiation 548

21.1 Introduction and Background 549 21.2 High-Accuracy Differentiation Formulas 552 21.3 Richardson Extrapolation 555 21.4 Derivatives of Unequally Spaced Data 557 21.5 Derivatives and Integrals for Data with Errors 558 21.6 Partial Derivatives 559 21.7 Numerical Differentiation with MATLAB 560 21.8 Case Study: Visualizing Fields 565 Problems 567

cha97962_fm_i-xvi.indd 10 07/11/16 8:25 pm

CONTENTS xi

Part six Ordinary Differential Equations 573

6.1 Overview 573 6.2 Part Organization 577

CHAPTER 22

Initial-Value Problems 579

22.1 Overview 581 22.2 Euler’s Method 581 22.3 Improvements of Euler’s Method 587 22.4 Runge-Kutta Methods 593 22.5 Systems of Equations 598 22.6 Case Study: Predator-Prey Models and Chaos 604 Problems 609

CHAPTER 23

Adaptive Methods and Stiff Systems 615

23.1 Adaptive Runge-Kutta Methods 615 23.2 Multistep Methods 624 23.3 Stiffness 628 23.4 MATLAB Application: Bungee Jumper with Cord 634 23.5 Case Study: Pliny’s Intermittent Fountain 635 Problems 640

CHAPTER 24

Boundary-Value Problems 646

24.1 Introduction and Background 647 24.2 The Shooting Method 651 24.3 Finite-Difference Methods 658 24.4 MATLAB Function: bvp4c 665 Problems 668

APPENDIX A: MATLAB BUILT-IN FUNCTIONS 674

APPENDIX B: MATLAB M-FILE FUNCTIONS 676

APPENDIX C: INTRODUCTION TO SIMULINK 677

BIBLIOGRAPHY 685

INDEX 687

cha97962_fm_i-xvi.indd 11 07/11/16 8:25 pm

Required=Results

Mobile

Connect Insight® Connect Insight is Connect’s new one-of-a-kind visual analytics dashboard—now available for both instructors and students—that provides at-a- glance information regarding student performance, which is immediately actionable. By presenting assignment, assessment, and topical performance results together with a time metric that is easily visible for aggregate or individual results, Connect Insight gives the user the ability to take a just-in-time approach to teaching and learning, which was never before available. Connect Insight presents data that empowers students and helps instructors improve class performance in a way that is efficient and effective.

73% of instructors who use Connect require it; instructor satisfaction increases by 28%

when Connect is required.

Students can view their results for any

Connect course.

Analytics

Connect’s new, intuitive mobile interface gives students and instructors flexible and convenient, anytime–anywhere access to all components of the Connect platform.

©Getty Images/iStockphoto

Using Connect improves retention rates by 19.8%, passing rates by

12.7%, and exam scores by 9.1%.

McGraw-Hill Connect® Learn Without Limits Connect is a teaching and learning platform that is proven to deliver better results for students and instructors.

Connect empowers students by continually adapting to deliver precisely what they need, when they need it, and how they need it, so your class time is more engaging and effective.

cha97962_fm_i-xvi.indd 12 07/11/16 8:25 pm

SmartBook® Proven to help students improve grades and study more efficiently, SmartBook contains the same content within the print book, but actively tailors that content to the needs of the individual. SmartBook’s adaptive technology provides precise, personalized instruction on what the student should do next, guiding the student to master and remember key concepts, targeting gaps in knowledge and offering customized feedback, and driving the student toward comprehension and retention of the subject matter. Available on tablets, SmartBook puts learning at the student’s fingertips—anywhere, anytime.

Adaptive

Over 8 billion questions have been answered, making McGraw-Hill

Education products more intelligent, reliable, and precise.

THE ADAPTIVE READING EXPERIENCE DESIGNED TO TRANSFORM THE WAY STUDENTS READ

More students earn A’s and B’s when they use McGraw-Hill Education Adaptive products.

www.mheducation.com

cha97962_fm_i-xvi.indd 13 07/11/16 8:25 pm

xiv

PREFACE

This book is designed to support a one-semester course in numerical methods. It has been written for students who want to learn and apply numerical methods in order to solve prob- lems in engineering and science. As such, the methods are motivated by problems rather than by mathematics. That said, sufficient theory is provided so that students come away with insight into the techniques and their shortcomings.

MATLAB® provides a great environment for such a course. Although other en- vironments (e.g., Excel/VBA, Mathcad) or languages (e.g., Fortran 90, C++) could have been chosen, MATLAB presently offers a nice combination of handy program- ming features with powerful built-in numerical capabilities. On the one hand, its M-file programming environment allows students to implement moderately compli- cated algorithms in a structured and coherent fashion. On the other hand, its built-in, numerical capabilities empower students to solve more difficult problems without try- ing to “reinvent the wheel.”

The basic content, organization, and pedagogy of the third edition are essentially pre- served in the fourth edition. In particular, the conversational writing style is intentionally maintained in order to make the book easier to read. This book tries to speak directly to the reader and is designed in part to be a tool for self-teaching.

That said, this edition differs from the past edition in three major ways: (1) new material, (2) new and revised homework problems, and (3) an appendix introducing Simulink.

1. New Content. I have included new and enhanced sections on a number of topics. The primary additions include material on some MATLAB functions not included in previ- ous editions (e.g., fsolve, integrate, bvp4c), some new applications of Monte Carlo for problems such as integration and optimization, and MATLAB’s new way to pass parameters to function functions.

2. New Homework Problems. Most of the end-of-chapter problems have been modified, and a variety of new problems have been added. In particular, an effort has been made to include several new problems for each chapter that are more challenging and dif- ficult than the problems in the previous edition.

3. I have developed a short primer on Simulink which I have my students read prior to covering that topic. Although I recognize that some professors may not choose to cover Simulink, I included it as a teaching aid for those that do.

cha97962_fm_i-xvi.indd 14 07/11/16 8:25 pm

Aside from the new material and problems, the fourth edition is very similar to the third. In particular, I have endeavored to maintain most of the features contributing to its pedagogical effectiveness including extensive use of worked examples and engineering and scien tific applications. As with the previous edition, I have made a concerted effort to make this book as “student-friendly” as possible. Thus, I’ve tried to keep my explanations straight- forward and practical.

Although my primary intent is to empower students by providing them with a sound introduction to numerical problem solving, I have the ancillary objective of making this introduction exciting and pleasurable. I believe that motivated students who enjoy engi- neering and science, problem solving, mathematics—and yes—programming, will ulti- mately make better professionals. If my book fosters enthusiasm and appreciation for these subjects, I will consider the effort a success.

Acknowledgments. Several members of the McGraw-Hill team have contributed to this project. Special thanks are due to Jolynn Kilburg, Thomas Scaife, Ph.D., Chelsea Haupt, Ph.D., and Jeni McAtee for their encouragement, support, and direction.

During the course of this project, the folks at The MathWorks, Inc., have truly dem- onstrated their overall excellence as well as their strong commitment to engineering and science education. In particular, Naomi Fernandes of The MathWorks, Inc., Book Program has been especially helpful and Jared Wasserman of the MathWorks Technical Support Department was of great help with technical questions.

The generosity of the Berger family has provided me with the opportunity to work on creative projects such as this book dealing with computing and engineering. In addition, my colleagues in the School of Engineering at Tufts, notably Masoud Sanayei, Babak Moaveni, Luis Dorfmann, Rob White, Linda Abriola, and Laurie Baise, have been very supportive and helpful.

Significant suggestions were also given by a number of colleagues. In particular, Dave Clough (University of Colorado–Boulder), and Mike Gustafson (Duke University) pro- vided valuable ideas and suggestions. In addition, a number of reviewers provided use- ful feedback and advice including Karen Dow Ambtman (University of Alberta), Jalal Behzadi (Shahid Chamran University), Eric Cochran (Iowa State University), Frederic Gibou (University of California at Santa Barbara), Jane Grande-Allen (Rice University), Raphael Haftka (University of Florida), Scott Hendricks (Virginia Tech University), Ming Huang (University of San Diego), Oleg Igoshin (Rice University), David Jack (Baylor Uni- versity), Se Won Lee (Sungkyunkwan University), Clare McCabe (Vanderbilt University), Eckart Meiburg (University of California at Santa Barbara), Luis Ricardez (University of Waterloo), James Rottman (University of California, San Diego), Bingjing Su (University of Cincinnati), Chin-An Tan (Wayne State University), Joseph Tipton (The University of Evansville), Marion W. Vance (Arizona State University), Jonathan Vande Geest (University of Arizona), Leah J. Walker (Arkansas State University), Qiang Hu (University of Alabama, Huntsville), Yukinobu Tanimoto (Tufts University), Henning T. Søgaard (Aarhus University), and Jimmy Feng (University of British Columbia).

It should be stressed that although I received useful advice from the aforementioned individuals, I am responsible for any inaccuracies or mistakes you may find in this book. Please contact me via e-mail if you should detect any errors.

PREFACE xv

cha97962_fm_i-xvi.indd 15 07/11/16 8:25 pm

xvi PREFACE

Finally, I want to thank my family, and in particular my wife, Cynthia, for the love, patience, and support they have provided through the time I’ve spent on this project.

Steven C. Chapra Tufts University

Medford, Massachusetts steven.chapra@tufts.edu

PEDAGOGICAL TOOLS

Theory Presented as It Informs Key Concepts. The text is intended for Numerical Methods users, not developers. Therefore, theory is not included for “theory’s sake,” for ex- ample no proofs. Theory is included as it informs key concepts such as the Taylor series, con- vergence, condition, etc. Hence, the student is shown how the theory connects with practical issues in problem solving.

Introductory MATLAB Material. The text in cludes two introductory chapters on how to use MATLAB. Chapter 2 shows students how to per form computations and create graphs in MATLAB’s standard command mode. Chapter 3 provides a primer on developing numerical programs via MATLAB M-file functions. Thus, the text provides students with the means to develop their own nu merical algorithms as well as to tap into MATLAB’s powerful built-in routines.

Algorithms Presented Using MATLAB M-files. Instead of using pseudocode, this book presents algorithms as well-structured MATLAB M-files. Aside from being useful com- puter programs, these provide students with models for their own M-files that they will develop as homework exercises.

Worked Examples and Case Studies. Extensive worked examples are laid out in detail so that students can clearly follow the steps in each numerical computation. The case stud- ies consist of engineering and science applications which are more complex and richer than the worked examples. They are placed at the ends of selected chapters with the intention of (1) illustrating the nuances of the methods and (2) showing more realistically how the methods along with MATLAB are applied for problem solving.

Problem Sets. The text includes a wide variety of problems. Many are drawn from en- gineering and scientific disciplines. Others are used to illustrate numerical techniques and theoretical concepts. Problems include those that can be solved with a pocket calculator as well as others that require computer solution with MATLAB.

Useful Appendices and Indexes. Appendix A contains MATLAB commands, Appendix B contains M-file functions, and new Appendix C contains a brief Simulink primer.

Instructor Resources. Solutions Manual, Lecture PowerPoints, Text images in Power- Point, M-files and additional MATLAB resources are available through Connect®.

cha97962_fm_i-xvi.indd 16 07/11/16 8:25 pm

1

Part One

Modeling, Computers, and Error Analysis

1.1 MOTIVATION

What are numerical methods and why should you study them? Numerical methods are techniques by which mathematical problems are formulated

so that they can be solved with arithmetic and logical operations. Because digital comput- ers excel at performing such operations, numerical methods are sometimes referred to as computer mathematics.

In the pre–computer era, the time and drudgery of implementing such calculations seriously limited their practical use. However, with the advent of fast, inexpensive digital computers, the role of numerical methods in engineering and scientific problem solving has exploded. Because they figure so prominently in much of our work, I believe that numerical methods should be a part of every engineer’s and scientist’s basic education. Just as we all must have solid foundations in the other areas of mathematics and science, we should also have a fundamental understanding of numerical methods. In particular, we

should have a solid appreciation of both their capabilities and their limitations.

Beyond contributing to your overall education, there are several additional reasons why you should study numerical methods:

1. Numerical methods greatly expand the types of problems you can address. They are capable of handling large sys- tems of equations, nonlinearities, and complicated geometries that are not uncommon in engineering and science and that are often impossible to solve analytically with standard calculus. As such, they greatly enhance your prob- lem-solving skills.

2. Numerical methods allow you to use “canned” software with insight. During

cha97962_ch01_001-026.indd 1 08/11/16 12:51 pm

2 PART 1 ModEling, CoMPuTERs, And ERRoR AnAlysis

your career, you will invariably have occasion to use commercially available prepack- aged computer programs that involve numerical methods. The intelligent use of these programs is greatly enhanced by an understanding of the basic theory underlying the methods. In the absence of such understanding, you will be left to treat such packages as “black boxes” with little critical insight into their inner workings or the validity of the results they produce.

3. Many problems cannot be approached using canned programs. If you are conversant with numerical methods, and are adept at computer programming, you can design your own programs to solve problems without having to buy or commission expensive software.

4. Numerical methods are an efficient vehicle for learning to use computers. Because nu- merical methods are expressly designed for computer implementation, they are ideal for illustrating the computer’s powers and limitations. When you successfully implement numerical methods on a computer, and then apply them to solve otherwise intractable problems, you will be provided with a dramatic demonstration of how computers can serve your professional development. At the same time, you will also learn to acknowl- edge and control the errors of approximation that are part and parcel of large-scale numerical calculations.

5. Numerical methods provide a vehicle for you to reinforce your understanding of math- ematics. Because one function of numerical methods is to reduce higher mathematics to basic arithmetic operations, they get at the “nuts and bolts” of some otherwise obscure topics. Enhanced understanding and insight can result from this alternative perspective.

With these reasons as motivation, we can now set out to understand how numerical methods and digital computers work in tandem to generate reliable solutions to mathemati- cal problems. The remainder of this book is devoted to this task.

1.2 PART ORGANIZATION

This book is divided into six parts. The latter five parts focus on the major areas of nu- merical methods. Although it might be tempting to jump right into this material, Part One consists of four chapters dealing with essential background material.

Chapter 1 provides a concrete example of how a numerical method can be employed to solve a real problem. To do this, we develop a mathematical model of a free-falling bungee jumper. The model, which is based on Newton’s second law, results in an ordinary differential equation. After first using calculus to develop a closed-form solution, we then show how a comparable solution can be generated with a simple numerical method. We end the chapter with an overview of the major areas of numerical methods that we cover in Parts Two through Six.

Chapters 2 and 3 provide an introduction to the MATLAB® software environment. Chapter 2 deals with the standard way of operating MATLAB by entering commands one at a time in the so-called calculator, or command, mode. This interactive mode provides a straightforward means to orient you to the environment and illustrates how it is used for common operations such as performing calculations and creating plots.

cha97962_ch01_001-026.indd 2 08/11/16 12:51 pm

1.2 PART oRgAniZATion 3

Chapter 3 shows how MATLAB’s programming mode provides a vehicle for assem- bling individual commands into algorithms. Thus, our intent is to illustrate how MATLAB serves as a convenient programming environment to develop your own software.

Chapter 4 deals with the important topic of error analysis, which must be understood for the effective use of numerical methods. The first part of the chapter focuses on the roundoff errors that result because digital computers cannot represent some quantities exactly. The latter part addresses truncation errors that arise from using an approximation in place of an exact mathematical procedure.

cha97962_ch01_001-026.indd 3 08/11/16 12:51 pm

4

Mathematical Modeling, numerical Methods, and Problem solving

1

CHAPTER OBJECTIVES

The primary objective of this chapter is to provide you with a concrete idea of what numerical methods are and how they relate to engineering and scientific problem solving. Specific objectives and topics covered are

• Learning how mathematical models can be formulated on the basis of scientific principles to simulate the behavior of a simple physical system.

• Understanding how numerical methods afford a means to generate solutions in a manner that can be implemented on a digital computer.

• Understanding the different types of conservation laws that lie beneath the models used in the various engineering disciplines and appreciating the difference between steady-state and dynamic solutions of these models.

• Learning about the different types of numerical methods we will cover in this book.

YOU’VE GOT A PROBLEM

suppose that a bungee-jumping company hires you. You’re given the task of predicting the velocity of a jumper (Fig. 1.1) as a function of time during the free-fall part of the jump. This information will be used as part of a larger analysis to determine the length and required strength of the bungee cord for jumpers of different mass.

You know from your studies of physics that the acceleration should be equal to the ratio of the force to the mass (Newton’s second law). Based on this insight and your knowledge

cha97962_ch01_001-026.indd 4 08/11/16 12:51 pm

1.1 A siMPlE MAThEMATiCAl ModEl 5

of physics and fluid mechanics, you develop the following mathematical model for the rate of change of velocity with respect to time,

dυ ___ dt = g − cd ___ m υ

2

where υ = downward vertical velocity (m/s), t = time (s), g = the acceleration due to gravity (≅ 9.81 m/s2), cd = a lumped drag coefficient (kg/m), and m = the jumper’s mass (kg). The drag coefficient is called “lumped” because its magnitude depends on fac- tors such as the jumper’s area and the fluid density (see Sec. 1.4).

Because this is a differential equation, you know that calculus might be used to obtain an analytical or exact solution for υ as a function of t. However, in the following pages, we will illustrate an alternative solution approach. This will involve developing a computer- oriented numerical or approximate solution.

Aside from showing you how the computer can be used to solve this particular prob- lem, our more general objective will be to illustrate (a) what numerical methods are and (b) how they figure in engineering and scientific problem solving. In so doing, we will also show how mathematical models figure prominently in the way engineers and scientists use numerical methods in their work.

1.1 A SIMPLE MATHEMATICAL MODEL

A mathematical model can be broadly defined as a formulation or equation that expresses the essential features of a physical system or process in mathematical terms. In a very gen- eral sense, it can be represented as a functional relationship of the form

Dependent variable = f ( independent variables , parameters, forcing functions ) (1.1)

where the dependent variable is a characteristic that typically reflects the behavior or state of the system; the independent variables are usually dimensions, such as time and space, along which the system’s behavior is being determined; the parameters are reflective of the system’s properties or composition; and the forcing functions are external influences acting upon it.

The actual mathematical expression of Eq. (1.1) can range from a simple algebraic relationship to large complicated sets of differential equations. For example, on the basis of his observations, Newton formulated his second law of motion, which states that the time rate of change of momentum of a body is equal to the resultant force acting on it. The mathematical expression, or model, of the second law is the well-known equation

F = ma (1.2)

where F is the net force acting on the body (N, or kg m/s2), m is the mass of the object (kg), and a is its acceleration (m/s2).

Upward force due to air resistance

Downward force due to gravity

FIGURE 1.1 Forces acting on a free-falling bungee jumper.

cha97962_ch01_001-026.indd 5 08/11/16 12:51 pm

6 MAThEMATiCAl ModEling, nuMERiCAl METhods, And PRoblEM solving

The second law can be recast in the format of Eq. (1.1) by merely dividing both sides by m to give

a = F __ m (1.3)

where a is the dependent variable reflecting the system’s behavior, F is the forcing func- tion, and m is a parameter. Note that for this simple case there is no independent variable because we are not yet predicting how acceleration varies in time or space.

Equation (1.3) has a number of characteristics that are typical of mathematical models of the physical world.

• It describes a natural process or system in mathematical terms. • It represents an idealization and simplification of reality. That is, the model ignores

negligible details of the natural process and focuses on its essential manifestations. Thus, the second law does not include the effects of relativity that are of minimal importance when applied to objects and forces that interact on or about the earth’s surface at velocities and on scales visible to humans.

• Finally, it yields reproducible results and, consequently, can be used for predictive purposes. For example, if the force on an object and its mass are known, Eq. (1.3) can be used to compute acceleration.

Because of its simple algebraic form, the solution of Eq. (1.2) was obtained easily. However, other mathematical models of physical phenomena may be much more complex, and either cannot be solved exactly or require more sophisticated mathematical techniques than simple algebra for their solution. To illustrate a more complex model of this kind, Newton’s second law can be used to determine the terminal velocity of a free-falling body near the earth’s surface. Our falling body will be a bungee jumper (Fig. 1.1). For this case, a model can be derived by expressing the acceleration as the time rate of change of the velocity (dυ/dt) and substituting it into Eq. (1.3) to yield

dυ ___ dt = F __ m (1.4)

where υ is velocity (in meters per second). Thus, the rate of change of the velocity is equal to the net force acting on the body normalized to its mass. If the net force is positive, the object will accelerate. If it is negative, the object will decelerate. If the net force is zero, the object’s velocity will remain at a constant level.

Next, we will express the net force in terms of measurable variables and parameters. For a body falling within the vicinity of the earth, the net force is composed of two opposing forces: the downward pull of gravity FD and the upward force of air resistance FU (Fig. 1.1):

F = FD + FU (1.5)

If force in the downward direction is assigned a positive sign, the second law can be used to formulate the force due to gravity as

FD = mg (1.6)

where g is the acceleration due to gravity (9.81 m/s2).

cha97962_ch01_001-026.indd 6 08/11/16 12:51 pm

1.1 A siMPlE MAThEMATiCAl ModEl 7

Air resistance can be formulated in a variety of ways. Knowledge from the science of fluid mechanics suggests that a good first approximation would be to assume that it is proportional to the square of the velocity,

FU = −cd υ2 (1.7)

where cd is a proportionality constant called the lumped drag coefficient (kg/m). Thus, the greater the fall velocity, the greater the upward force due to air resistance. The parameter cd accounts for properties of the falling object, such as shape or surface roughness, that af- fect air resistance. For the present case, cd might be a function of the type of clothing or the orientation used by the jumper during free fall.

The net force is the difference between the downward and upward force. Therefore, Eqs. (1.4) through (1.7) can be combined to yield

dυ ___ dt = g − cd ___ m υ

2 (1.8)

Equation (1.8) is a model that relates the acceleration of a falling object to the forces acting on it. It is a differential equation because it is written in terms of the differential rate of change (dυ/dt) of the variable that we are interested in predicting. However, in contrast to the solution of Newton’s second law in Eq. (1.3), the exact solution of Eq. (1.8) for the velocity of the jumper cannot be obtained using simple algebraic manipulation. Rather, more advanced techniques such as those of calculus must be applied to obtain an exact or analytical solution. For example, if the jumper is initially at rest (υ = 0 at t = 0), calculus can be used to solve Eq. (1.8) for

υ(t) = √ ___

gm ___ cd tanh ( √ ___

gcd ___ m t ) (1.9)

where tanh is the hyperbolic tangent that can be either computed directly1 or via the more elementary exponential function as in

tanh x = e x − e−x _______ ex + e−x (1.10)

Note that Eq. (1.9) is cast in the general form of Eq. (1.1) where υ(t) is the dependent variable, t is the independent variable, cd and m are parameters, and g is the forcing function.

EXAMPlE 1.1 Analytical solution to the bungee Jumper Problem

Problem statement. A bungee jumper with a mass of 68.1 kg leaps from a stationary hot air balloon. Use Eq. (1.9) to compute velocity for the first 12 s of free fall. Also deter- mine the terminal velocity that will be attained for an infinitely long cord (or alternatively, the jumpmaster is having a particularly bad day!). Use a drag coefficient of 0.25 kg/m.

1 MATLAB allows direct calculation of the hyperbolic tangent via the built-in function tanh(x).

cha97962_ch01_001-026.indd 7 08/11/16 12:51 pm

8 MAThEMATiCAl ModEling, nuMERiCAl METhods, And PRoblEM solving

solution. Inserting the parameters into Eq. (1.9) yields

υ(t) = √ _________

9.81(68.1) _________ 0.25 tanh ( √ _________

9.81(0.25) _________ 68.1 t ) = 51.6938 tanh(0.18977t) which can be used to compute

t, s υ, m/s

0 0 2 18.7292 4 33.1118 6 42.0762 8 46.9575 10 49.4214 12 50.6175 ∞ 51.6938

According to the model, the jumper accelerates rapidly (Fig. 1.2). A velocity of 49.4214 m/s (about 110 mi/hr) is attained after 10 s. Note also that after a sufficiently

0

20

40

60

0 4 8 12 t, s

Terminal velocity

υ, m

/s

FIGURE 1.2 The analytical solution for the bungee jumper problem as computed in Example 1.1. Velocity increases with time and asymptotically approaches a terminal velocity.

cha97962_ch01_001-026.indd 8 08/11/16 12:51 pm

1.1 A siMPlE MAThEMATiCAl ModEl 9

long time, a constant velocity, called the terminal velocity, of 51.6983 m/s (115.6 mi/hr) is reached. This velocity is constant because, eventually, the force of gravity will be in balance with the air resistance. Thus, the net force is zero and acceleration has ceased.

Equation (1.9) is called an analytical or closed-form solution because it exactly satis- fies the original differential equation. Unfortunately, there are many mathematical models that cannot be solved exactly. In many of these cases, the only alternative is to develop a numerical solution that approximates the exact solution.

Numerical methods are those in which the mathematical problem is reformulated so it can be solved by arithmetic operations. This can be illustrated for Eq. (1.8) by realizing that the time rate of change of velocity can be approximated by (Fig. 1.3):

dυ ___ dt ≅ Δυ ___ Δt =

υ(ti+1) − υ(ti ) ___________ ti+1 − ti (1.11)

where Δυ and Δt are differences in velocity and time computed over finite intervals, υ(ti) is velocity at an initial time ti, and υ(ti+1) is velocity at some later time ti+1. Note that dυ/dt ≅ Δυ/Δt is approximate because Δt is finite. Remember from calculus that

dυ ___ dt = lim Δt→0

Δυ ___ Δt

Equation (1.11) represents the reverse process.

υ(ti+1)

ti+1 t

∆t

υ(ti)

∆υ

ti

True slope dυ/dt

Approximate slope ∆υ υ(ti+1) − υ(ti)

ti+1 − ti∆t =

FIGURE 1.3 The use of a finite difference to approximate the first derivative of υ with respect to t.

cha97962_ch01_001-026.indd 9 08/11/16 12:51 pm

10 MAThEMATiCAl ModEling, nuMERiCAl METhods, And PRoblEM solving

Equation (1.11) is called a finite-difference approximation of the derivative at time ti. It can be substituted into Eq. (1.8) to give

υ(ti+1) − υ(ti ) ___________ ti+1 − ti

= g − cd ___ m υ(ti )

2

This equation can then be rearranged to yield

υ(ti+1) = υ(ti ) + [ g − cd ___ m υ(ti )2 ] (ti+1 − ti ) (1.12) Notice that the term in brackets is the right-hand side of the differential equation itself

[Eq. (1.8)]. That is, it provides a means to compute the rate of change or slope of υ. Thus, the equation can be rewritten more concisely as

υi+1 = υi + dυi ___ dt Δt (1.13)

where the nomenclature υi designates velocity at time ti, and Δt = ti+1 − ti. We can now see that the differential equation has been transformed into an equation that

can be used to determine the velocity algebraically at ti+1 using the slope and previous values of υ and t. If you are given an initial value for velocity at some time ti, you can easily compute velocity at a later time ti+1. This new value of velocity at ti+1 can in turn be employed to extend the computation to velocity at ti+2 and so on. Thus at any time along the way,

New value = old value + slope × step size

This approach is formally called Euler’s method. We’ll discuss it in more detail when we turn to differential equations later in this book.

EXAMPlE 1.2 numerical solution to the bungee Jumper Problem

Problem statement. Perform the same computation as in Example 1.1 but use Eq. (1.12) to compute velocity with Euler’s method. Employ a step size of 2 s for the calculation.

solution. At the start of the computation (t0 = 0), the velocity of the jumper is zero. Using this information and the parameter values from Example 1.1, Eq. (1.12) can be used to compute velocity at t1 = 2 s:

υ = 0 + [ 9.81 − 0.25 ____ 68.1 (0)2 ] × 2 = 19.62 m/s For the next interval (from t = 2 to 4 s), the computation is repeated, with the result

υ = 19.62 + [ 9.81 − 0.25 ____ 68.1 (19.62)2 ] × 2 = 36.4137 m/s

cha97962_ch01_001-026.indd 10 08/11/16 12:51 pm

1.1 A siMPlE MAThEMATiCAl ModEl 11

The calculation is continued in a similar fashion to obtain additional values:

t, s υ, m/s

 0 0 2 19.6200 4 36.4137 6 46.2983 8 50.1802 10 51.3123 12 51.6008 ∞ 51.6938

The results are plotted in Fig. 1.4 along with the exact solution. We can see that the numeri- cal method captures the essential features of the exact solution. However, because we have employed straight-line segments to approximate a continuously curving function, there is some discrepancy between the two results. One way to minimize such discrepancies is to use a smaller step size. For example, applying Eq. (1.12) at 1-s intervals results in a smaller error, as the straight-line segments track closer to the true solution. Using hand calcula- tions, the effort associated with using smaller and smaller step sizes would make such numerical solutions impractical. However, with the aid of the computer, large numbers of calculations can be performed easily. Thus, you can accurately model the velocity of the jumper without having to solve the differential equation exactly.

0

20

40

60

0 4 8 12 t, s

Terminal velocity

υ, m

/s

Approximate, numerical solution

Exact, analytical solution

FIGURE 1.4 Comparison of the numerical and analytical solutions for the bungee jumper problem.

cha97962_ch01_001-026.indd 11 08/11/16 12:51 pm

12 MAThEMATiCAl ModEling, nuMERiCAl METhods, And PRoblEM solving

As in Example 1.2, a computational price must be paid for a more accurate numeri- cal result. Each halving of the step size to attain more accuracy leads to a doubling of the number of computations. Thus, we see that there is a trade-off between accuracy and com- putational effort. Such trade-offs figure prominently in numerical methods and constitute an important theme of this book.

1.2 CONSERVATION LAWS IN ENGINEERING AND SCIENCE

Aside from Newton’s second law, there are other major organizing principles in science and engineering. Among the most important of these are the conservation laws. Although they form the basis for a variety of complicated and powerful mathematical models, the great conservation laws of science and engineering are conceptually easy to understand. They all boil down to

Change = increases − decreases (1.14)

This is precisely the format that we employed when using Newton’s law to develop a force balance for the bungee jumper [Eq. (1.8)].

Although simple, Eq. (1.14) embodies one of the most fundamental ways in which conservation laws are used in engineering and science—that is, to predict changes with respect to time. We will give it a special name—the time-variable (or transient) computation.

Aside from predicting changes, another way in which conservation laws are applied is for cases where change is nonexistent. If change is zero, Eq. (1.14) becomes

Change = 0 = increases − decreases or

Increases = decreases (1.15) Thus, if no change occurs, the increases and decreases must be in balance. This case, which is also given a special name—the steady-state calculation—has many applications in engi- neering and science. For example, for steady-state incompressible fluid flow in pipes, the flow into a junction must be balanced by flow going out, as in

Flow in = flow out For the junction in Fig. 1.5, the balance that can be used to compute that the flow out of the fourth pipe must be 60.

For the bungee jumper, the steady-state condition would correspond to the case where the net force was zero or [Eq. (1.8) with dυ/dt = 0]

mg = cd υ2 (1.16) Thus, at steady state, the downward and upward forces are in balance and Eq. (1.16) can be solved for the terminal velocity

υ = √ ___

gm ___ cd

Although Eqs. (1.14) and (1.15) might appear trivially simple, they embody the two funda- mental ways that conservation laws are employed in engineering and science. As such, they will form an important part of our efforts in subsequent chapters to illustrate the connection between numerical methods and engineering and science.

cha97962_ch01_001-026.indd 12 08/11/16 12:51 pm

1.3 nuMERiCAl METhods CovEREd in This book 13

Table 1.1 summarizes some models and associated conservation laws that figure prominently in engineering. Many chemical engineering problems involve mass balances for reactors. The mass balance is derived from the conservation of mass. It specifies that the change of mass of a chemical in the reactor depends on the amount of mass flowing in minus the mass flowing out.

Civil and mechanical engineers often focus on models developed from the conserva- tion of momentum. For civil engineering, force balances are utilized to analyze structures such as the simple truss in Table 1.1. The same principles are employed for the mechanical engineering case studies to analyze the transient up-and-down motion or vibrations of an automobile.

Finally, electrical engineering studies employ both current and energy balances to model electric circuits. The current balance, which results from the conservation of charge, is simi- lar in spirit to the flow balance depicted in Fig. 1.5. Just as flow must balance at the junction of pipes, electric current must balance at the junction of electric wires. The energy balance specifies that the changes of voltage around any loop of the circuit must add up to zero.

It should be noted that there are many other branches of engineering beyond chemical, civil, electrical, and mechanical. Many of these are related to the Big Four. For example, chem- ical engineering skills are used extensively in areas such as environmental, petroleum, and bio- medical engineering. Similarly, aerospace engineering has much in common with mechanical engineering. I will endeavor to include examples from these areas in the coming pages.

1.3 NUMERICAL METHODS COVERED IN THIS BOOK

Euler’s method was chosen for this introductory chapter because it is typical of many other classes of numerical methods. In essence, most consist of recasting mathematical operations into the simple kind of algebraic and logical operations compatible with digital computers. Figure 1.6 summarizes the major areas covered in this text.

Pipe 2 Flow in = 80

Pipe 3 Flow out = 120

Pipe 4 Flow out = ?

Pipe 1 Flow in = 100

FIGURE 1.5 A flow balance for steady incompressible fluid flow at the junction of pipes.

cha97962_ch01_001-026.indd 13 08/11/16 12:51 pm

14 MAThEMATiCAl ModEling, nuMERiCAl METhods, And PRoblEM solving

TABLE 1.1 devices and types of balances that are commonly used in the four major areas of engineering. For each case, the conservation law on which the balance is based is specified.

Field Organizing Principle Mathematical ExpressionDevice

Force balance:Mechanical engineering

Conservation of momentum

Upward force

Downward force

x = 0

= downward force − upward forcem d 2x

dt2

Machine

Current balance:

Voltage balance:

Around each loop Σ emf’s − Σ voltage drops for resistors = 0 Σ ξ − Σ iR = 0

For each node Σ current (i) = 0

Electrical engineering

Conservation of energy

Conservation of charge

+i2