MATLAB Mini-Project 4

Solving Nonlinear Equations

When the astronomer Johannes Kepler was attempting to determine the orbits of planets, he had measured quantities he called M and ALPHA, and he needed to solve the following equation to determine the quantity X: 

        M = X - ALPHA * sin X
At that time, there was no way known to solve this equation exactly; nor did anyone know a way to systematically determine a good approximation. That left a long and painful trial and error process.

In this project, we'll learn a little MATLAB, and also look at one way to determine an approximate answer to Kepler's problem, or, indeed, to almost any equation.

To solve the equation means to find a value of X that makes it true. Let's move all the terms to the right hand side, and call the result "F(X)": 
        F(X) = X - ALPHA * sin ( X ) - M
Now our problem has become finding a value X for which F(X)=0.

Rewriting the equation as a problem involving a function F(X) helps a lot. It may be hard to think about solving an equation, but it's easy to imagine the graph of the function F(X) as a curve, and to see that F(X)=0 when the curve crosses the X axis. We can also see that the sign of F(X), the size of F(X), and the derivative of F(X) might be useful pieces of information.

Exercise A: Write a MATLAB M file kepler.m which accepts a value for X, and returns the value of F(X). Suppose that M = 1 and ALPHA = 2. Your file should have the form: 

        function f = kepler ( x )
        f = ???
and you'll need to fill in the formula for F. Check that your function is working correctly by typing 
        kepler ( 0 )
which should give you a value of -1.

Get a feeling for what this function looks like, by plotting it between -10 and 10. If you look at this plot and at the formula, you should be able to guess an answer to the following questions:

MATLAB Syntax: you might find it very helpful to have the X axis displayed on your plots. To do so, you can use the LINE(X,Y) command, which plots a line connecting the points (X1,Y1) to (X2,Y2). If my X values go from -5 to 5, then I might use the following commands: 

        plot ( x, y )
        hold on
        line ( [-5,5], [0,0] )
        hold off

Well, now we've got our function F(X) set up so that MATLAB can evaluate it, and we have to face the more serious problem of figuring out how to locate a root of F(X), and how to use MATLAB to carry out our plan.

Here are three possible approaches to finding a solution to the problem. All of them are iterative. That is, you begin with an estimate of the solution, and you carry out some operation. Now you check to see if |F(X)| is small enough. If not, you carry out the operation again ...and again and again.

  1. The Table Method: Evaluate the function at a point. Based on that value, guess a better value of X, and evaluate F there.
  2. The Picture Method: Plot the function F(X). Estimate a range in which the root lies. Plot the function again, using the smaller range.
  3. The Bisection Method: Find a point where the function is positive, and one where it is negative. This is your starting interval. Look at the average of the two endpoints. This is your test point. Evaluate the function at the test point. Your new interval is the testpoint, plus whichever of the old endpoints had a function value of opposite sign.

Exercise B: Using the "table method", the "picture method", and bisection, try to find an approximate solution X. For the table and bisection methods, you should easily be able to get 10 digits of accuracy. For the picture method, just try to get 4 digits. For the picture and bisection methods, start with the range [0,3]. Try to keep track of the number of steps you took, your final value of X, and the corresponding value of |F(X)|.

For the table method, it might help a lot to have an M file do some of the work for you. I'm thinking of a command that looks something like this: 

        function picture ( a, b )
        x = ?
        y = ?
        plot ( x, y )
        hold on
        line ( [?,?], [?,?] )
        hold off
This would automate the process, and all you'd have to do is look at the picture and estimate new values for A and B.

Method...Steps......Final X......|F(X)|...
Table   
Picture   
Bisection   

Exercise C: You probably carried out the bisection method "by hand". Now you need to write down your work symbolically. Here's a first draft of an M file, called bisect_step.m, that takes a single step of the bisection method: 

        function [ a2, b2 ] = bisect_step ( a, b )
See if you can fill in the rest of the file, and figure out what it's doing. 
        fa = ?;
        fb = ?;
        c = ?;
        fc = ?;

        if ( sign ( fa ) == sign ( fc ) ) then
          a2 = ?;
          b2 = ?;
        else
          a2 = ?;
          b2 = ?;
        end
Test out your code on the Kepler problem and make sure you get the same results.

Once you're happy with bisect_step, I have a suggestion for you. Let's try to improve it! The big goal here is to make it so easy to use that all we have to do to take another step is hit the up arrow key. So we're going to change the code to look like this: 

        function new_interval = bisect_step ( interval )
Now we add two more lines immediately after the first one, to get the values of A and B out of INTERVAL: 
        a = interval(1);
        b = interval(2);
Then the rest of the routine is the same, but at the very end, we have to take A2 and B2 and pack them up into a single variable called NEW_INTERVAL: 
        new_interval = [ a2, b2 ];

Ready for the payoff? Type 

        interval = [ 0, 3 ];
        interval = bisect_step ( interval )
and then just keep hitting the up arrow key and ENTER to take the next step of bisection. Can you see what's happening here?

Exercise D: Now it's time to be a little more ambitious. We'd like to make an M file that takes as many steps of the bisection method as necessary to get the accuracy we want. Also, we'd like to be able to look for the roots of other functions. So we're going to write a fancier M file, called bisect.m, that looks like this: 

        function root = bisect ( f, a, b, tol )

        fa = f(a);
        fb = f(b);

        while ( abs ( fa ) > tol & abs ( fb ) > tol ) 

          c = ( a + b ) / 2;
          fc = f(c);

          if ( sign ( fa ) == sign ( fc ) )
            a = c;
            fa = fc;
          else
            b = c;
            fb = fc;
          end

        end

        if ( abs ( fa ) <= tol )
          root = a
        else
          root = b
        end
Can you explain what the pieces of this routine are doing? Read it carefully, and ask your mentor if you're confused.

Now, using this code, you should be able to type a command like 

        root = bisect ( 'kepler', 0, 3, 0.000001 )
and get an approximate root immediately.

Exercise E: Write three function files, f1.m, f2.m, and f3.m for the following functions, and then use your new bisection code to seek roots with an accuracy of 6 digits.
NameFormulaInterval...Root...
F1(X)X2-9[0,10] 
F2(X)COS(X)-X[0,3] 
F3(X)2*X-e-X[-10,100] 

If you'd like to work a little more on this project, ask your mentor about the secant method. This is a way to approximate the answer that is much faster than bisection.

Reference 1:
Kahaner, Moler and Nash,
Numerical Methods and Software,
Prentice Hall, 1988,
Library Call Number: TA345 K34.
Reference 2:
Charles Van Loan,
Introduction to Scientific Computing,
Prentice Hall, 1997,
Library Call Number: QA76.9 M35 V375.

Back to the Mini Projects page.

Last revised on 12 June 2001.