Homework 3. Assigned: Thu Feb 8. Due: Tue Feb 20 at midnight.

1. 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.
   
   
2. Implement a MATLAB function findzero for finding a root of a continuous function f on an interval [a,b], with the following requirements:
   • 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 right-hand side
   • 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 regula falsi.
   • 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
   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.

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