Symbolic Computing 
in Mathematica

Introduction

The previous lesson concerned the use of Mathematica as a flexible and powerful
scientific calculator.  If this were all Mathematica could do, it wouldn't be a
particularly interesting piece of software.  In this lesson, we'll go beyond
numerical computing and look at how Mathematica supports  symbolic computing.

When you first studied math in elementary school, numerical computing was all
that you did.  Math meant arithmetic in those days, and you learned to add,
subtract, multiply, and divide.  As you continued in school, you began to mix
symbols into your arithmetic by solving for unknowns in word problems.
Eventually, you began doing mathematics that didn't involve finding numerical
solutions at all: factoring polynomials, doing geometry, computing limits,
and integrating and differentiating functions are all examples.

In the remainder of this lesson we'll look at a number of Mathematica's facilities
for symbolic computing.  By the end you will know enough to begin exploiting
Mathematica in the other science, math, and engineering classes that you may be
currently taking.

Solving Equations

Introduction

One major activity at most any level of mathematics is solving equations, a
task at which Mathematica excels.

For example, let's solve a simple equation in one unknown:

Solve[2*x + 5 == 0, x]

Notice that Solve[ , ] takes two arguments: the equation to be solved and
the unknown to be solved for.  Since none of the numbers  in the equation 
are floating point numbers, Mathematica
prints out the solution as an exact rational number.  There are two ways to get
the answer into floating point form.  You can write at least one of the numbers
in floating point form, or you can use NSolve instead of Solve:

NSolve[2*x + 5 == 0, x]

Whereas Solve tries to solve equations symbolically, NSolve tries
to solve them numerically (there are very many other Mathematica functions which follow this same pattern, such as Integrate and NIntegrate for integrating functions analytically or 
numerically).  The distinction may be vague to you at the moment,
but it should be clarified in future lessons.

You can save away the result of solving an equation, just as you can save the
result of anything else, by using assignment as you learned in the last lesson:

result = Solve[2*x + 5==0,x]

The careless use of assignment can cause problems when you are doing symbolic
computations, and now is a good time to point that out.  Consider, for
example, the following computation:

Solve[result + 5 ==0, result]

Can you explain what happens?

Even though result looks like a variable in the call to solve, it really
isn't because we've previously given it a value.  So what you're really asking
Mathematica to solve is the equation 

 -5/2 + 5 = 0
 
  for the ``variable'' -5/2, which of course doesn't make any sense.

This is one of the most confusing aspects of Mathematica.  The best policy is to use
different symbols for the variables to which you assign values than you do for
the variables that are true mathematical unknowns.

It is possible to ``unassign'' a variable that contains a value:

result=.

Having ``unassigned'' result, we can now solve the equation that
couldn't be solved before:

Solve[result + 5 == 0, result]

Multiple Solutions

It is not unusual for an equation to have more than one solution.  For example,

\send{solve(x\^2 + x = 7);}{MathematicaV}{solve(x^2 + x = 7);[Return]}

gives us two answers in complete detail.  You can use assignment to save the
answer even when there is more than one solution.  For example,

answers = Solve[x^2 + x == 7, x]

Now the variable  answers contains  both solutions, as you can verify
by asking Mathematica for its value:

answers

If you want to extract, say, the second solution, you can do so by 
indexing the variable:

answers[[2]]

It is not hard to come up with an equation that has an infinite number of
solutions:

Solve[Cos[x] == 0, x]

Mathematica recognizes that fact and gives only one of the solutions.  It is not hard
to come up with an equation that Mathematica can't solve.  We'll shortly see one,
and discuss what to do about it.

As you have just seen, Mathematica can find none, some, or all of the solutions to
an equation, but it doesn't tell you how many it didn't find.  What this
means is that there's no substitute for having a bit of mathematical sense.
You need to have an idea of how many answers you're expecting so that you can
properly interpret Mathematica's solutions.

Compound Interest

Now that we've shown you some of the symbolic computing capabilities of Mathematica,
we're going to consider a series of problems from economics involving compound
interest.  We have several goals in doing this.

First, we want to use some of what you're just learned about Mathematica and at the
same time show you some more of its capabilities.  You should also pick up on
the fact that there's no way we can show you all that Mathematica has to offer.
Mathematica  (your textbook) will give you an idea of the breadth of
Mathematica.  As you already know, there is also substantial online help available.  
Most any new concept that you are likely to learn in an undergraduate math class  is
supported in some form by Mathematica.

Second, we want to you to come to appreciate that Mathematica will not do your math
homework for you.  Mathematica is not magic---you have to understand a lot about the
underlying mathematics in order to make progress with Mathematica.  If you are doing
a physics word problem, for example, you have to understand how to transform
the problem into a purely mathematical form.  At that point, if you choose,
Mathematica can help you arrive at a solution.  Think of Mathematica as something that will
help you with the uninteresting details of a problem, in much the same way that
a calculator does.

Fixed-Rate Interest

One of the most familiar economic models is the fixed interest rate model.  In
its various forms, it can be used to understand such phenomena as the annual
growth of an investment, the effects of inflation, the result of unchecked
population growth, and the spread of infectious diseases.

Here is the compound interest problem in its simplest form.  Suppose that you
invest P dollars in a bank account that pays annual interest at a rate R.
(For example, 5% interest would be an interest rate of .05.)  The interest is
credited to your account at the end of each year, at which point you begin
earning interest on the interest.  How much money will you have in your account
after n years?  Try to come up with a formula, involving P, R, and n,
that gives the answer.

Click here for the answer

At the beginning (year 0) you have P dollars.  After one year, you'll still
have the P dollars, plus P*R more, for a total of P*(1+R).  By the same
reasoning, after two years you'll have P*(1+R)*(1+R), Generalizing, at the
end of n years you'll have P*(1+R)^n dollars.  Notice that this formula
works even for n = 0, when it reduces to simply P.

Let's enter the formula into Mathematica so that you can see it printed out in an
easier-to-read form:

P*(1+R)^n

The formula is now sitting in the Mathematica window, but it isn't very useful for
computational purposes.  What would you have to do to figure out what $100
invested at a rate of .05 would be worth after 50 years?

Click here for the answer

You'd have to reenter the formula, only this time using 100 instead of P, .05
instead of R, and 50 instead of n.  Try it out for a few different numbers.

Because the formula's pretty short, this isn't too much of a pain to do.  But a
more complicated formulaÐor even this one, if you try to use it
tomorrowÐmight be difficult to remember.  It would be handy if you could
calculate compound interest by entering something like Compound[P,R,n].
Unfortunately, Mathematica has no built-in compound interest calculation function.

User-Defined Functions

Fortunately, it is possible for a Mathematica user to define new functions.
This is one of the key things that makes Mathematica less like a calculator and more
like a programming language.  It is also what will make you less like a
calculator user and more like a programmer.

Let's begin with a simple example.  Suppose that you'd like to have a function
called square that takes a number as an argument and returns its square
as its result (by convention, internal Mathematica functions start with a capital letter, 
while user defined functions start with a lower case letter.
Here's how to do that in Mathematica:

square[x_]:= x*x;

Now you can use the function square just as if it were built in.  Here's
one example of its use; you should experiment with others.

square[5]

The definition of the square function may look funny to you.  We are using
assignment to name the function square.  The x_ symbol is what
tells Mathematica that we are defining a function of the variable x.

When you use a user-defined function, Mathematica takes the argument that you provide
(5 in the example above), substitutes it for the parameter (x) throughout the
body of the function (x*x), yielding an expression (5*5), which it then
evaluates as usual (producing 25 in the example above).

Try writing a function called cube that takes one argument and returns
its cube.  Try it out to make sure it works.  When you're done,
congratulations! Ð you've just written a program.

Now let's return to our original problem.  We'd like to have a function called
compound that will calculate compound interest.  So let's define one:

compound[p_, r_, n_] := p*(1+r)^n;

There's one extra complication here.  The function compound takes three
arguments instead of one, which requires listing all three of the parameters
as shown.  Now to calculate the consequences of investing $100 at .05 interest
for 50 years, we need only do:

compound[100, .05, 50]

See if you can use the compound function to answer the following
questions involving compound interest.

(1) If you invest $1000 at 6% interest right now, how much will it be worth
in 40 years?

(2) What if you invest it at 12% for 20 years?

Exponential Growth

Let's put our formula for compound interest into a slightly different form.
Right now it involves raising 1+R to the power n.  Let's see if we can
express is instead as raising E (Mathematica's symbol for e, the base of the
natural logarithm) to some other power.  This transformation will lead us to
further insights later in the lesson.

If you are comfortable with logarithms and exponentials, you can easily do this
transformation with pencil and paper.  If you're rusty, Mathematica can help you out,
though you'll still have to provide a fair amount of guidance.  Let's begin
by asking Mathematica to solve an equation:

answer = Solve[(1+R)^n == E^x, x]

In other words, let's ask Mathematica to symbolically compute the proper power of
E.  (We're saving the answer in a variable so that we can use it again
below.)  There's something different about the equation that we're asking Mathematica
to solve here.  What is it?

Click here for the answer

The equation involves three unknowns: R, n, and x.  But we want a
solution for x in terms of R and n.  And that is exactly what Mathematica gave
us, since we specified that we wanted it to solve for x.

The answer that Mathematica gives us doesn't look as simple as it could be, although
you won't realize that unless you are familiar with logarithms.  Let's see if
Mathematica can simplify it further:

Simplify[answer]

There are any number of ways for Mathematica to give an answer, especially as answers
start becoming more complicated.  Sometimes asking Mathematica to simplify a solution
will give you an expression more to your liking.  At other times you might ask
Mathematica to factor an answer---and Factor is another built-in function.

Now you know all that you need to know to redefine the compound function
using a power of E.  Can you do it?

Click here for the answer

compound[P_,R_,n_] := P*E^(n*Log[1+R]);

In this formula, the Log[1+R] part of the exponent is called the  growth
rate.  In general, the exponent of E in a compound-interest calculation,
when divided through by n, yields the growth rate.  The reason that we
bothered to express our compound interest formula using a power of E was just
so we could discuss growth rates.

For any interest rate, the growth rate is never larger than the interest rate.
For small interest rates (such as 5%) though, the growth rate is only a little
bit smaller than the interest rate.  You can easily verify this using Mathematica.

Experiment with the new function for compound interest.  It should give you the
same answers as the old formula with one exception.  Just as it does with Pi,
Mathematica will not convert E into a floating point number unless you insist by
using NSolve.

Doubling an Investment

If we invest $1000 at 8% interest, compounded annually, how many years will
it be before our investment is worth $2000?  Using our definition of  
compound, it is not hard to set up an equation for Mathematica to solve that will
give you the answer to this question.  The answer is a little over 9, but be
sure that you can make Mathematica give you the answer before continuing.

Now let's generalize the question a bit.  If we invest $P at 8% interest,
compounded annually, how many years will it be before our investment is worth
$2*P?  You can use virtually the same equation to solve this problem, and
once again the answer is a little over 9.  Again, be sure that you can make
Mathematica give you the answer before continuing.

Apparently, the doubling time doesn't depend upon the initial investment, but
it does, of course, depend upon the interest rate.

Let's generalize the problem one step further.  If we invest $P at an interest
rate of R, compounded annually, how many years will be it before our investment
is worth $2*P?  This is exactly like the last question, except that we must use
an interest rate of R instead of .08.

Click here for the answer

Solve[compound[P,R,n] == 2*P, n]

Not surprisingly, the answer involves R.  Now take this answer and turn it
into a function called doubling that will take R as an argument and
return the doubling time for an investment at interest rate R.

Click here for the answer

doubling[R_] :=  Log[2] / Log[1+R];

Try the function out.  Notice that the answers are all given as products 
involving the natural log of 2.  See if you can doctor the definition of
doubling so that it will always give out a floating point answer.

Click here for the answer

The simplest thing is to change the ``2'' and the ``1'' in the body of the
function to ``2.0'' and ``1.0'' respectively.  Alternatively, you could
surround the body of the function with a call to N.

Approximate Doubling

Our formula for calculating doubling times is fairly simple, but there is an
even simpler approximate formula that does reasonably well for the small
interest rates common to everyday life.  To calculate the approximate time to
doubling, simply divide 0.69 by the interest rate.

The is called the rule of 69, and it is commonly used for quick calculations.
(You may have heard of it as the rule of 70, since it is usually easier to
divide something into 70 than to divide it into 69.)

Define a function called approx that takes an interest rate as its
argument and returns the approximate time to doubling using the rule of 69.
Experiment with doubling and approx to get a feel for how good the
approximation is.  You'll find that it's within 10% up to an interest rate
of around .20.  This is larger than the interest rates common to such
investments as bank accounts, mutual funds, mortgages, and credit cards.

As a good Mathematica user, you might what to know at exactly what interest rate the
error in the approximate doubling formula becomes 10% of the true answer.
It's not hard to ask the question using Solve, so let's do it:

Solve[doubling[R]-approx[R] / doubling[R] == .10, R]

Unfortunately, Mathematica is unable to solve the equation.  (Sometimes Mathematica will
spend a long time working on a problem that it can't solve.  You can tell that
it is working on such a problem because the cursor turns into a clock face.  If
you get tired of waiting for an answer, you can interrupt the computation by
typing Command. during the calculation.)

If we ask the same question using NSolve, though, we do get an
answer:

NSolve[doubling[R]-approx[R] / doubling[R] == .10, R]

What's going on?  The answer is that Solve tries to solve the problem
symbolically, which essentially means that Mathematica is pushing symbols
around to try to find an exact answer.  And it gets stuck.  By contrast, 
NSolve tries to solve the problem numerically, which means that it is
internally doing number crunching to find an approximate (but very close)
solution.

The approximate doubling formula is handy, but you may be wondering where
it came from.  If you have studied Taylor series expansions in your calculus
class, you can easily derive it.  Here's how.

Recall that the body of our approximation formula divides ln(2) by ln(1+R),
where R is the interest rate.  As you may know, an expression such as
ln(1+R) can be approximated by taking the first few terms of its Taylor
series expansion.

What?  You don't happen to know the Taylor series expansion of ln(1+R), and
you don't remember how to derive it?  That's OK---Mathematica can do it for you.  All
you need is enough mathematical knowledge to ask the question and understand
the answer.

Series[Log[1+R], {R, 0, 3}]

This gives us the first three terms in the Taylor series expansion.  (The last
term that Mathematica gives you is simply an indication that there are more terms,
all powers of R or higher, that are not displayed.  You can ask Mathematica to give
you more terms if you wish---use the help facility if you're interested how.)

Let's use the first term of the expansion, which is R, as an approximation
to ln(1+R).  The body of our function now reduces approximately to the
quotient of ln(2) and R.  And since ln(2) is approximately 0.69, we end
up with the rule of 69.

Compounding Intervals

So far, both our interest rates and our compounding interval have been annual.
In the real world, interest rates are indeed generally expressed as annual
rates, but interest is usually compounded more frequently.  A checking account,
for example, may earn interest that is compounded daily, while on a mortgage
interest is typically compounded monthly.

Let's see what would happen if interest were compounded bi-annually.  For
example, if instead of getting 6% interest at the end of the year, what would
happen if we got 3% interest in the middle of the year and then another 3%
at the end of the year?

Let's repeat the analysis that we did when we first considered compound
interest.  If we invest an amount P at an annual rate of R, then after six
months we would have (1+R/2)P.  At the end of year, then, we would have
(1+R/2)(1+R/2)P.  Try to come up with a formula for the amount of money in
the account after n full years.  Once you have it, define a function 
compoundtwice, along the lines of  compound.

Click here for the answer

Each compounding is at half the usual rate, but we compound twice as often per
year:

compoundtwice[P_, r_, n_] := P*(1+R/2)^(2*n)

Now that we have an expression for bi-annual compound interest, let's see if
we can express it as a power of E.  As we did before, let's start by
solving an equation:

answer = Solve[E^x == (1+R/2)^(2*n), x]

As before, this gives us the power to which we should raise E.  Let's see if
we can obtain a simpler formula, as we did last time:

Simplify[answer]

What we get via simplification doesn't look very simple.  What went wrong?

Click here for an answer

As we mentioned before, Simplify doesn't always give us exactly what
we want.  Here, it has expanded things out a bit more than we might like.
Fortunately, Mathematica allows you to specify what kind of simplification you'd
like it to perform.

Of course, given that you probably know something about simplifying logarithms
of powers, you could have done this by hand.  The flexible simplification
capabilities of Mathematica come in hand when the expression in question is too
complicated to deal with easily by hand.

We can now redefine compoundtwice in the new form:

compoundtwice[P_,R_,n_) := P*E^(2*n*Log[1+R/2];

What is the growth rate for an investment that is compounded bi-annually?  How
does this compare with the growth rate for an investment that is compounded
annually?

We now have the functions  compound and  compoundtwice.  Let's take
a look at them in Mathematica:

compound[P,R,n]

compoundtwice[P,R,n]

Compare them, and see if you can guess what the function for doing quarterly
compounding would look like.  What about monthly compounding?

Click here for the answer

All that changes is the exponent.  For quarterly compounding it would be 
4*n*Log[1+R/4], and for monthly compounding it would be  12*n*Log[1+R/12].
In other words, the exponent is always n times the growth rate.  This can be
easily verified using Mathematica.

Based upon this observation, see if you can come up with a function that gives
the amount of interest earned on an initial investment of P, at an interest
rate of R, over a period of n years, with interest compounded m times per
year.  This is just our original compound interest problem with the extra
variable m added.  Call your function  intervals.

Click here for the answer

intervals[P_,R_,n_,m_]:= P*E^(m*n*Log[1+R/m])

Back in the early 1970s, the interest rates that banks could pay on regular
savings accounts were strictly regulated.  Typically, a bank could pay no more
than 5% on such accounts.  After a time, banking deregulation and double-digit
inflation led to this limit being revoked and interest rates rose higher, at
least until the 1990s.

Since the amount of interest that they could pay depositors was regulated by
the government, one way that banks competed in the 1970s was by offering
favorable compounding intervals.  If one bank compounded monthly while another
compounded daily, the second bank would enjoy a competitive advantage even
though both were paying 5% annual interest.

Of course, there was no reason to stop at daily compounding.  A bank could, if
it chose, compound every hour or every minute or every second.  And in fact,
many banks took it to the limit and began compounding  continuously.

But what is continuous compounding?  Just as monthly compounding means
paying out interest 12 times per year at one-twelfth of the annual rate,
continuous compounding (at least conceptually) means paying out interest an
infinite number of times per year at an infinitesimal faction of the annual
rate.

Here's the key idea: we can derive an expression for the growth of a
continuously compounded investment by finding the limit of the compounding
expression  P*E^(m*n*ln(1+R/m)) as the number of compounding intervals
m goes to infinity.

It is probably not immediately obvious to you what the limit at infinity of
this expression is, and you may not even know how to go about calculating it.
But as long as you know that you need to calculate the limit, Mathematica can help
you with the calculation.

Let's warm up by using Mathematica to calculate a few interesting limits.  For
example, how about the limit of 1/x^2 as x goes to zero:

Limit[1/x^2, x->0]

Or the limit of 1/x as x goes to infinity:

Limit[1/x, x-> Infinity]

Or the limit of  1/ x as x goes to zero:

Limit[1/x, x-> 0]

Now let's return to our original problem: coming up with an expression for
continuously compounded interest.  To interpret the answer, you need to know
that  Exp[x] is another way of writing E^x.  This notation is
common in both Mathematica and in mathematics at large.

Limit[P*E^(m*n*Log[1+R/m]), m -> Infinity]

Here's a simpler way to compute the same limit.  Notice that we use the
function intervals that we previously defined, but we give it 
variables as arguments.

Limit[intervals[P,R,n,m], m -> Infinity]

What is the growth rate of a continuously-compounded investment?  What is
interesting about it?

Click here for the answer

As the number of compounding intervals increase, the growth rate comes ever
closer to the interest rate.  When an investment is compounded continuously,
the growth rate actually equals the interest rate.  In effect, what the
government was regulating in the 1970s was not the interest rate but the growth
rate!

Define a function continuous, along the same lines as compound and
compoundtwice, that calculates continuously compounded interest.

Mathematica Summary

In this section we will summarize the new Mathematica commands that you've
encountered in this lesson, and point you to further reading in the Mathematica handbook.


Solve allows you to symbolically solve an equation for an unknown.  You
must specify the equation and the unknown.  Actually, Solve is more
general than we have let on, as you can solve more than one equation for more
than one unknown.  For further information on Solve, see the Mathematica handbook.

Nsolve works much like Solve, but it finds numerical
solutions to equations.  You can learn more about NSolve in
the Mathematica book as well, or use the on-line documentation.

``Unassignment'', such as 

		a=.

refers to using assignment to remove the value that has
been assigned to a variable.  

Function definitions are used to add user-defined functions to Mathematica's
repertoire.  They are created as noted above, and can take any number of variables.

Simplify allows you to simplify expressions in various ways. 
There are other, more specialized, simplification commands as well, 
and we will come to them in time. 
  
Series will give you the Taylor series expansion of a function.  

Limit will calculate the limit of a function.

Homework

In this section we will pose a number of short problems involving
compound interest and the Mathematica functions that you've learned
about so far.  You are to use a Mathematica notbook to prepare 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.

When you are done with your report, you should save it to a file
called hw2.ma and submit if for grading.


(1) Between 1970 and 1990, the population of the United States
increased from approximately 200 million to approximately 250
million.  What has been the average annual rate of population
increase during that period?  What will the population of the United
States be in 2090 if this rate of increase continues?

(2) Derive a function that takes an interest rate R as its argument and
returns the number of years that it would take an investment at that rate to
double in value under continuous compounding.

(3) Derive a function triple that takes an interest rate R as its
argument and returns the number of years it would take an investment
to triple in value under interest that is compounded annually.

(4) Derive a good approximation to triple (called approx).  Like
the approximation to double earlier in this lesson, it should involve
only a single division.

(5) Consider a bank account that pays interest R, compounded m times per
year.  We say that this account is effectively paying interest at an
annual rate of S if an investment P in the account will be worth P(1+S)
at the end of the year.  Define a function effective that takes an
interest rate and a compounding interval as arguments and returns the effective
annual interest rate.

(6) By taking the limit of your effective function as the number of
compounding intervals approaches infinity, determine the effective interest
rate of a continuously compounded investment.