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