Numerical Computing in Mathematica

Introduction

With this lesson we will begin our study of Mathematica.  There are a number of
reasons why Mathematica is a good place to start when learning about computer-based
engineering problem solving.

First of all, Mathematica is both extremely powerful and relatively easy to use.  You
can do interesting things with Mathematica even if you are a complete beginner.  Its
support for numerical computation is far more comprehensive and flexible than
is provided by a scientific calculator.  And even though Mathematica is much more
powerful than a calculator, it is not really much harder to use.  Before
systems like Mathematica were available, the only way to do computer-based problem
solving was to write programs (usually in Fortran or C).  But with Mathematica you
can solve complex problems quickly and without doing much (if any) programming.
In order to solve the same kinds of problems in Fortran or C, you would need
much more training.

Second, Mathematica supports symbolic computation in addition to numerical
computation.  You are doubtless familiar with numerical computation, which is
what you do with a pocket calculator, and is what many scientific Fortran and C
programs do.  Symbolic computation involves manipulating mathematical symbols
much like you do when factoring polynomials or differentiating functions.

Third, Mathematica provides convenient facilities for visualizing your work.  With
Mathematica, you can easily produce two- and three-dimensional plots.  This will add
another dimension to your understanding of scientific and engineering
principles.

Finally, Mathematica is a form of programming.  We will use Mathematica to illuminate a
number of concepts that will arise later when we study C and Fortran.  The
background that you gain in using Mathematica will make it easier to learn these
languages.

In this and the lessons that follow, we will illuminate each of these aspects
of Mathematica.

Arithmetic

To use Mathematica, you type in a command (which you should think of as a question)
and Mathematica prints out the answer.  The thing that makes Mathematica powerful is the
kinds of questions you can ask and the kinds of answers it can supply.

Let's begin with simple arithmetic.  If you type in an arithmetic expression,
followed by the Enter key Ð exactly which key Mathematica takes as ``Enter'' will vary depending on your system.  There may be a key explicitly labeled Enter, or, quite possibly, some other key combination will have to be used.  A good possibility on many computer is to hit the Shift and Return keys are the same time Ð Mathematica will print back the value of the Mathematical expression.  For example, if you'd like to add 2 and 2 click somewhere in the box below (as with the previous  chapter, all Mathematica input cells in this  Tutorial have boxes put around them to give you a visual aid.  We will drop this convention starting in the next chapter).

2 + 2

The question you have just asked is ``what is two plus two,'' and the answer
that Mathematica prints back is ``four.''  Using Mathematica is that simple.

When you type a command to Mathematica, don't forget the Enter key ---it is a crucial
piece of punctuation.  Mathematica will not answer your command until it knows you
are through asking it.  Also, you have the option of ending your Mathematica code with a semicolon.  A semicolon tells Mathematica to execute the given line of input, but not to print out it's resul.  For example, if you type

a= 2 + 2;

and hit the enter key nothing will be printed out.  However, you can see that the above line of input was execute by typing

a

and the enter key, which returns the value assigned to a, 4.

You don't have to type a command on a single line, which is good because eventually we will have need to type some pretty long lines.  Instead, you can use the ``backslash'' character "\" which Mathematica recognizes to mean ``continue this line with what follows.''  For example, you could have typed:

a=2+2

Of course, there's not much point to using more than one line for such a simple
command, but this will be a handy feature when you start getting more ambitious
with Mathematica.

Naturally, Mathematica knows about more than addition.  Among other things, it knows
the following five binary arithmetic operations:

+		= addition 
- 		= subtraction 
*		= multiplication 
/		= division 
^		= exponentiation

You can do all kinds of arithmetic calculations using these operators.  Here
are a few.  Try to predict what the answer will be before trying each one out
in Mathematica.  Let's begin with one involving both addition and multiplication:

5+4*3

Did the answer surprise you?  Were you expecting an answer of 27 instead of 17?
Mathematica does multiplications and divisions before it does additions and
subtractions, an idea that should be familiar to you from algebra.  Using
another idea from algebra, how do you suppose you can force Mathematica to do the
addition first?

Click here for the answer

Simply use parentheses to group the expression differently, like this:

(5+4)*3

Here's an example using division:

6/21

Does anything about Mathematica's answer surprise you?

Click here for the answer

Mathematica gives the exact answer---the rational number two-sevenths---instead of
giving a decimal approximation as a calculator would do.  Mathematica tends not to
give approximate answers, and it goes to surprising lengths to avoid doing so.
You will see a number of examples of this behavior in this lesson.

Here's an example of taking a number to a power:

2^4

that raises 2 to the 4th power.  As another example,

2^1000

we'll raise 2 to the 1000th power.  This is another example of Mathematica's
reluctance to give approximate answers---it gives every single digit of the
answer! (Note also that Mathematica itself uses the \ character to denote line continuation).

Mathematica also knows about the factorial operator.  Recall that if n is a positive
integer, then n! is the product of the first n integers.  Perhaps it won't
surprise you at this point that

100!

produces every last digit of 100 factorial.

Saving Typing

Mathematica has a handy feature that will save you some typing, especially if you
tend to make typing mistakes.  Suppose that we want to take the average of ten
numbers:

1+2+3+4+5+6+7+8+9+10/10

The answer that we get looks wrong.  On close inspection, we notice that we
forgot to use parentheses around the additions.  At this point we could simply
reenter the corrected question, but there's an easier way to correct our
mistake.  We can actually edit the original question.

If you move the mouse cursor into the bad question and click the mouse, the
Mathematica edit cursor will appear.  You can then edit the question by typing
from the keyboard.  Place the cursor before the beginning of the question and
insert a left parenthesis.  Then insert a right parenthesis in front of the
division symbol.  Finally, type the  Enter key.  What happens?

Click here for the answer

The computation is performed again, the old answer is erased, and a new answer
is displayed.  The Mathematica ccursor is no longer active in that cell (as you can test by typign Enter again).

Assignment

When Mathematica computes an answer, it prints it out on the screen and then goes on
to the next question.  Often we would like to save an answer and use it later.
Mathematica provides {\it assignment} as a way of doing this.

For example, suppose that we'd like to average two numbers and save the result:

average = (12 + 83) / 2;

Now the symbol  average stands for the average of those two numbers.  We
can use it later just like any other number.  For example:

average*5

Numbers

Floating Points Numbers and Round Off Error

What if you really want to know a decimal approximation to a quotient?  The
simplest way is to end your numbers with decimal points.  Mathematica will take the
hint and use decimal points when giving its answers.  For example

6./21.

will give you the decimal approximation (to six places) that you might have
expected in the first place.  Similarly,

2.^1000.

will give you a six-digit approximation expressed in scientific notation.

What's going on here?  Why does Mathematica make such a big deal about expressing
answers exactly?

Mathematica is trying to avoid a common problem in computing: roundoff error.  If
you've experimented very much with a pocket calculator, you've probably
discovered that operations that should be mathematical inverses often fail to
be.  For example, the operation

10000000000000000. +1. - 10000000000000000.

should be exactly 1, but  Mathematica prints out zero!!  This  is because Mathematica stores numbers with decimal points only to a certain accuracy, and one is so small compared to 10^17 that it ``disappears'' in round-off error.  On the other hand,

10000000000000000  +1  - 10000000000000000

gives an exact answer, because Mathematica represents the integers exactly.

Imagine that you are doing arithmetic by hand.  If you restrict your
computations to integers and rational numbers (expressed as ratios of
integers), you can get an exact answer to any problem so long as you have
enough time and paper.  The same applies to Mathematica, although it uses computer
memory instead of paper.

Of course, your computations will take a tremendous amount of time if you end
up dealing with 100-digit integers.  You could save yourself a lot of time by
using decimal approximations to your fractions and by using scientific notation
to approximate your large integers.  You would, of course, be giving up a
certain amount of accuracy in the process.

In effect, Mathematica gives you a choice.  If you want answers that are correct to
the last digit, then stick to integers and ratios.  If you are satisfied with
an answer that is correct to a fixed number of significant digits, then use
decimal points when writing your numbers.  Numbers that are written with
decimal points are known in the computer world as floating-point numbers.
We'll use that term from here on out.

Evaluate each of the following floating point expressions and take a careful
look at the results.

1./9.

182347.1/234.4

.08383/38383.2

What do all of the results have in common?

Click here for the answer

Each of the results contains a six-digit mantissa, multiplied (if necessary) by
a power of ten.  By default, Mathematica computes floating point numbers to six
significant digits.

Controlling Precision

Sometimes six digits is too many, and sometimes it is not enough.  You can
easily control the number of digits that Mathematica uses in its floating point
calculations.  For example, evaluating

SetPrecision[1/9., 10]

will tell Mathematica to use a ten -digit floating point number to represent the real number version of 1/9.   Try reevaluating the three floating-point divisions above and you'll see the
difference.  You can change the length of Mathematica's floating point numbers any
time you please during a Mathematica session.  For the time being, though, let's
change it back to 10.

\send{Digits := 10

There is a tradeoff between the accuracy of numbers and the speed with which
computers can manipulate them.  This tradeoff isn't so important for the kinds
of problems that you've been doing, but it becomes extremely important in large
programs that deal with millions or billions of numbers.  Understanding the
computational differences between integers and floating point numbers is
crucially important in C and Fortran.  Not surprisingly, we'll have much more
to say about this topic later in the course.

Built-In Functions

Just like your pocket calculator, Mathematica can do more than simple arithmetic.
For example, if you'd like to determine the square root of 4, you can simply
do

Sqrt[4]

This is an example of the use of a Mathematica built-in function.  How is a function
used differently from a binary operator?

Click here for the answer

Binary operators (such as ``+'') always appear between their two arguments,
as in ``2+3''.  Built-in functions appear first, followed by their arguments.
The arguments appear in square brackets, and are separated by commas if there is
more than one.

Let's play around with Sqrt a bit more.  What do you suppose will happen
if you try to take the square root of 2?  Keep in mind what you already know
about Mathematica's treatment of integers.

Sqrt[2]

Since you asked a question involving an integer, Mathematica refuses to compromise
accuracy when reporting an answer.  Unfortunately, the square root of two is
neither an integer or a rational number.  The best that Mathematica can do is to tell
you, in effect, that the square root of two is the square root of two.  If what
you really want is a floating-point approximation, pose the question with a
floating point number.

Sqrt[2.0]

As we explained in the last section, one of the benefits of complete accuracy
is that inverse functions work as you might expect.  Thus

Sqrt[2.0] *  Sqrt[2.0]

yields something that may not be quite 2  due to round-off error (so Mathematica returns the answer as a floating point number), but

Sqrt[2] *  Sqrt[2]

yields the answer that you would expect.  What do you suppose we will get if we
take the square root of the integer 8?  Think about it before trying it out.

Sqrt[8]

As you should have expected we got an exact answer, but it is expressed in
lowest terms as two times the square root of two.

The more that you use Mathematica, the more you'll find it has the capacity to
surprise.  What do you suppose we'll get if we take the square root of -4?
Think about it before trying it out.

Sqrt[-4]

We get a complex number as a result!  Mathematica uses I to stand for the
square root of -1, as opposed to i (which is commonly used by mathematicians)
and j (which is commonly used by engineers).  You can also use complex
numbers in your questions.  For example, let's square I:

I*I

Trigonometry

Of course, Mathematica has more built-in functions than simply Sqrt.  For
example, it has all of the usual trigonometric functions and their inverses:

Sin		ArcSin  	Sec 	ArcSec 
Cos		ArcCos 		Csc		ArcCsc 
Tan		ArcTan		Cot		ArcCot

These functions deal with angles measured in radians, not in degrees.
Mathematica also knows about the constant Pi, which makes it easier to pose
questions that deal with radians.  Notice that Pi is spelled with a
capital ``P'' followed by a lower-case ``i''.  This is important to remember.

Let's try to find out the value of Pi by asking Mathematica about it:

Pi

Notice that Mathematica refuses to tell us anything about the value of Pi.  In
effect, it tells us that Pi is Pi.  Why is that?

Click here for the answer

Any number with a finite number
of digits that it printed out as the value of Pi would be an approximation,
since Pi has an infinite number of digits.  Instead, it just tells us that Pi
is Pi.  (It would be really slick if there were Greek letters on the computer
keyboard so that we could type the Greek symbol for Pi as input.)

The important thing to realize is that Mathematica treats Pi as an exact number by
representing its value with a symbol instead of with a number.  Let's
experiment with asking Mathematica some questions using the trig functions.  For
example, let's start with the cosine of zero.

Cos[ 0 ]

Mathematica obliges with an exact answer, since the cosine of zero happens to be an
integer.  Now let's try to find out the cosine of 1 radian.

Cos[1]

Because we are taking the cosine of an integer, Mathematica insists on giving us an
exact answer.  The best it can do is tell us that the cosine of one is the
cosine of one, much as it did with the square root of two.  Of course, we can
get a floating point approximation if we use a decimal point in the question.

Cos[ 1. ]

Since Pi is an exact number, we should expect to get back exact answers when we
ask questions involving Pi.

Cos[ Pi/2 ]

And voila, we get back the exact answer.  (Recall that Pi/2 radians is 90
degrees.)  How about the cosine of 60 degrees?

Cos[ Pi/3 ]

Again, the expected exact answer.  How about the cosine of 45 degrees?

Cos[ Pi/4 ]

Things are really getting interesting!  The answer isn't an integer or a
rational number, but Mathematica gives us an exact answer anyway: the square root of
two over two, which is equal to one over the sqaure root of two.  We can also work backwards.

ArcCos[ Sqrt[2]/2 ]

This time we get an exact answer that involves the exact number Pi.  Notice
what happens, though, if we use floating point numbers instead of integers.

ArcCos[ Sqrt[2. ]/2. ]

Forcing Evaluation

We have already seen that one way to convince Mathematica to give you a
floating-point answer is to use decimal points in all of your constants.  But
that is not always convenient.  For example, consider the following
calculation:

ArcCos[ 1/2 ] + ArcSin[ 1/2 ]

Another way to force Mathematica to give a floating-point answer is to use the N
function, which asks Mathematica to evaluate an expression to obtain a
floating-point number.  For example:

N[ ArcCos[ 1/2 ] + ArcSin[ 1/2 ]]

N takes an optional second argument, which specifies the number
of digits of precision you'd like your answer to have.  So we can easily
determine the first 100 digits of Pi:

N[ Pi, 100 ]

N also works with complex numbers by putting the real and imaginary
parts into floating point form.  Thus, compare

Sqrt[ -127 ]

and

N[ Sqrt[ -127 ]]

Learning More

Mathematica provides hundreds of built-in functions and constants.  No good purpose
would be served by cataloging them here---it is better to learn about them as
the need to use them arises.  You are probably wondering, though, how to at
least figure out what's available.

One source of information is printed documentation, such as The Mathematica Book,
but it is not comprehensive.  Whether you are learning about Mathematica, Fortran,
Unix, or any other software system, printed documentation tends to have a
number of limitations:

(1) It is almost always out-of-date.  Software changes almost continuously, so
printed documentation rapidly becomes dated.

(2) It is almost never there when you need it.  It is not likely that you'll
have your copy of the documentation with you every time you sit down at a
workstation.

(3) It is usually fairly expensive.

The best way to learn about a system such as Mathematica, once you know the basics,
is via online documentation.  Most systems---Mathematica and Unix are good
examples---provide extensive online documentation.  This material is generally
up-to-date, it is by definition available whenever you are using the software,
and it is free.  You should learn to use it.

For example, if you type a question mark to Mathematica it will respond by showing
you a list of available information.

?F*

Mathematica will display its answer to each help request in a separate window.
You'll need to position and place it with the mouse.

One of the topics is  Factorial, so let's try that out.

?Factorial

What you get is adescription of Mathematica's Factorial function.  Perhaps the FactorInteger function intrigues you.  We can learn more about it with

?FactorInteger

The information that is printed out may not make complete sense to you because
you are not yet a Mathematica expert.  But you should at least be able to figure out
that FactorInteger can be used to factor integers.  Once you know this much,
it's a good idea to try it out yourself.  When it comes right down to it, the
best way to learn about just about anything having to do with a computer is to
experiment.  For example

FactorInteger[ 1000 ]

shows us the prime factorization of 1000.  Perhaps you've always wondered what
the factorization of 100!\ is:

FactorInteger[ 100! ]

One useful Mathematica built-in feature is the %  symbol.  See if you
can figure out what it does by using the Mathematica help facility and by
experimenting.

Saving Your Work

You will be using Mathematica to do assignments throughout the quarter.  We will
usually ask you to prepare a report, and we will be asking you to
incorporate the output of Mathematica into your report.  To do this, you will need
to be able to save the results of a Mathematica session to a file.  Here's how.

In the Mathematica window, pull down the Notebook menu and select the Save As Special option.  One of the options will be to save the session (notebook) as plain text

All of the text that currently appears in your Mathematica window will now be saved
to the file that you specified.  To prove that to yourself, read the file that
you just saved into a text editor.  Notice that it is not exactly the same as what you
saved.  What is different?

Click here for the answer

The saved version of the Mathematica window is not as ``pretty'' as the original.
For example, square roots are displayed differently and the Greek letter for Pi
doesn't appear anywhere.

Homework

The inexactness that results from doing numerical computations with floating
point numbers often results in much more than a mistake at the last decimal
point.  If you're not careful, large errors can occur at crucial points in the
calculation that result in your answer being completely wrong.

If you are careful, though, most engineering problems can be adequately solved
with the 6 digit floating numbers that Mathematica provides by default.  In fact,
most engineering problems can be solved with smaller floating point numbers.
The reason is that very few physical measurements can be made with such
accuracy.  In view of this, Mathematica will allow you to change the size of the
floating point numbers that it uses.

In this set of exercises, you will use Mathematica to investigate the points that we
made in the last two paragraphs.

The exercises contain five questions, numbered in sequence and displayed in
boldface type.  You are to prepare a short report giving the answers to these
questions.  For each question, your report should contain the text of the
question as given below, your answer to the question, and the output from Mathematica
that shows how you arrived at your answer.

To prepare this report, you will need to save your Mathematica session to a file as
described earlier in this lesson.  Save it to a file called hw.txt.  Then
read the file into a text editor and edit it by removing the extraneous parts of the
Mathematica session and adding our questions and your commentary.  You will doubtless
want to use the mouse to copy and paste the questions from the notebook into the
text file.

Significant Digits

Let's begin by helping you get a feel for what it means to have ten-digit
floating point numbers.  Keep in mind that when we say a floating point number
has ten digits, we mean that, when written in scientific notation, its mantissa
has ten digits.

Suppose you are doing a calculation that depends upon knowing the surface area
of the earth.  Suppose further that your calculation must be accurate to ten
significant digits.  If your measure of the earth's surface area is not
accurate to ten digits, then there's no way your calculation can be accurate to
ten digits.  In other words, you're calculations can be no more accurate than
your data.

Is it reasonable to expect to be able to measure the surface area of the earth
to ten digits of accuracy?  Let's see exactly what that would entail.  First of
all, use Mathematica to calculate the surface area of the earth subject to the
following two assumptions:

Recall that the surface area of a sphere of radius $r$ is {\it 4 Pi r\^2}.

Thanks to Brigham Young, a city block in Salt Lake City is a square one-seventh
of a mile on a side.  What is the surface area of the earth expressed in city
blocks?  (In case you've forgotten, there are 5280 feet in a mile, 12 inches in
a foot, 2.54 centimeters in an inch, and 100,000 centimeters in a kilometer.)

Of course, the earth is not a perfect sphere and the radius that we have given
is not entirely precise either.  Furthermore, the surface of the earth is
extremely wrinkled.  All of this makes your estimate of the earth's surface
area, with its ten digits of ``accuracy,'' rather suspect.  Suppose that your
estimate is accurate to one-tenth of one percent.

Accumulating Error

The Outer Banks is a sand bar extending for hundreds of miles along the east
coast of the United States, mostly along the coast of North Carolina.  The sand
bar is broken by shallow channels that connect the Atlantic Ocean on the east
to the shallow salt water sounds on the west.  The result is a chain of long,
thin islands only a few feet above sea level.

The Outer Banks were the site of Sir Walter Raleigh's lost colony and the
Wright Brothers' first powered flight.  It is also known for surf fishing,
hazardous navigation, and impressive riptides.

The Wright Brothers' first powered flight was made from the highest natural
point along the Outer Banks at Kitty Hawk.  Kitty Hawk is a broad sand dune
that peaks at 38 feet above sea level.  The essence of this problem is to
determine how far out in Albemarle Sound one can see from the summit of Kitty
Hawk.

For a figure of the problem please refer to the main menue of the tutor.

You are sitting on the top of Kitty Hawk on a calm day with no waves.  You are
watching a friend swimming across the sound toward the mainland.  Because of
the curvature of the earth, your friend will eventually disappear across the
horizon.  How far away (in feet) can your friend get before he disappears from
view?  We're interested in finding the maximum possible straight line distance
from your eye to your friend.

See the figure for a sketch of what computations are required to solve this
problem.  Notice that the solution involves using the Pythagorean Theorem to
solve for the unknown side of a triangle.  The sketch is not to scale, so be
careful to notice that the unknown side is extremely small compared to the two
known sides, which are almost identical in length.

Now you should experiment with calculating the straight line distance using
floating point computations.  Begin by setting the floating point precision to
6digits.  Compute the value of the formula that you derived.  Be sure that
all of the numbers that you include in the formula are floating point
constants---not integers.  (The correct answer is approximately 39,857 feet.)

Repeat this process using successively fewer digits of accuracy.  Is what
happens what you expected?  Can you explain what's going on?  Think about these
questions before continuing.

Avoiding The Problem

Experienced computational scientists would not have approached this problem the
way that we encouraged you to do.  Instead, they would have first written the
formula to be solved in symbolic form (without any numbers) and looked for a
way to simplify the formula algebraically.

You should try that out now.  Using r to stand for the radius of the
earth and h to stand for the height of Kitty Hawk, come up with a formula
for the maximum line of sight distance.

Next, simplify the formula using your knowledge of algebra.

Finally, plug in the appropriate values for r and h and use Mathematica to
solve the problem using a number of digits that decreases from 10 all the way
down to 1.  What do you notice now?