Homework 3. Assigned: Thu Feb 8. Due: Tue Feb 20 at midnight.
Email your m-file implenting findzero as an attachment directly to the grader,
email@example.com, with a cc to me. Your grade will depend on how completely findzero
satisfies the requirements, so be sure you test it thoroughly: the grader will
test it on a variety of input. Answer question 1 in plain text in the body of the email or, if you prefer, in a separate attachment.
- Show that the square root algorithm you implemented in Homework 1 is actually Newton's
method: what is the function f(x) to which, when you apply Newton's method, you
get the iteration that you implemented? First write down a function f(x)
whose root x*, satisfying f(x*)=0, is the square root of 2,
and figure out the Newton iteration for this function: do you get the same
thing you implemented in homework 1? Show the details.
More generally define such an f that works for computing the square root of any number, say a.
Now recall the definition of quadratic convergence
on p.55 of Ascher-Greif, and recall that in class we showed that the
quadratic convergence constant M
for Newton's method is f''(x*)/(2 f'(x*)). To answer the question:
is the square root iteration
always quadratically convergent no matter what a you apply it to,
you need to find the value of M, and make sure it is finite --
in particular make sure f'(x*) is not 0: is this the case for all
a? By the way, this method is used to implement
the square root function in practice.
- Implement a MATLAB function findzero for finding a root of a continuous
function f on an interval [a,b], with the following requirements:
Test your program on lots of examples, for example constructing functions from combining
basic math functions such as exp(x)+sin(x) or log(x)+cos(x). Graph the
functions using plot so you have some idea what input values to use for a and b.
Compare the results you get with findzero to the built-in function fzero,
but bear in mind that if f has more than one root in
[a,b] the two routines may return different answers.
- it has 4 input arguments: an "autonomous function" f, and scalars
a, b and tol, and two output arguments, x and iters;
thus the first line of the m-file should be
function [x, iters] = findzero(f, a, b, tol)
- it should be invoked with a call something like [x, iters] = findzero(f, a, b, tol), where
the variable f has been assigned with a statement
something like f = @(x)fun(x), where fun is the
name of the MATLAB function that evaluates f(x) --- unfortunately
this only works in MATLAB Version 7, so if you are using an earlier version
(type "version" to find out), simply set f = 'fun' (a string) and,
inside findroot, call feval to evaluate f (type
help feval to learn how).
- it should set x to nan (or generate an error return) if
f(a) and f(b) are both positive, or both negative, or if b &le a
- otherwise it should implement a hybrid of the bisection method and the secant method that
attempts to return x with abs(f(x)) &le tol.
- the output iters should be the number of times f is called
- it should not go into an infinite loop if tol is too small. Instead, it should
terminate when machine precision has been reached; you can try the test
(b-a)/max(abs(a),abs(b)) &le eps, where eps is
the built-in machine epsilon, but you might need a factor of 2 on the
- one possible implementation is to always set the new point to the x-intercept of the
line passing through (a,f(a)) and (b,f(b)). This is the same as the secant
step using these 2 points. However, this will sometimes be as slow, or maybe even slower,
than bisection, measuring in terms of iters. This method is usually called
- a better idea is always use the 2 most recent points
to define the secant step, rejecting it if the x-intercept is outside the current interval [a,b], and taking a bisection step instead.
The most recent point will always be either a or b, but the second most recent
point might not the other of a or b. Here
a and b always bracket the root as in bisection.
This should take less iterations on some examples. But if you can't get this to work
reliably, just use regula falsi.
- it should not use derivatives (therefore it can't use Newton's method)
- it should not fail as long as f is continuous and it should generally be faster
(measured in terms of iters) than simple bisection as long as f is differentiable.
- it must have comments that clearly explain how it works
Note: the reason we care about iters is that in real applications, the evaluation
of f might take a long time.
Reading: Chapters 3-4 of Ascher-Greif (skip Section 3.2).