Numerical Methods - Assignment
Robustness of numerical methods
Problem.
Devise a robust and efficient scheme for finding all simple roots of a given (real nonlinear scalar) equation f (x) = 0 in a given interval x ∈ [a, b].
Background. Important qualities of numerical methods, in addition to convergence, are
(a) efficiency - requires a small number of function evaluations;
(b) robustness - fails rarely, if ever;
(c) minimality - uses a minimal amount of user input, e.g. does not require the derivative in addition to the function.
No single method meets all criteria. For instance, the Newton-Raphson map g[f ](x) is guaranteed to converge quickly if the initial iterate x0 is "close enough" to a simple root so that gxj [f](x0) < 1.
However, convergence is only local and the neighbourhood of convergence is difficult to assess ahead of time for a given f(x). In contrast, the bisection method is slow but it is guaranteed to converge provided only that a bracketing interval is found on which f(x) changes sign. A natural idea, then, is to combine the two in a fast and robust hybrid method.
Question 1. Write a MATLAB function [ai,ci,bi] = bisect(f, ai, bi, niter) to per- form the interval bisection method for a given function f on an interval given by its left and right ends ai and bi. The MATLAB function must check if f (x) changes sign in the interval and if so perform a specified number of iterations niter and then return the estimated root in ci and the left and right ends of the bracket where it is located in ai and bi; else return in ci a NaN MATLAB object.
Question 2. Write a MATLAB function ci = newtonraphson(f, x0, tol) to perform the Newton-Raphson method for the function f starting from initial iterate x0. The function must check if the absolute value of the derivative of the Newton-Raphson map at x0 is less than unity, and if so perform iterations and return in ci the estimated root; else return in ci a NaN value. To avoid non-essential input all derivatives needed must be computed using the finite-difference formula gprim = (g(x + eps)-g(x -eps))/(2*eps), where eps=1e-5. Use xk xk-1 < tol(1 + xk ) as a stopping condition for iterations; what is your interpretation of this condition?
Question 3. Write a MATLAB function rts=findroots(f,a,b,n,tol) to combine the bisection and the Newton-Raphson methods as follows. Uniformly partition the given inter- val [a, b] into n intervals of equal length. For each subinterval [ai, bi] over which the function changes sign use 3 bisection steps to approach the root more closely calling your bisect. Next, call your newtonraphson to perform Newton-Raphson iterations if the method con- verges, otherwise revert back to the last bisection bracket and apply 3 further bisection steps. Loop until the stopping condition f (xk) < tol is met. Return in the array rts all roots found.
Question 4. Write a MATLAB script Q04.m where the function findroots is used to find all of the roots of the Chebyshev polinomials Tn(x) of degrees n from 1 to 20 with tolerance tol=1e-12. In one figure, plot all roots found on the abscissa axis versus the polynomial degree on the ordinate axis. In a second figure, plot together the graphs of all Chebyshev polinomials Tn(x). For the polynomial T15 compare the roots found numerically with the known exact solutions and print the absolute error in each root in the format of the following example.
Roots and their Errors for the Chebyshev polinomial T[5]:
-9.51056516295153531182e-01 +0.00000000000000000000e+00
-5.87785252292473137103e-01 +0.00000000000000000000e+00
+3.21265594383589607162e-17 +9.33588993957266210749e-17
+5.87785252292473137103e-01 +1.11022302462515654042e-16
+9.51056516295153531182e-01 +0.00000000000000000000e+00
[Hint: you may use fprintf with format +30.20e'.]
Further hints and directions.
1. Since this Assignment is largely a programming exercise, it is essential to include all MATLAB codes in the pdf files of your report. Explain the steps your codes perform, list the variables used.
2. MATLAB codes for the interval bisection method and the Newton-Raphson method appear in our MATLAB lecture notes and in past Assignments. Feel free to modify these.
3. Do not use MATLAB inbuilt root-finding functions such as roots and fzero.
4. The definition and the exact roots of the Chebyshev polinomials are given in section 2.3.2 of the Lecture Notes.
5. Formatted output can be printed to screen or file e.g. by fprintf(' +30.20e +30.20e \n', [array1; array2]).
6. Write code in a"bottom-up" approach - implement the building blocks first; then test and verify for one equation; then loop over many equations.
7. While writing code it is useful to include MATLAB statements to compute and print out errors, orders of convergence and otherwise monitor what you code does.