# Lesson 3: Newton's Method and other Iterations # ================================= # # One of the most important numerical methods of finding # roots of functions is # Newton's Method. Maple's "fsolve" uses essentially this # method (with a few # added wrinkles). We'll look into how this works. # # First, a bit of a detour: defining functions in Maple. # # The simplest way to define a function in Maple is by a # mapping rule such as # > f:= x -> x^3 - 2*x^2 - 2; 3 2 f := x -> x - 2 x - 2 # # You can then evaluate the function at a certain value of the # independent # variable, putting the independent variable in parentheses. # > f(3), f(p); 3 2 7, p - 2 p - 2 # Note the difference between the function f and the expression # you'd get by > g:= x^3 - 2*x^2 - 2; 3 2 g := x - 2 x - 2 # The expression is tied to the particular name x for the # independent # variable. If you give x a certain value, the value of g changes # to # match it. On the other hand, the x in the definition of f is just # a "dummy variable'', a name that stands for "the independent # variable # of the function f''. The definition of f is not affected by any # values you assign to this name. > x:= 5: > f(t); 3 2 t - 2 t - 2 > g; 73 # Note: don't fall into the error of defining f(x) instead of f. # Maple won't complain, and f(x) will be what you defined it # as, but it won't use this formula for any other inputs. > x:= 'x'; x := x # (This undoes the assignment of 5 to x) > badf(x):= x^3 - 2*x^2 - 2; 3 2 badf(x) := x - 2 x - 2 > badf(y), badf(3); badf(y), badf(3) # You can also make an expression into a function with the # "unapply" command. If x is a variable and expr is an # expression depending on x, then unapply(expr, x) is the # function taking x to expr. # It's particularly useful when you have obtained a complicated # expression as a result of some complicated process and don't # want to type it in again. > h:= unapply(g, x); 3 2 h := x -> x - 2 x - 2 # # You can differentiate a function with the "D" operator, # obtaining another function. > D(f); 2 x -> 3 x - 4 x > D(f)(y); 2 3 y - 4 y # Newton's Method # ============ # # Newton's method is the following: # Suppose we want a solution to an equation f(x) = 0. # We start with an initial guess x0, and Newton's method # produces a sequence # x1, x2, ... that (when it works well) converges very rapidly to # a solution of the equation. As we will see, however, the # method does not always work well. The next item in the # sequence is given by the formula > x[n+1] = x[n] - f(x[n])/D(f)(x[n]); 3 2 x[n] - 2 x[n] - 2 x[n + 1] = x[n] - ------------------- 2 3 x[n] - 4 x[n] # It's convenient to write this as x[n+1] = newt(x[n]) where # "newt" is a function. Note the use of "evalf": we don't want # the exact answer (which could be # a fraction with large numerator and denominator if x is a # rational # number). > newt:= x -> evalf(x - f(x)/D(f)(x)); f(x) newt := x -> evalf(x - -------) D(f)(x) # Let's start with a typical guess, and see what happens in a # few iterations of the method. > x[0]:= 2; x[0] := 2 # This is the first entry in what Maple calls a "table''. If # you've done # some computer programming, you may be pleasantly # surprised to see that # you don't need to tell Maple that you're making a table, or # how big it # will be: the table is created automatically whenever you use # "[...]'' following a name, and will hold whatever entries you # put in it. The index inside the brackets is shown as a # subscript in output. > x[1]:= newt(x[0]); x[1] := 2.500000000 # Instead of calculating successive entries one-by-one this # way, we can use a "do loop" to get several at once. # # The do loop has a side effect: it leaves the index variable with # a value (one more # than the last value in the loop). This can cause trouble if you # try to use the same # variable as a symbolic variable later. To avoid this, I've # adopted the practice of always using double letters ii, jj, nn, # etc. for index variables of loops. > for nn from 2 to 6 do x[nn]:= newt(x[nn-1]) od; x[2] := 2.371428571 x[3] := 2.359405644 x[4] := 2.359304093 x[5] := 2.359304086 x[6] := 2.359304086 # Well, Maple seems to have settled on a value. Is this the # solution? > f(x[6]); 0 # # This was an example where Newton's method worked very # well. # But with other starting points we might not be so lucky. # To save the x values we just calculated, I'll use a different # name for the next table. > w[0]:= 1; w[0] := 1 > for nn from 1 to 6 do w[nn]:= newt(w[nn-1]) od; w[1] := -2. w[2] := -1.100000000 w[3] := -.3838107098 w[4] := .8053409552 w[5] := -1.369897554 w[6] := -.6206240787 # No sign of settling down yet. Actually it turns out that this # one will eventually # approach the root as well. # # Let's see how Newton's method works, using some graphs. # The idea is # to approximate the graph of f(x) for x near x[n] by its tangent # line at x[n]. # The next guess x[n+1] is where this tangent line hits the x # axis. # # We'll combine the graph of f(x) with the tangent line from # [x[n], f(x[n])] to # [x[n+1], 0]. As an added touch, we can plot the straight lines # from [x[n], 0] # to [x[n], f(x[n])] and from [x[n+1], 0] to [x[n+1], f(x[n+1])]. # # Besides plotting an expression as we've used it before, the # "plot" command # can take a list of points and join them with straight lines. A # list in Maple is # enclosed in square brackets, with the items separated by # commas. Each point # itself is a list (consisting of the x and y coordinates). It's # convenient # to define a function that takes x[n] as input and plots the # straight lines we # want. # > newtlines:= x -> plot([[x,0],[x,f(x)], [newt(x),0], [newt(x), f(newt(x))]]); newtlines := x -> plot([[x, 0], [x, f(x)], [newt(x), 0], [newt(x), f(newt(x))]]) # # We can combine the outputs of several plot commands in one # picture using # the "display" command. This is contained in a package called # "plots", together # with many other useful plotting commands. In order to use # it, we first have to # load the package, as follows: # > with(plots): # The display command itself can take as input a set of plot # outputs. A set in Maple is enclosed in braces {} (the main # difference between a set and a list is that order counts in a list # but not in a set). > display({ newtlines(x[0]), newtlines(x[1]),plot(f(x), x=1.9 .. 2.6) });\ ** Maple V Graphics ** # We could do more steps, but the next step from x[2] to x[3] # would be # practically invisible on this graph. Let's try the one that # didn't seem to converge. > display({ newtlines(w[0]), newtlines(w[1]), > newtlines(w[2]), newtlines(w[3]), > plot(f(x), x= -2.5 .. 1.5)});\ ** Maple V Graphics ** # Why the big difference between the two cases? The tangent # line at # x=x[0](=2), which of course is a very good approximation to # f(x) near # x = 2, is still not too bad an approximation near the root r = # 2.359304086. # It intersects the x axis at a point x[1] even closer to r, and the # tangent line there is an even better approximation to f(x) near # x=r. On the other hand, the tangent line at x=w[0](=1), # although a good approximation to f(x) near x=1, bears no # resemblance to f(x) near x=r. Because of this there is no # reason to expect # w[1] to be close to r, and therefore no reason to expect # Newton's method starting at w[0] to converge. In fact, the # slope of the tangent line at w[0] # is in the wrong direction, so that w[1] is even farther from r # than w[0] is. The w[n]'s jump around, and might never # come near r.