### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# NUMERICAL ANLYS I MATH 464

UW

GPA 3.76

### View Full Document

## 40

## 0

## Popular in Course

## Popular in Mathematics (M)

This 52 page Class Notes was uploaded by Addison Beer on Wednesday September 9, 2015. The Class Notes belongs to MATH 464 at University of Washington taught by Staff in Fall. Since its upload, it has received 40 views. For similar materials see /class/192112/math-464-university-of-washington in Mathematics (M) at University of Washington.

## Reviews for NUMERICAL ANLYS I

### 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: 09/09/15

Numerical Methods Design Analysis and Computer Implementation of Algorithms Anne Greenbaum l Timothy P Chartier September 277 2007 1These are course notes for Math 46456 given Autumn 2006 through Spring 2007 at the University of Washington by Anne Greenbaumi Professor Greenbaum and Professor Chartier retain the copyright to these notes They do not give anyone permission to copy computer les related to themi Send email to greenbau mathiwashingtoniedui Contents H 0 09 Basic Operations with MATLAB 6 11 Launching MATLAB 7 12 Vectors 7 13 Getting help 10 14 Matrices 11 15 Creating and running m les 12 16 Comments 12 17 Plotting 13 18 Creating your own functions 15 19 Printing 15 110 More loops and conditionals 17 111 Clearing variables 18 112 Logging your session 18 113 More advanced commands 19 Solution of a Single Nonlinear Equation in One Unknown 27 21 Bisection 29 22 Taylor s Theorem 33 23 Newton s Method 35 24 Quasi Newton Methods 41 241 Avoiding Derivatives 41 242 Constant Slope Method 42 243 Secant Method 42 25 Analysis of Fixed Point Methods 45 26 Fractals Julia Sets and Mandelbrot Setsquot 50 Floating Point Arithmetic 60 Costly Disasters Caused by Rounding Errors 60 32 Binary Representation and Base 2 Arithmetic 63 33 Floating Point Representation 64 34 IEEE Floating Point Arithmetic 66 35 Rounding 68 u U c 36 Correctly Rounded Floating Point Operations 69 37 Exceptions 70 Conditioning of Problems Stability of Algorithms 74 41 Conditioning of Problems 74 42 Stability of Algorithms 75 Solving Linear Systems and Least Squares Problems 80 51 Review of Matrix Algebra 81 52 Gaussian Elimination 81 521 Operation Counts 85 522 LU Factorization 87 523 Pivoting 88 524 Banded Matrices and Matrices for which Pivoting is Not Required 91 525 Implementation Considerations for High Performancequot 93 53 Other Methods for Solving Aw b 96 54 Conditioning of Linear Systems 98 541 Norms 99 542 Sensitivity of Solutions of Linear Systems 102 55 Stability of Gaussian Elimination with Partial Pivoting 105 56 Least Squares Problems 106 561 The Normal Equations 107 562 QR Decomposition 108 563 Fitting Polynomials to Data 110 Polynomial and Piecewise Polynomial Interpolation 117 The Vandermonde System 117 62 The Lagrange Form of the Interpolation Polynomial 117 63 The Newton Form of the Interpolation Polynomial 119 631 Divided Differences 121 64 The Error in Polynomial Interpolation 123 65 Piecewise Polynomial Interpolation 126 651 Piecewise Cubic Hermite Interpolation 128 652 Cubic Spline Interpolation 129 Numerical Differentiation and Richardson Extrapolation 134 71 Numerical Differentiation 134 72 Richardson Extrapolation 139 Numerical Integration 144 81 Newton Cotes Formulas 144 82 Formulas Based on Piecewise Polynomial Interpolation 148 83 Gaussian Quadrature 150 831 Orthogonal Polynomials 151 84 Romberg Integration 155 85 Periodic Functions and the Euler Maclaurin Formula 156 86 Singularities 160 9 Numerical Solution of the Initial Value Problem for Ordinary Differential Equa tions 164 91 Existence and Uniqueness of Solutions 165 92 OneStep Methods 168 921 Euler s Method 168 922 Higher Order Methods Based on Taylor Series 172 923 Midpoint Method 173 924 Methods Based on Quadrature Formulas 174 925 Classical Fourth Order Runge Kutta and Runge Kutta Fehlberg Methods 175 926 Analysis of One Step Methods 176 927 Practical Implementation Considerations 179 928 Systems of Equations 180 93 Multistep Methods 181 931 Adams Bashforth and AdamsMoulton Methods 181 932 General m Step Methods 183 933 Linear Difference Equations 186 934 The Dahlquist Equivalence Theorem 189 94 Stiff Equations 189 941 Absolute Stability 191 942 Backward Differentiation Formulae BDF Methods 195 943 Implicit Runge Kutta IRK Methods 196 95 Solving Systems of Nonlinear Equations in Implicit Methods 196 951 Fixed Point Iteration 197 952 Newton s Method 198 10 More Numerical Linear Algebra Eigenvalues Singular Values and Iterative Methods for Solving Linear Systems 101 Eigenvalue Problems 205 1011 The Power Method for Computing the Largest Eigenpair 212 1012 Inverse Iteration 215 1013 Rayleigh Quotient Iteration 216 1014 The QR Algorithm 217 1015 Google s PageRank 219 11 Numerical Solution of TwoPoint Boundary Value Problems 221 111 An Application Steady State Temperature Distribution 221 112 Finite Difference Methods 222 1121 Accuracy 224 1122 More General Equations and Boundary Conditions 229 113 Finite Element Methods 233 1131 Accuracy 238 12 Numerical Solution of Partial Differential Equations 243 121 Elliptic Equations 244 1211 Finite Difference Methods 245 1212 Finite Element Methods 249 122 Parabolic Equations 250 1221 Semidiscretization and the Method of Lines 251 1222 Discretization in Time 251 123 Separation of Variables 259 1231 Separation of Variables for Difference Equations 262 124 Hyperbolic Equations 263 1241 Characteristics 263 1242 Systems of Hyperbolic Equations 265 1243 Boundary Conditions 266 1244 Finite Difference Methods 266 125 Fast Methods for Poisson s Equation 270 1251 The Fast Fourier Transform 272 126 Multigrid Methods 275 Chapter 1 Basic Operations with MATLAB This book is concerned with the understanding of algorithms for problems of continuous mathe matics Part of this understanding includes the ability to implement such algorithms To avoid distracting implementation details however we would like to accomplish this implementation in the simplest way possible even if it is not necessarily the very most ef cient One system in which algorithm implementation is especially easy is called MATLAB 19 short for MATriX LABoratory While a helpful academic and instructive tool MATLAB is used in industry and government as well For instance systems engineers at the NASA Jet Propulsion Laboratory used MATLAB to understand system behavior before launching the Mars Exploration Rover MER spacecraft into space Figure 11 Left An artist s concept of the Mars rover Right Custom statistical MATLAB visualizations that were used to predict how the onboard systems would respond under various atmospheric conditions during descent to the Mars surface This chapter contains a short description of basic MATLAB commands and more commands are described in programs throughout the book For further information about using MATLAB we recommend 11 11 Launching MATLAB MATLAB is a high level programming language that is especially well suited to linear algebra type computations but it can be used for almost any numerical problem Following are some of the basic features of MATLAB that you will need to carry out the programming exercises in this book Depending on what system you are using you will start MATLAB either by double clicking on a MATLAB icon or by typing matlab or by some similar means When MATLAB is ready for input from you it will give a prompt such as gtgt You can use MATLAB like a calculator For instance if you type at the prompt gtgt 123 then MATLAB returns with the answer ans 7 Since you did not give a name to your result MATLAB stores the result in a variable called ans You can do further arithmetic using the result in ans gtgt ans4 and MATLAB will return with the result ans 1 7500 12 Vectors MATLAB can store row or column vectors The following commands gtgt V 1 2 3 4 gtme gtgt w 5 6 7 8 5 6 7 8 create a column vector v of length 4 and a row vector W of length 4 In general when de ning a matrix or vector semicolons are used to separate rows While commas or spaces are used to separate the entries Within a row You can refer to an entry of a vector by giving its index gtgt v2 ans 2 gtgt w3 ans 7 MATLAB can add two vectors of the same dimension but it cannot add v and W because v is 4 X 1 and W is 1 X 4 If you try to do this MATLAB will give an error message gtgt VW Error using gt Matrix dimensions must agree The transpose of W is denoted W gtgt w ans 0040301 You can add v and W using ordinary vector addition gtgtvw ans Suppose you wish to compute the sum of the entries in v One way to do this is the following gtgt v1 v2 v3 v4 ans 10 Another way is to use a for loop gtgt sumv 0 gtgt for i1z4 sumv sumv vi end gtgt sumv sumv 10 This code initializes the variable sumv to 0 It then loops through each value 2 1234 and replaces the current value of sumv with that value plus vi The line with the for statement actually contains three separate MATLAB commands It could have been written in the form for i1z4 sumv sumv Vi end MATLAB allows one line to contain multiple commands provided they are separated by commas or semicolons Hence in the oneline version of the for loop we had to put a comma or a semicolon after the statement for i1 4 This could have been included in the threeline version as well but it is not necessary Note also that in the threeline version we have indented the statements inside the for loop This is not necessary but it is a good programming practice It makes it easy to see which statements are inside and which are outside the for loop Note that the statement sumv 0 is followed by a semicolon as is the statement sumv sumv vi inside the for loop Following a statement by a semicolon suppresses printing of the result Had we not put the semicolon at the end of the rst statement MATLAB would have printed out the result sumv 0 Had we not put a semicolon after the statement sumv sumv vi then each time through the for loop MATLAB would have printed out the current value of sumv In a loop of length 4 this might be acceptable in a loop of length 4 million it probably would not be To see the value of sumv at the end we simply type sumv without a semicolon and MATLAB prints out its value 13 Getting help Actually the entries of a vector are most easily summed using a builtin MATLAB function called sum If you are unsure of how to use a MATLAB function or command you can always type help followed by the command name and MATLAB will provide an explanation of how the command works gtgt help sum SUM Sum of elements For vectors SUMX is the sum of the elements of X For matrices SUMX is a row vector with the sum over each column For ND arrays SUMX operates along the first nonsingleton dimension SUMXDIM sums along the dimension DIM Example If X 0 1 2 3 4 5 then sumX1 is 3 5 7 and sumX2 is 3 12 See also PROD CUMSUM DIFF In general you can type help in MATLAB and receive a summary of the classes of commands for which help is available If you are not sure of the command name for which you are looking there are two other helpful commands First typing helpdesk displays the help browser which is a very helpful tool Second if you are interested in commands related to summing you can type lookfor sum With this command MATLAB searches for the speci ed keyword sum in all help entries This command results in the following gtgt lookfor sum TRACE Sum of diagonal elements CUMSUM Cumulative sum of elements SUM Sum of elements SUMMER Shades of green and yellow colormap UIRESUME Resume execution of blocked Mfile UIWAIT Block execution and wait for resume RESUME Resumes paused playback RESUME Resumes paused recording It may take some searching to nd precisely the topic andor command that you are looking for 10 14 Matrices MATLAB also works with matrices gtgt A 1 2 3 4 5 6 7 8 0 A 1 2 3 4 5 6 7 8 0 There are many builtin functions for solving matrix problems For example7 to solve the linear system Aw b type Ab gtgt x Ab X 06667 03333 00000 Note that the solution is printed out to only four decimal places It is actually stored to about sixteen decimal places See Chapter 7 To see more decimal places7 you can type gtgt format long gtgtX 0 66666666666667 0 33333333333333 0 00000000000000 Other options include format short e and format long e to display numbers using scienti c notation You can check this answer by typing b Ax The notation Ax denotes standard matrix vector multiplication and standard vector subtraction is used when you subtract the result from b This should give a vector of 0 s if x solves the system exactly Since the machine carries only about 16 decimal digits we do not expect it to be exactly zero but as we will see later it should be just a moderate size multiple of 10 16 gtgt format short gtgt b Ax ans 10e15 00740 02220 0 This is a good result 15 Creating and running m les Typing MATLAB commands at the keyboard is ne if you are doing a computation once and will never need to make modi cations and run it again Once you exit MATLAB however all of the commands that you typed may be lost To save the MATLAB commands that you type so that they can be executed again you must enter them into a le called lenamem Then in MATLAB if you type lename it will run the commands from that le The m le can be produced using any text editor such as the one that comes up as part of the MATLAB window Once you save this le it will be available for future use 16 Comments Adding documentation in your MATLAB code allows you and others to maintain your code for future use Many a programmer has coded what appears to be a crystal clear implementation of an algorithm and later returned to be lost in the listing of commands Comments can help to alleviate this problem Adding comments to Matlab code is easy Simply adding a makes the remaining portion of that line a comment In this way you can have an entire line be a comment as in Z Solve Axb or you can append a comment after a MATLAB command as in x Ab 397 This solves the linear system Axb and stores the result in x The text following the is simply a comment for the programmer and is ignored by MATLAB during runtime 17 Plotting Tables and gures are usually more helpful than long strings of numerical results Suppose you are interested in viewing a plot of cos50w for O S x S 1 You can create two vectors one consisting of w values and the other consisting of corresponding 3 cos50w values and use the MATLAB plot command To plot the values of cos50w at w 0 01 02 1 type gtgt x 0011 397 Form the row vector of x values gtgt y cos50x 397 Evaluate cos50x at each of the x values gtgt plotxy 397 Plot the result Note that the statement x 00 1 1 behaves just like the for loop gtgt for i111 xi 0 1i1 end Note also that unless otherwise speci ed each of these statements produces a row vector In order to produce a column vector one could replace xi in the above for loop by xi 1 or one could replace the statement in the original code by x 0 0 1 1 Clearly this small number of evaluation points will not produce very high resolution as is seen in the plot produced by MATLAB depicted below n2 DA DE DE The plot command will graph the points w1y1 w2y2 which also implies that w and 3 must be vectors of the same length We did not use enough points to achieve a high resolution Therefore we will change an to include more points and include a title and labels for the axes in the plot with the commands below gtgt x 00011 397 Create a vector of 101 x values gtgt plotxcos50x 397 Plot x vs cos50x without explicitly storing cos50x gtgt title Plot of x versus cos50x gtgt ylabel cos50x gtgt xlabel x The plot resulting from these commands is Plat urx versus us5EIgtlt as5nx E E E E a E E E E n2 DA DE as X It is also possible to plot more than one function on the same graph To plot the two functions x cosw and gw w on the same graph type gtgt plotxcos50xxx The result is shown below us5EIgtlt gtlt n2 DA DE DE Notice the small boxed legend on the plot This was added to the plot with the following command legend cos50x x Another way to plot two functions on the same graph is to rst plot one then type hold on and then plot the other If you do not type hold on then the second plot will replace the rst but this command tells MATLAB to keep plots on the screen after they are created To go back to the default of removing old plots before new ones are added type hold off Again type help plot for more information on plotting or type doc plot which posts the Helpdesk documentation on the plot command Other useful commands are axis and plot3 For a bit of fun type the following commands gtgt x 0000110 gtgt cometxcos3x For more information on options available for plots type the commands hndlgraf hndlaxis and ardemo 18 Creating your own functions You can create your own functions to use in MATLAB There are two options for this If your function is simple eg x 32 253 then you may enter it using the command inline gtgt f inline xquot2 2x Inline function fx xquot2 2 Note the quot2 notation The expression xquot2 produces the square of x if x is a scalar but it gives an error message if x is a vector since standard vector multiplication is de ned only if the inner dimensions of the vectors are the same ie the rst vector is 1 by n and the second is n by 1 or the rst is n by 1 and the second is 1 by The operation quot applied to a vector however squares each entry individually Since we may wish to evaluate the function at each entry of a vector of x values we must use the quot operation To evaluate f at the integers between 0 and 5 type gtgt f 05 ans 0 3 8 15 24 35 If the function is more complicated you may create a le whose name ends in m which tells MATLAB how to compute the function Type help function to see the format of the function le Our function here could be computed using the following le called f m function output fx output xquot2 2x This function is called from MATLAB in the same way as above ie f x where x can be a scalar or vector 19 Printing While graphical output can be an important visualization tool numerical results are often presented in tables There are several ways to print output to the screen in MATLAB First to simply display a variable s contents use the command display 15 gtgt x 05z2 gtgt displayx X 0 05000 10000 15000 20000 In many cases the extra carriage return imposed by the display command clutters printed results Therefore another helpful command is disp which is similar to the display command gtgt dispx 0 05000 10000 15000 20000 You may still Wish to have the variable name printed Concatenating an array of text for output accomplishes this purpose gtgt disp x num2strx x 0 05 1 15 2 For more information type help num2str Tables can be created such as the following gtgt disp Score 1 Score 2 Score 3 disprand53 Score 1 Score 2 Score 3 04514 03840 06085 00439 06831 00158 00272 00928 00164 03127 00353 01901 00129 06124 05869 You may nd the fprintf command easier to use for tables of results For instance consider the simple loop gtgt fprintf x sqrt x nn for i15 fprintfquot7f fn isqrti end x sqrt x 1000000 1000000 2000000 1414214 3000000 1732051 4000000 2000000 5000000 2236068 The fprintf command takes format speci ers and variables to be printed in those formats The f format indicates that a number will be printed in xed point format in that location of the line and the n forces a carriage return after the two quantities i and sqrti are printed You can specify the total eld Width and the number of places to be printed after the decimal point by replacing f by say 784f to indicate that the entire number is to be printed in 8 spaces with 4 places printed after the decimal point You can send your output to a le instead of the screen by typing fid fopen sqrttxt w fprintf fid a X sqrtxnn for i115 fprintf fid quot740f 784fn isqrti end which prints the following table in a le called sqrttxt x sqrtx 1 1 0000 2 1 4142 3 1 7321 4 20000 5 22361 Again for more information refer to MATLAB documentation 110 More loops and conditionals We have already seen how for loops can be used in MATLAB to execute a set of commands a given number of times Suppose instead that one wishes to execute the commands until some condition is satis ed For example one might approximate a root of a given function f x by rst plotting the function on a coarse scale where one can see the approximate root then plotting appropriate sections on ner and ner scales until one can identify the root to the precision needed This can be accomplished with the following MATLAB code xmin input Enter initial xmin xmax input Enter initial xmax tol input Enter tolerance while xmaxxmin gt tol x xminzxmaxxmin100xmax y fx plotxy xmin input Enter new value for xmin xmax input Enter new value for xmax end The user looks at each plot to determine a value xmin that is just left of the root and a value xmax that is just right of the root The next plot then contains only this section so that closer values xmin and xmax can be determined In the next chapter we discuss more e icient ways of nding a root of fx Another important statement is the conditional if statement In the above code segment one might wish to let the user know if the code happens to nd a point at which the absolute value of f is less than some other tolerance say delta This could be accomplished by inserting the following lines after the statement y f x yminindex minabsy 397 This finds the minimum absolute value of y o and its index if ymin lt delta fprintf f f fn xindex yindex 397 This prints the x and y 397 values at this index end The if statement may also contain an else clause For example to additionally write a message when ymin is greater than or equal to delta one could modify the above if statement to say if ymin lt delta fprintf f f fn xindex yindex Z This prints the x and y K values at this index else fprintf No points found where fx lt fn delta end 111 Clearing variables You may clear a particular variable by typing gtgt clear x or all variables with gtgt clear all This is important when you want to be sure that all variable names have been erased from memory 112 Logging your session You can keep a record of your MATLAB session by typing gtgt diary hw1txt some other commands gtgt diary off This command records all subsequent commands that you type and all responses that MATLAB returns in a le named hw1txt You will want to name the le by replacing hw1txt with a more descriptive name related to your work Note however that you cannot then run the le from MATLAB this is simply a device for recording what happened during your keyboard session Note if you execute an M le in a session logged with the diary command you may want to type echo on77 before executing the M le In this way the commands in the M le are echoed along with MATLAB s response Otherwise the diary le will contain only the responses not the commands 113 More advanced commands It is perhaps apparent from the reference at the beginning of this chapter to the work on the Mars Exploration Rover7 that MATLAB has a large number of commands As we close this chapter7 let us demonstrate some of the graphical capabilities through an example The commands XY meshgrid3 125z3 Z peaksXY meshcXYZ axis 3 3 3 3 10 5 produce the plot In order to understand this plot7 search MATLAB documentation for the commands meshgrid peaks meschc and axis While this chapter will get you started using MATLAB7 effective use of the MATLAB documentation will be the key to proceeding to more complicated programs MATLAB S omle c E ow 15 mm and chef menus at The Mthorks h che rmd on me 19m he wee one of che author of LENPACK and MS PACK Roman hhmzee for nvmeucel compucmg Th gve easy wees 7 ck mmelcxehze che MATLAB em eve Mole lecewed hm beLhelol 5 degae 1m Celcech m 1961 and hD mum Molezwaeapie f enh d w p m e est 315 a che Unxveln The MeLhWolk5 m1 xme m 1939 he 91w walked 1m anel Hypelcube and Alden D1 Mole way eleeced ra che Namonel Academy of Engnaenng m 199w He hay lewwed honolely degeee m Lmkopmg Umvelsxcy Sweden che Umvemcy of Wazezloo and che Ieehmeel Umvezncy of Denmelk Exercises 1 2 chh the fullowmg mamaee and Veer315 1073 A 4 2 d AT h The veerat y m whxch By u 2 AB 1 The veerat 2 fat whxch A2 71 emery 3 0 2w for k 4 Uee MATLAB no pun acehle oivemee z 31111 and 352 for z o g 27 Label che came of you cable 20 5 Download the le plotfunctionl m for the web page for this class and execute it This should produce the two plots shown below The top plot shows the function fw 2cosw 7 63 for 76 S x S 37 and from this plot it appears that fw has three roots in this interval The bottom plot is a zoomed view near one of these roots7 showing that fw has a root near an 71454 Note the different vertical scale on this plot Note also that when we zoom in on this function it looks nearly linear over this short interval This will be important when we study numerical methods for approximating roots zuumedww a Modify this script so that the bottom plot shows a zoomed view near the leftmost root Write an estimate of the value of this root to at least 3 decimal places You may nd it useful to rst use the zoom feature in MATLAB to see approximately where the root is and then to choose your axis command for the second plot appropriately A U v Edit the script from part a to plot the function 4 i 73 x wsgir w over the range 0 S x S 4 and also plot a zoomed view near the leftmost root Write an estimate of the value of the root from the plots that is accurate to 3 decimal places Once you have de ned the vector x properly you will need to use appropriate componentwise multiplication and division to evaluate this expression y 4xsinx 3 2 xquot2 6 Plot each of the functions below over the range speci ed Produce 4 plots on the same page using the subplot command a fw Inc 7 1 for 73 S x S 3 Use abs in MATLAB fw w for 74 S x S 4 Use sqrt in MATLAB fw a expiw2 for 74 g x g 4 Use exp in MATLAB AA GU Vv 21 7 90 d W 1032 1 Use MATLAB to plot the circles for72 w 2 06223712 2 w 7 252 32 35 and zoom in on the plot to determine approximately where the circles intersect Hint One way to plot the rst circle is the following theta linspace0 2pi 1000 sqrt2 2 rcostheta y 1 rsintheta plotxy axis equal I X 397 so the circles look circular Could also use axis square Use the command hold on after this to keep this circle on the screen while you plot the second circle in a similar manner In this exercise you will plot initial stages of a process that creates a fractal known as Koch s Snowflake which is depicted below Stage 0 Stage 1 Stage 2 This exercise uses the MATLAB M le kochm which you will nd on the course s webpage The M le contains all the necessary commands to create the fractal except for the necessary plotting commands Edit this M le so that each stage of the fractal is plotted Hint 7 you will need to add a plot command before the for loop and also in the for loop Add the following commands to keep consistency between plots in the animation axis075 075 sqrt36 1 axis equal Note the cla command clears the axes Finally add the command pause 0 5 in appropriate places to slow the animation The fill command as opposed to plot produced the lled fractals depicted above We will create fractals using Newton s Method in the next chapter 22 9 0 H A magic square is an arrangement of the numbers from 1 to n2 in an n X 71 matrix where each number occurs exactly once and the sum of the entries of any row any column or any main diagonal is the same The MATLAB command magicn creates an n X n where n gt 2 magic square Create a 5 X 5 magic square and verify using the sum command in MATLAB that the sum of the columns rows and diagonals are equal Create a log of your session that records your work Hint To nd the sum of the diagonals read documentation for the diag and the flipud commands More advanced plotting commands can be useful in MATLAB programming a In the MATLAB command window type XYZ peaks30 surfXYZ b Give this plot the title 3 D shaded surface plot c Type colormap hot and observe the change in the plot d Print the resulting plot with the given title Computer graphics make extensive use of matrix operations For example rotating an object is a simple matrix vector operation In two dimensions a curve can be rotated through an angle 6 about the origin by multiplying every point that lies on the curve by the rotation matrix 7 sin 6 R 7 lt cosO gt 39 As an example let us rotate the square with vertices 1 75 2 75 2 5 and 1 5 through the angle 6 7r4 about the origin In MATLAB type the following code to generate the original and rotated squares plotted one on top of the other cos 6 sin 6 397 Create matrix whose first column contains coordinates of the first vertex 397 whose second column contains coordinates of the second vertex etc U 1 22 1 5 5 5 5 397 Create a red unit square 397 Note U1 creates a vector of the first row of U fillU1U2 r Z Use axis equal to make square look square then plot between 4 and 4 Z so that rotated squares can be seen Retain current plot so that Z subsequent graphing commands add to the existing graph axis equal axis4 4 4 4 hold on 10 397 Create rotation matrix theta pi4 R costheta sintheta sintheta costheta 397 Rotate each of the vertices and color rotated square blue U RU fillU1 U2 b gt Note that the fill command in MATLAB plots a lled polygon determined by two vectors of vertices a Continue rotating the square by applying the rotation matrix R to U over and over again and plotting the new squares in different colors You should nd that after applying R eight times the new square lands on top of the original Adapt this code to plot a triangle with vertices 50 62 and 4 1 Plot the triangles resulting from rotating this triangle 7r2 7r and 37r2 radians about the origin Plot all four triangles on the same set of axes b A O v Rotating by 6 radians and then again rotating by 76 radians leaves the gure unchanged This corresponds to multiplying rst by R6 and then by R76 Using MATLAB verify that these matrices are inverses of one another that is their product is the identity for 0 7r 3 and 7r 4 Using the trigonometric identities cos cos7 and isin sin7 prove that R and R7 are inverses of one another for any 6 A 1 V In Figure 12 we see a welliknown model in computer graphics known as the Stanford Bunny77 This triangulation is formed using 3851 triangles with 1889 vertices Suppose V is a matrix with three columns where row 2 contains the w y and z coordinate of the ith vertex in the model The image can be rotated by 6 radians about the y axis by multiplying V by Ry where cos 6 O 7 sin 6 Ry 0 1 0 sin 6 0 cos 6 Download the le bunnyzip from the web page for this course Unzip the le and run bunnym in MATLAB This will plot the triangulation of the Stanford Bunny a Edit the code to rotate the image by 7r6 radians about the y axis b The matrix containing the vertices is an 1889 X 3 matrix How many scalar multiplica tions are performed when you use matrix multiplication with Ry to rotate the image by 7r6 radians as you just did Keep this in mind as you watch how fast MATLAB performs the calculations and plots the results 24 k ail AK QNNE SAZAEEW Z waw 717 V lt 4 seswame v i Jr 44 V 2 a V sigma 7 3977 Vivi a g 5 V v w39 m lt nsas7 uzl LLIQLXN A V 452 39 W MEMVrM mi ANN E L mm gynmmuvuamw Am r sumsmax 39r V a gggmmh nsr g VA 4 V min k ngngmm unsAV mggn wnmk vav n A V s NEW m gsmu gnmi n A has AAL Vnwaummmmmn gnnmzmmwmmw nuanmi immvi m 39 A 1 ggg gg 4 gtxk EAEASEZE Figure 12 Triangulations are used in computer graphics to create models of objects in a computer Here we see the triangulation of the welliknown Stanford Bunny77 7 Need to attain permission Taken from httpgraphicsstanfordedudata3Dscanrep c Further edit the code to continuously rotate the image by 7r24 radians until the image has rotated one full rotation about the y axis Before drawing a new rendering of the rotated image you will want to clear the current plot and possibly include a short pause These two steps can be performed with the following code cla 397 clear current plot pause05 397 pause for 05 seconds You may want to look for the Rotate 3D button in the Figure containing the triangulation You can rotate the image using your mouse Of special note look at the bottom of the bunny and notice how you can look inside the model as if it were a porcelain gurine H 9quot Type helpdesk and nd documentation on the movie command or alternatively type doc movie At the bottom of the documentation on this command you will nd the following code Z peaks surf Z axis tight setgca nextplot replacechildren for j 1 20 surfsin2pij20ZZ Fj getframe end movie F 5 Cut and paste this code into the MATLAB command Window and describe the results Pro ciency in MATLAB programming can increase dramatically through effective use of the available documentation Chapter 2 Solution of a Single Nonlinear Equation in One Unknown Consider the problem facing an artillery o icer during battle As the sun rises he sees the position of the advancing enemy army Overnight their artillery has dug into a position d meters away The officer wishes to aim his cannon to hit the entrenched enemy He knows that the cannonball will leave the muzzle of the gun at an initial speed of 20 meters per second Neglecting air resistance the only force acting on the cannonball after it leaves the cannon is gravity and the height yt of the cannonball at any time t gt 0 satisfies the di erential equation y t 7g s 792 msecz To determine the height uniquely one can specify the initial height and vertical velocity of the cannonball 30 0 m y 0 v0 sin6 msec where 6 is the angle at which the cannon is aimed Note that the initial value 30 0 is an approximation the tip of the cannon muzzle surely isn t at ground level but it is not far o quot One can solve this di erential equation by integrating 20 gt 01 27 and since 30 110 sin 6 the constant cl must be 110 sin 6 Integrating again gives 1 2 yt figt vos1n6t 02 where 02 0 since 310 0 Now that we have a formula for the height at any time t gt 0 we can determine the amount of time before the cannonball hits the ground 1 2 6 07igt2vosin t t0 ort Finally neglecting air resistance the cannonball travels at a constant horizontal velocity 110 cos 6 so that it hits the ground a distance 2113 sin 6 cos 6 m 9 away from where it was red The problem thus becomes one of solving a nonlinear equation 22 6 6 vosm cos d g 21gt for the angle 6 at which to aim the cannon While one could solve this problem numerically perhaps using an available software package for solving nonlinear equations there are several things that one should take into account before returning to the cannoneer with a de nitive value for 6 1 The model developed above is only an approximation Consequently even if one returns with the exact solution 6 to 21 the red cannonball may end up missing the target A strong tailwind might cause the o icer s shot to sail over the approaching enemy Likewise the amount of gunpowder may vary between shots which corresponds to 110 in 21 changing for each shot Such factors may make adjustments to the model necessary 10 Equation 21 may not have a solution The maximum value of sinO cos 6 occurs at 6 7r4 2 where sinO cosO If 01 gt 70 then the problem has no solution In this case the enemy is simply out of range 9quot If a solution does exist it may not be unique lf 6 solves problem 21 then so does 6 i2k7r k 12 Also 7 6 solves the problem Either 6 or 721 7 6 might be preferable as a solution depending on what lies in front of the cannon q Actually this problem can be solved analytically Using the trigonometric identity 2 sin 6 cos 6 sin2 one nds 1 d9 6 7 arcs1n 7 2 110 28 In light of item 1 above however it may be prudent to go ahead and write a code to handle the problem numerically since more realistic models may require approximate solution by numerical methods By a good approximate solution we mean one that is acceptably close for the application under consideration A computed angle that results in the cannonball landing within a foot of its target is probably acceptable in this situation but if the equation arose from an application in laser surgery considerably more accuracy would likely be required 01 To solve a nonlinear equation numerically we usually start by writing it in the form f 0 Equation 21 can be written as 7 2112 sin 6 cos 6 7 d 9 This function is simple enough that it can be di ferentiated analytically We will see in the following sections that such functions can then be solved using Newton s method If we made the model more realistic however this might not be the case For example if we included the effect of air resistance in our model then we might be unable to solve the differential equation for yt analytically We could approximate its solution numerically for a given value of 6 In this case the problem could be solved using bisection the secant method or quasi Newton methods to be described in this chapter f6 0 The point of this discussion is that if one blindly tries to solve a physical problem using a numerical software package the results may be less than satisfactory Consideration of the model being used whether or not the problem has a solution and if so how many solutions and whether one solution is to be preferred to another are important preprocessing steps Likewise consideration of the level of accuracy required in the result will be of importance 21 Bisection Faced with an oncoming army the artillery o icer might use the following procedure to solve the problem He would guess the angle at which he should aim the cannon re it and see where the cannonball lands If it fails to hit the target he would adjust the angle and try again If he nds that at one angle 61 the cannonball goes too far and at another 62 it lands short then he might try some angle between 61 and 62 say their average 9192 This simple procedure is the method of bisection and it can be used to solve any equation at 0 provided f is a continuous function of w and we have starting values an and avg where f has different signs The basis for the bisection method is the Intermediate Value Theorem Theorem 211 ff is continuous on 11 Indy lies between fa and fb then there is a point x 6 11 where 3 If one of fa and fb is positive and the other negative then there is a solution to x O in 1 b In the bisection algorithm we start with an interval 11 where signfa 7 signfb We then compute the midpoint 424 and evaluate f We then replace the endpoint a or b where f has the same sign as f by TH and repeat 29 BISECTING A PHONE BOOK Have a friend choose a name at random from your local phone book Bet your friend that in less than thirty yes or no questions you will be able to guess the name that was chosen This task is simple Place about half of the pages in your left hand or right if you prefer and ask your friend Is the name in this half of the phone book77 In one question you have eliminated approximately half of the names Holding half of the non eliminated pages in your hand ask the question again Repeat this process until you have discovered the page which contains the name Continue to ask questions which halve the number of possible names until you have discovered the name Like the Bisection Algorithm you use the speed of exponential decay to your advantage Since 1230 is about 93 X 10 10 and the population of the United States is approximately 280 million you can see that any city phone book will suf ce and should leave you a comfortable margin for error Following is a simple MATLAB code that implements the bisection algorithm This routine uses bisection to find a zero of a usersupplied o o 397 continuous function f The user must supply two points a and b o o such that fa and fb have different signs The user also 397 supplies a convergence tolerance delta fa fa if fa root a break end 397 Check to see if a is a root fb fb if fb root b break end 397 Check to see if b is a root if signfasignfb Z Check for an error in input display Error fa and fb have same sign break end while absba gt 2delta 397 Iterate until interval width lt 2delta C ab2 397 Bisect interval fc fc 397 Evaluate f at midpoint if fc root c break end 397 Check to see if c is a root if signfcsignfa 397 Replace a by c fa by fc a c fa fc else b c fb fc end end root ab2 397 Replace b by c fb by fc 397 Take midpoint as approximation to root There are several things to notice about this code First it checks to see if the user has made an error in the input and supplied two points where the sign of f is the same When writing code for others to use it is always good to check for errors in input Of course typing in such things is timeconsuming for the programmer so when we write codes for ourselves we sometimes omit these things And sometimes live to regret it Next it checks to see if it has happened upon an exact root that is a point where f evaluates to exactly zero If this happens it returns with that point as the computed root This is not necessary it could continue to bisect the interval until the interval width reaches 2delta The MATLAB sign function returns 1 if the argument is positive 1 if it is negative and 0 if it is 0 so if fc were 0 the code would end up replacing b by C Example 1 Use bisection to nd a root of ne 33 7 2 7 1 andb2 in the interval 11 where a O D an bn Tn 1 0 2 1 1 2 1 2 15 0125 3 1 15 125 0609 4 125 15 1375 0291 5 1375 15 14375 0096 6 14375 15 146875 0011 7 14375 146875 1453125 0043 8 1453125 146875 14609375 0016 A graphical depiction of these results is given below 31 aisecuun7m x37x271 Rate of Convergence What can be said about the rate of convergence of the bisection algo rithm Since the interval size is reduced by a factor of 2 at each step the interval size after k steps is lb 7 ll2k which converges to zero as k 7 00 To obtain an interval of size 26 we need lial b b7 b7 7326 cgt 2k12 cgt k210g2lt gt7l 6 If the approximate root is taken to be the midpoint of this interval then it differs from an actual root by at most 6 Example 2 If a 1 and b 2 and we want to guarantee an error less than or equal to 104 how many iterations do we need to take 5 7 lt2104 gt 2k1gt7lll cgt kgtlog2lt gt71 b 7 1 b 7 2k 7 7 10 10 4 Since I 7 a 1 then k 2 1229 iterations are required to guarantee an error less than or equal to 104 After 13 iterations your answer could have an error much less than 10 4 but will not have an error greater than 104 Since each step reduces the error ie half the bracket size by a constant factor ie a factor of 2 the bisection algorithm is said to converge linearly The most dif cult part of using the bisection algorithm may be in nding an initial interval where the sign of f is different at the endpoints Once this is found the algorithm is guaranteed to converge Note that the bisection algorithm cannot be used for a problem like 32 0 pictured below since the sign of the function does not change 22 Taylor s Theorem Pinbany the mostrmai themem m numeucal ab91515 Tayluz y themem especxalb Taylm y themem thh lemamdel t 15 Wmth levxewmg chm xmpultant zemu A43 BROOK Tamouwsmwo Theorem 221 Taylor s theorem with remainder Let f w n be commmm on ab and 22 WWW mt m allz e ltab Thmthere u anmberfe ltab suchth b7 1 b7 Nb flta ltbr awn lt 2 M M Ma lt1 0 o WW lt5 lt22 Rulmula lt2 2 can alga be wutten u mg the numnatxun natatxun Mb i lt We L V xW N 7 n1 the mm 74 i a 7x 74 WW W 1W 7 M r i xltwgtlt x 30 is sometimes written as 0537 171 since there is a nite constant C namely f 1gtan1l such that R gt n w i l af a Sometimes the notation is used to mean in addition that there is no higher power 9 gt n 1 such that limzna 7 17 is nite This will be the case if the constant C above is nonzero If f is in nitely di ferentiable then the Taylor series expansion of f about the point a is w 7 a j fowl Mg V H 0 Depending on the properties of f the Taylor series may converge to f everywhere or throughout some interval about a or only at the point a itself For very well behaved functions like 63 which we will also denote as expw sinw and cosw the Taylor series converges everywhere The Taylor series about an 0 occurs often and is called a Maclaurin series The Taylor series contains an in nite number of terms If we are interested only in the behavior of the function in a neighborhood around an a then only the rst few terms of the Taylor series may be needed j fmmmmZji w 70 The Taylor polynomial Pnw may serve as an approximation to the function Example 3 Find the Taylor series of 63 expw about 1 First nd the value of the function and its derivatives at w 1 f 1 eXp1 f 1 eXp1 eXp1 Next construct the series using these values and a 1 7 1 2 7 1 3 expw 1 w71 f 1 flt3gt1 expm expo explt1gt explt1gt wily 7003 21 31 Example 4 Find the Maclaurin series for cosz First nd the value of the function cosz and its derivatives at z 0 f0 cos0 1 fWO sin0 0 no 75110 0 Ma cos0 1 f 0 7cos0 71 Next construct the series using these values and a 0 2 3 4 cosz f0z How f 0 1 Wow lt4gt0 Vo 2 4 31 4 4 71 0 1 7 z2 z4 7 W flycsz 1 E Z 2k Example 5 Find the fourth order Taylor polynomial P4 of cosz about I 0 Using the work from the previous exercise we nd 1 1 P4z 17 5902 if 1n the following section Taylor s Theorem will be used to derive a numerical method At the end of the chapter exercise 18 discusses how John Machin a contemporary of Brook Taylor used Taylor s Theorem to nd 100 decimal places of 7r 23 Newton s Method Isaac Newton a founder of Calculus described a method which became known as Newton s method in his Method of Fluzions and In nite Series 6 16 ISAAC NEWTON 16431727 35 We start with a geometric derivation of Newton s method Suppose we wish to solve fw O and we are given an initial guess 30 for the solution We evaluate n and then construct the tangent line to f at 30 and determine where it hits the w axis If f were linear this would be the root of f In general it might give a better approximation to a root than am does as illustrated in Figures 21 a and The slope of this tangent line is f w0 and it goes through the point wofw0 so its equation is y 7 fw0 f w0w 7 30 It hits the w axis y 0 at w 30 7 fw0f w0 assuming f w0 7 O Letting an be the point where this tangent line hits the w axis this process can be repeated by constructing the tangent line to f at 31 determining where it hits the w axis and calling that point 32 etc This process is known as Newton s method Newton s Method For solving x 0 Given an initial guess 30 for k O12 f 00k Set wk1 wk 7 f M a fxx272 b fgtlt smgtlt B u 5 4 2 1 X2 Xn 7n 5 Xlxl gtltEl 72 u 2 71 n 1 cfxx 3 dfx37x71 Figure 21 Newton s Method As a historical note Newton explained his method in 1669 using numerical examples and did not use the geometric motivation given here involving the approximation of a curve with its tangent line His work also did not develop the recurrence relation given above which was developed by the English mathematician Joseph Raphson in 1690 It is for this reason that the method is often called the NewtonRaphson Method Unfortunately Newton s method does not always converge Figures 21 c and d illustrate cases in which Newton s method fails In the example of c Newton s method diverges while in d it cycles between points never approaching the actual zero of the function One might combine Newton s method with bisection to obtain a more reliable algorithm Once an interval is determined where the function changes sign we know by the Intermediate Value Theorem that that interval contains a root Hence if the next Newton step would land outside that interval then do not accept it but do a bisection step instead With this modi cation the method is guaranteed to converge provided it nds an interval where the function has different signs at the endpoints This combination would solve example c in Figure 21 because once 0 and ml were determined the next Newton iterate on would be rejected and replaced by 0 wl2 Example d in Figure 21 would still fail because the algorithm does not locate any pair of points where the function has opposite signs Newton s method can also be derived from Taylor s theorem According to that theorem 3 7 wo2 f0 fw0 wiwof wo 2 f for some between am and x If an 0 then 7 2 o mo w 7 wowwe WM 23 and if am is close to the root xi then the last term involving xi 7 am2 can be expected to be small in comparison to the others Neglecting this last term and de ning an instead of xi to be the point that satis es the equation when this term is omitted we have 0 fwo 001 7 wof7007 or fw0 7 1 0 fw0gt Repeating this process for k 012 gives ack wk wk 7 7 24 1 gt This leads to the following theorem about the convergence of Newton s method Theorem 231 ff E 02 if 0 is su iciently close to a root as of f and if f w 7 0 then Newton s method converges to 304 and ultimately the convergence rate is quadratic that is there exists a constant C f w2fw such that lim wk1 I C 16700 wk 7 wilz 25 37 Before proving this theorem we rst make a few observations The quadratic convergence rate described in 25 means that for k large enough convergence will be very rapid If C is any constant greater than 0 it follows from 25 that there exists K such that for all k 2 K wk1 7 Jul S Clwk 7 5642 lf say C 1 and le 7 cell 01 then we will have leH 7 cell 3 10 2 wK2 7 cell 3 10 4 wK37m 3 10 s etc There are some important limitations however First the theorem requires that the initial guess avg be su iciently close to the desired root am but it does not spell out exactly what su iciently close means nor does it provide a way to check whether a given initial guess is su iciently close to a desired root The proof will make it clearer why it is di cult to say just how close is close enough Second it shows that the ultimate rate of convergence is quadratic but it does not say how many steps of the iteration might be required before this quadratic convergence rate is achieved This too will depend on how close the initial guess is to the solution Proof of Theorem It follows from equation 23 with avg replaced by 5 that NM 70 7 06102 f k f 2 WWW for some k between wk and cm Subtracting this from equation 24 for wk gives f kgt WW Now since f is continuous and f w 7 0 if C f w2f wlt then for any C gt 0 there is an interval about at in which f w2f w S C If for some K wK lies in this interval and if also le 7 cell lt 10 which is sure to hold even for K 0 if am is su iciently close to am then 26 implies that leH 7 ml S Cle 7 w2 lt le 7 w so that wKH also lies in this interval and satis es wK1 7 ml lt lC Proceeding by induction we nd that all iterates 5 k 2 K lie in this interval and so by 26 satisfy wk1 7 w wk 7 an 26 k1 S Ck w2 Clwkiwlllmiwl clacrwiwwwki wilwIwkilwl S Clwk7wnCwK7w IwK7w4lt 3 lt0le 7 MWMIW 7 ml Since Cle 7 cell lt 1 it follows that CIwK 7 wk1 K 7 O as k 7 00 and this implies that wk 7 w as k 7 00 Finally since wk and hence Q in 26 converge to am it follows from 26 that H k1 Ml f le lwk MIZ WWW 7G As noted previously the hypothesis of 30 being suf ciently close77 to a root at where suf ciently close77 can be taken to mean that 30 lies in a neighborhood of x where f 2f is less than or equal to some constant C and lam 7 ml lt 1C is usually impossible to check since one does not know the root xi There are a number of other theorems giving different suf cient conditions on the initial guess to guarantee convergence of Newton s method but as with this theorem the conditions are usually dif cult or impossible to check or they may be far more restrictive than necessary so we do not include these results here Early mathematicians looking at conditions under which Newton s method converges include Mourraille in 1768 and later Lagrange In 1818 Fourier looked at the question of the rate of convergence and in 1821 and again in 1829 Cauchy addressed this question For more information about the historical development of the method see 6 Examples 1 Compute 2 using Newton s method ie solve 32 7 2 0 Since fw 32 7 2 we have f w 2x and Newton s method becomes 2 wk 7 2 k01 2M 77 wk1 wk 7 Starting with am 2 and letting 6k denote the error wk 7 2 we nd 0 2 60 31 15 61 0086 02 14167 62 0025 33 14142157 63 216 7 6 z 356 The constant f wlt2f w for this problem is 12x2 x 03536 Note that for this problem Newton s method fails if 30 O For 30 gt 0 it converges to x2 and for 30 lt 0 it converges to 7x2 2 Suppose f w O Newton s method may still converge but only linearly Consider the problem fw E 2 0 Since f w 2w Newton s method becomes 3 w Li 1 k1 k 2wk 2 k Starting for example with 30 1 we nd that 60 1 31 61 32 7 62 Z wk 6k Instead of squaring the error in the previous step each step reduces the error by a factor of 2 3 Consider the problem fw E 33 0 Since f w 332 Newton s method becomes wk 2 w w 7 7 7w k1 k 3 k Clearly the error the difference between wk and the true root at 0 is multiplied by the 2 factor g at each step 4 Suppose fw E 7 What is the error reduction factor for Newton s method To understand example 2 recall formula 26 from the proof of Theorem 21 f k WW Since f w O in example 2 the term multiplying wk 7 my becomes larger and larger as wk approaches all so we cannot expect quadratic convergence However if we expand f about 4 using a Taylor series we nd fwk fM 0 MVHW wk Mf nk7 for some 71k between 4 and 5 Substituting this expression for flak above gives f k 2f nkgt M M lf f wi 7 0 then this establishes linear convergence under the assumption that 0 is su iciently close to 4 with the convergence factor approaching By retaining more terms in the Taylor series expansion of f wk about 4 this same approach can be used to explain examples 3 and 4 2 wk1 7 w wk 7 53 wk1 7 4lt Theorem 232 If f E Cp1 for p 2 1 if 30 is su iciently close to a root 3 of f and if fwi fltpgtw 0 butfltp1gtw 7 0 then Newton s method converges linearly to 4 with the error ultimately being reduced by about the factor pp 1 at each step that is lim wk1 7 I 7p 16700 Iwk7wil p1 Proof Taking formula 24 and subtracting w from each side gives 27 wk17m wk7wi7m Now expanding fwk and f wk about all gives wk CMerl 131 fwk mg Iwi lp pHRW fm 7 rem for some points Q and 71k between wk and 4 since the rst 13 1 terms in the rst Taylor series expansion and the rst 1 terms in the second are zero Making these substitutions in 27 gives 1 fltP1gtlt kgtgt 28gt r m 7 1 mm 40 Since fltp1gtw 7 O and fwd is continuous for any 5 gt 0 there is an interval about at in which the factor 1 7 1p 1fp1gt fp1gtn is within 5 of its value at g n 3 namely 1 7 1p 1 pp 1 whenever g and 7 lie in this interval In particular for e suf ciently small this factor will be bounded above by a constant C lt 1 It follows from 28 that if wk lies in this interval then wk1 7 ml 3 Clack 7 ml so wk also lies in this interval and if 30 lies in this interval then by induction wk1 7 ml 3 19 le 7 ml so wk 7 w as k 7 00 The limiting value of the error reduction factor is pp 1 As another example consider the functions f1 sinw and f2w sin2w both of which have x 7r as a root Note that cosw 7 O at w 7r while 2sinw cosw O at w 7r Note however that f 7r 7 0 Therefore we expect quadratic and linear convergence with a convergence factor of about 12 for Newton s method applied to f1w and f2w respec tively This is seen graphically in Figure 22 For f1 the root at 7r is said to have multiplicity 1 or to be a simple root while for f2 it is a root of multiplicity 2 Newton s method converges only linearly for a root of multiplicity greater than 1 Cunvemence umamcms memuu my mum and mum mm Ahsnmle ennr m sum 2 w 5m 5 m 15 m 25 112mmquot numbev Figure 22 Convergence of Newton s method for a function with root of multiplicity 1 versus a function with root of multiplicity 2 24 QuasiNewton Methods 241 Avoiding Derivatives Newton s method has the very nice property of quadratic convergence when suf ciently close to a root This means that one can hone in on a root from a nearby point very quickly much faster than with the bisection method for example The price is that one must evaluate both f and its derivative at each step For the simple examples presented here this is not dif cult In many problems however the function f is quite complicated It might not be one that can be written down analytically Instead one might have to run a program to evaluate In such cases di ferentiating f analytically is dif cult to impossible and even if a formula can be found it may 41 be very expensive to evaluate For such problems one would like to avoid computing derivatives or at least compute as few of them as possible while maintaining something close to quadratic convergence lterations of the form ack n1 00k 9k 7 Where 9k fka 29 are sometimes called quasiNewton methods 242 Constant Slope Method One idea is to evaluate f once at am and then to set gk in 29 equal to f w0 for all k If the slope off does not change much as one iterates then one might expect this method to mimic the behavior of Newton s method This method is called the constant slope method wk wk 7 where g fwo 210 To analyze this iteration we once again use Taylor s theorem Expanding f about the root wt we nd ack wk 7 mm 0wk 7 My using the notation instead of writing out the remainder term explicitly Subtracting w from each side of 210 and again letting ek E wk 7 w denote the error in 5 we have k1 6k i y 6191 gt 00 If ll 7 f wf w0 lt 1 then for 0 su iciently close to x close enough so that the Oei term above is negligible the method converges and convergence is linear In general we cannot expect better than linear convergence from the constant slope method A variation on this idea is to update the derivative occasionally Instead of taking gk in 29 to be f w0 for all k one might monitor the convergence of the iteration and when it seems to be slowing down compute a new derivative f This requires more derivative evaluations than the constant slope method but it may converge in fewer iterations Choosing between methods for solving nonlinear equations usually involves a tradeoff between the cost of an iteration and the number of iterations required to reach a desired level of accuracy 243 Secant Method The secant method is de ned by taking gk in 29 to be 9 M 211 wk 7 1971 Note that in the limit as 1 approaches 5 gk approaches f wk so that gk can be expected to be a reasonable approximation to f wk when 1 is close to 5 It is actually the slope of the 42 line through the points wk1fwk1 and wkfwk called a secant line to the curve f since it intersects the curve at two points The secant method is then de ned by w w 7fkgtk k71gt 1911 k famifmfl k 12 212 To begin the secant method7 one needs two starting points am and an An illustration of the secant method for nding a root of at 33 7 32 7 1 with 30 05 and 31 2 is given in Figure 23 xx 7y 71 Figure 23 Finding a root of at 3 7 2 71 using am 05 and 31 2 with the secant method Example Consider again the problem at E 32 7 2 0 Taking 30 1 and 1 27 the secant method generates the following approximations and errors x2 13333 62 700809 303 14 63 700142 304 14146 e4 426 7 4 305 14142114 e5 7216 7 6 From the table7 it is di cult to determine exactly what the convergence rate of the secant method might be It seems faster than linear the ratio ek1ek is getting smaller with k But it seems slower than quadratic ek1 eZ seems to be growing larger with k Can you guess what the order of convergence is It turns out that the order of convergence of the secant method is 1435 x 162 I bet you didn t guess that This result is obtained from the following lemma Lemma 241 ff 6 CZ 139wa and 31 are su ciehtly close to a root 34 of f and if fwi 7 0 then the error ek in the secant method satis es lim 6k 01 213 19700 61961971 where C1 f w2fw Statement 213 means that for k large enough ek1 z 061961971 214 where the approximate equality can be made as close as we like by choosing k large enough Before proving the lemma let us demonstrate why the stated convergence rate follows from this result Assume that there exist constants a and 04 such that for all k suf ciently large Ic1g a Ic03 where again we can make this as close as we like to an actual equality by choosing k large enough In a formal proof this assumption would need justi cation but we will just make the assumption in our argument Since also lekl m aIek1D we can write k71 ka1D r Now assuming that 214 holds substituting these expressions for 6191 and ek1 into 214 gives aleklu C k ka1 Clekll D a lD r Moving the constant terms to one side this becomes a11a071 g k11a7ozr Since the lefthand side is constant independent of k the right hand side must be as well which means that the exponent of lekI must be zero 1 1 a 7 04 0 Solving for 04 we nd 1 i xE 2 In order for the method to converge 04 must be positive so the order of convergence must be the positive value 04 1 Proof of Lemma Subtracting w from each side of 212 gives fwk k 6k71gt fltwkgt fk717 where we have substituted 6k 7 6191 acc 7 53 7 53194 7 53 for wk 7 1 Combining terms this becomes k1 6k 7 w ack 7 fwk71 fwk k fltk71gt k71 00k i 001971 00k i 001971 ack fltk71gt 7 k1 6196194 215 44 where the last factor inside the brackets in 215 converges to 1f wlt as wk and 194 approach 53 The rst factor can be written in the form fwkfw 7 fwk71fz zk w Ik71 I 9 M71 If we de ne 9w 2 fw 576 216 then the rst factor inside the brackets in 215 is 7 g4ck1acc 7 53194 which converges to g w as wk and 1 approach 53 Differentiating expression 216 we nd that w MV W f aw 9 w win 7 and taking the limit as go 7 xi using L Hopital s rule gives a ltwgt 132 w2wwigt wgt if 2 Thus assuming that 6k and 6191 converge to O as k 7 00 it follows from 215 that lim 6k PM 16700 61961971 2 f w Finally the assumption that 30 and 31 are su iciently close to 3 means that for any C gt f wf w we can assume that am and an lie in an interval about at in which 2 S C61 eo and that say 1S120and 0 12Cgt It follows that 2 S min6160 S 14C so avg also lies in this interval By induction it follows that all iterates wk lie in this interval and that lekl S ek1 S S flgleol Hence 6k and 6191 converge to O as k 7 00 D 25 Analysis of Fixed Point Methods Sometimes instead of writing a nonlinear equation in the form fw 0 one writes it in the form of a xed point problem an Mac that is the problem is to nd a point x that remains xed under the mapping w A natural approach to solving such a problem is to start with an initial guess am and then iterate according to n1 MM 217 Both Newton s method and the constant slope method can be thought of in this way For Newton s method since wk wk 7fwkf wk the function u whose xed point is being sought is 4753 w 7 fwf w Note that 4753 x if and only if fw 0 assuming that f w 7 O For the constant slope method since wk wk 7 fwkg we have 4753 w 7 fwg where g fw0 Note that for this function as well 4753 x if and only if fw O The secant 45 method wk wk 7 7 wk1fwk 7 fwk1 does not t this general pattern since the right hand side is a function of both wk and 194 Figure 24 illustrates cases where xed point iteration converges and diverges In particular the iteration can produce monotone convergence oscillatory convergence periodic behavior or chaoticquasi periodic behavior pm VZHanMxy Xx X as 1 pm 423 lanhx pm 1H3nh127x 42 Figure 24 Fixed Point lteration The iteration may display monotonic convergence upper left oscillatory convergence upper right monotonic divergence lower left or oscillatory divergence lower right Geometrically xed points lie at the intersection of the line 3 w and the curve 3 Mac as seen in Figure 24 Mathematically such points occur when the output of a function is equal to its input 4753 3 Fixed points correspond to equilibria of dynamical systems that is to points where the solution neither grows nor decays Example Find the xed points of 4753 32 7 6 First plot 3 32 7 6 and y w on the same axes to get an idea of where the intersection points are To determine the xed points mathematically solve the quadratic equation 2 7 6 w or 2 7 w 7 6 w 7 3w 2 O The solutions are ac 3 and w 72 What determines whether a xed point iteration will converge The following theorem gives a sa icient condition for convergence Theorem 251 Assume that Lp E C1 and ILpgt lt 1 in some interval about a xed point at If 30 is in this interval then the xed point iteration 217 converges to 53 Proof Expanding Mack in a Taylor series about xi gives n1 WW WW M xvilt62 M 0 mm for some k between wk and 53 Subtracting w from each side this becomes k1 wow51c and taking absolute values on each side gives 6k1 S w kgt leklr Thus if Ivaam lt 1 for all x in an interval about 53 and if 30 is in this interval then future iterates remain in this interval and ek decreases at each step The error reduction factor approaches IW WQII hm k1 I 19700 Iekl w wr Example As an example let us consider three xed point iteration schemes that can be derived from the equation 33 632 7 8 0 which has a unique solution in the interval 1 2 The functions whose xed points we seek are 1 Lp1w w36w2w78 8 2 W2w w6and 8733 3 Lpgltwgt4 6 Simple algebra shows that each of these functions has a xed point in the interval 12 where fw 33 6x2 7 8 has a root for example 47253 x if x is in the interval and 32 8w 6 or 533476902 780 Table 21 shows the behavior of xed point iteration applied to each of these functions with initial guess at 15 Run the MATLAB code fixedpoint m to replicate this table or experiment with different initial guesses n we 7 w my 7 w my 7 w 0 15 15 15 1 10375 1032795559 08779711461 2 1764990234 1066549452 1104779810 3 5516973759 1063999177 1052898680 4 1679 X1029 1064191225 1067142690 5 4734 gtlt1087 1064176759 1063386479 6 1061 X10263 1064177849 1064388114 7 1064177767 1064121800 8 1064177773 1064192663 9 1064177772 1064173811 10 1064178826 11 1064177492 12 1064177847 13 1064177753 14 1064177778 15 1064177771 16 1064177773 17 1064177772 Table 21 Fixed point iteration applied to three different functions derived from the equation 33 632 7 8 0 These results are explained by the previous theorem 1 First consider 47153 33 6x2 w 7 8 a function for which xed point iteration quickly diverged Note that 3m2 1290 1 gt 1 for all x 6 12 so that one must expect divergence 8 2 Next consider 2 w In this case 90 x2 w6 32 lt 1 for alle 12 7 6 W2 Since 30 15 6 12 convergence is guaranteed 48 3 Proving the convergence of xed point iteration for 473cc cc is left as an exercise for the reader Actually omen need not exist in order to have convergence The iteration 217 will converge to a xed point cc if g is a contraction that is if there exists a constant L lt 1 such that for all cc and 3 My 7mm S Lw7yr 218 Theorem 252 If Lp is a contraction then it has a unique xed point cc and the iteration wk1 Lpck converges to ac from any avg Proof We will show that the sequence 6210 is a Cauchy sequence and hence converges To see this let 9 and 16 be positive integers with 16 gt 9 and use the triangle inequality to write lack 7 ch S lack 7 wk1wk1 7 wkigl cj1 7 ml 219 Each difference lwm7 5cm1 can be written as o7cm1 7 Lpwm2 which is bounded by LIwm1 7 5cm2 where L lt 1 is the constant in 218 Repeating this argument for layml 7 wmql we nd that Iwm7wm71 S L2wm72 7wm73l and continuing in this way gives Iwm7wm71 S Lm llwl 7w0 Making these substitutions in 219 we nd Iwk 7 mil 3 L1H Lk 2 L7w1 7 wol L7 L 7i w17w0 If 16 and j are both greater than or equal to some positive integer N we will have 1 17L lwlrwjl SLN w17w07 and this quantity goes to O as N 7 00 This shows that the sequence 5 16 01 is a Cauchy sequence and hence converges To show that wk converges to a xed point of ya we rst note that the fact that o is a contraction implies that it is continuous Hence o7lim1CHOO wk limCHOO Letting cc denote limCHOO 5 we have Mai E M lim 5 lim Mack lim wk 4 16700 16700 16700 so that the limit of the sequence wk is indeed a xed point Finally if y is also a xed point then since o is a contraction we must have 155 341 117550 31013 L1 341 where L lt 1 This can hold only if 31 4 so the xed point is unique D 49 JULIA OF THE JULIA SET GASTON JULIA 189371978 7 was only 25 years of age when he pub lished his 199 page Memolre sur l lt mtlon des fonetlons rationelles The publication led to fame in the 1920 s but was virtually forgotten until Benoit Mandelbrot s study of fractals The leather strap seen cov ering a portion of Julia s face was due to a severe injury in World War 1 that led to the loss of his nose 26 Fractals Julia Sets and Mandelbrot Sets Fixed point iteration and Newton s method can be used for problems de ned in the complex plane as well Consider for example the problem of nding a xed point of Lpz E 22 It is easy to predict the behavior of the xed point iteration zk1 1f zo lt 1 then the sequence zk converges to the xed point 0 1f zo gt 1 then the iterates grow in modulus and the method diverges lf zD 1 then 1 for all k If 20 1 or 20 71 the sequence quickly settles down to the xed point 1 On the other hand if say 20 egmg then 21 647W3 and 22 687W3 627W3 20 so the cycle repeats It can be shown that if 20 62 where 04 is irrational then the sequence of points zk never repeats but becomes dense on the unit circle The sequence of points 20 21 Lpzo zg LpLpzD is called the orbit of 20 under Lp lf Lpz is a polynomial function then the set of points 20 for which the orbit remains bounded is called the lled Julia set for Lp and its boundary is called the Julia set Thus the lled Julia set for 22 is the closed unit disk while the Julia set is the unit circle These sets are named after the French mathematician Gaston Julia In 1918 he and Pierre Fatou investigated the behavior of these sets which can be far more interesting than the simple example presented above Their work received renewed attention in the 1980 s when the advent of computer graphics made numerical experimentation and visualization of the results easy and fun Figure 25 displays the lled Julia set for the function Lpz 22 7 125 The boundary of this set is a fractal meaning that it has dimension neither one nor two but some fraction in between The gure is also selfsimilar meaning that if one repeatedly zooms in on a small subset of the picture the pattern still looks like the original The plot in Figure 25 was obtained by starting with many different points 20 throughout a rectangular region containing the two xed points of Lp and running the xed point iteration until one of three things occurred Either Z 2 in which case we conclude that the orbit is unbounded and 20 is not in the set or zk comes within 10 6 of a xed point and stays within that distance for 5 iterations in which case we conclude that the iterates are converging to the xed point and 20 is in the set or the number of iterations k reaches 100 before either of the previous two conditions 50 Fwnxm n w Emu HEB mmn men 6an a H mm cans 5 saint nwmm Sm nounr m 2 sum mmasmdam n a son n05ltmxw w 3 w mama v050 ran 92 9835 rosdmmm Eat rmdnm S m wwmsp 5 sum mmquot Zoom 393 Edam 433 W 7 H nw m I W P 99 F5 W A H mm m a F51 W n am H mm m P mnn m0 393 mam masmdnm amt drosdmmm rm 3 no you nm aim Eon m ESE wsxm n a gram or BEE n on m mdwrumm m n mmnmluim w m n 5929 max sum 01v 3 rm 57 3m a E mmsnwnmm mammm E a v 5 max saint sum 01v 5 re amt WEm a xwvmm 325 IE an n 5 30m Fm max sum mumquot m E r 0 9 mac drosdmmmuw 0x 393 sum wvv 0 m 355x 0 05m sum mama voinm 9 399 may no 9ltmxwm ltlt macaw saxoswr minrmx Eitwma 05 n o mmdnm 5 n 938 5355 3 355m 98 o a v3 0 m Hum pm on wmmsiiw n w mmasmunm m m m a 0 3 9835 Hmmm an n 30 95 mo 0 Sman 9 05m 35 Va 0 mano 0 8 95 a manidnw rm nomm 05m Em mo 3 w H03 592 m3 mo vmstouam Hm 535er aro mm on m mum win a 33m xmmi 33 mswwmm a sum 33 mm mmquot m nexxmnn iii 9 m n vwxm amnmu are m HamE mmxmdn mmnm sums 05m 33 nounusmm 93 n 5 mi 5 nm 3 m v Onrmx 33m 3 m mam 356539 5 sum noHanxiw SE ma mm nm 5 sum rvmanim Eat Svgiw sum 98v in xmmEnm 45 53 5 H90 Qatari wrosn 39mm mmquot Him 9383535 on 55man xmmEnm w m n 3 29 36333 5 70 v Evvrmm wormgwnnm 05m 5 sum 503 0 armnw Hiding min w w m mmm E xmmEnm 0 w mo 53 Hyman HEB m on H 05m 535 HW om x m m 5 no on gum 358 A m 31 5 0 735 w H H 95555 on 555 tnmm on n E Emm 3 55 mmquot on n quot5035quot mam 59833quot E35 Bum KSBBSW munmm v050 no a dam 00555qu Eam mx Swm Sam gt552on may mmv vmm 5 ywnxm n a 159m Ewnw voinm wH

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### 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'

## Why people love StudySoup

#### "Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."

#### "I used the money I made selling my notes & study guides to pay for spring break in Olympia, Washington...which was Sweet!"

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.