SUNY Geneseo Department of Computer Science

Lesson 9—The “If” Statement

Purpose

Background

“If” statements are mechanisms for making choices in a program: if such-and-such is true, then do this thing, but if something else is true instead do another thing, et cetera. This lesson teaches you how to use “if” statements and their variants in Matlab, and asks you to apply them to the problem of finding places where a function is equal to 0.

You should read or watch the following before starting the in-class part of this lesson:

• Sections 4.1 through 4.3 of Attaway’s text
• “The if-else statement” mini-lecture from Edinburgh. (This mini-lecture solves a problem described on page 51 of the course booklet that comes with Edinburgh’s Matlab course—you might want to open that booklet, also available at the Edinburgh site, while watching the lecture.)
• The “Conditional statements” sub-section of the “Programming” section of the “Logic and Flow Control” chapter of Mathworks’ tutorials.
• The series of four video lectures on the Bisection method at https://www.youtube.com/playlist?list=PL3C53B3B12414AF72

In addition, if you want to do the bonus part of this lesson, you may find the following video lectures helpful:

• The introduction to Newton’s method (also known as the Newton-Raphson method) at https://www.youtube.com/watch?annotation_id=annotation_134081&feature=iv&src_vid=lFYzdOemDj8&v=nRbsfnUftT4
• The example of Newton’s method at https://www.youtube.com/watch?annotation_id=annotation_171114&feature=iv&src_vid=nRbsfnUftT4&v=lFYzdOemDj8

Exercise

The endpoint of this lesson is a Matlab script that calculates the date of an upcoming closest approach between Earth and Mars. You will get to this endpoint through a number of steps: first, you will use equations for the orbits of Earth and Mars to devise an expression for the distance between them as a function of time. To find out where this distance is at its minimum, you will use the fact that minima (and maxima) of a function correspond to places where the function’s derivative is 0. So you will find the derivative of your distance function with respect to time, and use the bisection method (and, as a bonus, Newton’s method also) to find an appropriate zero of the derivative. But before you do any of this, a bit of warm-up with “if’ statements is in order….

Warm-Up

Write a Matlab script that reads a number from the user, and prints “Big Number” if the number is greater than 10, and “Small Number” if the number is 10 or less.

Then extend your script so that it prints “Big Number” for numbers between 10 and 100, “Really Big Number” for numbers greater than 100 but less than 1000, “Small Number” for numbers between 0 and 10, and “Out-of-Bounds Number” for numbers less than 0 or above 1000.

In order to test your script (and maybe to understand exactly what it’s supposed to do), come up with examples of numbers that should produce each possible message from your script; think particularly about “boundary cases,” i.e., examples that lie exactly at the cut-off between one case and another. In fact, the directions for the extended script are deliberately vague about exactly where some boundaries lie—you need to think about where you want them to be and write your script accordingly, although I don’t particularly care what choices you make.

Earth and Mars

Earth’s orbit around the sun can be approximated by a circle in the so-called ecliptic plane (the ecliptic plane being defined as the plane of Earth’s orbit—so it is no surprise at all that whatever shape the orbit has is in that plane). This circle can be described by the parametric equations

x(t) = cos( 2πt/365.25 )
y(t) = sin( 2πt/365.25 )

where x and y give the Earth’s position relative to an arbitrary coordinate system. Distances in this coordinate system are measured in astronomical units, the average distance between Earth and the sun. The t parameter can be thought of as time, measured in days. The origin of the coordinate system is the sun.

In the same coordinate system, Mars’s orbit can be approximated by the circle

x(t) = 1.64 cos( 2πt/687 + 0.19 )
y(t) = 1.64 sin( 2πt/687 + 0.19 )

While these equations aren’t perfect, I believe that they reasonably accurately describe the current relative positions of Earth and Mars. In particular, t = 0 corresponds to March 19, 2014; positive t values give positions on days after that date, and negative ones give positions on days before that date.

Step 1. Write Matlab code that draws a plot showing the orbits of Earth and Mars as given by the above equations. As a mini-bonus, make the script somehow mark the positions of Earth and Mars on the day you do this step (the first day of this lesson, March 28, would be t = 9). See if you can figure out from this extended plot where in the sky you should look to see Mars (it helps to know that the plot is basically a sketch of the solar system as seen from very high above the north pole).

Step 2. Using the orbit equations given above, derive an equation for the distance between Earth and Mars as a function of t. Recall that from the Pythagorean theorem, the distance between points (x₁,y₁) and (x₂,y₂) is √((x₁-x₂)²+(y₁-y₂)²). Write another Matlab script (or extend the first one) to plot this distance over the interval t = 0 to t = 687 (one full orbit of Mars around the sun).

While the problem “calculate the distance between Earth and Mars for the next Mars year” doesn’t lend itself to easy examples for design and testing, you can test your work in the spirit of example-driven testing, i.e., by thinking about some typical and special features of the distance and checking that your code replicates them. In particular your plot from Step 1 should show that Earth and Mars are fairly close at low t values, with Mars slightly ahead of Earth in its orbit. But because Earth moves faster and on a smaller circle than Mars does, Earth is catching up to Mars, i.e., the distance between them is decreasing. Once Earth passes Mars however, the distance will start increasing, and will continue doing so until some point in the future, when it will start decreasing again. Finally, the distance at the end of one Mars year will not be the same as at the beginning of the year, because Earth will be somewhere else in its orbit at that time. Inspect your plot of distance to see if it has these general characteristics—if it does, your distance function and the code based on it are probably correct.

Step 3. Find the derivative of your distance function with respect to t, and code that derivative as a Matlab function. There are at least 3 ways you could do this, each with advantages and disadvantages:

1. You can just work out the derivative by hand. This isn’t as hard to do as it might seem, it mostly just involves meticulous use of the chain rule. This requires some manual work, but doing it can expose some structure in the derivative that will help you code it in Matlab.
2. Learn about Matlab’s symbolic math features (i.e., features that allow Matlab to do algebra on expressions in order to calculate new expressions, instead of treating expressions just as formulas to follow in order to produce numbers) and use them to have Matlab calculate the derivative symbolically for you. This saves you the work of manually finding a derivative, but will probably take more time in total, since you will need to learn new Matlab features and apply them. But you will come out of the experience knowing a little more about Matlab, too.
3. Write a Matlab function that can evaluate the derivative by approximating the definition of “derivative.” The definition says that

   f′(x) = lim_h→0( f(x+h) - f(x) ) / h

   A typical way to approximate this is to set h to a very small value, evaluate the function in question at the x value in question and at that x value plus h, and divide the difference by h (i.e., do exactly what the right hand side of the definition says, except replace the idea of h approaching 0 with a single fixed h that is very close to 0). The advantage of this approach is that it is something you’ll probably need to do anyhow (albeit for a different function) if you do the Newton’s method bonus; the big disadvantage is that this technique can be very sensitive to round-off error in the computer, and so you will probably need to fine-tune the value of h in order to get a good approximation to the derivative.

However you calculate the derivative, extend your script from Step 2 to plot a scaled version of it on the same axes as the distance between Earth and Mars. I suggest that you plot a “scaled version” of the derivative because the actual derivative is very small compared to the distance, and probably won’t be very recognizable if plotted without scaling. But if you multiply the derivative by, say, 100, and plot the result, you’ll get a curve that you can see and relate to the distance curve.

Once again, you can test your derivative code by checking that the plot has features you would expect of the derivative: that it is positive where the distance between Earth and Mars is increasing, negative where the distance is decreasing, zero where the distance seems to reach its minimum and maximum, etc.

Step 4. Write Matlab code (either as another extension to the script from Step 3, or as a new function) that uses the bisection method to find the t value at which the derivative of distance is 0 for the next closest approach between Earth and Mars. Use your plot of distance and its derivative from Step 3 to get the initial interval for bisection.

One way to test your result from this step is to compare it to geometric intuition. Intuitively, the closest approach between Earth and Mars should happen when they lie on the same line from the sun (i.e., when the sun, at the origin of the coordinate system, Earth, and Mars form a straight line). This happens when the angle between Earth’s position and the X axis is the same as the angle between Mars’s position and the X axis. You can figure out when this happens by setting the arguments to the cosine (or sine) functions in the equations for Earth’s and Mars’s positions equal, and solving for t. However, while geometric intuition can be a good guide to what to expect from a program such as this one, it’s a good idea to base the official output on something more rigorous (such as the zero of the derivative).

When will the closest approach happen? What will the distance between the two planets be at this time?

Bonus

If you finish all of the above and want more to do, write Matlab code that uses Newton’s method to find the zero of the derivative of distance. Newton’s method is another popular algorithm for finding zeros of functions, and is well worth having some experience with. Compare the answers that Newton’s method and bisection produce for the date of Earth’s next closest approach to Mars. They should be the same.

Follow-Up

In your study group’s first face-to-face meeting with me following your “Complete By” date above, I will look over your solutions to the exercise, ask you any further questions I have, and answer any questions from you. I may also ask you to demonstrate your script(s) in Matlab. Please bring at least one copy of the necessary “.m” file(s) to your meeting. This will speed the meeting along.

In study group meetings before your “Complete By” date, I will ask for progress reports from you on this lesson, answer questions you have, look at any intermediate code you want to show me, etc.