Visualization with Mathematica

Introduction

Visualization is an extremely important engineering activity for at least two
different reasons:

1) The ability to visualize problems and their solutions (for example, by
graphing or animating them) makes it easier to develop and apply intuition or
``common sense.''

2) Visual aids are almost always necessary in order to effectively communicate
engineering results.

Collectively, the different types of visualization represent a significant
fraction of all engineering activity.

As a result, engineering has a rich tradition surrounding the use of simple
plots.  Plots represent represent the bulk of the visual aids used in
engineering reports.  Your ability to effectively incorporate plots into your
writing will dramatically affect how your writing is perceived.  Although
effective communication is no substitute for quality work, it is probably the
next most important factor in engineering success.

There are a wide variety of software packages available on workstations and
personal computers for producing plots.  Mathematica has a simple, though not
especially flexible, plotting capability.  In this lesson we will study the
problem of using Mathematica to produce two-dimensional plots.

As you advance in your career, you will doubtless use other plotting packages.
You should keep in mind from the outset that even the fanciest package will be
useless unless you know how to use plots to effectively convey information.

Simple 2D Graphs

Producing a two-dimensional graph in Mathematica is quite easy.  To do this you must
give Mathematica:

(1) A symbolic expression containing at most one variable.

(2) The variable that is contained within the expression.  The horizontal axis
of the graph will correspond to values of this variable; the vertical axis will
correspond to values of the expression.

(3) The two endpoints of the horizontal axis.

For example, the Mathematica command

Plot[2*x^3 - 3*x^2 + 1, {x, 0, 1}];

will plot the value of 2*x^3 - 3*x^2 + 1 as x ranges from 0 to 1.

Multiple Graphs

You can graph more than one function on the same plot.  Both must be functions
of the same variable.  For example, to graph the function from the last section
along with  sin(x) on the same plot:

Plot[{2*x^3 - 3*x^2 + 1, Sin[x]}, {x, 0, 1}];

Notice that we simply enclose the two expressions to be graphed inside curly
braces and separate them with commas.  Mathematica doesn't label the two curves, so
you have to exercise some care in telling them apart.

Command-Click the mouse inside the plot window.  You will see that the horizontal and
vertical coordinates of the mouse point are displayed in a field at the top or bottom of the notebook  window (exactly where depends on which version of Mathematica you have).
This is useful, for example, if you want a rough idea of where two curves
intersect.  Where do the two curves that you've just plotted intersect?

Click here for the answer

They intersect at approximately x=0.5087 and y=0.4870.  Try using Solve[]
to obtain the answer symbolically.  You'll discover that Mathematica isn't able to
come up with the solution.

What Can Go Wrong?

It is easy to make mistakes when asking Mathematica to produce a plot, and
unfortunately Mathematica does not always tell you what's gone wrong.  Let's look at
a few examples so you'll know what to expect.

All functions that you plot must be functions of a single variable.  Mathematica will
give you an error message if you violate this restriction.  For example, try
plotting x*y and see what happens.

Mathematica will only plot a function at the points at which it is real-valued.  Try
plotting Sqrt[x] as x ranges from -5 to 5.  What happens?

Click here for the answer

Mathematica displays a curve only for values of x between 0 and 5.  This makes
sense, since Sqrt[x] is undefined for negative numbers.

This sensible behavior can be confusing in a couple of contexts.  To see what
we mean, try plotting x*I as x ranges from 0 to 10.  What happens?

Click here for the answer

The graph is completely blank.  This makes sense if you remember that  I
is not a variable, but is Mathematica's symbol for the square root of -1.  As a
result, the expression that you tried is not real-valued for any value of
x.

Now try plotting z(x) as x ranges from 0 to 10.  Can you explain
what happens?

Click here for the answer

Again the plot is empty.  In this case, the expression is undefined for all
values of x simply because the function z is undefined.  The error
message that we get here clues us in to the problem.

Visualizing Oscillation

Engineering problems can rarely be completely described by a single real-valued
function of one real-valued variable.  (This is both good news and bad news.
It is what makes engineering interesting, and it is what makes engineering
hard!)  This means that most of the skill in visualization (at least at the
level of two-dimensional plotting) is the ability to figure out what to graph,
and what not to graph, in order to convey insight into the real problem.

In general, this is accomplished by graphing a simplified version of the more
general function that describes the solution to a problem.  In other words,
having solved a problem, an engineer must then find a function whose graph
conveys the essentials of the solution without overwhelming the reader with
extraneous detail.

The rest of this lesson illustrates this idea in action.  We will study the
problem of  forced oscillation, which shows up in all areas of
engineering.  That is, we will try to understand how a physical system (such as
a car) vibrates when subjected to a repetitive external force (such as bumps in
the road).  We will try to impress upon you the critical role that
two-dimensional plotting has in communicating the solutions to these kinds of
problems.

Oscillation

Understanding the behavior of many physical systems involves understanding how
they vibrate.  We are usually interested in the frequency and the 
amplitude of a vibration.

A car driving across a bumpy road will bounce up and down.  The spacing and
height of the bumps, the speed of the car, and the makeup of the car's
suspension system determines the frequency and amplitude at which the car
bounces.  Here, the frequency is the number of bounces per second, and the
amplitude is the maximum amount up or down by which the car bounces.  (Most
passengers would prefer a low amplitude bounce!)

Oscillations such as this are traditionally modeled by scientists and engineers
with functions of the form

A cos(2 Pi f t)

where A is the amplitude of the oscillation, f is the frequency of the
oscillation, and t is time.  Take some time now to use Mathematica to generate a
variety of plots for this function for different values of A and f.  What
will the two axes of the plot be?

Click here for the answer

The horizontal axis will be time measured in seconds, and the vertical axis
will be the vertical position of a fixed point (say, the driver's seat) of the
car.  A positive position means that the car has bounced up, and a negative
position means that the car has bounced down.

Do your first plot by using 1 for both the amplitude and the frequency.  You
should notice that the resulting curve oscillates between a maximum of 1 and
and a minimum of -1, which reflects the amplitude.  You should also notice that
the curve goes through an entire oscillation once every second, which reflects
the frequency.  Frequency is typically measured in units of Hz, where 1
Hz means 1 complete oscillation per second.

Try some more visualization.  What happens if you double A or cut it in half?
What happens if you double or halve f?  What happens if you double A and
f simultaneously?  Double one and halve the other?

Be sure that you obtain a good physical intuition for this basic
single-frequency oscillation function before you move on.  You will rapidly
become lost in the rest of the lesson if you don't develop a good intuition
here.  Before you continue, you should discover that the function is quite
simple to understand.  Hopefully, you will also conclude that visualizing this
function was instrumental in understanding it.

System Response

Things don't generally vibrate by themselves.  A car, for example, vibrates in
response to the bumps in the road.  Suppose that you are designing the
suspension system (springs and shock absorbers) of a car.  Ideally, you'd like
your car to have a smooth ride no matter the pattern of bumps on the road.
That is, you'd like the car to have low-amplitude (or even no-amplitude)
vibrations under all driving conditions.

A car's suspension is a ``linear system.''  This means that the frequency with
which the car bounces up and down will always be equal to the frequency at
which the bumps in the road occur.  For example, if you go over bumps at the
rate of five per second (in other words, at a frequency of 5 Hz), your car will
bounce up and down at a frequency of 5 Hz.  There's nothing you can do about
it.

You can design the car's suspension to minimize the amplitude of
the bounce.  But this job is not as straightforward as you might imagine.  Any
suspension system that you design will have different behaviors at different
frequencies.  You might, for example, design a suspension that produces moderate
bounces at low frequencies (bumps spaced far apart) and small bounces at
high frequencies (bumps spaced close together).

What all this means is that if you are designing a car suspension, you will be
interested in its gain.  The gain of a system is a function of frequency,
and is the ratio between the amplitude of the response of the system being
vibrated and the amplitude of the system causing the vibration.

For example, suppose that you are driving a car along a road with regularly
spaced 12-inch high bumps.  Suppose that the car has a gain of 1.0 at a
frequency of 1 Hz and a gain of 0.25 at a frequency of 10 Hz.  This means that
if you hit one bump per second, the car will bounce up and down 12 inches per
vibration.  (We obtain this by multiplying the size of the bumps by the gain.)
If you go ten times faster, the car will bounce up and down only three inches
per vibration.  (Of course, it will be bouncing ten times more frequently.)

If you take more advanced classes in mechanical or electrical engineering,
you'll learn how to derive the gains for mechanical and electrical systems.
Since it is not our intention to teach you how to derive gains in this course,
we're going to give you some that we have already derived.  You will then be
able to experiment with them.

Given below are functions, called transfer functions, that describe the
behavior of the suspensions of three kinds of cars: a full-size American car
from around 1980 (boat), a small off-road Japanese truck (truck), and an
expensive European touring sedan (sedan).

boat[s_]:=s -> (9*s + 53) / (27*s*s + 9*s + 53);

truck[s_]:=(160*s + 400) / (9*s*s + 160*s + 400);

sedan[s_]:=(2.143*(s*s+49))/((0.25*s*s+2.86*s+9.14)*
(0.25*s*s+2.104*s+11.488));

Each transfer function takes a complex number as its argument and computes a
complex number as its result.

Given a transfer function G, you can obtain the gain for a given frequency
f by taking the absolute value of G(f I).  (Don't forget that I is
Mathematica's symbol for the square root of -1.)  You may not be familiar with the
notion of the absolute value of a complex number.  The important things to know
are that it is a real number and that it can be computed with the
function Abs[ ].

Graph the gain of each kind of car as a function of frequency, for frequencies
ranging from 0 to 20 Hz.  Keep in mind that the gain shows you how the height
(amplitude) of the bumps will translate into the height of the bounces for
different frequencies of bumps.  Since Mathematica automatically adjusts the range
of the vertical axis, the curves will be hard to compare unless all three are
graphed on the same plot.  But since Mathematica doesn't label the curves, it'll be
hard to tell them apart unless you first plot them individually!

What can you conclude about the behavior of each kind of car by looking at its
gain curve?  Are any frequencies especially good or bad?  Do any of the
suspensions seem particularly exotic?  Think about what the graphs tell you
before you look at our commentary.

Click here for commentary

Notice that the ``boat'' bounces wildly at low frequencies, with the gain
peaking at around 4.3 at a frequency of two Hz.  This means that if you
drive over a series of three-inch bumps at a rate of two/second, the car will
bounce up and down by over twelve inches.  You could easily lose control of
such a car if you hit the wrong kind of bumpy pavement.  Fortunately, this
behavior disappears at higher frequencies.

The ``truck,'' with its stiff suspension, has a gain of close to one for all
frequencies.  This means that the truck will bounce up and down by about the
same height as the bumps, regardless of how fast you're driving.  Of course,
this is about the same behavior that a piece of lumber or a brick would have if
it were sliding down the road!

The ``sedan'' has quite an interesting curve.  The gain steadily goes to zero
at seven Hz, and then goes back up again to a small non-zero value.  There's a
reason for this.  The human body is particularly uncomfortable when it is
vibrated seven times per second.  This designers of the ``sedan'' have
evidently taken than into consideration.

Phase Delay

Now that you know about gain, you can easily calculate the oscillation of a car
as it bounces down the road given the gain of its suspension and the function
describing the rate at which we hit the bumps in the road.  Suppose that the
road is described by the expression

A cos(2 Pi f t)

Remember that A characterizes the height of the bumps and f characterizes
their spacing.  Since we are measuring f in units of Hz, which are
oscillations per second, the value of the variable t should be specified in
units of seconds.

If the transfer function of the car is given by the function G(f I), what
expression will describe the car's bouncing?

Click here for the answer

We simply multiply the expression that describes the road by the gain of the
car's suspension:

A Abs[G(f I)] Cos[2 Pi f t]

Use what you've learned so far to plot, on the same graph, the oscillation of a
road and the oscillation of a car.  Let A = 1 and produce several plots for
different values of f and for different kinds of cars.  You can pick some
interesting values for f by looking at the plot of the overall gain for the
car.  On the resulting plots, the horizontal axis will be time (in seconds) and
the vertical axis will be amplitude.

What's unrealistic about the graphs that you have just produced?  Hint: no
physical system can react instantly to an input!

Click here for the answer

When you graph the road and the car, the two curves are perfectly in phase.
That is, each peaks when the other peaks.  In a real system, the oscillation of
the car would lag the oscillation of the curve.

In fact, the amount by which the oscillation of a particular car lags the
oscillation of a road depends upon the characteristics of the car's suspension
and upon the frequency of the bumps.  (In other words, the time lag depends
upon the same two factors as the gain.)  This time lag is known as the 
phase delay.

The phase delay is just as easy to compute as the gain if you know the transfer
function.  The phase delay at a given frequency f is the angle of
G(f I).  Mathematica  uses the function Arg[z] to give the angle of the complex number
z.

Graph the phase delay of each car's suspension over frequencies ranging from 0
to 20 Hz, much as you graphed the gain.  As with the gain, you will want to
graph them all on the same plot.

In these plots, the horizontal axis will be in units of Hz, but what are the
units of phase delay on the vertical axis?  A negative phase delay means that
the response (the bounce) will follow the signal (the bump).  It gives the
number of radians by which the response follows the signal.  Since there are
2p radians in a complete oscillation, a phase delay of --p, for
example, would mean that the response follows the signal by half an
oscillation.

Now try this.  If you know that the phase delay of a system at a frequency of
f Hz is -p radians, by how many seconds will the response trail the
signal?

Click here for the answer

A phase delay of -p radians means that the response will trail the signal by
the following fractional part of a complete oscillation:  p/(2 p).  At a
frequency of f Hz, each oscillation lasts 1/f seconds.  So the delay in
seconds will be   p/(2 f p).

For example, a phase delay of --p radians at 10 Hz translates into a delay of
 Pi/(2 10 p) = 1/20 of a second.

Putting it All Together

You already know that the formula describing oscillation is of the form

A cos(2 p f t)

You have been graphing a lot of such curves during this lesson.  In all of
these graphs, what is true when t = 0?

Click here for the answer

In all of them, the curve is at its high point (i.e., it has the value A)
when t=0.  This is a consequence of the fact that cos(0) = 1.

The formula describing a phase-delayed oscillation is

 A cos(2 p f t -- p)
 
where p is the phase delay.  Experiment with plotting this formula, letting
A and f be one, for various values of p between 0 and --2p radians.
It'll help to plot more than one curve on the same graph.  What happens when
p is --2p?

Now you know all that you need to know to graph the response of a particular
car to a particular bumpy road.  If the road is described by the formula

A cos(2 p f t)

and you know the transfer function G of the car, what formula will
describe the phase-delayed oscillation of the car?

Click here for the answer

We already know how to how to take the gain into account:

A Abs[ G(f I)] Cos[2 Pi f t]

To account for the phase delay, we simply compute and subtract the phase delay:

A Abs[G(f I)] Cos[2 Pi f t -- Arg[G(f I)]]

Now you can pick a frequency and a car and plot its response to the bumps in
a road.

Role of Visualization

We have been studying car suspensions in order to illustrate the role of
visualization in understanding the mathematics that describe a physical
situation.  Let's step back and think about what we've been doing.

The problem has been to determine how a car bounces up an down in response to
regularly spaced bumps in a road.  We were able to simplify the problem
somewhat because a car suspension is a linear system.  (We asked you to take
the simplification on faith.)

In the end, our problem was that given

* the frequency at which bumps occur,

we wanted to determine, for a given car,

* the ratio of the height of each bounce to the height of each bump (the
gain), and

* the amount of time by which each bounce trails the bump that caused it (the
phase delay).

We did the mathematics, which may or may not have been very revealing to you.
We could have graphed the overall results for each kind of car, but this would
have required a three-dimensional graph (since there are two outputs and one
input), which would have been hard to interpret.

In the end, we took advantage of the fact that the gain and the phase delay are
independent of one another and used separate two-dimensional graphs to
understand each individually.  These graphs have the advantages that they are
easy to understand and communicate the solution completely.

Plots of gain and phase delay as a function of frequency crop up in all fields
of engineering.  For example, they are used to describe the behaviors of stereo
amplifiers and the way that bridges vibrate in response to traffic.  In the
checkoff section, you will see how they can be used to understand one aspect
of the Earth's climate!

Step Responses

We have been doing a steady-state analysis of car suspension.  That is,
we have looked at what happens as a car rolls over a series of bumps, ignoring
what happens when the bouncing is just getting started after the first bump or
two.

A car designer might also be interested in what happens when a car hits a
single bump, such as a break in the pavement.  We can use
transfer functions to analyze this problem also.

Let's look at what happens when a car, rolling along flat pavement, rolls onto
a ``step.''  A step is a place where the road suddenly becomes a few inches
higher, such as a curb.

The calculations that we need to do involve Laplace transforms, which you will
learn about as you continue in engineering.  In order for Mathematica to perform
Laplace transform operations, you must first load the Laplace transform library
into Mathematica:

Needs["Calculus`LaplaceTransform`"]

Not all of the power ofMathematica is available to you immediately upon starting
Mathematica.  Much of Mathematica's more advanced capabilities are contained in 
packages, which you must explicitly load when you need them.  So, for
example, if you need to do Laplace transforms, you have to ask Mathematica to load
its Laplace transform library.  There are many other libraries available from
within Mathematica.

The buttons below will graph the vertical position of each kind of car as a
function of time.  We assume that the car rolls over a step of height 1 at time
0.  Each curve starts out at zero and stabilizes at one.  This is because,
after rolling up onto the step, the car ends up higher than it was before.
Notice that the time scales are different for each plot.

Plot[InverseLaplaceTransform[boat(s)/s, s, t], 
{t, 0, 20}];

Next the truck:

Plot[InverseLaplaceTransform[truck(s)/s, s, t], 
{t, 0, 20}];

Last the sedan:

Plot[ InverseLaplaceTransform[sedan(s)/s, s, t], 
{t, 0, 20}] ;

Labeling Plots

If you look at one of the plots that you've been producing during this lesson,
you'll notice that they are not labeled very well, indeed not at all as a default.  Neither the horizontal nor the vertical axis not labeled, and
the graph as a whole has no title.  It is possible to put better labels onto
plots, and it is important to do so when you hand in plots as part of your
homework assignments.  In this section we'll see how.

Let's suppose that you need to plot the function 2t^2, where
t is time measured in seconds and the result is altitude measured in feet.
Our first attempt gives a relatively uninformative plot:

Plot[2*t^2, {t, 0, 1}];

Fortunately, Mathematica's Plot[  ]  function supports a wide range of options for things like adding labels to the plot (indeed, if you are willing to put in the time it is possible to do almost anything in a plot).  

The options of relevance are:

AxesLabel 

PlotLabel

and

Frame Label 

(if you choose to put a ``frame'' around your plot.  These options are all specified with the Mathematica ``goes to'' symbol -> obtained by typing a hyphen and a geater than sign.  For example, in the above case we could write:

Plot[2*t^2, {t, 0, 1}, AxesLabel-> {Time, Altitude}, 
PlotLabel-> Altitude versus Time];

Execute this command and see how the plot comes out now.

Why might it be desirable to do this?  Click here for the answer.

Labeling your plots makes them clearer to those who have to read them. By itself a plot is little more than a line inside a rectangle, it is the numbers and labels on the axes, and the overall description of what is being plotted, that give a plot the power to convey a lot of information in succinct fashion. 

To learn more about giving your plots labels and names, read the Mathematica Book, page 418 ff. where there are many descriptive examples.

Printing Plots

The description of how to print plots  is extremely system dependent.  This section will have to be added by the individual instructor.

Mathematica  Summary

In this section we will summarize the new Mathematica  commands that you've
encountered in this lesson.

Plot is a Mathematica function for producing two-dimensional plots.  It takes
as arguments the expression to plot, the name of the variable contained by the
expression, and the low and high bounds of the horizontal axis.  Infinity and
negative infinity may be used as bounds.

Needs is used to read a package into Mathematica  so that you can have access
to advanced Mathematica capabilities..

Homework

Suppose that you had access to a table showing the average temperature in Salt
Lake City for every day for the last one hundred years.  You could easily
average over those hundred years to produce an average daily temperature
for each day of the year.  If you then plotted these averages by using the
horizontal axis for days of the year and the vertical axis for temperatures,
you would come out with an oscillation curve.  The high point would occur
some time in the summer, and the low point would occur some time in the winter.

Now suppose that you had access to a table showing the average length of each
day of the year in Salt Lake City.  If you graphed this, you would see a
similar curve.  The high point would occur on June 21, and the low point would
occur on December 21.

There is, of course, a relationship between these two curves.  To make an
analogy with what we've been studying in this lesson, the average day lengths
are like the bumps in a road, and the average daily high temperatures are like
the bouncing of a car.  We can model the climate of Salt Lake City as a linear
system that is responding to the oscillation in the amount of sunlight.  Notice
that the amount of sunlight and the average daily temperature oscillate
at the same frequency: one year.  This is what we'd expect from a linear
system.

You now know that to predict how a car will vibrate as it rolls over the bumps
in a road, you only need to know the car's transfer function.  Similarly, to
predict how the average daily high temperature will change over the course of a
typical year, you need to know a transfer function that characterizes the
climate of Salt Lake City.

Salt Lake City and San Francisco are at roughly the same latitude, which means
that they see roughly the same annual variation in day length.  But the two
cities have very different climates.  What this means is that Salt Lake City
and San Francisco have very different transfer functions!

A place that is far from the oceans, like Salt Lake City, is said to have a
continental climate.  A typical transfer function for a continental climate
might be

continental[s_]:=  1 / (4.4*s*10^6 + 1)

A place that is near an ocean, like San Francisco, is said to have a coastal
climate.  A typical transfer function for a coastal climate might be

coastal[s_]:= 1 / (3.7*s*10\^7 + 1)

Let's see how we can use these transfer functions to understand the difference
in the average high temperatures in Salt Lake City and San Francisco.

The transfer functions given above (just like all the transfer functions in
this lesson) assume that oscillations are given in units Hz.  What this means
is that if we were to calculate the gain or phase delay at a frequency of 1 Hz,
we would be looking at what the world would be like if the earth revolved
around the sun once per second.  Obviously, we are more interested in what
happens at frequency of 1 revolution/year rather than 1 revolution/second.

Write a short report answering the following questions.  When is
the hottest day of the year in the coastal climate?  What about in the
continental climate?  How much hotter is the hottest day of the year in the
continental climate than in the coastal climate?

Explain how you drew your conclusions by including pieces of your Mathematica session
into your report in appropriate places.  Your report should incorporate your
answers to questions 1 and 2 below, and should refer to the plots from
questions 3--6 below.  Be sure to label the plots sensibly.

(1) If a frequency of 1 Hz corresponds to one revolution/second, what
frequency corresponds to one revolution/year?  The answer will be a tiny
fraction of 1.  (Hint: there are 60 seconds in a minute, 60 minutes in an hour,
24 hours in a day, and 365.257 days in a year.)

(2) At a frequency of one revolution/year, as you calculated above, how
many radians corresponds to a phase delay of a single day?

(3) Graph the gain of the continental climate transfer function over
frequencies ranging from 0.25 revolutions/year to 4 revolutions/year.  (You'll
need to convert these frequencies back to revolutions/sec in order to determine
the proper ranges to plot.)

(4) Graph the gain of the coastal climate transfer function over frequencies
ranging from 0.25 revolutions/year to 4 revolutions/year.

(5) Graph the phase delay of the continental climate transfer function over
the same frequency range as in question 3.

(6) Graph the phase delay of the coastal climate transfer function over the
same frequency range as in question 4.