Getting Started
Create a new script (m-file) in MATLAB and add a clear command to the top of your script.
Problem 1: Use Newton's Method to find the roots of f(x) = x2 − 2 = 0
Newton's formula is used to find the next approximation of the root of an equation. Had Newton only invented this method during his career, he would have been famous. Instead he also invented calculus, classical mechanics, and many other things. Newton's formula is:
Use MATLAB to find successive numeric approximations for the root of the equation x2 − 2 = 0 by implementing the algorithm given below. As you do this:
An algorithm for using Newton's formula find a root of a function to within ε = 1.0 × 10−15
Save initial guess (x_curr)
Compute f(x_curr)
while abs(f(x_curr)-0) is more than epsilon
Compute f'(x_curr)
Use Newton's method to compute next approximation (x_next)
Copy x_next to x_curr
Recompute f(x_curr)
(Note: an algorithm like the above would not be given to you on an exam. We would expect that you could develop such an algorithm given a successive numeric approximation formula.)
How long (how many approximations) does it take to reach the limit of MATLAB's precision?
What happens when you try different initial guess values? (Do you find the same or a different root or does it not converge?) Try several different initial guesses of varying distances from the roots.
How could you modify your code so that it found where x2 − 2 = R (where R is a value input by the user)?
Problem 2: Use Newton's Method to find the roots of this more complex equation.
In your script:
Problem 3: Order of Accuracy — Simpson's Method
We will now program Simpson's Method for approximating the definite integral of functions. Simpson's Method uses the concept that the integral is the area bounded by the function and the x-axis and the integral limits on the left and right. Consider the function
and the integral 
Plot the function y = f(x) over x = 0 to 2π to get a feeling for the area the integral is computing.
We have written an implementation of Simpson's Method for you. We have chosen to break Simpson's Method into two distinct functional components to help you. Of course, there are other ways to implement Simpson's method for computing definite integrals. Download simpson.m to your current working directory.
Notice that there are four input parameters to the simpson function: a function handle representing the function to be integrated, the lower limit of the integral, the upper limit of the integral, and the number of subintervals to approximate the solution. As required by Simpson's Method, the number of subintervals must be an even integer.
The result should be equal to 4π. Next integrate the same function over the same limits using 8 subintervals.
You should get the same result. Canceling errors allow Simpson's Method to exactly compute this periodic function between periodic limits.
Repeat to set s(2) equal to the results for N=6 and again to set s(3) for N=8 and again to set s(4) for N=12 and once more to set s(5) for N=16.
Does the value of the integral appear to be converging to a single value? If so, how quickly?
To get a really good estimate for the exact value for the integral of 2 + sin(x) between 0 and 3π/2, call the simpson function using 64 subintervals and name this result exact_s.
What does semilogy do?
How much more accurate is the approximation for 16 subintervals than the approximation for 4 subintervals? Is this consistent with the stated accuracy of O(h4) ? Hint: Assume exact_s is the exact value for your computations.
Which line(s) of your script would you have to change to compute a different integral?
Problem 4: Turning Newton's Method into a Function
We now want to turn the Newton's Method algorithm you implemented in Problems 1 and 2 into a more general function (like simpson from Problem 3 or fzero) that will allow a user to specify the function whose roots they are interested in.
Create a new function in a file named newton.m. The beginning of the function should be:
function [root, iters] = newton(fctn, f_prime, x0) % newton % % Uses Newton's method to find a root of fctn(x) near x0. % % Newton's method is used to find a root for fctn(x) starting with the % initial guess, x0. Successive approximations are done until either % fctn(x) is within epsilon of 0 (where epsilon = 10^-15) % or % 100 iterations have been done (i.e., 100 new guesses have been % computed) % % Input: % fctn a function handle for function whose root is being determined % f_prime a function handle for function representing the first % derivative of fctn(x) % x0 the initial guess for the location of the root % % Output: [root, iters] where % root the final estimate for the root % iters the number of iterations (until either fctn(x) is within % epsilon of 0 or until 100 iterations have been done) %
Using your code from Problems 1 and 2 as a guide, implement the newton function so that it behaves as described in its function comment.
Test your newton function using the function from Problem 1 (i.e., find the roots of x2 − 2 = 0) Note that in order to be able to hang on to both values returned, you'll need to use code like:
[ xvalue, numIters ] = newton(func, func_prime, 0)
Next, test your newton function using the function from Problem 2.
Now let's investigate how the initial guess affects whether Newton's Method converges, which root it converges upon, and how many iterations it takes:
Discuss your results with your team. Which roots did you find? Did Newton's always converge? How did the number of iterations compare to how close the guess was to the root it found?