Numerical methods for engineers 7th edition chapra pdf download

Citation preview

Numerical Methods for Engineers S E VEN TH ED IT I ON

Steven C. Chapra Berger Chair in Computing and Engineering Tufts University

Raymond P. Canale Professor Emeritus of Civil Engineering University of Michigan

NUMERICAL METHODS FOR ENGINEERS, SEVENTH EDITION Published by McGraw-Hill Education, 2 Penn Plaza, New York, NY 10121. Copyright © 2015 by McGraw-Hill Education. All rights reserved. Printed in the United States of America. Previous editions © 2010, 2006, and 2002. 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 0 DOC/DOC 1 0 9 8 7 6 5 4 ISBN 978–0–07–339792–4 MHID 0–07–339792–x Senior Vice President, Products & Markets: Kurt L. Strand Vice President, General Manager, Products & Markets: Marty Lange Vice President, Content Production & Technology Services: Kimberly Meriwether David Executive Brand Manager: Bill Stenquist Managing Director: Thomas Timp Global Publisher: Raghothaman Srinivasan Developmental Editor: Lorraine Buczek Marketing Manager: Heather Wagner Director, Content Production: Terri Schiesl Senior Content Project Manager: Melissa M. Leick Buyer: Jennifer Pickel Cover Designer: Studio Montage, St. Louis, MO Cover Image: Peak towering above clouds: Royalty-Free/CORBIS; Skysurfers: Getty Images/Digital Vision/RF Media Project Manager: Sandra M. Schnee Compositor: Aptara®, Inc. Typeface: 10/12 Time Roman Printer: R. R. Donnelley 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. Numerical methods for engineers / Steven C. Chapra, Berger chair in computing and engineering, Tufts University, Raymond P. Canale, professor emeritus of civil engineering, University of Michigan. — Seventh edition. pages cm Includes bibliographical references and index. ISBN 978-0-07-339792-4 (alk. paper) — ISBN 0-07-339792-X (alk. paper) 1. Engineering mathematics—Data processing. 2. Numerical calculations—Data processing 3. Microcomputers—Programming. I. Canale, Raymond P. II. Title. TA345.C47 2015 518.024’62—dc23

2013041704

To Margaret and Gabriel Chapra Helen and Chester Canale

CONTENTS PREFACE

xiv

ABOUT THE AUTHORS

xvi

PART ONE MODELING, COMPUTERS, AND ERROR ANALYSIS 3

PT1.1 Motivation 3 PT1.2 Mathematical Background PT1.3 Orientation 8

5

CHAPTER 1 Mathematical Modeling and Engineering Problem Solving 1.1 A Simple Mathematical Model 11 1.2 Conservation Laws and Engineering Problems 21 CHAPTER 2 Programming and Software

18

27

2.1 Packages and Programming 27 2.2 Structured Programming 28 2.3 Modular Programming 37 2.4 Excel 39 2.5 MATLAB 43 2.6 Mathcad 47 2.7 Other Languages and Libraries 48 Problems 49 CHAPTER 3 Approximations and Round-Off Errors 3.1 Significant Figures 56 3.2 Accuracy and Precision 3.3 Error Definitions 59 3.4 Round-Off Errors 65 Problems 79

iv

58

55

11

CONTENTS

v

CHAPTER 4 Truncation Errors and the Taylor Series

81

4.1 The Taylor Series 81 4.2 Error Propagation 97 4.3 Total Numerical Error 101 4.4 Blunders, Formulation Errors, and Data Uncertainty Problems 108 EPILOGUE: PART ONE 110 PT1.4 Trade-Offs 110 PT1.5 Important Relationships and Formulas 113 PT1.6 Advanced Methods and Additional References

106

113

PART TWO ROOTS OF EQUATIONS 117

PT2.1 Motivation 117 PT2.2 Mathematical Background PT2.3 Orientation 120 CHAPTER 5 Bracketing Methods

119

123

5.1 Graphical Methods 123 5.2 The Bisection Method 127 5.3 The False-Position Method 135 5.4 Incremental Searches and Determining Initial Guesses Problems 142 CHAPTER 6 Open Methods

145

6.1 Simple Fixed-Point Iteration 146 6.2 The Newton-Raphson Method 151 6.3 The Secant Method 157 6.4 Brent’s Method 162 6.5 Multiple Roots 166 6.6 Systems of Nonlinear Equations 169 Problems 173 CHAPTER 7 Roots of Polynomials

176

7.1 Polynomials in Engineering and Science 7.2 Computing with Polynomials 179 7.3 Conventional Methods 182

176

141

vi

CONTENTS 7.4 Müller’s Method 183 7.5 Bairstow’s Method 187 7.6 Other Methods 192 7.7 Root Location with Software Packages Problems 202 CHAPTER 8 Case Studies: Roots of Equations

192

204

8.1 Ideal and Nonideal Gas Laws (Chemical/Bio Engineering) 204 8.2 Greenhouse Gases and Rainwater (Civil/Environmental Engineering) 8.3 Design of an Electric Circuit (Electrical Engineering) 209 8.4 Pipe Friction (Mechanical/Aerospace Engineering) 212 Problems 215 EPILOGUE: PART TWO 226 PT2.4 Trade-Offs 226 PT2.5 Important Relationships and Formulas 227 PT2.6 Advanced Methods and Additional References

PART THREE LINEAR ALGEBRAIC EQUATIONS 231

PT3.1 Motivation 231 PT3.2 Mathematical Background PT3.3 Orientation 241 CHAPTER 9 Gauss Elimination

233

245

9.1 Solving Small Numbers of Equations 245 9.2 Naive Gauss Elimination 252 9.3 Pitfalls of Elimination Methods 258 9.4 Techniques for Improving Solutions 264 9.5 Complex Systems 271 9.6 Nonlinear Systems of Equations 271 9.7 Gauss-Jordan 273 9.8 Summary 275 Problems 275 CHAPTER 10 LU Decomposition and Matrix Inversion 10.1 LU Decomposition 278 10.2 The Matrix Inverse 287 10.3 Error Analysis and System Condition Problems 297

291

278

227

207

CONTENTS

vii

CHAPTER 11 Special Matrices and Gauss-Seidel

300

11.1 Special Matrices 300 11.2 Gauss-Seidel 304 11.3 Linear Algebraic Equations with Software Packages Problems 316 CHAPTER 12 Case Studies: Linear Algebraic Equations

311

319

12.1 Steady-State Analysis of a System of Reactors (Chemical/Bio Engineering) 319 12.2 Analysis of a Statically Determinate Truss (Civil/Environmental Engineering) 322 12.3 Currents and Voltages in Resistor Circuits (Electrical Engineering) 326 12.4 Spring-Mass Systems (Mechanical/Aerospace Engineering) 328 Problems 331 EPILOGUE: PART THREE 341 PT3.4 Trade-Offs 341 PT3.5 Important Relationships and Formulas 342 PT3.6 Advanced Methods and Additional References

342

PART FOUR OPTIMIZATION

345

PT4.1 Motivation 345 PT4.2 Mathematical Background PT4.3 Orientation 351

350

CHAPTER 13 One-Dimensional Unconstrained Optimization

355

13.1 Golden-Section Search 356 13.2 Parabolic Interpolation 363 13.3 Newton’s Method 365 13.4 Brent’s Method 366 Problems 368 CHAPTER 14 Multidimensional Unconstrained Optimization 14.1 Direct Methods 371 14.2 Gradient Methods 375 Problems 388

370

viii

CONTENTS CHAPTER 15 Constrained Optimization

390

15.1 Linear Programming 390 15.2 Nonlinear Constrained Optimization 401 15.3 Optimization with Software Packages 402 Problems 413 CHAPTER 16 Case Studies: Optimization

416

16.1 Least-Cost Design of a Tank (Chemical/Bio Engineering) 416 16.2 Least-Cost Treatment of Wastewater (Civil/Environmental Engineering) 421 16.3 Maximum Power Transfer for a Circuit (Electrical Engineering) 425 16.4 Equilibrium and Minimum Potential Energy (Mechanical/Aerospace Engineering) 429 Problems 431 EPILOGUE: PART FOUR 438 PT4.4 Trade-Offs 438 PT4.5 Additional References 439

PART FIVE CURVE FITTING 441

PT5.1 Motivation 441 PT5.2 Mathematical Background PT5.3 Orientation 452 CHAPTER 17 Least-Squares Regression

443

456

17.1 Linear Regression 456 17.2 Polynomial Regression 472 17.3 Multiple Linear Regression 476 17.4 General Linear Least Squares 479 17.5 Nonlinear Regression 483 Problems 487 CHAPTER 18 Interpolation

490

18.1 Newton’s Divided-Difference Interpolating Polynomials 18.2 Lagrange Interpolating Polynomials 502 18.3 Coefficients of an Interpolating Polynomial 507 18.4 Inverse Interpolation 507 18.5 Additional Comments 508 18.6 Spline Interpolation 511 18.7 Multidimensional Interpolation 521 Problems 524

491

CONTENTS CHAPTER 19 Fourier Approximation

ix

526

19.1 Curve Fitting with Sinusoidal Functions 527 19.2 Continuous Fourier Series 533 19.3 Frequency and Time Domains 536 19.4 Fourier Integral and Transform 540 19.5 Discrete Fourier Transform (DFT) 542 19.6 Fast Fourier Transform (FFT) 544 19.7 The Power Spectrum 551 19.8 Curve Fitting with Software Packages 552 Problems 561 CHAPTER 20 Case Studies: Curve Fitting

563

20.1 Linear Regression and Population Models (Chemical/Bio Engineering) 563 20.2 Use of Splines to Estimate Heat Transfer (Civil/Environmental Engineering) 567 20.3 Fourier Analysis (Electrical Engineering) 569 20.4 Analysis of Experimental Data (Mechanical/Aerospace Engineering) 570 Problems 572 EPILOGUE: PART FIVE 582 PT5.4 Trade-Offs 582 PT5.5 Important Relationships and Formulas 583 PT5.6 Advanced Methods and Additional References

PART SIX NUMERICAL DIFFERENTIATION AND INTEGRATION 587

PT6.1 Motivation 587 PT6.2 Mathematical Background PT6.3 Orientation 599

597

CHAPTER 21 Newton-Cotes Integration Formulas 21.1 The Trapezoidal Rule 605 21.2 Simpson’s Rules 615 21.3 Integration with Unequal Segments 21.4 Open Integration Formulas 627 21.5 Multiple Integrals 627 Problems 629

603

624

584

x

CONTENTS CHAPTER 22 Integration of Equations

633

22.1 Newton-Cotes Algorithms for Equations 22.2 Romberg Integration 634 22.3 Adaptive Quadrature 640 22.4 Gauss Quadrature 642 22.5 Improper Integrals 650 Problems 653 CHAPTER 23 Numerical Differentiation

633

655

23.1 High-Accuracy Differentiation Formulas 655 23.2 Richardson Extrapolation 658 23.3 Derivatives of Unequally Spaced Data 660 23.4 Derivatives and Integrals for Data with Errors 661 23.5 Partial Derivatives 662 23.6 Numerical Integration/Differentiation with Software Packages Problems 670 CHAPTER 24 Case Studies: Numerical Integration and Differentiation

663

673

24.1 Integration to Determine the Total Quantity of Heat (Chemical/Bio Engineering) 673 24.2 Effective Force on the Mast of a Racing Sailboat (Civil/Environmental Engineering) 675 24.3 Root-Mean-Square Current by Numerical Integration (Electrical Engineering) 677 24.4 Numerical Integration to Compute Work (Mechanical/Aerospace Engineering) 680 Problems 684 EPILOGUE: PART SIX 694 PT6.4 Trade-Offs 694 PT6.5 Important Relationships and Formulas 695 PT6.6 Advanced Methods and Additional References

PART SEVEN ORDINARY DIFFERENTIAL EQUATIONS 699

PT7.1 Motivation 699 PT7.2 Mathematical Background PT7.3 Orientation 705

703

695

CONTENTS CHAPTER 25 Runge-Kutta Methods

xi

709

25.1 Euler’s Method 710 25.2 Improvements of Euler’s Method 721 25.3 Runge-Kutta Methods 729 25.4 Systems of Equations 739 25.5 Adaptive Runge-Kutta Methods 744 Problems 752 CHAPTER 26 Stiffness and Multistep Methods 26.1 Stiffness 755 26.2 Multistep Methods Problems 779

755

759

CHAPTER 27 Boundary-Value and Eigenvalue Problems

781

27.1 General Methods for Boundary-Value Problems 782 27.2 Eigenvalue Problems 789 27.3 Odes and Eigenvalues with Software Packages 801 Problems 808 CHAPTER 28 Case Studies: Ordinary Differential Equations

811

28.1 Using ODEs to Analyze the Transient Response of a Reactor (Chemical/Bio Engineering) 811 28.2 Predator-Prey Models and Chaos (Civil/Environmental Engineering) 818 28.3 Simulating Transient Current for an Electric Circuit (Electrical Engineering) 822 28.4 The Swinging Pendulum (Mechanical/Aerospace Engineering) 827 Problems 831 EPILOGUE: PART SEVEN 841 PT7.4 Trade-Offs 841 PT7.5 Important Relationships and Formulas 842 PT7.6 Advanced Methods and Additional References

PART EIGHT PARTIAL DIFFERENTIAL EQUATIONS 845

PT8.1 Motivation 845 PT8.2 Orientation 848

842

xii

CONTENTS CHAPTER 29 Finite Difference: Elliptic Equations

852

29.1 The Laplace Equation 852 29.2 Solution Technique 854 29.3 Boundary Conditions 860 29.4 The Control-Volume Approach 866 29.5 Software to Solve Elliptic Equations 869 Problems 870 CHAPTER 30 Finite Difference: Parabolic Equations

873

30.1 The Heat-Conduction Equation 873 30.2 Explicit Methods 874 30.3 A Simple Implicit Method 878 30.4 The Crank-Nicolson Method 882 30.5 Parabolic Equations in Two Spatial Dimensions Problems 888 CHAPTER 31 Finite-Element Method

885

890

31.1 The General Approach 891 31.2 Finite-Element Application in One Dimension 31.3 Two-Dimensional Problems 904 31.4 Solving PDEs with Software Packages 908 Problems 912

895

CHAPTER 32 Case Studies: Partial Differential Equations

915

32.1 One-Dimensional Mass Balance of a Reactor (Chemical/Bio Engineering) 915 32.2 Deflections of a Plate (Civil/Environmental Engineering) 919 32.3 Two-Dimensional Electrostatic Field Problems (Electrical Engineering) 921 32.4 Finite-Element Solution of a Series of Springs (Mechanical/Aerospace Engineering) 924 Problems 928 EPILOGUE: PART EIGHT 931 PT8.3 Trade-Offs 931 PT8.4 Important Relationships and Formulas 931 PT8.5 Advanced Methods and Additional References

932

CONTENTS

xiii

APPENDIX A: THE FOURIER SERIES

933

APPENDIX B: GETTING STARTED WITH MATLAB

935

APPENDIX C: GETTING STARTED WITH MATHCAD BIBLIOGRAPHY INDEX

957

954

943

PREFACE It has been over twenty years since we published the first edition of this book. Over that period, our original contention that numerical methods and computers would figure more prominently in the engineering curriculum—particularly in the early parts—has been dramatically borne out. Many universities now offer freshman, sophomore, and junior courses in both introductory computing and numerical methods. In addition, many of our colleagues are integrating computer-oriented problems into other courses at all levels of the curriculum. Thus, this new edition is still founded on the basic premise that student engineers should be provided with a strong and early introduction to numerical methods. Consequently, although we have expanded our coverage in the new edition, we have tried to maintain many of the features that made the first edition accessible to both lower- and upper-level undergraduates. These include: • Problem Orientation. Engineering students learn best when they are motivated by problems. This is particularly true for mathematics and computing. Consequently, we have approached numerical methods from a problem-solving perspective. • Student-Oriented Pedagogy. We have developed a number of features to make this book as student-friendly as possible. These include the overall organization, the use of introductions and epilogues to consolidate major topics and the extensive use of worked examples and case studies from all areas of engineering. We have also endeavored to keep our explanations straightforward and oriented practically. • Computational Tools. We empower our students by helping them utilize the standard “point-and-shoot” numerical problem-solving capabilities of packages like Excel, MATLAB, and Mathcad software. However, students are also shown how to develop simple, well-structured programs to extend the base capabilities of those environments. This knowledge carries over to standard programming languages such as Visual Basic, Fortran 90, and C/C11. We believe that the current flight from computer programming represents something of a “dumbing down” of the engineering curriculum. The bottom line is that as long as engineers are not content to be tool limited, they will have to write code. Only now they may be called “macros” or “M-files.” This book is designed to empower them to do that. Beyond these five original principles, the seventh edition has new and expanded problem sets. Most of the problems have been modified so that they yield different numerical solutions from previous editions. In addition, a variety of new problems have been included. The seventh edition also includes McGraw-Hill’s Connect® Engineering. This online homework management tool allows assignment of algorithmic problems for homework, quizzes, and tests. It connects students with the tools and resources they’ll need to achieve success. To learn more, visit www.mcgrawhillconnect.com. McGraw-Hill LearnSmart™ is also available as an integrated feature of McGraw-Hill Connect® Engineering. It is an adaptive learning system designed to help students learn faster, study more efficiently, and retain more knowledge for greater success. LearnSmart assesses xiv

PREFACE

xv

a student’s knowledge of course content through a series of adaptive questions. It pinpoints concepts the student does not understand and maps out a personalized study plan for success. Visit the following site for a demonstration. www.mhlearnsmart.com As always, our primary intent in writing this book is to provide students with a sound introduction to numerical methods. We believe that motivated students who enjoy numerical methods, computers, and mathematics will, in the end, make better engineers. If our book fosters an enthusiasm for these subjects, we will consider our efforts a success. Acknowledgments. We would like to thank our friends at McGraw-Hill. In particular, Lorraine Buczek and Bill Stenquist, who provided a positive and supportive atmosphere for creating this edition. As usual, Beatrice Sussman did a masterful job of copyediting the manuscript and Arpana Kumari of Aptara also did an outstanding job in the book’s final production phase. As in past editions, David Clough (University of Colorado), Mike Gustafson (Duke), and Jerry Stedinger (Cornell University) generously shared their insights and suggestions. Useful suggestions were also made by Bill Philpot (Cornell University), Jim Guilkey (University of Utah), Dong-Il Seo (Chungnam National University, Korea), Niall Broekhuizen (NIWA, New Zealand), and Raymundo Cordero and Karim Muci (ITESM, Mexico). The present edition has also benefited from the reviews and suggestions by the following colleagues: Betty Barr, University of Houston Jalal Behzadi, Shahid Chamran University Jordan Berg, Texas Tech University Jacob Bishop, Utah State University Estelle M. Eke, California State University, Sacramento Yazan A. Hussain, Jordan University of Science & Technology Yogesh Jaluria, Rutgers University S. Graham Kelly, The University of Akron Subha Kumpaty, Milwaukee School of Engineering Eckart Meiburg, University of California-Santa Barbara Prashant Mhaskar, McMaster University Luke Olson, University of Illinois at Urbana-Champaign Richard Pates Jr., Old Dominion University Joseph H. Pierluissi, University of Texas at El Paso Juan Perán, Universidad Nacional de Educación a Distancia (UNED) Scott A. Socolofsky, Texas A&M University It should be stressed that although we received useful advice from the aforementioned individuals, we are responsible for any inaccuracies or mistakes you may detect in this edition. Please contact Steve Chapra via e-mail if you should detect any errors in this edition. Finally, we would like to thank our family, friends, and students for their enduring patience and support. In particular, Cynthia Chapra, Danielle Husley, and Claire Canale are always there providing understanding, perspective, and love. Steven C. Chapra Medford, Massachusetts [email protected] Raymond P. Canale Lake Leelanau, Michigan

ABOUT THE AUTHORS 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 Surface Water-Quality Modeling and Applied Numerical Methods with MATLAB. Dr. Chapra 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 is a Fellow of the ASCE, and has received a number of awards for his scholarly contributions, including the Rudolph Hering Medal (ASCE), and the Meriam-Wiley Distinguished Author Award (American Society for Engineering Education). He has also been recognized as the outstanding teacher among the engineering faculties at Texas A&M University, the University of Colorado, and Tufts University. Raymond P. Canale is an emeritus professor at the University of Michigan. During his over 20-year career at the university, he taught numerous courses in the area of computers, numerical methods, and environmental engineering. He also directed extensive research programs in the area of mathematical and computer modeling of aquatic ecosystems. He has authored or coauthored several books and has published over 100 scientific papers and reports. He has also designed and developed personal computer software to facilitate engineering education and the solution of engineering problems. He has been given the Meriam-Wiley Distinguished Author Award by the American Society for Engineering Education for his books and software and several awards for his technical publications. Professor Canale is now devoting his energies to applied problems, where he works with engineering firms and industry and governmental agencies as a consultant and expert witness.

xvi

Numerical Methods for Engineers

PART ONE

MODELING, COMPUTERS, AND ERROR ANALYSIS PT1.1

MOTIVATION Numerical methods are techniques by which mathematical problems are formulated so that they can be solved with arithmetic operations. Although there are many kinds of numerical methods, they have one common characteristic: they invariably involve large numbers of tedious arithmetic calculations. It is little wonder that with the development of fast, efficient digital computers, the role of numerical methods in engineering problem solving has increased dramatically in recent years. PT1.1.1 Noncomputer Methods Beyond providing increased computational firepower, the widespread availability of computers (especially personal computers) and their partnership with numerical methods has had a significant influence on the actual engineering problem-solving process. In the precomputer era there were generally three different ways in which engineers approached problem solving: 1.

2.

3.

Solutions were derived for some problems using analytical, or exact, methods. These solutions were often useful and provided excellent insight into the behavior of some systems. However, analytical solutions can be derived for only a limited class of problems. These include those that can be approximated with linear models and those that have simple geometry and low dimensionality. Consequently, analytical solutions are of limited practical value because most real problems are nonlinear and involve complex shapes and processes. Graphical solutions were used to characterize the behavior of systems. These graphical solutions usually took the form of plots or nomographs. Although graphical techniques can often be used to solve complex problems, the results are not very precise. Furthermore, graphical solutions (without the aid of computers) are extremely tedious and awkward to implement. Finally, graphical techniques are often limited to problems that can be described using three or fewer dimensions. Calculators and slide rules were used to implement numerical methods manually. Although in theory such approaches should be perfectly adequate for solving complex problems, in actuality several difficulties are encountered. Manual calculations are slow and tedious. Furthermore, consistent results are elusive because of simple blunders that arise when numerous manual tasks are performed.

During the precomputer era, significant amounts of energy were expended on the solution technique itself, rather than on problem definition and interpretation (Fig. PT1.1a). This unfortunate situation existed because so much time and drudgery were required to obtain numerical answers using precomputer techniques. 3

4

MODELING, COMPUTERS, AND ERROR ANALYSIS

FORMULATION Fundamental laws explained briefly

SOLUTION

FIGURE PT1.1 The three phases of engineering problem solving in (a) the precomputer and (b) the computer era. The sizes of the boxes indicate the level of emphasis directed toward each phase. Computers facilitate the implementation of solution techniques and thus allow more emphasis to be placed on the creative aspects of problem formulation and interpretation of results.

Elaborate and often complicated method to make problem tractable

FORMULATION In-depth exposition of relationship of problem to fundamental laws

SOLUTION Easy-to-use computer method

INTERPRETATION

INTERPRETATION

In-depth analysis limited by timeconsuming solution

Ease of calculation allows holistic thoughts and intuition to develop; system sensitivity and behavior can be studied

(a)

(b)

Today, computers and numerical methods provide an alternative for such complicated calculations. Using computer power to obtain solutions directly, you can approach these calculations without recourse to simplifying assumptions or time-intensive techniques. Although analytical solutions are still extremely valuable both for problem solving and for providing insight, numerical methods represent alternatives that greatly enlarge your capabilities to confront and solve problems. As a result, more time is available for the use of your creative skills. Thus, more emphasis can be placed on problem formulation and solution interpretation and the incorporation of total system, or “holistic,” awareness (Fig. PT1.1b). PT1.1.2 Numerical Methods and Engineering Practice Since the late 1940s the widespread availability of digital computers has led to a veritable explosion in the use and development of numerical methods. At first, this growth was somewhat limited by the cost of access to large mainframe computers, and, consequently, many engineers continued to use simple analytical approaches in a significant portion of their work. Needless to say, the recent evolution of inexpensive personal

PT1.2

MATHEMATICAL BACKGROUND

5

computers has given us ready access to powerful computational capabilities. There are several additional reasons why you should study numerical methods: 1.

2.

3.

4.

5.

PT1.2

Numerical methods are extremely powerful problem-solving tools. They are capable of handling large systems of equations, nonlinearities, and complicated geometries that are not uncommon in engineering practice and that are often impossible to solve analytically. As such, they greatly enhance your problem-solving skills. During your careers, you may often have occasion to use commercially available prepackaged, or “canned,” computer programs that involve numerical methods. The intelligent use of these programs is often predicated on knowledge of the basic theory underlying the methods. 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. Numerical methods are an efficient vehicle for learning to use computers. It is well known that an effective way to learn programming is to actually write computer programs. Because numerical methods are for the most part designed for implementation on computers, they are ideal for this purpose. Further, they are especially well-suited to illustrate the power and the limitations of computers. 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 acknowledge and control the errors of approximation that are part and parcel of large-scale numerical calculations. Numerical methods provide a vehicle for you to reinforce your understanding of mathematics. 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.

MATHEMATICAL BACKGROUND Every part in this book requires some mathematical background. Consequently, the introductory material for each part includes a section, such as the one you are reading, on mathematical background. Because Part One itself is devoted to background material on mathematics and computers, this section does not involve a review of a specific mathematical topic. Rather, we take this opportunity to introduce you to the types of mathematical subject areas covered in this book. As summarized in Fig. PT1.2, these are 1.

2.

Roots of Equations (Fig. PT1.2a). These problems are concerned with the value of a variable or a parameter that satisfies a single nonlinear equation. These problems are especially valuable in engineering design contexts where it is often impossible to explicitly solve design equations for parameters. Systems of Linear Algebraic Equations (Fig. PT1.2b). These problems are similar in spirit to roots of equations in the sense that they are concerned with values that

6 FIGURE PT1.2 Summary of the numerical methods covered in this book.

MODELING, COMPUTERS, AND ERROR ANALYSIS

(a) Part 2: Roots of equations

f(x)

Solve f(x) = 0 for x.

Root x

(b) Part 3: Linear algebraic equations Given the a’s and the c’s, solve a11x1 + a12x2 = c1 a21x1 + a22x2 = c2

x2

Solution

for the x’s.

x1

(c) Part 4: Optimization Determine x that gives optimum f(x).

f(x)

Minimum x

(d) Part 5: Curve fitting f(x)

f(x)

Interpolation

Regression

(e) Part 6: Integration I = 兰ab f (x) dx Find the area under the curve.

x

x f(x)

I x

PT1.2 FIGURE PT1.2 (concluded)

MATHEMATICAL BACKGROUND

7

( f ) Part 7: Ordinary differential equations Given dy ⌬y ⯝ = f (t, y) dt ⌬t

y

solve for y as a function of t. yi + 1 = yi + f (ti , yi ) ⌬t

Slope = f(t i , y i )

⌬t ti

(g) Part 8: Partial differential equations

ti + 1

t

Given ⭸2u + ⭸2u = f(x, y) ⭸x2 ⭸y2

y

solve for u as a function of x and y

x

3.

4.

5.

satisfy equations. However, in contrast to satisfying a single equation, a set of values is sought that simultaneously satisfies a set of linear algebraic equations. Such equations arise in a variety of problem contexts and in all disciplines of engineering. In particular, they originate in the mathematical modeling of large systems of interconnected elements such as structures, electric circuits, and fluid networks. However, they are also encountered in other areas of numerical methods such as curve fitting and differential equations. Optimization (Fig. PT1.2c). These problems involve determining a value or values of an independent variable that correspond to a “best” or optimal value of a function. Thus, as in Fig. PT1.2c, optimization involves identifying maxima and minima. Such problems occur routinely in engineering design contexts. They also arise in a number of other numerical methods. We address both single- and multi-variable unconstrained optimization. We also describe constrained optimization with particular emphasis on linear programming. Curve Fitting (Fig. PT1.2d). You will often have occasion to fit curves to data points. The techniques developed for this purpose can be divided into two general categories: regression and interpolation. Regression is employed where there is a significant degree of error associated with the data. Experimental results are often of this kind. For these situations, the strategy is to derive a single curve that represents the general trend of the data without necessarily matching any individual points. In contrast, interpolation is used where the objective is to determine intermediate values between relatively error-free data points. Such is usually the case for tabulated information. For these situations, the strategy is to fit a curve directly through the data points and use the curve to predict the intermediate values. Integration (Fig. PT1.2e). As depicted, a physical interpretation of numerical integration is the determination of the area under a curve. Integration has many

8

MODELING, COMPUTERS, AND ERROR ANALYSIS

6.

7.

PT1.3

applications in engineering practice, ranging from the determination of the centroids of oddly shaped objects to the calculation of total quantities based on sets of discrete measurements. In addition, numerical integration formulas play an important role in the solution of differential equations. Ordinary Differential Equations (Fig. PT1.2f ). Ordinary differential equations are of great significance in engineering practice. This is because many physical laws are couched in terms of the rate of change of a quantity rather than the magnitude of the quantity itself. Examples range from population-forecasting models (rate of change of population) to the acceleration of a falling body (rate of change of velocity). Two types of problems are addressed: initial-value and boundary-value problems. In addition, the computation of eigenvalues is covered. Partial Differential Equations (Fig. PT1.2g). Partial differential equations are used to characterize engineering systems where the behavior of a physical quantity is couched in terms of its rate of change with respect to two or more independent variables. Examples include the steady-state distribution of temperature on a heated plate (two spatial dimensions) or the time-variable temperature of a heated rod (time and one spatial dimension). Two fundamentally different approaches are employed to solve partial differential equations numerically. In the present text, we will emphasize finite-difference methods that approximate the solution in a pointwise fashion (Fig. PT1.2g). However, we will also present an introduction to finite-element methods, which use a piecewise approach.

ORIENTATION Some orientation might be helpful before proceeding with our introduction to numerical methods. The following is intended as an overview of the material in Part One. In addition, some objectives have been included to focus your efforts when studying the material. PT1.3.1 Scope and Preview Figure PT1.3 is a schematic representation of the material in Part One. We have designed this diagram to provide you with a global overview of this part of the book. We believe that a sense of the “big picture” is critical to developing insight into numerical methods. When reading a text, it is often possible to become lost in technical details. Whenever you feel that you are losing the big picture, refer back to Fig. PT1.3 to reorient yourself. Every part of this book includes a similar figure. Figure PT1.3 also serves as a brief preview of the material covered in Part One. Chapter 1 is designed to orient you to numerical methods and to provide motivation by demonstrating how these techniques can be used in the engineering modeling process. Chapter 2 is an introduction and review of computer-related aspects of numerical methods and suggests the level of computer skills you should acquire to efficiently apply succeeding information. Chapters 3 and 4 deal with the important topic of error analysis, which must be understood for the effective use of numerical methods. In addition, an epilogue is included that introduces the trade-offs that have such great significance for the effective implementation of numerical methods.

PT1.3

ORIENTATION

PT 1.1 Motivation

9

PT 1.2 Mathematical background

PT 1.3 Orientation

PART 1 Modeling, Computers, and Error Analysis

PT 1.6 Advanced methods

1.1 A simple model

PT 1.5 Important formulas

CHAPTER 1 Mathematical Modeling and Engineering Problem Solving

EPILOGUE PT 1.4 Trade-offs

1.2 Conservation laws

2.1 Packages and programming

4.4 Miscellaneous errors

4.3 Total numerical error

2.2 Structured programming

CHAPTER 4 Truncation Errors and the Taylor Series

4.2 Error propagation

CHAPTER 2 Programming and Software

2.4 Excel

CHAPTER 3 Approximations and Round-Off Errors

4.1 Taylor series

2.5 MATLAB

2.7 Languages and libraries 2.6 Mathcad

3.4 Round-off errors

3.1 Significant figures 3.3 Error definitions

2.3 Modular programming

3.2 Accuracy and precision

FIGURE PT1.3 Schematic of the organization of the material in Part One: Modeling, Computers, and Error Analysis.

10

MODELING, COMPUTERS, AND ERROR ANALYSIS TABLE PT1.1 Specific study objectives for Part One. 1. Recognize the difference between analytical and numerical solutions. 2. Understand how conservation laws are employed to develop mathematical models of physical systems. 3. Define top-down and modular design. 4. Delineate the rules that underlie structured programming. 5. Be capable of composing structured and modular programs in a high-level computer language. 6. Know how to translate structured flowcharts and pseudocode into code in a high-level language. 7. Start to familiarize yourself with any software packages that you will be using in conjunction with this text. 8. Recognize the distinction between truncation and round-off errors. 9. Understand the concepts of significant figures, accuracy, and precision. 10. Recognize the difference between true relative error et, approximate relative error ea, and acceptable error es, and understand how ea and es are used to terminate an iterative computation. 11. Understand how numbers are represented in digital computers and how this representation induces round-off error. In particular, know the difference between single and extended precision. 12. Recognize how computer arithmetic can introduce and amplify round-off errors in calculations. In particular, appreciate the problem of subtractive cancellation. 13. Understand how the Taylor series and its remainder are employed to represent continuous functions. 14. Know the relationship between finite divided differences and derivatives. 15. Be able to analyze how errors are propagated through functional relationships. 16. Be familiar with the concepts of stability and condition. 17. Familiarize yourself with the trade-offs outlined in the Epilogue of Part One.

PT1.3.2 Goals and Objectives Study Objectives. Upon completing Part One, you should be adequately prepared to embark on your studies of numerical methods. In general, you should have gained a fundamental understanding of the importance of computers and the role of approximations and errors in the implementation and development of numerical methods. In addition to these general goals, you should have mastered each of the specific study objectives listed in Table PT1.1. Computer Objectives. Upon completing Part One, you should have mastered sufficient computer skills to develop your own software for the numerical methods in this text. You should be able to develop well-structured and reliable computer programs on the basis of pseudocode, flowcharts, or other forms of algorithms. You should have developed the capability to document your programs so that they may be effectively employed by users. Finally, in addition to your own programs, you may be using software packages along with this book. Packages like Excel, Mathcad, or The MathWorks, Inc. MATLAB® program are examples of such software. You should become familiar with these packages, so that you will be comfortable using them to solve numerical problems later in the text.

1

CHAPTER

Mathematical Modeling and Engineering Problem Solving Knowledge and understanding are prerequisites for the effective implementation of any tool. No matter how impressive your tool chest, you will be hard-pressed to repair a car if you do not understand how it works. This is particularly true when using computers to solve engineering problems. Although they have great potential utility, computers are practically useless without a fundamental understanding of how engineering systems work. This understanding is initially gained by empirical means—that is, by observation and experiment. However, while such empirically derived information is essential, it is only half the story. Over years and years of observation and experiment, engineers and scientists have noticed that certain aspects of their empirical studies occur repeatedly. Such general behavior can then be expressed as fundamental laws that essentially embody the cumulative wisdom of past experience. Thus, most engineering problem solving employs the two-pronged approach of empiricism and theoretical analysis (Fig. 1.1). It must be stressed that the two prongs are closely coupled. As new measurements are taken, the generalizations may be modified or new ones developed. Similarly, the generalizations can have a strong influence on the experiments and observations. In particular, generalizations can serve as organizing principles that can be employed to synthesize observations and experimental results into a coherent and comprehensive framework from which conclusions can be drawn. From an engineering problem-solving perspective, such a framework is most useful when it is expressed in the form of a mathematical model. The primary objective of this chapter is to introduce you to mathematical modeling and its role in engineering problem solving. We will also illustrate how numerical methods figure in the process.

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 general sense, it can be represented as a functional relationship of the form forcing independent Dependent , parameters, b 5f a variables functions variable

(1.1) 11

12

MATHEMATICAL MODELING AND ENGINEERING PROBLEM SOLVING

Problem definition

THEORY

Mathematical model

DATA

Problem-solving tools: computers, statistics, numerical methods, graphics, etc.

Numeric or graphic results Societal interfaces: scheduling, optimization, communication, public interaction, etc.

Implementation

FIGURE 1.1 The engineering problemsolving process.

where the dependent variable is a characteristic that usually 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 the system. 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 wellknown equation F 5 ma

(1.2)

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

1.1 FU

A SIMPLE MATHEMATICAL MODEL

13

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

F m

(1.3)

where a 5 the dependent variable reflecting the system’s behavior, F 5 the forcing function, and m 5 a parameter representing a property of the system. 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 several characteristics that are typical of mathematical models of the physical world: 1. 2. FD

FIGURE 1.2 Schematic diagram of the forces acting on a falling parachutist. FD is the downward force due to gravity. FU is the upward force due to air resistance.

3.

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 the mass of an object are known, Eq. (1.3) can be used to compute acceleration.

Because of its simple algebraic form, the solution of Eq. (1.2) can be 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 freefalling body near the earth’s surface. Our falling body will be a parachutist (Fig. 1.2). A model for this case can be derived by expressing the acceleration as the time rate of change of the velocity (dy兾dt) and substituting it into Eq. (1.3) to yield dy F 5 m dt

(1.4)

where y is velocity (m/s) and t is time (s). Thus, the mass multiplied by the rate of change of the velocity is equal to the net force acting on the body. 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 (Fig. 1.2), the net force is composed of two opposing forces: the downward pull of gravity FD and the upward force of air resistance FU: F 5 FD 1 FU

(1.5)

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

(1.6)

where g 5 the gravitational constant, or the acceleration due to gravity, which is approximately equal to 9.81 m/s2.

14

MATHEMATICAL MODELING AND ENGINEERING PROBLEM SOLVING

Air resistance can be formulated in a variety of ways. A simple approach is to assume that it is linearly proportional to velocity1 and acts in an upward direction, as in FU 5 2cy

(1.7)

where c 5 a proportionality constant called the drag coefficient (kg/s). Thus, the greater the fall velocity, the greater the upward force due to air resistance. The parameter c accounts for properties of the falling object, such as shape or surface roughness, that affect air resistance. For the present case, c might be a function of the type of jumpsuit or the orientation used by the parachutist 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 mg 2 cy dy 5 m dt

(1.8)

or simplifying the right side, dy c 5g2 y m dt

(1.9)

Equation (1.9) 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 (dy兾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.9) for the velocity of the falling parachutist 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 parachutist is initially at rest (y 5 0 at t 5 0), calculus can be used to solve Eq. (1.9) for y(t) 5

gm ( 1 2 e 2 (cym)t ) c

(1.10)

Note that Eq. (1.10) is cast in the general form of Eq. (1.1), where y(t) 5 the dependent variable, t 5 the independent variable, c and m 5 parameters, and g 5 the forcing function. EXAMPLE 1.1

Analytical Solution to the Falling Parachutist Problem Problem Statement. A parachutist of mass 68.1 kg jumps out of a stationary hot air balloon. Use Eq. (1.10) to compute velocity prior to opening the chute. The drag coefficient is equal to 12.5 kg/s. Solution. Inserting the parameters into Eq. (1.10) yields y(t) 5

9.81(68.1) ( 1 2 e 2 (12.5y68.1)t ) 5 53.44 ( 1 2 e 20.18355t ) 12.5

which can be used to compute 1

In fact, the relationship is actually nonlinear and might better be represented by a power relationship such as FU 5 2cy2. We will explore how such nonlinearities affect the model in problems at the end of this chapter.

1.1

A SIMPLE MATHEMATICAL MODEL t, s

15

v, m/s

0 2 4 6 8 10 12 `

0.00 16.42 27.80 35.68 41.14 44.92 47.54 53.44

According to the model, the parachutist accelerates rapidly (Fig. 1.3). A velocity of 44.92 m/s is attained after 10 s. Note also that after a sufficiently long time, a constant velocity, called the terminal velocity, of 53.44 m/s 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.10) is called an analytical, or exact, solution because it exactly satisfies 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. As mentioned previously, numerical methods are those in which the mathematical problem is reformulated so it can be solved by arithmetic operations. This can be illustrated

FIGURE 1.3 The analytical solution to the falling parachutist problem as computed in Example 1.1. Velocity increases with time and asymptotically approaches a terminal velocity.

Terminal velocity

v, m/s

40

20

0

0

4

8 t, s

12

16

MATHEMATICAL MODELING AND ENGINEERING PROBLEM SOLVING

v(ti +1) True slope dv/dt v Approximate slope v v(ti +1) – v(ti ) = t –t i +1 i t

v(ti )

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

ti +1

ti

t

t

for Newton’s second law by realizing that the time rate of change of velocity can be approximated by (Fig. 1.4): y(ti11 ) 2 y(ti ) dy ¢y > 5 dt ¢t ti11 2 ti

(1.11)

where Dy and Dt 5 differences in velocity and time, respectively, computed over finite intervals, y(ti) 5 velocity at an initial time ti, and y(ti+1) 5 velocity at some later time ti+1. Note that dy/dt > ¢yy¢t is approximate because Dt is finite. Remember from calculus that dy ¢y 5 lim dt ¢tS0 ¢t Equation (1.11) represents the reverse process. Equation (1.11) is called a finite divided difference approximation of the derivative at time ti. It can be substituted into Eq. (1.9) to give y(ti11 ) 2 y(ti ) c 5 g 2 y(ti ) m ti11 2 ti This equation can then be rearranged to yield y(ti11 ) 5 y(ti ) 1 c g 2

c y(t ) d (ti11 2 ti ) m i

(1.12)

Notice that the term in brackets is the right-hand side of the differential equation itself [Eq. (1.9)]. That is, it provides a means to compute the rate of change or slope of y. Thus, the differential equation has been transformed into an equation that can be used to determine the velocity algebraically at ti11 using the slope and previous values of

1.1

A SIMPLE MATHEMATICAL MODEL

17

y and t. If you are given an initial value for velocity at some time ti, you can easily compute velocity at a later time ti11. This new value of velocity at ti11 can in turn be employed to extend the computation to velocity at ti12 and so on. Thus, at any time along the way, New value 5 old value 1 slope 3 step size Note that this approach is formally called Euler’s method. EXAMPLE 1.2

Numerical Solution to the Falling Parachutist Problem Problem Statement. Perform the same computation as in Example 1.1 but use Eq. (1.12) to compute the velocity. Employ a step size of 2 s for the calculation. Solution. At the start of the computation (ti 5 0), the velocity of the parachutist is zero. Using this information and the parameter values from Example 1.1, Eq. (1.12) can be used to compute velocity at ti11 5 2 s: y 5 0 1 c 9.81 2

12.5 (0) d 2 5 19.62 m/s 68.1

For the next interval (from t 5 2 to 4 s), the computation is repeated, with the result y 5 19.62 1 c 9.81 2

12.5 (19.62)d 2 5 32.04 m/s 68.1

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

v, m/s

0 2 4 6 8 10 12 `

0.00 19.62 32.04 39.90 44.87 48.02 50.01 53.44

The results are plotted in Fig. 1.5 along with the exact solution. It can be seen that the numerical 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 l-s intervals results in a smaller error, as the straight-line segments track closer to the true solution. Using hand calculations, 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 falling parachutist without having to solve the differential equation exactly. As in the previous example, a computational price must be paid for a more accurate numerical result. Each halving of the step size to attain more accuracy leads to a doubling

18

MATHEMATICAL MODELING AND ENGINEERING PROBLEM SOLVING

Terminal velocity Approximate, numerical solution

v, m/s

40

Exact, analytical solution 20

FIGURE 1.5 Comparison of the numerical and analytical solutions for the falling parachutist problem.

0

0

4

8

12

t, s

of the number of computations. Thus, we see that there is a trade-off between accuracy and computational effort. Such trade-offs figure prominently in numerical methods and constitute an important theme of this book. Consequently, we have devoted the Epilogue of Part One to an introduction to more of these trade-offs.

1.2

CONSERVATION LAWS AND ENGINEERING Aside from Newton’s second law, there are other major organizing principles in 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 5 increases 2 decreases

(1.13)

This is precisely the format that we employed when using Newton’s law to develop a force balance for the falling parachutist [Eq. (1.8)]. Although simple, Eq. (1.13) embodies one of the most fundamental ways in which conservation laws are used in engineering—that is, to predict changes with respect to time. We give Eq. (1.13) the special name 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.13) becomes Change 5 0 5 increases 2 decreases or Increases 5 decreases

(1.14)

1.2

CONSERVATION LAWS AND ENGINEERING

19

Pipe 2 Flow in = 80

Pipe 1 Flow in = 100

FIGURE 1.6 A flow balance for steady-state incompressible fluid flow at the junction of pipes.

Pipe 4 Flow out = ?

Pipe 3 Flow out = 120

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 computation—has many applications in engineering. 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 5 flow out For the junction in Fig. 1.6, the balance can be used to compute that the flow out of the fourth pipe must be 60. For the falling parachutist, steady-state conditions would correspond to the case where the net force was zero, or [Eq. (1.8) with dy兾dt 5 0] mg 5 cy

(1.15)

Thus, at steady state, the downward and upward forces are in balance, and Eq. (1.15) can be solved for the terminal velocity y5

mg c

Although Eqs. (1.13) and (1.14) might appear trivially simple, they embody the two fundamental ways that conservation laws are employed in engineering. As such, they will form an important part of our efforts in subsequent chapters to illustrate the connection between numerical methods and engineering. Our primary vehicles for making this connection are the engineering applications that appear at the end of each part of this book. Table 1.1 summarizes some of the simple engineering models and associated conservation laws that will form the basis for many of these engineering applications. Most of the chemical engineering applications will focus on 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. Both the civil and mechanical engineering applications will focus on models developed from the conservation 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 applications to analyze the transient up-and-down motion or vibrations of an automobile.

20

MATHEMATICAL MODELING AND ENGINEERING 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 upon which the balance is based is specified. Field

Device

Chemical engineering

Organizing Principle

Mathematical Expression

Conservation of mass

Mass balance: Input

Reactors

Output

Over a unit of time period mass = inputs – outputs Civil engineering

Conservation of momentum

Structure

Force balance: + FV – FH

+ FH – FV

At each node  horizontal forces (FH ) = 0  vertical forces (FV ) = 0 Mechanical engineering

Machine

Conservation of momentum

Force balance:

Upward force x=0 Downward force

2 m d x2 = downward force – upward force dt

Electrical engineering

Conservation of charge Current balance: + –

For each node  current (i ) = 0

+ i1

– i3 + i2

Circuit i1R1

Conservation of energy Voltage balance:

␰

i2R2 i3R3

Around each loop  emf’s –  voltage drops for resistors = 0   –  iR = 0

PROBLEMS

21

TABLE 1.2 Some practical issues that will be explored in the engineering applications at the end of each part of this book. 1. Nonlinear versus linear. Much of classical engineering depends on linearization to permit analytical solutions. Although this is often appropriate, expanded insight can often be gained if nonlinear problems are examined. 2. Large versus small systems. Without a computer, it is often not feasible to examine systems with over three interacting components. With computers and numerical methods, more realistic multicomponent systems can be examined. 3. Nonideal versus ideal. Idealized laws abound in engineering. Often there are nonidealized alternatives that are more realistic but more computationally demanding. Approximate numerical approaches can facilitate the application of these nonideal relationships. 4. Sensitivity analysis. Because they are so involved, many manual calculations require a great deal of time and effort for successful implementation. This sometimes discourages the analyst from implementing the multiple computations that are necessary to examine how a system responds under different conditions. Such sensitivity analyses are facilitated when numerical methods allow the computer to assume the computational burden. 5. Design. It is often a straightforward proposition to determine the performance of a system as a function of its parameters. It is usually more difficult to solve the inverse problem—that is, determining the parameters when the required performance is specified. Numerical methods and computers often permit this task to be implemented in an efficient manner.

Finally, the electrical engineering applications employ both current and energy balances to model electric circuits. The current balance, which results from the conservation of charge, is similar in spirit to the flow balance depicted in Fig. 1.6. 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. The engineering applications are designed to illustrate how numerical methods are actually employed in the engineering problem-solving process. As such, they will permit us to explore practical issues (Table 1.2) that arise in real-world applications. Making these connections between mathematical techniques such as numerical methods and engineering practice is a critical step in tapping their true potential. Careful examination of the engineering applications will help you to take this step.

PROBLEMS 1.1 Use calculus to solve Eq. (1.9) for the case where the initial velocity, y(0) is nonzero. 1.2 Repeat Example 1.2. Compute the velocity to t 5 8 s, with a step size of (a) 1 and (b) 0.5 s. Can you make any statement regarding the errors of the calculation based on the results? 1.3 Rather than the linear relationship of Eq. (1.7), you might choose to model the upward force on the parachutist as a secondorder relationship, FU 5 2c¿y2 where c9 5 a bulk second-order drag coefficient (kg/m). (a) Using calculus, obtain the closed-form solution for the case where the jumper is initially at rest (y 5 0 at t 5 0).

(b) Repeat the numerical calculation in Example 1.2 with the same initial condition and parameter values, but with second-order drag. Use a value of 0.22 kg/m for c9. 1.4 For the free-falling parachutist with linear drag, assume a first jumper is 70 kg and has a drag coefficient of 12 kg/s. If a second jumper has a drag coefficient of 15 kg/s and a mass of 80 kg, how long will it take him to reach the same velocity the first jumper reached in 9 s? 1.5 Compute the velocity of a free-falling parachutist using Euler’s method for the case where m 5 80 kg and c 5 10 kg/s. Perform the calculation from t 5 0 to 20 s with a step size of 1 s. Use an initial condition that the parachutist has an upward velocity of 20 m/s at t 5 0. At t 5 10 s, assume that the chute is instantaneously deployed so that the drag coefficient jumps to 60 kg/s.

22

MATHEMATICAL MODELING AND ENGINEERING PROBLEM SOLVING

1.6 The following information is available for a bank account: Date

Deposits

Withdrawals

220.13

327.26

216.80

378.51

450.35

106.80

127.31

350.61

5/1

Interest

Balance 1522.33

6/1 7/1 8/1

(b) Plot the solution on a semilog graph (i.e., ln c versus t) and determine the slope. Interpret your results. 1.8 A group of 35 students attend a class in a room that measures 11 m by 8 m by 3 m. Each student takes up about 0.075 m3 and gives out about 80 W of heat (1 W 5 1 J/s). Calculate the air temperature rise during the first 20 minutes of the class if the room is completely sealed and insulated. Assume the heat capacity, Cy, for air is 0.718 kJ/(kg K). Assume air is an ideal gas at 208C and 101.325 kPa. Note that the heat absorbed by the air Q is related to the mass of the air m, the heat capacity, and the change in temperature by the following relationship:

9/1

Q5m Note that the money earns interest which is computed as

#

T2

CydT 5 mCy (T2 2 T1 )

T1

The mass of air can be obtained from the ideal gas law:

Interest 5 i Bi where i 5 the interest rate expressed as a fraction per month, and Bi the initial balance at the beginning of the month. (a) Use the conservation of cash to compute the balance on 6/1, 7/1, 8/1, and 9/1 if the interest rate is 1% per month (i 5 0.01/month). Show each step in the computation. (b) Write a differential equation for the cash balance in the form dB 5 f (D(t), W(t), i) dt where t 5 time (months), D(t) 5 deposits as a function of time ($/month), W(t) 5 withdrawals as a function of time ($/month). For this case, assume that interest is compounded continuously; that is, interest 5 iB. (c) Use Euler’s method with a time step of 0.5 month to simulate the balance. Assume that the deposits and withdrawals are applied uniformly over the month. (d) Develop a plot of balance versus time for (a) and (c). 1.7 The amount of a uniformly distributed radioactive contaminant contained in a closed reactor is measured by its concentration c (becquerel/liter or Bq/L). The contaminant decreases at a decay rate proportional to its concentration—that is,

PV 5

m RT MwT

where P is the gas pressure, V is the volume of the gas, Mwt is the molecular weight of the gas (for air, 28.97 kg/kmol), and R is the ideal gas constant [8.314 kPa m3/(kmol K)]. 1.9 A storage tank contains a liquid at depth y, where y 5 0 when the tank is half full. Liquid is withdrawn at a constant flow rate Q to meet demands. The contents are resupplied at a sinusoidal rate 3Q sin2(t). y

0

FIGURE P1.9

decay rate 5 2kc where k is a constant with units of day21. Therefore, according to Eq. (1.13), a mass balance for the reactor can be written as dc dt

5

Equation (1.13) can be written for this system as d(Ay) dt

2kc

change decrease a b5a b in mass by decay (a) Use Euler’s method to solve this equation from t 5 0 to 1 d with k 5 0.175d21. Employ a step size of Dt 5 0.1. The concentration at t 5 0 is 100 Bq/L.

a

5 3Q sin 2 (t) 2

Q

change in b 5 (inflow)  2 (outflow) volume

or, since the surface area A is constant dy Q Q 5 3 sin 2 (t) 2 dt A A

PROBLEMS

23

Use Euler’s method to solve for the depth y from t 5 0 to 10 d with a step size of 0.5 d. The parameter values are A 5 1250 m2 and Q 5 450 m3/d. Assume that the initial condition is y 5 0. 1.10 For the same storage tank described in Prob. 1.9, suppose that the outflow is not constant but rather depends on the depth. For this case, the differential equation for depth can be written as dy a(1 1 y) 1.5 Q 5 3 sin 2 (t) 2 dt A A Use Euler’s method to solve for the depth y from t 5 0 to 10 d with a step size of 0.5 d. The parameter values are A 5 1250 m2, Q 5 450 m3/d, and a 5 150. Assume that the initial condition is y 5 0. 1.11 Apply the conservation of volume (see Prob. 1.9) to simulate the level of liquid in a conical storage tank (Fig. P1.11). The liquid flows in at a sinusoidal rate of Qin 5 3 sin2(t) and flows out according to Qout 5 3(y 2 yout ) 1.5

y . yout

Qout 5 0

y # yout

where flow has units of m3/d and y 5 the elevation of the water surface above the bottom of the tank (m). Use Euler’s method to solve for the depth y from t 5 0 to 10 d with a step size of 0.5 d. The parameter values are rtop 5 2.5 m, ytop 5 4 m, and yout 5 1 m. Assume that the level is initially below the outlet pipe with y(0) 5 0.8 m.

rtop

Qin

ytop

y yout

s

Qout 1

0

FIGURE P1.11

1.12 In our example of the free-falling parachutist, we assumed that the acceleration due to gravity was a constant value. Although this is a decent approximation when we are examining falling objects near the surface of the earth, the gravitational force decreases as we move above sea level. A more general representation based on Newton’s inverse square law of gravitational attraction can be written as g(x) 5 g(0)

R2 (R 1 x) 2

where g(x) 5 gravitational acceleration at altitude x (in m) measured upward from the earth’s surface (m/s2), g(0) 5 gravitational acceleration at the earth’s surface (> 9.81 m/s2), and R 5 the earth’s radius (> 6.37 3 106 m). (a) In a fashion similar to the derivation of Eq. (1.9) use a force balance to derive a differential equation for velocity as a function of time that utilizes this more complete representation of gravitation. However, for this derivation, assume that upward velocity is positive. (b) For the case where drag is negligible, use the chain rule to express the differential equation as a function of altitude rather than time. Recall that the chain rule is dy dy dx 5 dt dx dt (c) Use calculus to obtain the closed form solution where y 5 y0 at x 5 0. (d) Use Euler’s method to obtain a numerical solution from x 5 0 to 100,000 m using a step of 10,000 m where the initial velocity is 1500 m/s upward. Compare your result with the analytical solution. 1.13 Suppose that a spherical droplet of liquid evaporates at a rate that is proportional to its surface area. dV 5 2kA dt where V 5 volume (mm3), t 5 time (min), k 5 the evaporation rate (mm/min), and A 5 surface area (mm2). Use Euler’s method to compute the volume of the droplet from t 5 0 to 10 min using a step size of 0.25 min. Assume that k 5 0.08 mm/min and that the droplet initially has a radius of 2.5 mm. Assess the validity of your results by determining the radius of your final computed volume and verifying that it is consistent with the evaporation rate. 1.14 Newton’s law of cooling says that the temperature of a body changes at a rate proportional to the difference between its temperature and that of the surrounding medium (the ambient temperature), dT 5 2k(T 2 Ta ) dt where T 5 the temperature of the body (8C), t 5 time (min), k 5 the proportionality constant (per minute), and Ta 5 the ambient temperature (8C). Suppose that a cup of coffee originally has a temperature of 708C. Use Euler’s method to compute the temperature from t 5 0 to 10 min using a step size of 2 min if Ta 5 208C and k 5 0.019/min. 1.15 As depicted in Fig. P1.15, an RLC circuit consists of three elements: a resistor (R), and inductor (L) and a capacitor (C). The flow of current across each element induces a voltage drop.

24

MATHEMATICAL MODELING AND ENGINEERING PROBLEM SOLVING

Kirchhoff’s second voltage law states that the algebraic sum of these voltage drops around a closed circuit is zero, iR 1 L

Q1

Q3

Q5

q di 1 50 dt C

where i 5 current, R 5 resistance, L 5 inductance, t 5 time, q 5 charge, and C 5 capacitance. In addition, the current is related to charge as in

Q2

Q4

Q6

Q7

dq 5i dt (a) If the initial values are i(0) 5 0 and q(0) 5 1 C, use Euler’s method to solve this pair of differential equations from t 5 0 to 0.1 s using a step size of Dt 5 0.01 s. Employ the following parameters for your calculation: R 5 200 V, L 5 5 H, and C 5 10–4 F. (b) Develop a plot of i and q versus t. Resistor Inductor Capacitor q di L c iR dt

i

FIGURE P1.15

1.16 Cancer cells grow exponentially with a doubling time of 20 h when they have an unlimited nutrient supply. However, as the cells start to form a solid spherical tumor without a blood supply, growth at the center of the tumor becomes limited, and eventually cells start to die. (a) Exponential growth of cell number N can be expressed as shown, where m is the growth rate of the cells. For cancer cells, find the value of m. dN 5 mN dt (b) Write an equation that will describe the rate of change of tumor volume during exponential growth given that the diameter of an individual cell is 20 microns. (c) After a particular type of tumor exceeds 500 microns in diameter, the cells at the center of the tumor die (but continue to take up space in the tumor). Determine how long it will take for the tumor to exceed this critical size. 1.17 A fluid is pumped into the network shown in Fig. P1.17. If Q2 5 0.6, Q3 5 0.4, Q7 5 0.2, and Q8 5 0.3 m3/s, determine the other flows.

Q10

Q9

Q8

FIGURE P1.17

1.18 The velocity is equal to the rate of change of distance x (m), dx 5 y(t) dt

(P1.18)

(a) Substitute Eq. (1.10) and develop an analytical solution for distance as a function of time. Assume that x(0) 5 0. (b) Use Euler’s method to numerically integrate Eqs. (P1.18) and (1.9) in order to determine both the velocity and distance fallen as a function of time for the first 10 s of free-fall using the same parameters as in Example 1.2. (c) Develop a plot of your numerical results together with the analytical solution. 1.19 You are working as a crime-scene investigator and must predict the temperature of a homicide victim over a 5-hr period. You know that the room where the victim was found was at 108C when the body was discovered. (a) Use Newton’s law of cooling (Prob. 1.14) and Euler’s method to compute the victim’s body temperature for the 5-hr period using values of k 5 0.12/hr and Dt 5 0.5 hr. Assume that the victim’s body temperature at the time of death was 378C, and that the room temperature was at a constant value of 108C over the 5-hr period. (b) Further investigation reveals that the room temperature had actually dropped linearly from 20 to 108C over the 5-hr period. Repeat the same calculation as in (a) but incorporate this new information. (c) Compare the results from (a) and (b) by plotting them on the same graph. 1.20 Suppose that a parachutist with linear drag (m 5 70 kg, c 5 12.5 kg/s) jumps from an airplane flying at an altitude of a kilometer with a horizontal velocity of 180 m/s relative to the ground. (a) Write a system of four differential equations for x, y, yx 5 dx/dt and yy 5 dy/dt.

PROBLEMS (b) If the initial horizontal position is defined as x 5 0, use Euler’s methods with Dt 5 1 s to compute the jumper’s position over the first 10 s. (c) Develop plots of y versus t and y versus x. Use the plot to graphically estimate when and where the jumper would hit the ground if the chute failed to open. 1.21 As noted in Prob. 1.3, drag is more accurately represented as depending on the square of velocity. A more fundamental representation of the drag force, which assumes turbulent conditions (i.e., a high Reynolds number), can be formulated as

25 (b) At steady-state, use this equation to solve for the particle’s terminal velocity. (c) Employ the result of (b) to compute the particle’s terminal velocity in m/s for a spherical silt particle settling in water: d 5 10 mm, r 5 1 g/cm3, rs 5 2.65 g/cm3, and m 5 0.014 g/(cm?s). (d) Check whether flow is laminar. (e) Use Euler’s method to compute the velocity from t 5 0 to 2215 s with Dt 5 2218 s given the parameters given previously along with the initial condition: y (0) 5 0. FD FB