Jeffrey R. Chasnov - Department of Mathematics, HKUST

1 numerical MethodsJeffrey R. ChasnovAdapted for: numerical Methods for EngineersClick to view a promotional videoThe Hong Kong University of Science and TechnologyDepartment of MathematicsClear Water Bay, KowloonHong KongCopyrightc 2012 by Jeffrey Robert ChasnovThis work is licensed under the creative Commons Attribution Hong Kong License. Toview a copy of this license, visit or senda letter to creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, follows were my lecture notes for Math 3311: introduction to numerical Meth-ods, taught at the Hong Kong University of Science and Technology. Math 3311,with two lecture hours per week, was primarily for non-mathematics majors andwas required by several engineering also have some free online courses on Coursera. A lot of time and effort has goneinto their production, and the video lectures for these courses are of high can click on the links below to explore these you want to learn differential equations, have a look atDifferential Equations for EngineersIf your interests are matrices and elementary linear algebra, tryMatrix Algebra for EngineersIf you want to learn vector calculus (also known as multivariable calculus, or calcu-lus three), you can sign up forVector Calculus for EngineersAnd if your interest is numerical methods, have a go atNumerical Methods for EngineersJeffreyR.

2 ChasnovHong KongFebruary 2021iiiContents1 IEEE .. with a decimal or binary point .. of binary numbers .. numbers .. unsigned integers as hex numbers .. single precision format: .. numbers .. of computer numbers .. numbers .. smallest positive integer that is not exact in single Machine epsilon .. IEEE double precision format .. Roundoff error example ..52 Root Method .. s Method .. Method .. 2= using Newton s Method .. of fractals using Newton s Method .. of convergence .. s Method .. Method .. 103 Systems of Elimination .. pivoting .. counts .. of nonlinear equations .. 204 Least-squares a straight line .. to a linear combination of functions .. 245 interpolation .. polynomial .. polynomial .. polynomial .. linear interpolation.

3 Spline interpolation .. interpolation .. 33vCONTENTS6 formulas .. rule .. rule .. s rule .. rules .. rule .. s rule .. versus global error .. integration .. 397 Ordinary differential of analytical solutions .. value problem .. value problems .. problem .. methods: initial value problem .. method .. Euler method .. Runge-Kutta methods .. Runge-Kutta methods .. Runge-Kutta Methods .. of differential equations .. methods: boundary value problem .. difference method .. method .. methods: eigenvalue problem .. difference method .. method .. 53viCONTENTSC hapter 1 IEEE DefinitionsBit=0 or 1 Byte=8 bitsWord=Reals:4 bytes (single precision)8 bytes (double precision)=Integers:1, 2, 4, or 8 byte signed1, 2, 4, or 8 byte Numbers with a decimal or binary point Decimal:10310210110010 110 210 310 4 Binary:232221202 12 22 32 Examples of binary Hex numbers{0, 1, 2, 3.}

4 , 9, 10, 11, 12, 13, 14, 15}={0, 1, 2, , a,b,c,d,e,f} 4-bit unsigned integers as hex IEEE SINGLE PRECISION IEEE single precision format:s 0e 1 2 3 4 5 6 7 8f 9 31# =( 1)s 2e 127 = signe = biased exponentp=e-127 = = significand (use binary point) Special numbersSmallest exponent:e = 0000 0000, represents denormal numbers ( )Largest exponent:e = 1111 1111, represents , if f = 0e = 1111 1111, represents NaN, if f6=0 Number Range:e = 1111 1111 = 28- 1 = 255reservede = 0000 0000 = 0reservedso,p = e - 127 is1 - 127 p 254-127-126 p 127 Smallest positive normal number= 0000 0000 2 126' 10 38bin: 0000 0000 1000 0000 0000 0000 0000 0000hex: 00800000 MATLAB: realmin( single )Largest positive number= 1111 1111 2127=(1+ (1 2 23)) 2127'2128' 1038bin: 0111 1111 0111 1111 1111 1111 1111 1111hex: 7f7fffffMATLAB: realmax( single )Zerobin: 0000 0000 0000 0000 0000 0000 0000 0000hex: 00000000 Subnormal numbersAllow (in software)Smallest positive number = 0000 0001 2 126= 2 23 2 126' 10 452 CHAPTER 1.

5 IEEE EXAMPLES OF COMPUTER Examples of computer numbersWhat is , & 1/2 in hex ? ( 1)0 2(127 127) ,s=0,e=0111 1111,f=0000 0000 0000 0000 0000 000bin: 0011 1111 1000 0000 0000 0000 0000 0000hex: 3f80 ( 1)0 2(128 127) ,s=0,e=1000 0000,f=0000 0000 0000 0000 0000 000bin: 0100 00000 1000 0000 0000 0000 0000 0000hex: 4000 00001/2= ( 1)0 2(126 127) ,s=0,e=0111 1110,f=0000 0000 0000 0000 0000 000bin: 0011 1111 0000 0000 0000 0000 0000 0000hex: 3f00 Inexact numbersExample:13= ( 1)0 14 (1+13),so thatp=e 127= 2 ande=125=128 3, or in binary,e=0111 1101. How isf=1/3 represented in binary? To compute binary number, multiply successivelyby 2 as .. that 1/3 exactly in binary is .. With only 23 bits to representf, thenumber is inexact and we havef=01010101010101010101011,where we have rounded to the nearest binary number (here, rounded up).

6 Themachine number 1/3 is then represented as00111110101010101010101010101011or in 1. IEEE MACHINE Find smallest positive integer that is not exact in single pre-cisionLetNbe the smallest positive integer that is not exact. Now, I claim thatN 2=223 .. 1,andN 1=224 .. integerNwould then require a one-bit in the 2 24position, which is not avail-able. Therefore, the smallest positive integer that is not exact is 224+1=16 777 MATLAB, single(224) has the same value as single(224+1). Since single(224+1)is exactly halfway between the two consecutive machine numbers 224and 224+2,MATLAB rounds to the number with a final zero-bit in f, which is Machine epsilonMachine epsilon(emach)is the distance between 1 and the next largest number. If0 <emach/2, then 1+ =1 in computer math. Also sincex+y=x(1+y/x),if 0 y/x<emach/2, thenx+y=xin computer number 1 in the IEEE format is written as1=20.

7 0,with 23 0 s following the binary point. The number just larger than 1 has a 1 in the23rd position after the decimal point. Therefore,emach=2 23 10 is the distance between 1 and the number just smaller than 1? Here, thenumber just smaller than one can be written as2 1 .. 1=2 1(1+ (1 2 23)) =1 2 24 Therefore, this distance is 2 24= spacing between numbers is uniform between powers of 2, with logarithmicspacing of the powers of 2. That is, the spacing of numbers between 1 and 2 is 2 23,between 2 and 4 is 2 22, between 4 and 8 is 2 21, etc. This spacing changes fordenormal numbers, where the spacing is uniform all the way down to the machine number just greater than5A rough estimate would be 5(1+emach) =5+5emach, but this is not exact. Theexact answer can be found by writing5=22(1+14),so that the next largest number is22(1+14+2 23) =5+2 21=5+ 1.

8 IEEE IEEE DOUBLE PRECISION IEEE double precision formatMost computations take place in double precision, where round-off error is re-duced, and all of the above calculations in single precision can be repeated fordouble precision. The format iss 0e 1 2 3 4 5 6 7 8 9 10 11f 12 63# =( 1)s 2e 1023 = signe = biased exponentp=e-1023 = = significand (use binary point) Roundoff error exampleConsider solving the quadratic equationx2+2bx 1=0,wherebis a parameter. The quadratic formula yields the two solutionsx = b b2+ the solution withb>0 andx>0 (thex+solution) given byx= b+ b2+1.( )Asb ,x= b+ b2+1= b+b 1+1/b2=b( 1+1/b2 1) b(1+12b2 1)= in double precision, realmin 10 308and we would likexto be accurateto this value before it goes to 0 via denormal numbers. Therefore,xshould becomputed accurately tob 1/(2 realmin) 2 10307.

9 What happens if wecompute ( ) directly? Thenx=0 whenb2+1=b2, or 1+1/b2=1. That is1/b2=emach/2, orb= 2/ emach 1. IEEE ROUNDOFF ERROR EXAMPLEFor a subroutine written to compute the solution of a quadratic for a generaluser, this is not good enough. The way for a software designer to solve this problemis to compute the solution forxasx=1b(1+ 1+1/b2).In this form, if 1+1/b2=1, thenx=1/2bwhich is the correct asymptotic 1. IEEE ARITHMETICC hapter 2 Root FindingSolvef(x) =0 forx, when an explicit analytical solution is Bisection MethodThe bisection method is the easiest to numerically implement and almost alwaysworks. The main disadvantage is that convergence is slow. If the bisection methodresults in a computer program that runs too slow, then other faster methods maybe chosen; otherwise it is a good choice of want to construct a sequencex0,x1,x2.

10 That converges to the rootx=rthat solvesf(x) =0. We choosex0andx1such thatx0<r<x1. We say thatx0andx1bracket the root. Withf(r) =0, we wantf(x0)andf(x1)to be of oppositesign, so thatf(x0)f(x1)<0. We then assignx2to be the midpoint ofx0andx1,that isx2= (x0+x1)/2, orx2=x0+x1 sign off(x2)can then be determined. The value ofx3is then chosen as eitherthe midpoint ofx0andx2or as the midpoint ofx2andx1, depending on whetherx0andx2bracket the root, orx2andx1bracket the root. The root, therefore,stays bracketed at all times. The algorithm proceeds in this fashion and is typicallystopped when the increment to the left side of the bracket (above, given by(x1 x0)/2) is smaller than some required Newton s MethodThis is the fastest method, but requires analytical computation of the derivative off(x). Also, the method may not always converge to the desired can derive Newton s Method graphically, or by a Taylor series.