Newton's Method for Finding Roots
Objective: This module shows how to implement Newton’s Method using a
Maple procedure
in order to find a real root of an equation f(x) = 0. This iteration method
leads naturally
to a do loop that can terminate when we are close to the root. This loop is
placed inside the
procedure. The Maple commands in this module are
plot, abs, proc
The geometric idea behind Newton’s method is that the x-intercept of a tangent
line to
the function f(x) is close to the root of the equation f(x) = 0.
The method begins with
a guess, x0, that is based on knowing where f is positive and negative. For
example, if
f(x) = x3 + 2x2 − 3, then the Intermediate Value Theorem tells us that a root is
between
0 and 2, since f(0) = −3 < 0 and f(2) = 13 > 0. Suppose x0 = 0.5. Then we
construct
the tangent line at x0 and find its x-intercept, x1, as the next guess at the
root. The method
continues until the guesses converge (that is, show no change) to the root.
Recall that the tangent line at x0 has slope f'(x0). The tangent line is then y
= f'(x0) +
D(f)(x0)(x − x0) where we make use of the D(f) command. We plot the function
together
with the tangent line to show that the intercept of the line comes close to the
root of f(x):
> f:=x->xˆ3+2*xˆ2-3;
> x[0]:=0.5;
> tanline[0]:=x->f(x[0])+D(f)(x[0])*(x-x[0]);
> plot({f(x),tanline[0](x)},x=0..2,y=-5..5);
We can compute the intercept of the tangent line by solving for x:
f(x0) + D(f)(x0)(x − x0) = 0 implies that x − x0 = −f(x0)/D(f)(x0), or
x = x0 − f(x0)/D(f)(x0)
We can call this x value tanzero(x[0]). From the formula above, we see that
tanzero
can be written as a function:
> tanzero:= x->evalf(x-f(x)/D(f)(x));
Note that the Newton’s method approximation is put inside a evalf(). The purpose
of
evalf() is to force Maple to evaluate the answer in the decimal form, instead of
an algebraic
expression. Moreover, if we let x[1]:=tanzero(x[0]), we can see that the tangent
line
at x1 has an intercept that is a better approximation to the root than that at
x0. To see this,
execute the following Maple commands:
> x[1]:=tanzero(x[0]);
> tanline[1]:=x->f(x[1])+D(f)(x[1])*(x-x[1]);
> plot({f(x),tanline[0](x),tanline[1](x)},x=0..2,y=-5..5);
4.1 Execute the commands x[1]:=tanzero(x[0]); x[2]:=tanzero(x[1]); repeatedly
until the output converges to 1. For each execution of the command, what is the
number of decimal places that are correct in each approximation?
4.2 Implement the sequence of commands in Exercise 4.1 using a for loop. Before
the loop
begins, define the initial guess using the label x[0]:=0.3;.
4.3 Terminate the loop in Exercise 4.2 if two successive approximations are very
close; that is,
if abs(x[i-1]-x[i])<0.001. Use the initial guess x[0]:=0.2; Your loop should
have enough iterations that the if .. then .. statement was for sure invoked.
The proc()...end; Command
Another programming feature, similar to the for...do...od; and the if...then...fi;
command structures, is the Maple procedure. A procedure is a group of commands
that performs
a particular task. The proc command typically has input variables (in order to
complete
a task) and returns one or more answers (depending on the task). For example,
the following is
a procedure that computes x2 + 3x, given an input value x:
> f:=proc(x)
local output;
output:=xˆ2+3*x;
RETURN(output);
end;
Notice that the proc ... end command begins with stating the function name (f),
and we
define f as a procedure. Inside the procedure, we calculate the output value
using the variable
output, and then RETURN the output value. To execute this procedure, say, we
would like to
evaluate f(2), simply type
> f(2);
Within the procedure, you need to define any labels you use as local. In the
previous
example, the variable output is a local variable inside the procedure only.
Notice that we can also define f(x) by f:=x->x∧2+3*x. In fact, Maple turns this
definition
to a procedure, as shown by the output (if you turn the output into Maple
notation mode).
Procedures are, however, a lot more powerful in performing complex tasks. For
example, the
following is a procedure that computes the first and second derivatives of a
function f. Notice
that f has to be defined as an expression:
> Derivs:=proc(f)
local df,ddf;
df:=diff(f,x);
ddf:=diff(diff(f,x), x);
RETURN(df,ddf);
end;
To execute this procedure, say, we want to find the first and second derivative
of x3 + 2x − 5,
the user simply enters
> Derivs(xˆ3+2*x-5);
4.4 Implement the above procedure Derivs in your worksheet. Explain what the
code for
Derivs did in each line.
Here’s another example that determines if a quadratic function f(x) = ax2+bx+c
has two real
roots, one real root, or two complex roots. The procedure uses the extension of
the if ...
then ... else statement to if ... then ... elif ... else. Look up
the structure in the Maple help section.
> RootType:=proc(a,b,c)
local discrim;
discrim:=bˆ2-4*a*c;
if discrim > 0 then RETURN("two real roots")
elif discrim < 0 then RETURN("two complex roots")
else RETURN("one real root")
fi;
end;
4.5 Enter the above procedure into your worksheet. Execute the command
RootType(1,3,2).
What does this command do? Which polynomial is it checking?
4.6 Write a procedure called ispositive where the procedure returns “true” if
the input
value is positive, otherwise, it returns “false”. Test the accuracy of your
procedure by
evaluating ispositive(1), ispositive(-2), ispositive(0). (Hint: start with
ispositive:=proc(x) )
4.7 Incorporate Newton’s method that we developed in 4.2 into a procedure newton(f, s,
n) where f is the input function in terms of x; s represents the initial guess
x[0]; and n is the
number of iterations that you want to carry out in the loop. Test your procedure
by evaluating
newton(x->x∧3+2*x∧2-3, 0.2, 3), newton(x->x∧3+2*x∧2-3, 5, 10).
|
