SLIDEREAD

Numerical Methods for EngineersLecture Notes forJeffrey R. Chasnov

The Hong Kong University of Science and TechnologyDepartment of MathematicsClear Water Bay, KowloonHong Kongby Jeffrey Robert ChasnovThis work is licensed under the Creative Commons Attribution 3.0 Hong Kong License. To view a copy of thislicense, 171 SecondStreet, Suite 300, San Francisco, California, 94105, USA.

PrefaceView the promotional video on YouTubeThese are the lecture notes for my upcoming Coursera course , Numerical Methods for Engineers (forrelease in January). Before students take this course, they should have some basic knowledgeof single-variable calculus, vector calculus, differential equations and matrix algebra. Students shouldalso be familiar with at least one programming language. In this course, however, I will exclusivelyuse Matlab. I teach the basics of Matlab in the first week, but this may be too short an introductionfor students of limited programming ability and they may need to supplement their programminglessons elsewhere.I have divided these notes into chapters called Lectures, with each Lecture corresponding to avideo on Coursera. I have also uploaded all my Coursera videos to YouTube, and links are placed atthe top of each Lecture.There are problems at the end of each lecture, some that require analytical solutions and othersthat require Matlab programs. Solutions to the analytical questions and Learner Templates for theMatlab programs can be found in the Appendix.On the Coursera platform, at the end of each week there is also both an assessed multiple-choicequiz and a Matlab project. Details of the Matlab projects and their Learner Templates can also befound in these lecture notes.Jeffrey R. ChasnovHong KongNov iii

ContentsIScientific Computing11Binary numbers52Double precision73Matlab as a calculator94Scripts and functions115Vectors136Line plots177Matrices218Logicals259Conditionals2710 Loops2911 Project I: Logistic map (Part A)3112 Project I: Logistic map (Part B)33IIRoot Finding3513 Bisection method3914 Newton’s method4115 Secant method4316 Order of convergence4517 Convergence of Newton’s method4718 Fractals from Newton’s method49v

CONTENTSvi19 Coding the Newton fractal5120 Root finding in Matlab5521 Project II: Feigenbaum delta (Part A)5722 Project II: Feigenbaum delta (Part B)5923 Project II: Feigenbaum delta (Part C)61III63Matrix Algebra24 Gaussian elimination without pivoting6725 Gaussian elimination with partial pivoting6926 LU decomposition with partial pivoting7127 Operation counts7528 Operation counts for Gaussian elimination7729 Operation counts for forward and backward substitution7930 Eigenvalue power method8131 Eigenvalue power method (example)8332 Matrix algebra in Matlab8533 Systems of nonlinear equations8934 Systems of nonlinear equations (example)9135 Project III: Fractals from the Lorenz equations93IV95Quadrature and Interpolation36 Midpoint rule9937 Trapezoidal rule10138 Simpson’s rule10339 Composite quadrature rules10540 Gaussian quadrature10741 Adaptive quadrature10942 Quadrature in Matlab11143 Interpolation113

CONTENTSvii44 Cubic spline interpolation (Part A)11545 Cubic spline interpolation (Part B)11746 Interpolation in Matlab12147 Project IV: Bessel functions and their zeros123VOrdinary Differential Equations12548 Euler method12949 Modified Euler method13150 Runge-Kutta methods13351 Second-order Runge-Kutta methods13552 Higher-order Runge-Kutta methods13753 Higher-order odes and systems13954 Adaptive Runge-Kutta methods14155 Integrating odes in Matlab (Part A)14356 Integrating odes in Matlab (Part B)14557 Shooting method for boundary value problems14958 Project V: Two-body problem (Part A)15159 Project V: Two-body problem (Part B)153VIPartial Differential Equations60 Boundary and initial value problemsPractice quiz: Classify partial differential equations15515916161 Central difference approximation16362 Discrete Laplace equation16563 Natural ordering16764 Matrix formulation16965 Matlab solution of the Laplace equation (direct method)17166 Jacobi, Gauss-Seidel and SOR methods17567 Red-black ordering177

viiiCONTENTS68 Matlab solution of the Laplace equation (iterative method)17969 Explicit methods for solving the diffusion equation18170 Von Neumann stability analysis18371 Implicit methods for solving the diffusion equation18572 Crank-Nicolson method for the diffusion equation18773 Matlab solution of the diffusion equation18974 Project VI: Two-dimensional diffusion equation193AppendixA Problem solutions and Matlab Learner Templates195197

Week IScientific Computing1

3In this week’s lectures, we learn how to program using Matlab. We learn how real numbers arerepresented in double precision and how to do basic arithmetic with Matlab. We learn how to usescripts and functions, how to represent vectors and matrices, how to draw line plots, how to use logicalvariables, conditional statements, for loops and while loops. Your programming project will be towrite a Matlab code to compute the bifurcation diagram for the logistic map.

4

Lecture 1Binary numbersView this lecture on YouTubeWe do our arithmetic using decimals, which is a base-ten positional number system. For example,the meaning of the usual decimal notation is illustrated by524.503 5 102 2 101 4 100 5 10 1 0 10 2 3 10 3 .Each position in a decimal number corresponds to a power of 10. A digit in this number system isdefined as any of the numerals from 0 to 9, while digit is also the English word meaning a finger or atoe. Decimals probably arose from counting using ten fingers.Computers count using two states, and this has led to the use of binary numbers, which is a basetwo positional number system. Here, instead of digits, we use bits, defined as either a 0 or a 1. Themeaning of a binary number is illustrated by101.011 1 22 0 21 1 20 0 2 1 1 2 2 1 2 3 .In both decimal and binary, certain fractions can be represented only by infinitely repeated numerals. In decimal, only fully reduced fractions whose denominators are a product of powers of 2 and 5do not repeat. For example, we have1 0.5,2but we also have1 0.3,31 0.25,41 0.16,61 0.2,51 0.125;81 0.142857,71 0.1,9where we use the bar notation to indicate repeating numerals.In binary, only fully reduced fractions whose denominators are a product of powers of 2 do notrepeat. For example, using our more familiar digits to represent the fraction, the corresponding binarynumbers are given by1 0.1,2and1 0.01,31 0.0011,51 0.01,41 0.001,81 0.0010,61 0.001,71 0.000111.9Binary numbers with infinite repeats can not be represented exactly using a finite number of bits andare a potential source of round-off errors.5

6LECTURE 1. BINARY NUMBERSProblems for Lecture 11. Using binary, round the fractions 1/3, 1/5, 1/6, 1/7 and 1/9 to six places after the binary point.Which numbers round down and which numbers round up?Solutions to the Problems

Lecture 2Double precisionView this lecture on YouTubeEight bits make a byte. Most numerical computation is done in double precision, with numbersstored using eight bytes. The format of a double precision number issfez} { z} {z} { · · · · · · · · 01 2 3 4 5 6 7 8 9 10 11 1263# ( 1)s 2e 1023 1.fwhere s is the sign bit, e is the biased exponent, and 1.f is the significand. Here, e is a whole numberwritten in binary, 1023 is in decimal, and f is the fractional part of the binary number following abinary point.The 64 bits of a double precision number are distributed so that the sign uses one bit, the exponentuses 11 bits (in decimal, 0 e), and the significand uses 52 bits. The distribution of bitsreconciles two conflicting needs: that the numbers should range from the very large to the very small,and that numbers should be closely spaced.Both the largest and smallest exponents are reserved. When e is all ones (e in decimal),f 0 is used to represent infinity (written in Matlab as Inf) and f 6 0 is used to represent ‘nota number’ (written in Matlab as NaN). NaN typically results from a 0/0, / or operation.When e is all zeros, the double precision representation changes from 1. f to 0. f , allowing these denormal numbers to gracefully underflow towards zero. The largest positive double precision number isrealmax 1.7977e 308, and the smallest positive normal number is realmin 2.2251e-308.Another important number is called machine epsilon (written in Matlab as eps) and is defined asthe distance between one and the next largest machine number. If 0 δ eps/2, then 1 δ 1 incomputer arithmetic. And sincex y x (1 y/x ),when y/x eps/2, then x y x. In double precision, machine epsilon is equal to eps 2.2204e-16.Note that the spacing between numbers is uniform between powers of 2, but changes by a factor oftwo with each additional power of two. For example, the spacing of numbers between one and two iseps, between two and four is 2*eps, between four and eight is 4*eps, and so on.We have already learned that not all rational numbers can be represented exactly in double precision. For example, the numbers 1/3 and 1/5 are inexact. In most computations this doesn’t matter.But one needs to be aware that a calculation using inexact numbers that should result in an integer(such as zero) may not because of roundoff errors.7

8LECTURE 2. DOUBLE PRECISIONProblems for Lecture 21. Determine the double precision formats of the numbers 1, 1/2 and 1/3.2. Using the format of a double precision number, determine the largest machine number realmax.3. Using the format of a double precision number, determine the smallest positive normal machinenumber realmin.4. Using the format of a double precision number, determine the distance between one and the nextlargest machine number.Solutions to the Problems

Lecture 3Matlab as a calculatorView this lecture on YouTubeThe basic arithmetic operators are the usual , -, *, /, . The single quote ' transposes a matrix. Theorder of operations follow the standard mathematical rules, and parenthesis can be used. Commonmathematical functions are sin, cos, tan, exp, log. The base-10 logarithm is log10. The square-rootis sqrt and the absolute value is abs. There is a constant pi, though one will need to define e exp(1).Complex numbers use the imaginary unit i or j. Some examples are log(exp(1))ans 1 (sqrt(5) 1)/2ans 1.6180 Matlab also has two special numbers called Inf (for infinity) and NaN (for not-a-number). Thesespecial numbers commonly occur when your program either has a bug or some other run-time error.Simple occurrences of Inf and NaN are 1/0ans Inf 0/0ans NaN A semicolon at the end of a command suppresses output. Two commands can also be placed on oneline, using either a comma or a semicolon.: x 0; y 1; x,yx 0y 1 9

LECTURE 3. MATLAB AS A CALCULATOR10Problems for Lecture 31. Use Matlab to compute the following expressions. 5 1a)2 10 10 5 15 1 22 b)5c)25 125 d)11 52 1e) e3f ) ln (e3 )g) sin (π/6)h) cos (π/6)i) sin2 (π/6) cos2 (π/6)Solutions to the Problems

Lecture 4Scripts and functionsView this lecture on YouTubeFor a quick solution to a problem, one can type directly into the Command Window. But whenyou want to type a longer sequence of commands, or write a program that will be run several times,you are better off writing either a script or a function. These saved programs are called m-files becauseof their .m extension.A script is essentially a sequence of commands that can be typed directly into the CommandWindow, but are saved into a file. When working with a file, editing and debugging become easier,and commands can be quickly added or changed.You can easily start a new script in the HOME tab, and save it with a name and in a directoryof your choice. When developing code, many Matlab programmers like to start a script with thecommandsclear all; close all; clc;These three commands clear the work space, close all figures, and clear the command window. It is aquick way to reset Matlab to prevent previously work from interfering with your new script.It is often useful for a script to print data during execution. To print to the command window or afile, you will need to learn how to use the Matlab function fprintf.m.A function is different than a script because the variables defined in a function are removed fromthe Workspace after the function ends. You can also create a new function in the HOME tab. Whencreating a new function, my Editor helpfully shows me the following:function [outputArg1,outputArg2] untitled(inputArg1,inputArg2)%UNTITLED Summary of this function goes here%Detailed explanation goes hereoutputArg1 inputArg1;outputArg2 inputArg2;endA function can optionally have both input arguments and output arguments. The title of the functionshould match the name of the file it is saved under. Functions may also be written inside of scripts.Matlab calls these local functions.In this course, we will mainly write scripts. But we will also make use of many Matlab providedfunctions. After this course, you may sometimes want to download functions from the Matlab FileExchange. Maybe you will even write a function yourself for others to use.11

LECTURE 4. SCRIPTS AND FUNCTIONS12Problems for Lecture 41. The first few Fibonacci numbers Fn are given by 1, 1, 2, 3, 5, 8, 13, 21, 34, . . . , where starting from 2,each number is the sum of the preceding two numbers. Write a functionF Fibonacci(n)that returns the nth Fibonacci number. You should use Binet’s formula, given byFn where Φ Φn ( φ)n ,55 1,2 φ 5 1.2Because of round-off errors, your function should return a number rounded to the nearest integerusing the built-in Matlab function round.m.Solutions to the Problems

Lecture 5VectorsView this lecture on YouTubeMatlab