Numerical Methods for Engineers Second Edition Numerical Methods for Engineers Second Edition D. Griffiths Colorado School of Mines Golden, Colorado I. Smith University of Manchester Manchester, UK Front cover illustration: Stiffness matrix eigenmodes of a 20-node hexahedral finite element after Smith and Kidger (1991), computed using Program 4. Courtesy of Dr.
CRC Press Taylor & Francis Group 6000 Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487-2742 © 2006 by Taylor & Francis Group, LLC CRC Press is an imprint of Taylor & Francis Group, an Informa business No claim to original U. Government works Version Date: 20110713 International Standard Book Number-13: 978-1-4200-1024-4 (eBook - PDF) This book contains information obtained from authentic and highly regarded sources. Reasonable efforts have been made to publish reliable data and information, but the author and publisher cannot assume responsibility for the validity of all materials or the consequences of their use. The authors and publishers have attempted to trace the copyright holders of all material reproduced in this publication and apologize to copyright holders if permission to publish in this form has not been obtained.
If any copyright material has not been acknowledged please write and let us know so we may rectify in any future reprint. Except as permitted under U. Copyright Law, no part of this book may be reprinted, reproduced, transmitted, or utilized in any form by any electronic, mechanical, or other means, now known or hereafter invented, including photocopying, microfilming, and recording, or in any information stor- age or retrieval system, without written permission from the publishers. For permission to photocopy or use material electronically from this work, please access www.com (http://www.com/) or contact the Copyright Clearance Center, Inc.
(CCC), 222 Rosewood Drive, Danvers, MA 01923, 978-750-8400. CCC is a not-for-profit organization that pro- vides licenses and registration for a variety of users. For organizations that have been granted a pho- tocopy license by the CCC, a separate system of payment has been arranged. Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are used only for identification and explanation without intent to infringe.
Visit the Taylor & Francis Web site at http://www.com and the CRC Press Web site at http://www.com To Valerie, Will and James List of Programs 2.1 Gaussian elimination for linear simultaneous equations 17 2.2 Gaussian elimination using [L][U] factorization 23 T 2.3 Gaussian elimination using [L][D][L] factorization 28 T 2.4 Cholesky [L][L] factorization using banded storage 34 T 2.5 Cholesky [L][L] factorization using skyline storage 36 2.6 [L][U] factorization with pivoting 40 2.7 Cholesky [L][L]T factorization using skyline storage 43 prescribed solutions by penalty method 2.8 Jacobi iteration for linear simultaneous equations 49 2.9 Gauss-Seidel iteration for linear simultaneous equations 54 2.10 Successive overrelaxation for linear simultaneous equations 58 2.11 Steepest descent for linear simultaneous equations 62 2.12 Conjugate gradients for linear simultaneous equations 65 2.13 Stabilized bi-conjugate gradients for 70 linear simultaneous equations 2.14 Preconditioned conjugate gradients for 74 linear simultaneous equations 2.1 Iterative substitution for a single root 92 3.2 Bisection method for a single root 98 3.3 False position method for a single root 102 3.4 Newton-Raphson method for a single root 106 3.5 Modified Newton-Raphson method for a single root 109 3.6 Iterative substitution for systems of equations 114 3.7 Newton-Raphson for systems of equations 119 3.8 Modified Newton-Raphson for systems of equations 123 List of Programs (continued) 4.1 Vector iteration for “largest” eigenvalue 137 and its eigenvector 4.2 Shifted vector iteration for eigenvalue 140 and its eigenvector 4.3 Shifted inverse iteration for nearest eigenvalue 145 and its eigenvector 4.4 Vector iteration for [K]{x} = λ[M]{x} 151 4.5 Conversion of [K]{x} = λ[M]{x} 155 to symmetrical standard form 4.6 Jacobi diagonalization for eigenvalues 163 of symmetrical matrices 4.7 Householder reduction of symmetrical matrix 169 to tridiagonal form 4.8 Lanczos reduction of symmetrical matrix 173 to tridiagonal form 4.9 [L][R] transformation for eigenvalues 177 4.10 Characteristic polynomial method 184 for eigenvalues of symmetrical tridiagonal matrix 5.1 Interpolation by Lagrangian polynomials 196 5.2 Interpolation by forward differences 204 5.3 Interpolation by cubic spline functions 211 5.4 Curve fitting by least squares 232 6.1 Repeated Newton-Cotes rules 262 6.2 Repeated Gauss-Legendre rules 275 6.3 Adaptive Gauss-Legendre rules 280 6.4 Gauss-Laguerre rules 287 6.5 Multiple integrals by Gauss-Legendre rules 302 7.1 One-step methods for systems of ODEs 337 7.2 Theta-method for linear ODEs 345 7.3 Fourth order predictor-corrector methods 356 7.4 Shooting method for second order ODEs 371 8.1 Explicit finite differences in 1D 424 8.2 Simple FE analysis of Example 8.3 432 Contents 1 Introduction and Programming Preliminaries 1 1.4 External Fortran subprogram libraries .5 A simple Fortran program .6 Some simple Fortran constructs .8 User-supplied functions and subroutines .9 Errors and accuracy .4 Intrinsic and library-supplied precision routines. 13 2 Linear Algebraic Equations 15 2.1 Observations on the elimination process .3 Equation solution using factorization .1 Observations on the solution process by factorization .4 Equations with a symmetrical coefficient matrix .1 Quadratic form and positive definiteness .6 Compact storage for variable bandwidths .8 Equations with prescribed solutions .1 The iterative process .2 Very sparse systems .3 The Gauss-Seidel method .1 The method of ‘steepest descent’ .2 The method of ‘conjugate gradients’ .3 Convergence of iterative methods .13 Comparison of direct and iterative methods .3 Multiple roots and other difficulties .2 False position method .1 Newton-Raphson method .2 A modified Newton-Raphson method .6 Acceleration of convergence .7 Systems of nonlinear equations .1 Iterative substitution for systems .2 Newton-Raphson for systems .3 Modified Newton-Raphson method for systems .1 Orthogonality and normalization of eigenvectors .2 Properties of eigenvalues and eigenvectors .3 Solution methods for eigenvalue equations .1 Shifted vector iteration .2 Shifted inverse iteration .3 Intermediate eigenvalues by deflation .4 The generalized eigenvalue problem [K]{x} = λ[M]{x} .1 Conversion of generalized problem to symmetrical stan- dard form .1 Comments on Jacobi diagonalization .2 Householder’s transformation to tridiagonal form .3 Lanczos transformation to tridiagonal form .4 LR transformation for eigenvalues of tridiagonal matri- ces .6 Characteristic polynomial methods .1 Evaluating determinants of tridiagonal matrices .2 The Sturm sequence property .3 General symmetrical matrices, e. 188 5 Interpolation and Curve Fitting 193 5.2 Difference methods .3 Difference methods with equal intervals .3 Interpolation using cubic spline functions .4 Numerical differentiation .1 Interpolating polynomial method .2 Taylor series method .2 Linearization of data .2 Newton-Cotes rules .5 Higher order Newton-Cotes rules (n > 3) .6 Accuracy of Newton-Cotes rules .7 Summary of Newton-Cotes rules .8 Repeated Newton-Cotes rules .9 Remarks on Newton-Cotes rules .3 Gauss-Legendre rules .3 Two-point Gauss-Legendre rule, (n = 2) .4 Three-point Gauss-Legendre rule, (n = 3) .5 Changing the limits of integration .6 Accuracy of Gauss-Legendre rules .4 Adaptive integration rules .5 Special integration rules .1 Gauss-Chebyshev rules .2 Fixed weighting coefficients .4 Sampling points outside the range of integration .2 Integration over a general quadrilateral area. 307 7 Numerical Solution of Ordinary Differential Equations 317 7.2 Definitions and types of ODE .3 Initial value problems .1 One-step methods .2 Reduction of high order equations .3 Solution of simultaneous first order equations .4 θ-methods for linear equations .5 Predictor-corrector methods .7 Error propagation and numerical stability .8 Concluding remarks on initial value problems .4 Boundary value problems .1 Finite difference methods .3 Weighted residual methods.
386 8 Introduction to Partial Differential Equations 393 8.2 Definitions and types of PDE .3 First order equations .4 Second order equations .5 Finite difference method .6 Finite element method. 434 A Descriptions of Library Subprograms 445 B Fortran 95 Listings of Library Subprograms 447 C References and Additional Reading 469 Index 475 Chapter 1 Introduction and Programming Preliminaries 1.1 Introduction There are many existing texts aimed at introducing engineers to the use of the numerical methods which underpin much of modern engineering prac- tice. Some contain “pseudocode” to illustrate how algorithms work, while others rely on commercial “packaged” software such as “Mathematica” or “MATLAB”. But the vast majority of the large computer programs which lie behind the design of engineering systems are written in the Fortran language, hence the use of the latest dialect, Fortran 95, in this book.
Nearly fifty entire programs are listed, many being made more concise through the use of “li- braries” of approximately twenty five subprograms which are also described and listed in full in Appendices A and B. Free Fortran 95 compilers are also widely available and users are therefore encouraged to modify the programs and develop new ones to suit their particular needs.2 Running programs Chapters 2-8 of this textbook describe 49 Fortran 95 programs covering a wide range of numerical methods applications. Many of the programs make use of a subprogram library called nm_lib which holds 23 subroutines and functions. In addition, there is a module called precision which controls the precision of the calculations.
Detailed instructions for running the programs described in this textbook are to be found at the web site: www.edu/~vgriffit/NM After linking to this site, consult the readme.txt file for information on: 1) how to download all the main programs, subprograms and sample data files, 1 2 Numerical Methods for Engineers 2) how to download a free Fortran 95 compiler, e., from the Free Software Foundation at www.org, 3) how to compile and execute programs using a simple batch file In the interests of generality and portability, all programs in the book as- sume that the data and results files have the generic names nm95.dat and nm95. A batch file that can be downloaded from the web with the main programs and libraries called run2.bat copies the actual data file (say fred.dat) to nm95.dat before execution. Finally, after the program has run, the generic results file fe95.res is copied to the actual results file name (say fred. If users wish to use the run2.bat batch file, data files and results files must have the extensions *.
See the readme.txt file as described above for more details.3 Hardware The use of Fortran means that the numerical methods described can be “ported” or transferred, with minimal alteration, from one computer to an- other. For “small” problems and teaching purposes a PC will suffice whereas Smith and Griffiths (2004) show how the same philosophy can be used to solve “large” problems, involving many millions of unknowns, on parallel “clusters” of PCs or on “supercomputers”.4 External Fortran subprogram libraries Another advantage of the use of Fortran is the existence of extensive sub- program libraries already written in that language. These can therefore easily be inserted into programs such as those described herein, avoiding unneces- sary duplication of effort. Some libraries are provided by manufacturers of computers, while others are available commercially, although academic use is often free.
A good example is the NAG subroutine library (www.uk) which contains over 600 subroutines. These cover a wide range of numerical appli- cations, and are organized into “Chapters” which are listed in Table 1. The NAG “Chapters” can very broadly be classified into “deterministic” numerical analyses, statistical analyses and “utility” routines. The present book deals only with the first of these classes and even then with a subset.
In the first column of Table 1.1 the number in parentheses indicates the chapter of the Introduction and Programming Preliminaries 3 current text which forms an introduction to the same topic. Other Fortran libraries with which readers may be familiar include HSL, IMSL, LINPACK and EISPACK, the last two being sublibraries dealing with linear algebra (Chapter 2 in this book) and eigenvalue problems (Chapter 4) respectively. LINPACK and EISPACK have been largely superseded by the more modern LAPACK1. It can be seen that the majority of deterministic analysis methods will be dealt with in the following chapters.
The selection is governed by limitations of space and of teaching time in typical courses. Attention is directed specifically towards coverage of probably the most im- portant areas of numerical analysis concerning engineers, namely the solution of ordinary and partial differential equations.