Numerical Methods, MAT 350
Autumn 2001

Homework: Roots of Nonlinear Equations

Due Friday, September 14, 2000 in class

This homework will be more challenging than the first homework set. You should start working on it as soon as possible.

In this assignment, we will investigate the performance of several root finding routines acting on a number of test functions. We will examine the following test functions:

  1. In my research, I must find the root of the following function on the interval [0, 1]:
    0.8 -Pi x cot(Pi x)
  2. Here is a polynomial on the interval [0, 1]:
    - 64 + 336 x - 588 x^2 + 343 x^3
  3. This example is from solid state physics. This function describes the occupation of levels in a Fermi gas as a function of energy. The interval is [0,1].
    0.5 - 1/(exp(60 - 100 x) + 1)
  4. An oscillating function on the interval [0,1]. Newton's method had better watch out!
    1.2 - 3 x + 0.3 sin(24 x)
  1. In the file that will also contain your main program, write four subroutines containing the four test functions. First, I give the FORTRAN examples of the functions:
    1. pi = 3.1415926535897932D0
      f1 = 0.8D0 - Pi*x/tan(Pi*x)
    2. f2 = -64 + 336*x - 588*x**2 + 343*x**3
    3. f3 = 0.5D0 - 1.0D0/(1.0D0 + exp(60.0D0 - 100.0D0*x))
    4. f4 = 1.2D0 - 3.0D0*x + 0.3D0*sin(24.0D0*x)
    and then examples in C:
    1. pi = 3.1415926535897932;
      return 0.8 - pi*x/tan(pi*x);
    2. return - 64 + 336*x - 588*pow(x,2) + 343*pow(x,3);
    3. return 0.5 - 1.0/(1.0 + exp(60.0 - 100.0*x));
    4. return 1.2 - 3.0*x + 0.3*sin(24.0*x);
    You should use double precision for all floating point variables and constants. In your main program, if you are using FORTRAN, you will have to declare the names of the four test functions to be EXTERNAL. In C, you will have to include prototypes of your four test functions at the beginning of the file. In order to use the various math functions in C, you will also have to include the math header file:
    #include <math.h>
    at the beginning of the file.
  2. If you will be using Newton's method (see below), your subroutines will also have to return the first derivative. An example of how one might do this in FORTRAN is:
    C  This function returns 0.5 - 1/(exp(60 - 100 x) + 1) and its derivative
          double precision function f3(x,fprime)
          implicit none
          double precision x,fprime,temp
          temp =  exp(60.0D0 - 100.0D0*x)
    C  This is the derivative of the function:
          fprime = (-100.0D0*temp)/(1.0D0 + temp)**2
    C  This is the function itself:
          f3 = 0.5D0 - 1.0D0/(1.0D0 + temp)
    In C, the derivative fprime would be passed to the subroutine as a pointer.
  3. For each of the four test functions, explain why it would --- or would not --- be difficult to find the root: Is the function approximately linear near the root? Is the function oscillatory?
  4. Write a function that finds a root of a given function. Initial of students assigned to each method:
  5. Explain how your root-finding routine works. That is, start with a paragraph explaining, step by step, how the algorithm works. Then explain, step by step, the subroutine you wrote.
  6. The root-finding routine should be saved in a file separate from your main program. It should have the following arguments: If you are doing Newton's method, you only need to supply one initial value.
  7. Using the ends of the interval as starting values, calculate the root of each of the four test functions to an accuracy of 5.0E-16. Print out your answer showing 16 digits past the decimal point. In Fortran, you do this by adding a format statement:
     11   format('my results for f1=',F20.16)
    In C, you do this using:
      printf("my result for f1=%.16g\n",y);
    (I have no idea how you do this in C++.) Which cases converged?
  8. If you are using Newton's method, you must specify one starting point. Use the left hand side of the interval as your starting point. If the routine does not converge using this starting point, try to find a starting value that does converge.
  9. On the Internet, find a root-finding routine based on Brent's Algorithm. Use this routine to find the root for all your test functions. Is there agreement for the two methods in the cases that converged?
  10. Let delta be the difference between the exact value for the root, determined above, and the value after n function calls. Plot log10(abs(delta)) as a function of the number of function calls n (log10 is logarithm base 10). You should be ready to explain these results.
  11. You should be prepared to make a five minute presentation of your results to the class on the due date. If you will not attend this class, you must tell me as soon as possible.