COMPUTATIONAL METHODS MCEN 3030
Popular in Course
Popular in Mechanical Engineering
This 35 page Class Notes was uploaded by Deron Schmitt on Thursday October 29, 2015. The Class Notes belongs to MCEN 3030 at University of Colorado at Boulder taught by Staff in Fall. Since its upload, it has received 10 views. For similar materials see /class/232015/mcen-3030-university-of-colorado-at-boulder in Mechanical Engineering at University of Colorado at Boulder.
Reviews for COMPUTATIONAL METHODS
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/29/15
Lectures 1 and 2 Base 10 numbers decimal Numbers using b 10 as a basis are the usual numbers we are using in our every day life The dj are the digits of that number dN1dN2d2d1d0 Zdjioi 1 j0 where dj 6 0172374576777879 Example 144 210739 c10100 d1101 c12102 j0 4 40 100 144 with dj0f0rj gt2 On computers Things are quite different The smallest building block for coding informa tion on a computer is called the bit in017 they can be arranged in bytes 8 bits and words On a 32 bits machine a word has a 4 bytes length while on a 64 bits machine words are 8 bytes long The coderprogrammer doe not use but can directly bits he uses mostly integers and oating pointsln matlab7 computations are performed in double precision 64 bits Consequently since the word size is nite on a computer the default types in a programming language need some sort of unique representation and some choices need to be made To encode the sign of a number one bit is required Consequently7 if we wish to represent all the possible integers using 4 bytes word we have l71S id j 2 where dj in07 1 The biggest and smallest numbers representable are given respectively by dj 1 and s 0 and dj 1 and s 1 231 7 1 231 7 1 271 mam 11111111111111111111111111111111 214748364710 Notethe riite sum ofa geometric series earl be computed usirig Gn 20 bi bquot171 1271 39 For oating point numbers that can take values after the 77decimal77 point7 a fair representation is the following I7 f Awe Z djbii j1 The nite set generated by the above representation will be denoted M57107 5mm 51mm denorm C R Note that all the parameters of the set are integers except for denorm true or false Which will be explained later7 also7 e E Emmemm 17 emm The sum in 21 djb j is called the mantissa Single precision We have 1 bit for the sign 8 for the exponent and 23 for the mantissa totaling 32 bits m Double precision 1 bit for the sign 11 for the exponent and 52 for the mantissa totaling 64 bits lsl e m l Even with this de nition7 a problem still persists the problem of uriique representation Example Suppose p 6 arid b 10 lets try to represerit 01 01 7 100000 gtlt100 001000 gtlt 102 Which one of the above is the right representation Normalized oating point number EN In the previous example7 the ambiguity can be removed by adding the con straint d1 31 0 The only possible representation again in the previous ex ample is 100000 gtlt 100 The mantissa will cover the interval b 171 7 b p lb7 1 7 e and the exponent is allowed to vary from emm to em Denormalized oating point number FD When e emm then b1 0 is allowed Here7 the mantissa covers the interval b z b 1 7 b p 6 1b 7 6 The two different choice for our nite oating point sets are Fb7p7 Emir emazy true EN b71137 Emir 5min U FDULZL Emir 5min Fb7p7 emiru emam false FNULIL Emir 57mm the de nition for zero 0 is added to the set F and corresponds to e emm with m 0 Most computers nowadays are using the denormalized oating point num ber system This is why the picture about the binary representation in the book p62 g 37 is not completely true There is no gap around zero in the denormalized representation Hidden bit In the normalized oating point number system7 in the binary case b 27 the d1 bit does not need to be stored explicitly since d1 is always 1 The same remark goes for denormalized oating point number system using b 2 b1 0 always This enables computations with a 53 bits mantissa in double precision instead of 52 and7 in single precision7 24 bits can be employed The hidden bit plays an important role for performing rounding operations in M2710emimemahtruefalse Before7 touching the subject of roundo errors lets rst concentrate on how arithmetic operations are happening on a computer Arithmetic operations on a computer 0 Normalize the 2 numbers appearing in the operation 0 Denormalize the smallest 0 Perform the operation 7 7 7 gtlt o Normalize the result Example A Sum 11102 and 12 in base 2 b 2 with p 3 3 digits after the point 8 11102 12 5 111gtlt 24 100 gtlt 2 normalize s 111 gtlt 24 010 gtlt 22 denormalize smallest s 111 gtlt 24 001 gtlt 23 denormalize smallest s 111 gtlt 24 000 gtlt 24 denormalize smallest s 111 000 gtlt 24 perform operation 5 111 gtlt 24 normalize nothing to do here notice that the number 12 once denormalized is replaced by 02 This type of error is called an under ow Example B Sum 111 gtlt 25m and 001 gtlt 25m in base 2 with p 3 s 111 gtlt 25 001 gtlt 25 already normalizedl s 111 001 gtlt 25m perform addition 5 0100 gtlt 25mm r1 normalize s NaN the maximum value of the edponent after normalization was bigger than the the maximum allowed This type of problem is called an over ow error Rounding in l The results of arithmetic operations with operands numbers in the oat ing point number set 19b7p7emmemahtruefalse often needs more than p mantissa places and an exponent outside the range emimemm for their representation as oating points numbers The logical approach to overcome this problem is to roundo results to an exactly representable number in F Example C Sum 123456 and 987654 in 19107 67 79 97true note that the exact result is 12346587654 X 105 s 123456 gtlt 105 987654 gtlt 100 normalized s 1234560 gtlt 105 0000098 gtlt 105 denormalz39ze s 1234560 0000098 gtlt 105 add 5 1234658 gtlt 105 not rounded s 123466 gtlt 105 rounded Now7 the role played by the extra bit in our base 2 representation is clear It serves into doubling the numbers in l to enable rounding Another small example for a calculation in base 2 with rounding Example D Sum 1110 gtlt 23 and 1100 gtlt 20 in base 2 withp 4 s 11100 gtlt 23 11000 gtlt 20 already normalized s 11100 gtlt 23 01100 gtlt 21 denormalize smallest s 11100 gtlt 23 00110 gtlt 22 again 5 11100 gtlt 23 00011 gtlt 23 again avoiding under ow s 11100 00011 gtlt 23 add 5 11111 gtlt 23 no rounding s 10000 s 01000 0 gtlt 23 rounded 0 gtlt 24 normalize Representing mathematically arithmetic in l The operations described so far can be casted in a general framework Let 0 be the rounding operation It performs the following Vd E R there is a unique ox E F Simply7 given a number x E R it is possible to nd 2 numbers in F say 1 and x2 enclosing x Then7 the goal of o is to round z to its closest representation in F A very simple algorithm for doing so consist in taking the midpoint be tween the two numbers in l namely i 1 ln base two this is easily done by considering the extra hidden bit If the aforementioned z is in 17 i then x is rounded to 1 If x is in 7 2 then x is set to 2 A probable case is z i in this case7 a rule stating that we round away from zero would lead to z x2 because we supposed 1 3 2 The last step of this lecture consist in representing the rounding operation in a mathematically manageable way Let us assume that the rounding op eration on a real number induces an error representable by 0 S p 3 eps lt 1 then the rounding operation is represented by ox x1 p The quantity eps is called the machine epsilon homework It follows that for arty exact arithmetic operation Ly E R and rounding operator 0 we have a way to represent the iriecoaet operation 9 z yy1p Pseudo arit hmetic AssociatiVity does not hold DistributiVity does not hold CommutatiVity holds ExampleVerify that in gerieral the associativity for the multiplication does riot hold in l or in other words z 8 y 8 w 31 x 8 y 8 w Start with z 8 y 8 w y z Xy1p1 96in X21p1192 Oh the other harid we have for x y 8 w z y xgtlty1p3 zzgtltygtltz1p31p4 arid iri gerieral 1 p11 p2 31 1 p31p4 To keep in mind 966 of numbers similar in magnitude but having opposite signs7 an error called cancellation occurs This is7 by far7 the most frequent cause of problems in numerical algorithms Example Suppose we are using M107 67 799true Let a 764352 gtlt 101 a 764359 gtlt 101 and b 764286 gtlt 101 Compute the following ecopressions aibaibab andab a i b 660000 gtlt10 3 a i b 730000 gtlt10 3 a b 152863 gtlt102 a b 152864 gtlt102 Notice that here a small perturbation ofa 370 00007 generates an error in the rst digit of the mantissa for the subtraction while it generates only an error on the last digit in the ease of the addition Chapter 2 Roots of Nonlinear Equations The problem of nding roots of nonlinear equations is quite common in engineering and science We shall begin with the discussion of solving a single nonlinear equation 0 When f is linear or even quadratic equation nding the root is trivial For a general nonlinear equation one can use trial and error method by simply guessing the value of z and evaluating whether is zero or not If not the new guess is made and is evaluated again to see whether the new guess provides a better estimate of the root This process of guessing can be automated based on the information that can be obtained from the function Different techniques of organizing such iterations have different names Some are based on functional evaluations and some use the gradient informationi Three of these techniques namely bisection method secant method and Newton s method will be discussed nexti 21 Bisection Method Bisection method is the easiest way to organize convergent iterationsi The bisection method is based on simple subsequent subdivision ideal The method proceeds as follows Start with two initial guesses 11 and 12 evaluate the function and obtain two values n and zz such that fzl lt 0 lt The next guess 13 11 zz2 is simply a point that divides the interval 1112 in half Evaluate the function at 13 For the next iteration use 11 and 13 if fzg gt 0 or 13 and 12 if fzg lt 0 Continue this subdivision process until the distance between two points becomes less than a p39m39o39m39 de ned distance 239e lri1 7 g 6 The main drawback of the bisection method is slow convergence and dif culty generalizing it to multiple dimensions 22 Secant Method The secant method is an improvement of the bisection method The improvement is achieved by not only using the information about the signs of the guesses but also using the actual values of the function The algorithm proceeds as follows Suppose that using two initial guesses 11 and 12 we evaluate the function and obtain two values yl n and yz zz such that fzl lt 0 lt With the secant method we form the straight line between the points zlyl and zzy2i The equation for this line is yy1mzizl where iyiiyl 12711 m is the reciprocal of the slope of the line The next guess is the value for 13 at which the above straight line approximation to the function predicts y 0 That is the intersection of the horizontal line y 0 CHAPTER 2 ROOTS OF NONLINEAR EQUATIONS Flguxe 21 Illmtzatlon of secant method wtth the sttatght ltne Whlch ytelds 0y1m23711 In genezal the successlve ltexates axe obtamed fzom the fozmula tt 0 yset Ms71 12 23 ts the ttetatton tndex and Wter 7 y 2172 ate the teotbtooals of the s ope of the successwe sttatght ltnes seoants ltetattons ate contlnued unttl N21 ts sumotently close to 0 23 Newton7s Method A dteletent btooedute of obtatntng new values of 21 ts to teblaoe the seoant method by Newton s method fox a t of the equatton tatne ms at 2 21 and so on Thls yte ds the followlng equatton o y if 21 E m 2 7 The next guess 21 ts obtatned by solytng Eq 21 f0 2 when y o Whlch tesults tn 21 de ned by 2 2 Thls ttetatto method genetally oonyetges quadtattoally As an tllusttatton of the method agatn let s oonst he sam oblem befoze N o s method teoutte only one tntttal ess If 2 ts a guess m o a bettet guess 22 ts usually ob d by extxapolatlon to the 8x15 f t too of the tangent to y 2 1 y 21st 12 z 2512 5512 Obvlously211 21 only If mag 0 ln geneta thts wtll not ooout f0 any mte tquot tnstead we ttetate unttl law 7 21 g 5 f0 some sumotently small 5 Note that tn the case when analyttoal dettyattye dfdz ts n t known tt can be evaluated numettoally by eyaluattng functlon at addtttonal looatton close to 21 and use appzoxlmate yalue fo the dettyattye if N flt2t 5 ms Eta e 5 One may enoountet dtmoulty tn obtatntrg a conve ofthezoot Tn tht l and If thts dtstanoe tnoteases deotease tt aaootdtng to some a pmom de ned tule tged solutton If f2 ohanges outyatute at the looatton t t t m as NEWTON S METHOD 1N MULTIPLE DIMENSIONS ii ii Figuxe 2 2 lllustzation ofNewton s method 24 Newton s Method in Multiple Dimensions A genexahzatlon of the Newton s method to multiple dimensions is xelatwel e coozdmate is zeplaced by the Vectoz x e zn e is the dimensionality space the t eth guess ex is zeplaced by a Vectoz x en em the function f is zeplaced by a Vectoz function r f1 fn A L 39 39 39 L Let us i11ustiete the pxoblem of nding nonlineaz equation ystzalghtfozwazd the of e 9 r r 37 M ioots of nonlineaz multidimensional equations by soiying a single fx 0 As in one dimensional case multidimensional Newton s method ieouiies only one initial guess X1 A d ttei guess X2 is obtained by extzapolatlon of the 1ine tangent to y fx at x X en so on Note that the 1ine is de ned in 2 dimensional space e1 zmyj yn This yields the folloWing equation i e y y We Vftm 39 X xt wheie Vf is the Jeoobien J de ned by i 23 BIA BIA Vi F JVf The next guess q1 is obtained system of hneei equations ow de J m X Xe Solution of the system hneei equationspioVides the next guess x1 i T e Itelatlon ontinue unti1 ersi 7 e foi 1y sm o t aco J can e by so1Ving Eq 2 3 foi x when y 0 which iesu1ts in the o11oWing ned by ewe 12 2 4 s h c XeHz g sonne su oient all e Anal gous1y to one dimensional case he J blan b uated numeucally by evaluating function at e additional location close to on and use appxoxlmate Value oi the deiiyetiye wg who ems z rem 2m 5 sin A 39 di ouity in bt 39 39 im ohenges cuzvatuze at the location of the loot ln thi s ease one ieouiies to put some nrdimenslonal 1inniteis on the distance between two successive guesses and if this distance inoiesses deoiesse it accozdmg to some a mom de ned iu1e Matlab Basics M C E N 3 030 MCENRUBU Labs UKquotquot7 n Il Outline Data Structures Special Operator Colon z 0 Matrix Generation Arithmetic Operators Linear Algebra N ICENRUBU Labs 03 7 D llX Data Structures MCENRURU Labs 98 47 D 3l8 Matrix Maj or element of Matlab is the matrix Scalar or number 1 x 1matrix Vector row or column 1 x n or m x 1 matrices Matrix n x m matrix which is the general case Examples Scalar a 10 sizea Vector bl 2 3 sizeb c 1 2 3 sizec Matrix d1 2 3 45 6 sizedsized rowsized l 1 colsized12 MCENRURU Labs 98117 n wll Matrix Access Matrix elements using subscripts Ai j element in row i and column j of A A12A22 summation A22 2 modi cation Al nm rst n elements of mth column of A A m colon by itself refers to all elements Anl m rst m elements of nth row of A An row by itself refers to all elements MCENRURH Labs 08 47 n 5l8 Special Operator Colon z RICENHL Labs AKquotquot7 n olx Colon z One of the most useful operators in MATLAB It can create vectors subscript arrays and specify for iterations 1 Create vectors jzk is the same as jj1k with default increment 1 EX 110 jzk is empty ifj gt k EX 101 jzizk is the same as jjij2i k EX 100 750 TVICENRURH Labs WWW n 7IX Colon z 2 Subscript arrays the use of the colon to pick out selected rows columns and elements of matrices Aj is the jth column of A Ai is the ith row of A A is the same as A Anmj is nth to mth elements of the jth column 3 Specify for iterations for m 110 for loop am lm assign values end MCENRUBH Labs WWW n XIX Matrix Generation RICENHL Labs AKquotquot7 n Lle Matrix Generation Question How to generate the following matrix OOOOUH OONNJQQ 00000000 COrPNOO 0103000 Answer Many approaches MCENRUSU Labs 9807 I lolx Matrix Generation 1 Use an explicit list of elements Al900062800073700084600095 separate elements of a row with blanks or commas use a semicolon to indicate the end of each row surround entire list of elements with square brackets Note Impractical for a large matrix MCENRUSU Labs MK7 I llIX Matrix Generation 2 Use builtin functions ones zeros eye diag etc generating matrix A zeros5 modifying matrix A1 5 151 5 A251469 A1425916 a15 b69 c916 Adiaga0 diagb1 diagc 1 For the Project 3 it is easy to implement the matrix using the builtin functions 3 Others MCENRUSU Labs HR7 I lilx Arithmetic Operators MCENM Labs HR7 I lSl Arithmetic Operators MATLAB has two types of arithmetic operations Matrix arithmetic operations are de ned by the rules of linear algebra Array arithmetic operations are carried out element by element The period character is for the elementbyelement operations Note that the addition and subtraction are as same as the elementby element operations and from the above de nitions of operators MCENRUSU Labs 9807 I HaIX Arithmetic Operators Matrix arithmetic Operators Addition Subtraction Multiplication Division 0 A Power Complex conjugate transpose Specify evaluation order MCENBUSU Labs 9307 I lSlx Arithmetic Operators Array arithmetic Operators Addition Subtraction Elementbyelement Multiplication Elementbyelement Division 0 Element by element Power Elementbyelement Complex conjugate transpose Specify evaluation order MCENRUSU Labs 0807 I lotl8 Linear Algebra MCENMBU Labs 9307 I l7IX Linear Algebra Bmwm m mmmmmmAymz xlA b use Gauss elimination LU 1 x2 A use lu decomposition u Lb C I O P 1uA 6 use lu decomposition and p is the permutation matrix ULPb L H o C x3 00 x4invAb too expensive for matrix inverse MCENRUSU Labs 9307 I XiIX Lectures 3 to 4 Systems of linear equations We are interested in solving linear equations where a by hand elimination or application of Cramer7s rule is senseless Model problem Let A E Rmxm and b E E Find z E R such that Ax b is satis ed Although if this seems irrelevant at rst7 matrix problems can be found everywhere in the applied sciences o Discretization of PDEs Finite difference discretizations Finite element discretizations Finite volumes Flows around airplanes Elasticity beam problems7 eigenmodes of a bridge 7 Heat transfer Semi Conductors design the next generation of chips using the current one 7 Geophysical ows 0 Statistics eg least squares 0 Google the search engine is based on an eigenvalue problem In the above model problem7 it is pretty obvious that for z to exist the inverse of the matrix must exist We will say that a matrix A is regular if one the following criteria is satis ed o A 1 exist 0 rankA m o The rows are linearly independent 0 The columns are linearly independent o detA 73 0 0 Ad b has a unique solution Vb E R o z 0 is the only solution to Ax 0 o 0 is not an eigenvalue of A o A is symmetric positive de nite The last criteria involves the notion of positive de niteness SPD The latter is de ned as follows If A AT and pTAp gt 0 for all z E R then A is symmetric positive de nite or SPD This notion is quite useful in practice Eigenvalues can be found by solving the polynomial factorization problem detA 7 AI 0 For an m gtlt m matrip we get in eigenvalues For each eigenvalue and eigenvector can be found using Av Av Lets denote by M the matrip having as its columns the eigenvectors ofA then AM MA If the inverse ofM exists then M lAM A and thus A is similar to the diagonal matria A diag12m The inverse of the A matria is A 1 MA lM l Obviously if one of the eigenvalues is zero the inverse ofA cannot be found More generally any square matricc can always be but in its Jordan form which is simply an upper triangular matrip Else said we can nd matrices X and X 1 such that X lAX is an upper triangular matrip The values on the diagonal are the eigenvalues of the matrip More de nitions Norms are very useful in mathematics science and engineering They are tools used to prove theorems and more practically to compute distances and measure errors For a vector v xyzT in R3 it is possible to de ne the distance between the head of the vector and the origin 0 000T like dv 7 0 Mp2 y2 22 This distance function is in fact the Euclidian norm De nition of a norm Any function Q a R that satis es 1 HAH 2 0 2 Oif and only ifA0 3 Hamill lozHlAH for 04 E R 4 HABH HAHHBH forAB 69 For vectors7 Q is replaced by R and for matrices7 the space 9 is changed to Rmxm De nition of p norms for vectors Let 1 E R then for p 6 100 we have m A Hva Z W i1 while for p 00 the de nition is Hvlloo max vltlvztm lvml Notice that for p 2 and m 3 we get the usual Euclidian norm De nition of the p norm of a matrix For A E Rmxm llA llp HAM M p E willfg w HAM lp 3 Another norm useful in practice is called the Frobenius norm and is de ned as the sum of the squares of the elements of the matrix More precisely m m 22 law i1 j1 HAHF Some important facts about norms o HABH HAHHBH forAyBeRmm o S for A ERmxm and z ER The condition number of a linear system It was shown in class that even if the problem of the intersection of two lines has a unique solution that under the effect of a small perturbation of their slopes the resulting perturbed solution could be far away from from the unperturbed one In these notes the notion is generalized to matrices and a number called the condition number will be de ned First things rst let us de ne the perturbation of a regular matrix A E Rmxm by AA also 6 Rmxm and its relative error by pA The same can be done for the right hand side b7 0 of the model problem Pb HAbHHbH We are now in a position to quantify the sensitivity of z with respect to perturbation of the matrix A For t E R we consider the following problem A tAAxt b tAb with 0 z A lb 5 Notice that for t 0 we have the original model problem Am b Differentiating with respect to t the last equation one obtains d dx dx ampA tAAxt Ag AAth tAAE Ab A mag Mm Ab 6 if the above equation is evaluated at t 0 we obtain a de nition for the derivative of z around It 0 lt0 E A 1Ab 7 AA0 The Taylor expansion of t around It 0 is dx 2 9W 960 Elam 0it 7 de ning the relative error on the solution 7 xl and using the de nition for the derivative the Taylor expan sion can be rewritten like 1 7 ll ll ll cll A lAb A lAA 0 H H H zlt WHOM ll ll ll ll llA lllllAbll llA lllllAAllWONl ll ll WM 7 Ab HA WHHlMHWHO Z lx llAbH HAAH lA llHlAK7 HAHHzH HAM HAbH MM 2 77t0t gtk le HAM l llA lllllAllpbAltl0t2 t2 S llA lllllAlMpHJAW 0 S S W 002 W 0052 S llA lllllAlM OI HMH llwll In the above calculation the fact that g and thus W S m was employed to go from gtk to gtk Also division by is OK since l 7E 0 and A is regular When A and b are undergoing a rst order perturbation the relative error on the solution x is bounded by a rst order relation on the relative errors made on A and b ampli ed by the number llA llHlAll S llA lllllAlKJb 9AM 002 De nition of the condition number of a matrix If A E Rmxm and regular then its condition number is EUREKA E 0AEllA 1llllAll 8 Example Lets de ne Ab Ab and the solution to Ax b as 300 150 100 020 001 714 A 150100075 b 100 Ab 7001 x 16 100 075 060 100 001 20 Solving for i in Ai b Ab leads to the solution i 7 z Ax 0257 71367130 HAbHOOHblloo 001 and 068 Thus there is a factor of 68 in increase in the error Using the norm H00 gives the condition number 748 for A indicating that even more unfavorable cases can exist for this particular A On a computer7 if HAeps 2 1 then the numerical solution of the linear system may not even be correct in the rst digit Direct methods for solving Ax b
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'