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 XAt 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 ) - MNow 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.
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 offThis 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 = ?; endTest 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 endCan 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.
Name | Formula | Interval | ...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.
Back to the Mini Projects page.
Last revised on 12 June 2001.