Derive Lab Manual-1
Derive Lab Manual-1
Derive Lab Manual-1
R & D Publishing
The procedures and applications presented in this supplement have been
included for their instructional value. They have been tested with care but
are not guaranteed for any purpose. The authors do not offer any warranties
or representations, nor do they accept any liabilities with respect to the
programs or applications.
c Copyright 2000 by the Authors.
All rights reserved.
No part of this publication may be reproduced, stored in a retrieval system,
or transmitted, in any form or by any means, electronic, mechanical, photo-
copying, recording, or otherwise, without the prior written permission of the
authors.
Printed in the United States of America.
http://www.math.hawaii.edu/CalcLabBook/
and to see our class web page containing our syllabus, assignments and gen-
eral information for students see
http://www.math.hawaii.edu/lab/
Honolulu, Hawaii
May 7, 2001
Contents
Preface ix
Calculus Reform and Computers . . . . . . . . . . . . . . . . . . . ix
How To Teach From The Manual . . . . . . . . . . . . . . . . . . . x
Advice For The Students . . . . . . . . . . . . . . . . . . . . . . . xiii
Setting Up The Computer Lab . . . . . . . . . . . . . . . . . . . . xiii
World-Wide Web Site For Our Lab . . . . . . . . . . . . . . . . . . xv
iii
1 Curve Sketching 25
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.2 Working with Graphs . . . . . . . . . . . . . . . . . . . . . . . 26
1.3 Exponential vs Polynomial Growth . . . . . . . . . . . . . . . 29
1.4 Laboratory Exercises . . . . . . . . . . . . . . . . . . . . . . . 32
2 The Derivative 37
2.1 The Derivative as a Limit of Secant Lines . . . . . . . . . . . . 37
2.2 Local Linearity and Approximation . . . . . . . . . . . . . . . 41
2.3 Laboratory Exercises . . . . . . . . . . . . . . . . . . . . . . . 44
4 Curve Fitting 55
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.2 Fitting Polynomials to Data Points . . . . . . . . . . . . . . . 55
4.3 Exponential Functions and
Population Growth . . . . . . . . . . . . . . . . . . . . . . . . 60
∗
4.4 Approximation Using Spline Functions . . . . . . . . . . . . 61
4.5 Laboratory Exercises . . . . . . . . . . . . . . . . . . . . . . . 63
8 Series 131
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
8.2 Geometric Series . . . . . . . . . . . . . . . . . . . . . . . . . 132
8.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
8.4 Approximating Infinite Series . . . . . . . . . . . . . . . . . . 137
8.5 Laboratory Exercises . . . . . . . . . . . . . . . . . . . . . . . 147
vii
5.7 Basins of attraction for x3 − 3x = 0 . . . . . . . . . . . . . . . 86
5.8 Bisection method for finding roots . . . . . . . . . . . . . . . . 88
ix
studied by using the computer to explore, to visualize and to suggest further
directions to study. We also try to convey that studying calculus can be
fun to do and that it is very important in understanding other topics in
mathematics and other fields.
This manual uses Derive because of its ease of use. Students enter ex-
pressions by filling out a form which includes special mathematical symbols
such as π and ∞. Then, at the click of a button one can differentiate, inte-
grate or plot the expressions. This calculator type interface is very easy for
students to learn and after about 15 minutes of introduction during the first
visit to the lab the students are ready to start using the program. There
are also functions that are equivalent to the various buttons and menu com-
mands. Knowing these functions enables one to write programs that extend
the power of Derive. We will give numerous examples demonstrating both
the simple calculator mode and the powerful programming language.
• Chapter 5 covers the Newton method for finding real roots to equations.
Emphasis is given to both the speed at which the method converges and
the care which is required in the selection of the starting point. There
is also a brief review of the algebra of complex numbers and applica-
tions of Newton’s method to finding complex roots to polynomials. A
discussion of fractals and chaos is given toward the end of the chapter.
• Chapter 10 shows how to plot polar and parametric curves. There are
discussions of rotations of graphs, speed of particles traveling along a
parametric curve and some connections with complex numbers.
Ralph S. Freese
David A. Stegenga
Part of our web page
Chapter 0
0.1 Overview
In this course you will learn to use the computer mathematics program De-
rive. This program, along with others such as Maple and Mathematica,
are very powerful tools for doing calculus. They are capable of doing exact
computations with arbitrary precision. This means that you can work with
numbers of any size or number of decimal places (most spreadsheets only
use 10-20 significant digits). These programs can simplify mathematical ex-
pressions by canceling common factors and doing other algebraic operations.
They can do symbolic calculus such as differentiation and integration, solve
equations and factor polynomials. When possible these programs solve these
problems exactly and when exact solutions do not exist, such as factoring
high degree polynomials or integration of some non-polynomial expressions,
then numerical methods are applied to obtain approximate results.
Probably the most important numerical technique is to graph and com-
pare functions. This will be a key feature of the labs. Typically we will
explore a topic by first graphing the functions involved and then trying to
do symbolic calculus on them using the insight gained from the picture. If
the problem is too difficult algebraically we then try numerical techniques
to gain further insight into the problem. It is this combination of graphics,
algebra and numerical approximation that we want to emphasize in these
labs.
Calculus is a hard subject to learn because it involves many ideas such as
slopes of curves, areas under graphs, finding maximums and minimums, ana-
1
2 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
lyzing dynamic behavior and so on. On the other hand, many computations
involve algebraic manipulations, simplifying powers, dealing with basic trig
expressions, solving equations and other techniques. Our goal is to help you
understand calculus better by concentrating on the ideas and applications in
the labs and let the computer do the algebra, simplifying and graphing.
Another important goal of the lab is to teach you a tool which can used
from now on to help you understand advanced work, both in mathematics
and in subjects which use mathematics. There are many features such as
matrices and vector calculus which we will not discuss but can be learned
later as you continue with your studies in mathematics, physics, engineering,
economics or whatever. Any time you have a problem to analyze you can
use the computer to more thoroughly explore the fundamental concepts of
the problem, by looking at graphs and freeing you from tedious calculations.
This chapter contains a brief introduction on how to use Derive. We
suggest you sit down at the computer and experiment as you look over the
material. Derive is very easy to learn thanks to its system of menus. The
few special things you need to remember are discussed below and can also
be found using the help feature in Derive.
This manual assumes that you will be using Derive for Windows (DfW)
Version 4. The recently released DfW Version 5 (DfW5) can also be used but
there have been some changes to the menus and buttons that we try to point
out using footnotes. See also Appendix B for more information on Version 5.
you click the Author menu and then click Expression. In this manual we will
indicate that two step combination by simply Author/Expression.1
In this manual we use a typewriter like font, eg., a(b + c) to indicate
something you might type in. We use a sans serif font for special keys on the
keyboard like Enter (the return key) and Tab. Most of Derive has easy to
use menus and buttons which we will described below. It only takes a few
minutes of practice to become capable at using DfW. You can skim over this
section now or then come back to it later as the need arises.
After clicking the button, you enter a mathematical expression, i.e., you
type it in and then press the Enter key or else click OK. You enter an expres-
sion using the customary syntax: addition +-key, subtraction --key, division
/-key, powers ^-key and multiplication *-key (however; multiplication does
not require a *, i.e., 2x is the same as 2*x). Derive then displays it on the
screen in two-dimensional form with raised superscripts, displayed fractions,
and so forth. You should always check to make sure the two-dimensional
form agrees with what you thought you entered (see Editing below to see
how to correct typing errors). Table 1 gives some examples.
If you get a syntax error when you press enter (or click OK) the problem
is usually mismatched parentheses. Carefully check that each left parenthe-
sis is matched with a corresponding right one. Also be careful to use the
round parentheses and not the square brackets since they are used for vector
notation; see Section 0.14 on page 18.
Note from (3) and (4) and from (6) and (7) of Table 1 that it is sometimes
necessary to use parentheses. Also note in (8), that to get the fraction you
want, it is necessary to put parentheses around the numerator and denomi-
nator. See what happens if you enter (8) without the parentheses. Also try
entering some expressions of your own. There are two ways to enter square
roots. One way is using the 0.5 or 1/2 power as in (9) and the other is to
enter the special square root character as in (10).
1
In DfW5, authoring has been improved and simplified but the above methods still
works. See Appendix B for details.
4 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
Table 1:
2
Leonhard Euler (oi’lar) was a 18th century Swiss mathematician.
0.5. EDITING 5
0.5 Editing
Suppose that you author an expression, click OK and then observe that you
typed something wrong. To enter a correction you would click the
button again and then click the right mouse button. A menu opens up with
several options, one of which is Insert Expression.3 Clicking this option puts
the highlighted expression of the current algebra window into the author box.
You edit this expression as you would in any windows program. That is you
position the cursor either by clicking or using the arrow keys. By highlighting
or selecting a subexpression and typing you replace the selected text with the
new text. One can use the Edit/Copy Expressions4 menu or Cntl-C to place
a highlighed expression from any algebra window onto the clipboard and
then in the authoring form right click the mouse and click Paste to copy the
clipboard contents. The simpler method of just right clicking the mouse and
then Insert is the best way as long as your expressions are in a single window.
There is also an option for inserting the expression enclosed in parenthesis.
A key equivalent to these techniques is to press either the F3 or the F4 key.
You select or highlight expressions in the algebra window by clicking on
them. For more complicated expressions you can click several times until the
desired subexpression is selected. This requires a little practice but you can,
x2
for example, select the x + 2 part of the expression sin x(x+2) by clicking on
it 4 times (each click takes you deeper into the expression).
The displayed expressions are numbered. You can refer to them as #n.
So, for example, with the expressions in Table 1, you could get sin(x)/x2 by
Authoring #5/#2.
When you start Derive it is in a character mode. This means it treats
each single character as a variable, so if you type ax Derive takes this to
3
In DfW5 you click in the Author box first and then right click Insert. A new feature
is to right click the erroneous expression, choose Edit and press the Enter key after you
make the corrections.
4
This is simply Copy in DfW5.
6 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
be a times x. This mode is what is best for calculus. The exception to this
are the functions Derive knows about. If you type xsinx, Derive knows
you want x sin(x). Actually on the screen you will see x SIN(x): Derive
displays all variables in lower case and all functions in upper case.
(x + 1)(x − 2). This means that the roots of f (x) are x = −1, 2, i.e., these
are the only solutions to f (x) = 0.
We can also do this by using the SOLVE function. To do this we highlight
the equation, say x2 −x−2 = 0 and click the . If you forget the function
of a button just hold the cursor on it and a brief explanation will appear.
An alternate method is to choose Solve/Algebraically from the drop down
menu. The quadratic formula is used to solve for the roots so it is possible
the answer will involve square roots (and even complex solution, e.g., x2 + 1
has no real roots but it does have two complex ones, namely, x = ±i).
If f (x) is not a quadratic polynomial then Derive may not be able to
factor it; nevertheless, it may be able to solve the equation f (x) = 0. As
an example, sin x = 0 has infinitely many solution x = mπ where m is
any integer. If we use Derive to solve this equation it gives the 3 solutions
corresponding to m = −1, 0, 1 (these are the principle solutions and all others
are obtained by adding or subtracting multiples of 2π).
Finally, the simple equation sin x − x2 = 0 cannot be solved exactly in
Derive although it is obvious that x = 0 is one solution and by viewing
the graph we see another one with x ≈ 1. In order to approximate this
solution we need to choose Solve/Numerically.6 You will then be asked for a
range of x’s (initially it is the interval [−10, 10]). Since there are (at least)
2 solutions in [−10, 10], you should restrict the interval to say [.5, 1] which
seems reasonable based on the graphical evidence. The result is that Derive
gives the solution x = .876626. We will discuss how this computation is done
later in Chapter 5.
Note that Solve/Numerically will only give one solution (or none if there
are none) even if the interval you choose cantains several solutions. To find
additional solutions you need to use Solve/Numerically again but with an
interval avoiding the first solution.
0.8 Substituting
√
If you have an expression like x2 + 1/x and you want to evaluate this
with say x = 3 or if you solved an equation f (x) = 0 and you want to
substitute that value of x back into f (x), you start by highlighting the desired
6
In DfW5 you choose Solve/Expression and then select the Numerically solution
method option.
0.9. CALCULUS 9
expression. Next you click the button and the substitution form opens
up. You need to fill in the substitution value so you would just type 3 in the
first example. On the other hand, if the substitution value is large, say lots
of digits or some other complicated expression in the algebra window, the
easiest way is to move the form out of the way (just hold down the left mouse
button in the top strip and drag to another location) and select the desired
expression by clicking on it. Then, paste it into the form by right clicking
and choosing Insert. If there happen to be other variables in the expression
you may have to change the variable in the variable list box.
0.9 Calculus
This menu item is very important for us. After choosing the Calculus drop
down menu, you get a submenu with Limit, Differentiate, Taylor series,
Integrate, Sum, Product, and Vector. There are also buttons for most of
these operators. After you have authored an expression, you can differenti-
ate it by either clicking the button or choosing Calculus/Differentiate
from the drop down menu. The form will have entries for what variable to
use and how many times to differentiate, but it usually guesses right so you
can just click OK. Then simplify.
To integrate an expression, first author it or highlight it if it is already
in the algebra window, then either click the button or else choose
Calculus/Integrate. The form will have entries for what expression to inte-
grate; it will guess you want to integrate the highlighted expression. It will
have an entry for what variable you what to integrate with; again it will prob-
ably guess right. It will also have entries for the limits of integration. If you
want an indefinite integral, just click the appropriate button and click OK.
For a definite integral click the appropriate button, type in the upper/lower
limit, then click OK. See Figure 0.2 on the next page for several examples
using Differentiate and Integrate on the Calculus menu.
The options Calculus/Limits is similar to the above. To find
x2 − 1
lim
x→−1 x + 1
you enter the expression, then either click or choose Limits from the
Calculus menu. You fill in the variable (which is x) and the limit point which
10 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
Y
3
xi = x · x2 · x3 = x6 .
i=1
0.10 Plotting
Supposed you want to graph the function x sin x. You simply author the
expression, by clicking the pencil button , to be plotted and then
click the button. A plot window will then opens up and the icon-
bar will change to a new set of buttons. You then click the button
again (it’s position is different in the plot window) and the graph will be
drawn. There are several different ways to view the algebra and plot window
together. The one we used to produce the pictures in the manual is to first
select the algebra window (if you are currently in the plot window you can
go to the algebra window by clicking the button) and then choose
Window/Tile Vertically from the menus. This will split the screen into two
windows: an algebra window on the left and a plotting window on the right.
These windows each have a number in their upper left hand corner. You can
have several plot windows associated with a single algebra window but you
cannot plot together expressions from different algebra windows. You can
switch windows by either clicking the top strip of the window or clicking the
or buttons. Actually you can click anywhere in the window to
select it but the top strip avoids changing the highlighted expression in the
algebra window or moving the cross in the plot window.
You can plot several functions in the same plot window. Move to the
algebra window, highlight the expression you want to plot, switch to the plot
window and then click the in the plot window. Now both expressions
will be graphed. You can plot as many as you want this way. The plot
window also has a menu option, Edit/Delete Plot, for removing some or all
of the expressions to be plotted. Pressing the Delete-key also removes the
current plot.
When you plot, there is a small crosshair in the plot window, initially at
the (1, 1) position. You can move the cross around using the arrow keys or
by clicking at a new location. The coordinates of the cross are give at the
bottom of the screen. This is useful for such things as finding the coordinates
of a maximum or a minimum, or where two graphs meet. In order to center
the graph so that the cross is in the center of the window, click the
button. This is useful for zooming in and out to get a better view of the
graph. There are several buttons for doing this in the plot window. Take a
look and you will see a button for zooming in, namely , and for zooming
0.10. PLOTTING 13
out and various ways of changing just the x-scale or just the y-scale.
You should try clicking these buttons to see exactly what happens.
In general, these buttons change the scale of the plot window by either
doubling or halving it. You can customize these by using the button7
(that’s a picture of a balance scale). Just click this button and fill out the
form the appears with your own numbers. You can see the current scale at
the bottom of the screen.
We mentioned above how to plot any number of graphs simultaneously by
repeatedly switching between the algebra window and the graphics window.
Another technique for plotting three or more functions is to plot a vector
of functions. This just means authoring a collection of functions, separated
by commas and surrounded by square brackets. For example, plotting the
expression [x, x^2, x^3] will graph the three functions: x, x2 , and x3 .
In order to plot a collection of individual points one enters the points as a
matrix, for example authoring the expression [[-2,-2], [0,-3], [1,-1]]
and then plotting it will graph the 3 points: (−2, −2), (0, −3) and (1, −1).
A quick way of authoring a vector is to use the button and a quick
way to enter a matrix is to use the button. One then just fills out the
form that opens up. So for example with the 3 points above we would click
the matrix button and select 3 rows and 2 columns. The form will
open up and we then fill in the 6 numbers above in the obvious order. You
move between fields by either clicking or using the Tab-key.
When plotting points you have a choice of connecting the points with a
line segment so that it appear like the graph of a function. You do this by
choosing Options from the menu bar. There are lots of interesting items on
this menu that will allow you to customize plotting colors, the size points are
plotted, axes and so on. To connect points we choose Points and then check
the Yes button.8 We can also modify the size of the points by clicking the
appropriate button. See the Figure 0.4 on the next page where each of these
techniques is demonstrated. Choosing Option/Plot Color controls the color
of a plot and whether the colors should rotate with each new graph.
7
This functionality has been replaced in DfW5 with the Set/Plot Region and Set/Plot
Range menu options.
8
The default is for points to be connected but if data points are to be plotted then it
is usually best to select No.
14 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
n! = n · (n − 1) · · · 2 · 1 n = 1, 2, . . .
where the properties of the Derive function IF(test, true, false) should
be pretty obvious. The definition forces the function to circle back over and
16 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
0.14 Vectors
Vectors are quite useful in Derive, even for calculus. They are also useful
in plotting. To enter the 3 element vector with entries a, b, and c, we can
Author [a, b, c] directly. It is important to note the square brackets which
are used in Derive for vectors; commas are used to separate the elements.
An easier approach is just click the button and fill in the three values
on the vector input form.
Derive also provides a useful function for constructing vectors whose
elements follow a specific pattern. The vector function is a good way to
make lists and tables in Derive. For example, if you Author vector(n^2,
n, 1, 3), it will simplify to [1, 4, 9]. The form of the vector function is
0.14. VECTORS 19
We will have other application that will require us to refer to the individ-
ual expression inside of a vector. This is done with the Derive SUB function
(which is short for subscript). Thus, for example, [a,b,c] SUB 2 simplifies
to the second element b. Derive will display this as [a, b, c]2 which explains
the name. For a matrix or vector of vectors then double subscripting is used
so that, for example, if
1 2
y :=
3 4
then Authoring y SUB 2 SUB 1 will be displayed as y2,1 and simplify to 3
(because it’s on row 2 and column 1).
expressions to an existing window. If you forget the name of your files just
type either A: or H: and press the enter–key to see a listing of your files.
During the course of your session with the computer you will make lots of
typing and mathematical mistakes. Before saving your work to a file or before
printing and turning your lab in for grading you should erase the unneeded
entries and clean up the file. The three buttons can be used for
this purpose. For example, if you select several expressions by say dragging
the mouse pointer over them with the left button held down and then press
the button these expression will be removed. Clicking will undo
the last delete. You can move a block of highlights lines by holding down
the right mouse button and dragging the block to a new location. Of course,
when you delete or move some lines then the line numbers will no longer be
in a proper sequence of #1, #2, . . . . You can correct this by pressing the
renumbering button . You should practice these commands on some
scratch work to make certain you understand them.12
One way to use the move command is to write comments in the file and
placing them before computation. Many of the *.MTH files that we wrote
for this lab manual use this technique. To do it, just author a line of text
enclosed in double quotes, for example, "Now substitute x=0.". Then,
move this comment to the appropriate location.
You can print all the expressions in the algebra window (even the ones
you can’t see) by pressing the button. You do the same thing to
print a graph. Just activate the plot window and press . Typically,
students turn in the labs by printing out the algebra window and penciling in
remarks and simple graphs. More extensive graphs can be printed out. Some
combination of hand writing and printouts should be the most efficient.13
0.16 Help
You can obtain on-line help by choosing Help. This help feature provides
information on all Derive functions and symbols. Suppose that you want
to know how to enter the second derivative of a functions f (x) by typing.
12
DfW5 has changed the methods for moving expressions, correcting expressions and
adding comments. See Appendix B for details.
13
The new functionality of embedded graphics in DfW5 allows for more flexibility in
these matters.
22 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
Q1. I tried to plot the line ax + b and instead I got an error message about
“too many variables”. What did I do wrong? You must define a, b
to have numerical values, otherwise Derive treats your function as
f (a, b, x) which it cannot plot.
Q2. I tried to plot the family of parabolas x2 + c and I got an error about too
many variables. What did I do wrong? Same problem as above, except
now Derive is trying to plot a surface z = f (x, c)in 2-dimensional
space. You probably want to enter and Simplify a vector of functions
such as
VECTOR(x^2 + c, c, 0, 4).
Q4. I solved for the 3 roots of the cubic x3 −2x2 +x−2 and I got x = 2 which
I guessed from the graph but the other two solutions were x = î and −î.
Where do these last two come from? If you factor the cubic instead of
using Solve you would get (x − 2)(x2 + 1). The complex solutions come
from that quadratic term. In calculus, we just ignore those complex
solutions. For example, numerically solving the above cubic will give
only real solutions.
Q6. I tried to author the inverse tangent function arctan x and I got a · r ·
c · tan x instead. What’s wrong? Derive recognized the tan x part but
treated the other symbols as individual constants. Use atan x.
Q8. I tried to author x^n and I got a syntax error! How was this possible?
The problem here is that either x or n is previously defined as a function.
For example, maybe you had authored x(t) := sin t. You can check on
this by scrolling up to find the definition. If instead, you know the
problem is that x(t) is defined and you want to remove that definition,
then just author x:=. In extreme cases you might just open a new
window and copy over some of your expressions using the Copy and
Paste technique.
Q9. I entered and simplified sin(2π) and I got SIN([2π]) instead, what hap-
pened? You authored sin[2pi] instead of sin(2pi). Derive treats
square brackets not as parenthesis but as a device for defining vectors,
see Section 0.14.
24 CHAPTER 0. INTRODUCTION AND DERIVE BASICS
Curve Sketching
1.1 Introduction
Before the widespread use of computers and graphing calculators, graphing
a function f (x) was done by a combination of techniques including:
• Plotting some judiciously chosen points.
• Locating the intercepts on the x-axis by finding solutions to the equa-
tion: f (x) = 0.
• Using Calculus to finding the local minima and maxima and where f (x)
is increasing and decreasing.
• Finding the points where the graph is curving upward and where it
is curving downward. The transition point in between is called an
inflection point.
• Finding the horizontal and vertical asymptotes.
Graphing is easier with a computer algebra system. In order to under-
stand the behavior of a function f (x) we can plot it using Derive’s plot
window (see Section 0.10 for instructions on plotting.). Moreover, we can
also find the local minima and maxima and the other items above if we need
them. We do this in Chapter 3. It is also possible to make a small change in
the function f (x) and then plot that change to see how the graph is affected.
We can also see important aspects of the graph by zooming in or out and
moving around to different parts of the graph. In this lab you will develop
your skills at graphing with the computer.
25
26 CHAPTER 1. CURVE SKETCHING
a sin(b(x − x0 ))
where a, b and x0 are given numbers. This raises the question of how the
graph of a function, such as sin x, changes when subject to the above mod-
ifications. You should observe the changes by comparing with the original
function but you should also think about why the changes make sense, for
example, what does changing a do, what is the geometrical significance of
the point x = x0 on the x-axis.
As our first task, lets plot y = f (x) + c where f (x) = sin x and where c
varies over the range −2, −1, . . . , +2. We start by clicking the pencil button
, entering sin x + (-2) and pressing OK. Next, click the button
1.2. WORKING WITH GRAPHS 27
and the Plot Window appears. Notice that there are a different set of buttons
at the top. The most important ones are the plot button , the various
zoom buttons such as which doubles the scale in each direction and
the cross buttons such as which centers the plot on the cross. The
cross is the little reference point whose coordinates are shown in the lower
left corner of the DfW window. You move the cross by either pointing and
clicking or else using the arrow keys.
Clicking the button will plot expression that is highlighted in the
Algebra Window. We return to the Algebra Window by clicking the
button and then entering the next expression sin x + (-1). Pressing
puts us back in the Plot Window and by clicking again we plot the
second expression in the same window as the first so that we can compare
them. The second graph appears to be the same as the first except that it
is shifted up one unit. We proceed in a similar manner with the remaining
three expressions with each successive graph shifting up another unit.
We like to have the Algebra Window and the Plot Window visible at
the same time. We do this by splitting the screen using the Window/Tile
Vertically command. Of course, you can work with only one window at a
time but it is easy to change the active window by simply clicking on the
name strip. You can actually click anywhere in the window to activate it
but this is not a good idea in the Plot Window since this will also move the
cross. Typically, the cross is carefully located at some important feature of
the graph and you don’t want to move it unnecessarily.
There are a few shortcuts that we can do to make this task a little easier.
First of all it’s not necessary to retype each expression over and over again.
Instead, we click on the expression we want to modify (this causes it to be
highlighted) and then we press . We then right click in the author
form and a menu opens up. We choose the Insert expression option and the
highlighted expression appears in the author box. We then just change the
-2 to -1 and press OK. That does it and we can now plot the new expression.
A very fast way of plotting similar expressions all at once is to uses
vectors 1 . Select the Calculus/Vector menu and enter sin x + c in the form.
Click the Variable list box and select the variable c. For this example, set
the Starting value to -2 and the Ending value to 2. The Step Size can be left
1
See Section 0.14 on page 18 for more information about the vector function
28 CHAPTER 1. CURVE SKETCHING
at 1 although in other examples you might want to change this. Press OK.
Next we need to simplify the resulting expression
Warning You can change the function definition by simply Authoring f(x)
:= with a new expression. One problem with making definitions is that you
tend to forget them, especially as you move on to a new problem. Thus, if you
use the same letter f (x) in every problem you will definitely have problems
and get confused as to the current definition of the expression f (x). A good
strategy is to use definitions as little as possible and only when necessary.
Use different letters, say f1(x):= for the first problem and f2(x):= for the
second, and so on. You might also want to eliminate a definition as soon
as you are done with it. It is important to note that once you define a
function by this method it will not go away if you simply erase that line
from your algebra window because it is in the computer’s memory. The
way to completely remove a definition using the letter f is to author the
expression f:=. This gives f an empty definition.
1.3. EXPONENTIAL VS POLYNOMIAL GROWTH 29
You might note that after you have defined a function in Derive it
thereafter is displayed with uppercase letters. This is just a visual cue and
it is not necessary to type expressions using uppercase letters2
Unfortunately, Derive can not solve all equations (nor can anybody else!)
and so we need to also consider approximate numerical solutions. The tech-
nique for doing this uses calculus and will be discussed in Chapter 5.
For example, our equation x4 = ex cannot be solved exactly. If you
try then Derive just returns the empty solution vector []. We can get
approximate solutions if we use the Solve/Numerically3 menu. You enter
the equation x^4 = e^x and then choose an appropriate interval which you
believe contains the desired solution. This choice is usually made by looking
at the graph. For example, the negative solution in Figure 1.2 appears to
be in the interval [−1, 0]. Making the interval larger will not increase the
number of solutions since the answer, unlike exact equation solving, has at
most one answer. If there are no solutions in the interval, then you need to
specify a new interval. In our example, the other solution in Figure 1.2 can
be found if we choose the interval from 0 to 2.
3
In DfW5 you choose Solve/Expression and then select the Numerically solution
method option.
1.3. EXPONENTIAL VS POLYNOMIAL GROWTH 31
Are there any other solutions? It is pretty clear that there are no other
solutions for x < 0 but what about large x? From the graph it appears that
x4 grows much faster than ex . But of course ex has “exponential growth” so
perhaps ex > x4 for large enough x. To test this we try zooming out several
times in the Plot Window by clicking . But both functions grow very
quickly and it is a little difficult to be successful with this method. Back in
the Algebra Window we can try to numerically solve x4 = ex for larger x.
If we choose the interval to be 2 to 20, we get the solution x = 8.61316.
So the graphs cross at this point. To find the y value of this point, we
use the substitute button to substitute 8.61316 into x4 . The result
approximates to 5503.64.
To see this on the graph we need to zoom out once so that the x–scale
includes x = 8.61316. Then we need to zoom out on the y–axis without
zooming out on the x–axis. The need to do this is why there are those
extra zoom buttons. After zooming out several times we obtain the graph
of Figure 1.3. Another approach is to use the Set/Scale menu and enter an
x-scale of 5 and a y-scale of 5000.
There are a couple things this demonstration shows. First that in order
to see the important features of a graph it may take some skill at moving
around and manipulating the scale of the graph. Moreover, even though we
can clearly see the two graphs intersecting at x = 8.61316 in Figure 1.3, we
can no longer see the other two solutions. So it may not be possible to see
all the important features in one plot. In the exercises you will learn how to
move and scale in the plot window and to use the algebra window in order
to find all the important features of one or more graphs.
a. 81/2
b. sin( π4 )
c. sin( π4 )/51/2
3. Graph cos x, 2 cos x, and cos(2x) and explain what the transformations
f (x) → f (ax) and f (x) → af (x) do to the graph of f (x).
a. Graph g(t).
b. What sort of transformations should be applied to sin t to make
its graph look like the graph of g(t)?
c. Use Derive’s crosshair to find the maximum value of g(t) and to
find the first root of g(t) to the left of zero.
d. Use these numbers to find a and b so that g(t) = a sin(t + b), at
least approximately.
*e. Find exact values of a and b so that
a. Plot f (x) for a = −2, −1, 0, 1, and 2. You can use the vector
function if you like but plot each expression one at a time so you
can see its graph.
b. For each of the five values of a determine whether the graph is
curving upward, curving downward, no curving at all or curving
both ways. If the graph has both upward and downward curving
graphically approximate the points where there is a transition
between upward and downward curving. These are the inflection
points. Print out the Plot Window and mark it up with your
pencil indicating which graph corresponds to which value of a,
inflection points, etc.
7. Find the points where the curves ln x and x1/4 intersect. Make two (or
more) graphs with different scales showing the places where the curves
intersect.
your graphs show the main features such as the x and y–intercepts, the
critical points, and the inflection points.
3x
a. sin(x) cos(20x) b. √
4x2 + 1
1
c. d. x sin(1/x)
1 + 5000(x − 1)2
*10. Let
−2x3 + 6x2 − 3x + 5
g(x) =
4x2 − 6x − 7
a. Graph g(x) so that your graph shows the main features of this
function.
b. This graph has a slant asymptote, i.e., an asymptote which is a
line with nonzero slope. Zoom out a few times until you can see
this slant asymptote.
c. Find the formula for the slant asymptote by using Simplify/Expand.
*11. In reading this chapter you might have wondered if ex and x4 intersect
some place beyond c = 8.61316 · · · . You could use Derive to verify
1.4. LABORATORY EXERCISES 35
that there is no solution say between 8.7 and 100 and this would be
strong evidence that they don’t intersect beyond c, but not a proof. So
in this problem you are to find a proof that ex and x4 don’t intersect
beyond c (without using Derive). Hint: By taking 4th roots we must
show ex/4 > x for all large x. Now using Calculus, show that the slope
of ex/4 − x is positive for all x ≥ 8, i.e., the derivative is positive. Use
this to show ex > x4 for all x > c.
36 CHAPTER 1. CURVE SKETCHING
Chapter 2
The Derivative
F(x) := x^3/3
SL(a, h) := f(a) + (f(a+h) - f(a))/h (x - a)
The first step defines f (x) to be the function x3 /3. The second defines a
function SL(a, h) which gives the secant line through the points (a, f (a))
and (a + h, f (a + h)). For example, if we Simplify SL(1,1) we get 7x−63
so
that the equation of the secant line determined by the points (1, f (1)) and
(2, f (2)) is y = 7x−6
3
.
37
38 CHAPTER 2. THE DERIVATIVE
If the drawing is too quick to see the animation, try the following method
instead. Erase the 5 secant lines in the plot window by pressing the Del key
5 times. In the algebra window select an individual line in the vector by
repeatedly clicking on it. Then, activate the plot window and press .
Finally, repeat this process several times to see the pattern evolving in the
plot window; namely, that the lines are rotating about the point (1, 13 ) on
the graph. Moreover, the lines appear to be tending to the tangent line.
In Figure 2.1 on the preceding page this pattern is clearly shown: the
secant lines tend to the tangent line by rotating in a clockwise manner, i.e.,
with decreasing slope. We can use Derive to illustrate this effect using
calculus. That is, as h tends to 0, the secant line tends to the tangent line
by taking the limit: Author SL(1,h) and choose Calculus/Limit, taking the
variable to be h (not x). After Simplifying we get
3x − 2 2
lim SL(1, h) = =x− .
h→0 3 3
3
which, in fact, yields the tangent line to x /3 at x = 1. Check this out for
yourself by plotting this function on your previous graph. Since the slope
of the secant line is (f (a + h) − f (a))/h, this explains why we define the
derivative as
f (a + h) − f (a)
(1) f 0 (a) = lim
h→0 h
and why the derivative is the slope of the tangent line.
A quicker way of entering the 5 lines above is to observe that there is
a pattern to the values 1, 12 , 14 , . . . ; namely 21n for n = 0, 1, . . . , 4. We can
use this fact and the VECTOR function on the Calculus menu to simplify this
task. Select Calculus/Vector and enter SL(1,1/2^n) in the form. Note that
using uppercase letters is not necessary and that the highlighted expression
will be replaced with whatever you type. For the Variable, scroll down and
select n. Next we take the Starting value to be 0 since 20 = 1 and the Ending
value to be 4 since 24 = 16. Click OK and simplify the resulting expression
VECTOR(SL(1, 2^-n), n, 0, 4). The result is a vector of five secant lines
as above. You will find that this is a convenient method of producing a large
number of expressions without typing them individually.
We can later change the definition of f (x) to a different function and
use the SL(a, h) function to get secant lines to the new function. The file
F-SECANT.MTH contains the definitions of SL(a, h) and the tangent line
function, TL(a), discussed below.
40 CHAPTER 2. THE DERIVATIVE
Tangent Lines In order to get a function TL(a) for the tangent line at a
analogous to the secant line function SL(a,h), we need to be a little careful
since the most obvious definition; namely,
TANGENT(sin x, x, 0) .
The answer is y = x. You should verify this visually by plotting sin x and x
together on the same graph. Now center the cross at the origin and zoom in
several times. On the other hand, using the formulas
d
sin x = cos x and sin 0 = 0, cos 0 = 1
dx
we can easily derive the above formula for the tangent line.
The point is that we now have a powerful tool for approximating the
values of sin x as long as x is small; namely,
sin x ≈ x for x ≈ 0 .
Looking back at our sin( 13 ) example above, we see that this simple technique
has given a result which is accurate to roughly 2 decimal places.
We can also use the above approach to approximate the derivative of a
function and plot the result. For example, we know that the derivative of
f (x) = x3 is 3x2 by using the standard formulas. On the other hand, the
function of x, g(x, h) = f (x+h)−f
h
(x)
with h fixed at some value like h = .01
2.2. LOCAL LINEARITY AND APPROXIMATION 43
is a good approximation to 3x2 as one can see from Figure 2.2. The figure
actually shows both plots although they appear to be only one curve. In
Derive you should enter and Simplify the above expression (it sometimes
helps to Expand the result to further simplify it). Then compare the graph
with 3x2 by plotting both expressions together.
In this case, we can easily estimate how the close the two graphs are by
algebraically simplifying (enter in Derive and press ) to get:
(x + h)3 − x3
− 3x2 = 3hx + h2 .
h
Now for h = .01, the right hand side above no larger than .03|x| + .0001 and
as long as we take −1 ≤ x ≤ 1 then this error term is approximately .03 and
our accuracy is nearly 2 decimal places.
Let us try this technique with sin x. Maybe you haven’t learned yet that
the derivative of sin x is cos x. Here is way to strongly suggest this result
44 CHAPTER 2. THE DERIVATIVE
x 1
a. 1− b. 2x − x2
2 2
1 2
c. x −x d. sin x
2
2. a. Using the TANGENT(u,x,a)
√ function find the equation of the tan-
gent line for f (x) = 3 x (enter cube roots as x^(1/3)) at the point
a = 8 and plot it along with the graph of f (x).
√3
b. In
√ part a you found the tangent line to x at a = 8. Estimate
3
9 by finding the y–value of this line when x = 9. Compare your
answers with DfW’s own approximation to this quantity obtained
by clicking the button.
2.3. LABORATORY EXERCISES 45
3. The acceleration due to gravity, a, varies with the height above the
surface of the earth. If you go down below the surface of the earth, a
varies in a different way. It can be shown that, as a function of r, the
distance from the center of the earth, a is given by
GM r3 for r < RAD
RAD
a(r) =
GM for r ≥ RAD
r2
where RAD is the radius of the earth, M is the mass of the earth, and
G is the gravitational constant. All three of these are constants. In
order to define the function a(r) and examine its graph, we’ll use the
numerical values3 : GM := 4.002 × 1014 and RAD := 6.4 × 106 meters.
a. Define a(r) using the technique in Section 0.13 and plot its graph.
Rescale as necessary to give a good picture4 .
b. Is a a continuous function of r?
c. Is a a differentiable function of r? Explain your answer.
d. Refine your answer to Part a so that the domain of a(r) is r ≥ 0.
a. Graph Γ(x) and the four secant lines to Γ(x) through the points
(3, Γ(3)) and (3 + h, Γ(3 + h)), for h = 1/2, 1/4, 1/8, and 1/16.
[It is known that Γ(3) = 2, but you don’t really need this here.]
3
Simply Author the expression GM:=4.002 10^14 and similarly for RAD.
4
The horizontal scale should reflect the value of RAD and the vertical scale can be
determined from the numerical value of a(RAD)
46 CHAPTER 2. THE DERIVATIVE
*5. Let f (x) = x2 sin(1/x) for x 6= 0 and f (0) = 0. In this problem you will
show that f (x) is continuous and differentiable for all x but f 0 (x) is not
continuous at 0. This means to find f 0 (0) you must use the definition of
the derivative; you cannot just find f 0 (x) and take the limit as x → 0.
3.1 Introduction
Calculus is a beautiful and important subject. It derives its importance
from its ability to describe and model basic phenomena in so many fields.
Besides physics, chemistry and engineering, it is used in biology, economics,
and probability. In order for calculus to be useful to you, you will need
to understand calculus graphically and numerically as well as algebraically.
Algebraically you learn how to differentiate functions given as complicated
expressions. But you also need to understand the derivative visually as the
slope of the tangent line and physically as a rate of change of y with respect
to x.
With Derive it is easy to learn all three of these aspects and to see the
relations between them.
47
48 CHAPTER 3. BASIC ALGEBRA AND GRAPHICS
Looking further at the graph we can see that f (x) does not have a local
minimum or maximum at x = 0; in fact f 0 (x) < 0 on both sides of 0. The
graph also shows that x = 9/8 is where the minimum occurs and that f (x) is
decreasing on (−∞, 9/8] and increasing on [9/8, ∞). If we highlight the first
3.3. MAX-MIN PROBLEMS 49
• First we are only interested in those points x which lie on our interval
[a, b]. It probably will be necessary to obtain decimal approximations
to our critical points x to see if they satisfy a ≤ x ≤ b.
• The second problem is that the absolute maximum may not occur at
the critical points at all. Instead, it occurs at either the left or right
endpoint, i.e., x = a or x = b. This situation should be clear from the
graph.
• The third problem is that perhaps Derive can not solve the equation
f 0 (x) = 0. As you look at the graph, you should be able to guess the
critical point, at least approximately. Then you can ask Derive to
approximate the critical point by using the Solve/Numerically1 menu.
You will need to specify an interval on which Derive will search for
a solution. But this should be obtained by again looking the graph
and choosing a convenient interval containing the approximate critical
point. Choosing too big of an interval can get you into trouble because
Derive may end up finding a solution which is not the one you want.
You would then need to refine your interval to exclude the undesired
critical points.
Now plot this. Zoom out by clicking and see if it appears that g(x) has
a horizontal asymptote. A nice technique to do this is to leave the vertical
scale alone and zoom out in the horizontal direction. There are several ways
to do this; one way is use the zoom buttons on the menu bar. Can you guess
1
In DfW5 you choose Solve/Expression and then select the Numerically solution
method option.
3.4. ZOOMING AND ASYMPTOTES 51
which button does this? Another way is change the scale, say x = 100 (click
the button with the balance scale on it). Use the cross-hair to estimate the
value of y that g(x) is tending to for large x. You should get y = 3 (see
Figure 3.2).
(The answer above is too wide for the window to display completely so you
will need to scroll to see the 3. Use the horizontal scroll bar at the bottom
of the window to do this.)
x 1
a. 1− b. 2x − x2
2 2
1 2
c. x −x d. sin x
2
a. Enter the expression for f (x) and plot it. Using the crosshair as
described in Section 0.10, find an approximation for the both the x
and y–coordinates of the local minimum of f (x).
b. Using the same method find an approximation for the unique x
satisfying the equation f (x) = 0, i.e., the value of x where the
graph crosses the x–axis.
c. In the algebra window click or use Solve to find the root
exactly. Approximate this as a decimal and see how it compares
with your answer in Part (b).
3.5. LABORATORY EXERCISES 53
d. In the algebra window find the derivative f 0 (x) and use it to find
the exact coordinates of the local minimum you approximated in
Part (a). (This will give you the critical value x; to get the y–
coordinate substitute the value of the x–coordinate into f (x) by
clicking or by using the Simplify/Substitute/Variable menu;
see Section 0.8.)
3. Consider the class of all function of the form f (x) = x3 + bx2 + cx.
a. Author the expression form x3 + bx2 + cx. Click the
substitution button and enter some specific values for b,
2
c, then plot the result . Do this for several different choices for b
and c and observe the critical points and inflection points of the
different graphs.
b. Using Derive’s calculus facilities in the algebra window, show
that the function f (x) = x3 + bx2 + cx always has exactly one
inflection point, regardless of the values of b, c.
c. Again using Derive’s calculus facilities, show that f (x) can have
either 0, 1 or 2 critical points. Determine for what values of b and
c does f (x) have no critical point?
d. Choose values b, c which demonstrate that f (x) may have either
0, 1 or 2 critical points and plot their graphs.
x 7
4. Graph the function f (x) = ( 1+x 2 ) . At first there appears to be no
part of the graph showing in the graphics window but this can not be
since f (0) = 0. Try replotting the graph in another color by either just
clicking several times or using the Options/Plot Color menu.
Now the graph appears to the horizontal line y = 0 but this can not be
since clearly f (x) = 0 only for x = 0.
c. In the Graphics window, use the Zoom buttons or else the Set/Range
menu in such a way that both the local maximum and local min-
imum points are visible. Furthermore, make the y-scale compara-
ble to the y-coordinate of the local maximum.
d. After you get a good looking graph, print out the result.
*6. Suppose we have the situation of the previous problem except that now
the metal for the top and the bottom of the can costs 1.5 times as much
as the metal for the side. What is h/r for the can of minimun cost?
4.1 Introduction
Consider the population of a certain country, P (t), as a function of time. We
may not know exactly what P (t) is but instead just have a table of data,
for example, the population at the beginning of each year for the last few
years. We are interested in finding an appropriate curve for the data. We
might try comparing the data against a linear function ax + b, a quadratic
function ax2 + bx + c or say an exponential curve of the form P (t) = aert .
Under a certain model of population growth P (t) will have this last form.
Our problem is to determine the parameters a, b, . . . , from the data. Once
we do this then we can use P (t) to estimate the population at times between
the data and predict the population in the future.
This kind of problem of fitting a function from a family of functions to
numerical data arises frequently in many applied areas including statistics.
In this lab we use the computer to help visualize data and fit the data to a
function from a class of functions. We begin with the class of all polynomial
functions.
55
56 CHAPTER 4. CURVE FITTING
let’s consider the problem of finding a polynomial function f (x) which goes
through these points. That is, we want f (x) to satisfy f (xi ) = yi for i =
1, . . . , n. If xi 6= xj , for i 6= j, then it turns out that we can always find a
unique polynomial of degree at most n − 1 going through these points. This
is quite obvious when n = 2 since there is a unique line passing through any
2 distinct points.
If we are given 3 points, (x1 , y1), (x2 , y2), and (x3 , y3 ), and want to find a
quadratic polynomial passing through these points, we let f (x) = ax2 +bx+c
be an arbitrary quadratic. Since f (xi ) = yi for i = 1, 2, and 3, we obtain
three (linear) equations
ax21 + bx1 + c = y1
(1) ax22 + bx2 + c = y2
ax23 + bx3 + c = y3
in the unknowns a, b, and c. (Remember, we are given the points (xi , yi) so
they are known and we want to find the unknowns a, b, and c.) We then
solve this system of 3 equations for the 3 unknowns a, b, and c.
For example, suppose we want to find a quadratic polynomial f (x) =
ax2 + bx + c passing through (0, 0), (1, 2), and (2, 8). The way to do this with
DfW is to first author f(x) := ax^2 + bx + c then choose Solve/System
from the menu bar, set the number of equation to be 3, and then enter the
three equations (you can either use the Tab-key after entering an equation
or click the next equation box)
f(0) = 0 f(1) = 2 f(2) = 8
Click on the Equation Variables box and select the variables to solve for as a,
b, and c. Click OK and then simplify the resulting expression SOLVE( [F(0)
= 0, F(1) = 2, F(2) = 8], [a, b, c]) (see Section 0.7 on page 6). De-
rive returns
[a = 2 b = 0 c = 0]
So in this case f (x) = 2x2 .
We can double check this result by plotting the function 2x2 determined
above along with the 3 × 2 data matrix
0 0
1 2
2 8
4.2. FITTING POLYNOMIALS TO DATA POINTS 57
use single characters for our variables. But to solve (1) we need the variables
x1, x2, etc. When we enter x1, Derive will think of this as x times 1, which
is not what we want. So we need to declare that x1, x2, etc., are variables.
To declare x1 as a variable you can author x1 :=. You need to do this for
all three x’s and y’s. Be sure to use the assignment characters := and not
the equation character =.
As a convenience for the exercises the ADD-HEAD file has provided dec-
larations for x0, . . . , x4 and y0, . . . , y4. But if you need more variables or
if you use a different letter other than x or y then you will need to declare
them.
Now enter the data matrix by clicking the button and selecting 3
rows and 2 columns. Then, enter x1 press Tab and enter y1. Continuing,
enter all the remaining xi ’s and yi ’s. Click the OK button and the result
should be
x1 y1
x2 y2
x3 y3
f (x1 ) = y1
(2) f (x2 ) = y2
f 0 (x3 ) = y3 (That’s the derivative!)
In other words we specify that f (x) must pass through (x1 , y1) and (x2 , y2 )
and that its slope at x = x3 is y3 . We define f(x) as before and define g(x)
4.2. FITTING POLYNOMIALS TO DATA POINTS 59
as the derivative1
as before and then factor the answer, we see that the denominator of each of
the 3 fractions is
(You can do this just as in Figure 4.2 except you need to define g(x) and
use g(x3) = y3 in place of f(x3) = y3.) Of course we are assuming that
1
Note that you can’t simply define g(x):=DIF(f(x),x). We ran into this problem
earlier on page 41 (see Section 0.12 on page 16 for a complete explanation). One solution
to this problem is to use the utility file and define g(x):= SUBST( DIF(f(u),u), u, x).
60 CHAPTER 4. CURVE FITTING
CURVEFIT(x,[[x1,y1],[x2,y2]],[[x3,y3]])
Population models are studied more thoroughly in Chapter 11 using the the-
ory of differential equations but for now we will just consider the exponential
model. Here P (t) is the population at time t and a is the population at the
starting time, t0 . Problem 7 uses this model.
There are two parameters in (4), a and r. These parameters can be
determined if we know the population at two different times, t1 and t2 , i.e.,
if we know P (t1 ) = y1 and P (t2 ) = y2 . This gives the equations
aer(t1 −t0 ) = y1
aer(t2 −t0 ) = y2
∗
4.4. APPROXIMATION USING SPLINE FUNCTIONS 61
but solving for a and r is a little more difficult since this is not a linear
system of equations. Hence, we can not solve this system using the function.
Instead, the way to do it is to use the first equation to solve for a and then
substitute that value into the second equation and then solve the resulting
equation for r. You can use SOLVE on any single equation but not systems.
Another approach is to observe that the equations are linear in the quan-
tities ln a and r because, if we let c = ln a, they are equivalent to:
c + r(t1 − t0 ) = ln y1
c + r(t2 − t0 ) = ln y2
Of course, once we find c then a = ec , so you’re done. Problem 6 will require
solving these equations.
This solution is
where we note that each data point can be referred to as data SUB i or
alternately, using the symbol bar as data↓i.
Now to find our second quadratic f2 (x) connecting the second pair of data
points and making sure that the graph of the two functions is smooth at x = 1
we simply solve for f2 (x) using the equation: f10 (1) = f20 (1). Continuing in
this way we get quadratics f1 (x), f2 (x), . . . , fn−1 (x) corresponding to each of
the n−1 intervals: [data1,1 , data2,1 ], [data2,1 , data3,1 ], . . . , [datan−1,1 , datan,1].
Note that we have used the double subscript notation to get the x-values in
the first column of the data matrix.
We combine these functions into a single function using the CHI function.
Here CHI(a,x,b) is 0 unless a ≤ x ≤ b in which case it is 1. Thus, the
combined function is
X
3
(5) f (x) = fk (x) · CHI(datak,1, x, datak+1,1 )
k=1
In our example we can solve for the three quadratics and get
The definition behind the SPLINE function (see the file ADD-UTIL) is
fairly straightforward. The function fk (x) depends on the previous function
0
fk−1 (x) and more specifically on the quantity fk−1 (xk ), where the k th interval
is [xk , xk+1 ]. It turns out that it is more efficient to make a vector out of the
n − 1 slopes m = [m1 , m2 , . . . ] using the formula
yk+1 − yk
(6) mk = 2 − mk−1 k = 2, . . . , n − 1 .
xk+1 − xk
which can be derived using DfW(see the file F-SPLINE.MTH). Using this
formula one produces the vector of slopes using the ITERATES function. The
formula looks a little complicated at first but should look straigtforward
after some careful examination (see either the file ADD-UTIL or the SLOPE
function in the file F-SPLINE.MTH).
Using these quantities one then computes fk (x) using
.
See Figure 4.3 where we use this method of approximation to approximate
the function y = sin x using n = 7, which is the smallest integer greater than
2π. Thus, based on the numbers sin 1, . . . , sin 7 plus the derivative at x = 0,
i.e., m0 = 1, we get a good approximation to the sine function.
f (1) = 1
f (3) = 4
f 0 (2) = 2
(see equation 2). Use the second form of the CURVEFIT function to
find the solution. Note that the ddata is a 1×2 matrix in this case.
What does Derive tell you? What if you change the derivative
to f 0 (2) = 3/2? Can you explain what the answer means?
4.5. LABORATORY EXERCISES 65
For the following problems you will need to use the variables x1, x2, x3, x4,
y0, y1, y2, y3, and y4. As long as you have loaded the ADD-HEAD file these
variable will be declared. On the other hand, if you use a variable such as
z0 then you will need to declare it by authoring z0:=. See the discussion on
page 57.
3. Let (x1 , y1) and (x2 , y2 ) be two points in the plane with x1 6= x2 . Let
m = xy22 −y
−x1
1
be the slope of the line through these points. The Mean
Value Theorem says that if f (x) is a differentiable function which passes
through these points then f 0 (x3 ) = m for some x3 between x1 and x2 .
Show that if f (x) has the form ax2 + bx+ c then we can take x3 = (x1 +
x2 )/2, i.e., show that f 0 ((x1 + x2 )/2) = m if f (x) has this form. Hint:
Solve the system f (x1 ) = y1 , f (x2 ) = y2 for b and c and substitute
those values back into ax2 + bx + c. Then show that the derivative of
the resulting expression is m when x = (x1 + x2 )/2. Of course this
means that all quadratic functions through (x1 , y1 ) and (x2 , y2 ) have
the same slope at (x1 + x2 )/2.
4. Use the CURVEFIT function to find the quadratic function f (x) = ax2 +
bx + c that satisfies
f (x0 ) = y0
f (x1 ) = y1
f (x2 ) = y2
Integrate the resulting function over the interval [x0, x2]. Observe that
your answer is a pretty big expression that requires scrolling to view.
Now substitute in this expression x1 = (x0+x2)/2 using the but-
ton and simplify. Note that x1 is the midpoint of the interval [x0, x2].
The answer should be a very simple formula in terms of x0 , x2 , y0 , y1
and y2 . In the next chapter this calculation will be the basis for the
Simpson Method of numerical integration.
such that
f (x1 ) = y1
f (x2 ) = y2
f (x3 ) = y3
f 00 (x4 ) = y4
Show that this is always possible if x1 , x2 , and x3 are all distinct and
x4 6= (x1 + x2 + x3 )/3. The algebra in this problem gets fairly messy.
6. Let (x1 , y1 ) and (x2 , y2 ) be two points in the plane with y1 > 0, y2 > 0,
and x1 6= x2 . Let f (x) = aerx be an exponential function. Show that
it always possible to find a and r so that f (x) passes through these
points. Hint: you need to solve the equations
aerx1 = y1
aerx2 = y2
To do this first solve for a in one of these and substitute the answer
into the other and then solve for r. Make sure that your formulas for
a and r only involve the data points.
7. Table 4.1 on page 68 shows the population of the US for every decade
from 1800–1900. Consider two models for the data: an exponential
model P (t) = aert and a linear model L(t) = bt + c.
a. Use the data for 1800 and 1810 to determine 2 linear equations in
the 2 unknowns b and c. Use to find b and c.
b. Do the same thing for the a and r but this time you will need to
apply the previous exercise to solve for a and r.
c. Show that your model in Part (b) is the same as 5.3er(t−1800) where
r is the value you obtained in Part (b). Use properties of the
exponential function to explain why this is true.
d. What does each model predict for 1830?
e. How do the models compare during the first 50 years? 100 years?
Do this by graphing both functions and the population data. Ad-
just the scale to get a good picture. Note: The data in Table 4.1
can be copied from the file F-population.
4.5. LABORATORY EXERCISES 67
*9. Let (x1 , y1), (x2 , y2), and (x3 , y3 ) be three points in the plane
with x1 6= x2 , x1 6= x3 , and x2 6= x3 . Show that all cubic functions,
f (x) = ax3 + bx2 + cx + d which go through all three of these points
have the same second derivative at x1 +x32 +x3 . (Hint: Just solve the 3
equations in the 3 unknowns b, c, and d in terms of the 4th unknown a.
Differentiate twice and substitute in the above value of x. Check that
the answer does not depend on a.)
68 CHAPTER 4. CURVE FITTING
5.1 Introduction
This lab explains two techniques for numerically solving equations, Newton’s
famous method and the bisection method. If we have any equation we want
to solve for x, we can subtract one side from the other to get an equation of
the form f (x) = 0. Of course, in case f (x) is a polynomial then solving this
equation means finding the roots of f (x).Thus, for quadratic polynomials we
would ordinarily use the quadratic formula. However, we will be considering
very general functions which typically involve trigonometric functions, loga-
rithms and exponentials and hence algebraic methods are usually hopeless.
Newton’s method is called a dynamic process and is related to interesting
topics such as chaos and fractals. We will explore these concepts later in this
chapter.
69
70 CHAPTER 5. FINDING ROOTS USING COMPUTERS
to r.
Since y − y0 = m(x − x0 ) is the equation of the line through (x0 , y0 ) with
slope m, the equation for the tangent line of f (x) through (x0 , f (x0 )) is
Solving for x when y = 0 gives x = x0 − f (x0 )/f 0(x0 ). Thus we get the
(n + 1)st approximation from the nth by the formula:
f (xn )
(1) xn+1 = xn −
f 0 (xn )
(x1 , 0) ≈ (2.36363, 0). The process is now repeated, starting with the guess
x1 .
It is convenient to view the computations as an iteration process:
f (x)
(2) NG(x) = x −
f 0 (x)
which changes a guess x into a (hopefully) better guess NG(x). (Note that
xn+1 = NG(xn ).) You can think of NG as standing for ‘Newton guess’ or
for ‘next guess’. To make a Derive function to do this for our function
f (x) = x2 + x − 1 we define
(3) NG(x) := x-(x^2+x-1)/(2x+1)
Now starting with x0 and successively applying this function to the previous
result produces a sequence of approximations:
x0
x1 = NG(x0 )
x2 = NG(x1 ) = NG(NG(x0 ))
x3 = NG(x2 ) = NG(NG(x1 )) = NG(NG(NG(x0 )))
..
.
which we hope get closer and closer to the exact answer. In the limit we
want this sequence of approximations to converge to the root.
We can compute several approximates by first Authoring NG(5), and then
approximating. Now we can author NG, press the right mouse button and then
click (Insert expression) or press F4. This will bring down the highlighted
expression in parentheses giving NG(2.36363) which we approximate (just
press Simplify instead of OK) again and then repeat this process.
A somewhat fancier method is to use the Derive’s ITERATES func-
tion. ITERATES(u,x,a,n) simplifies to an (n + 1)-vector whose first en-
try is a and each subsequent entry is obtained by substituting the pre-
vious entry for x in u. Thus, ITERATES(x^2,x,2, 4) returns the vector
[2, 4, 16, 256, 65536]. (The function ITERATE is similar, but just gives the
last value, so ITERATE(x^2,x,2,4) gives 65536.) We can get the first 4
approximates by authoring ITERATES(NG(x), x, 5, 4) and approximating
the result.
72 CHAPTER 5. FINDING ROOTS USING COMPUTERS
Loading the utility file ADD-HEAD adds two functions to the system that
make computing the Newton iterations easier. The function NEWT(u,x,a)
computes the Newton guess of the expression u, in the variable x, starting
with an initial guess at x0 = a. In our previous example of f (x) = x2 + x − 1
with starting point x0 = 5 we would enter NEWT(x^2+x-1,x,5). To get a
vector containing the starting point and the first 4 Newton iterates you author
and simplify NEWT(x^2+x-1,x,5,4). The general syntax is NEWT(u,x,a,k).
Looking at the algebra window in Figure 5.1 we see the above function
along with the first 4 iterates starting at x0 = 5. The graphic demonstration
shows the Newton method in action by plotting a part of the tangent line
until it crosses the x-axis. The picture clearly shows how well the Newton
method works since one has to zoom-in several times near the actual root in
order to see the last two iterations. The utility function DRAW NEWT(u,x,a,k)
simplifies to a matrix which plots the figure shown in Figure 5.1.
Alternately, that file contains the necessary definition for doing the graph-
ics directly. The basic idea is to make a vector out of several triples of points
which have the form (x, 0), the initial guess on the x–axis, (x, f (x)), the cor-
responding point on the curve, and (NG(x), 0), the place where the tangent
to the curve at (x, f (x)) intersects the x–axis. When we graph these points
we want the lines connecting them to be drawn. If this is not the case then
adjust the Options/Points menu.
You might note that a little trick is used in the above of DRAW N in
the file F-NEWT. The special form of the VECTOR(u,x,v) function sets x
equal to each value in the vector v = [v1 , . . . , vn ] and makes the new vector
[u(v1 ), . . . , u(vn )].
√
Example. Suppose we want a numerical approximation of 2. We think
of it as a solution to the equation x2 − 2 = 0. Then formula (2) gives the
very simple expression:
f (x) x2 − 2 x2 + 2 x 1
NG(x) = x − 0
= x − = = +
f (x) 2x 2x 2 x
which shows that the xn ’s are getting closer and closer to r. To estimate the
number of decimals in common take the special case that λ = 0.1 = 10−1.
Then, |xn − r| ≤ 10−n |x0 − r|, i.e., mn = n, and we see that we gain one
decimal place with each computation.
5.3. WHEN DO THESE METHODS WORK 75
In general, the same sort of thing will occur with any 0 < λ < 1. If λ is
greater than 0.1 then the number of decimals mn will increase at a slower rate
but nevertheless it can be shown to increase without bound. For example,
in Derive if we approximate VECTOR(1/2^n,n,1,16,3) then we get
which suggests that mn ≈ n/3. In fact, if your numerically solve the equation
2−n = 10−m for m using Derive you will get m = 0.301029n.
If we assume that f 00 (x) exists and that f 0 (x) 6= 0 on I we can compute
NG0 (x) as follows:
0 d f (x)
NG (x) = x− 0
dx f (x)
0 2
f (x) − f (x)f 00 (x) f (x)f 00 (x)
=1− = .
f 0 (x)2 f 0 (x)2
Since f (r) = 0 we see that NG(r) = 0! Thus, not only is | NG0 (x)| less than
one near r but it can be made as small as you like. This follows from the
continuity of the function NG(x) which in turn follows from the continuity
of f 0 (x) and the fact that it is never zero near r. The implication of this fact
is a much faster rate of convergence, i.e., the mn increases much faster than
the above computations suggest.
f (x)
NG(x) = x −
f 0 (x)
Before we prove this fact, let’s see what it really means. Suppose for
example that M = 50 and we pick x in I that is within 2 decimals places of
76 CHAPTER 5. FINDING ROOTS USING COMPUTERS
the exact value r, i.e., |x − r| < .01. The above relationship then says that
our Newton guess starting at x, NG(x), satisfies:
and since .005 is much smaller than .01 we have a big improvement over our
initial guess of x. Moreover, since our new guess is closer to r than our initial
guess we see that it is also in the interval I and hence we can repeat the same
calculation using NG(x) in place of x.
We could make further calculations, as above, to see how rapidly the
sequence of approximates converges to the exact root r or we could have
Derive do it for us. If we author the expression
which clearly shows the very rapid rate of convergence to the exact answer r.
Proof. By making the interval I a little smaller if necessary we may assume
that both f 0 , f 00 are continuous and bounded on I. Let us assume that
where c1 , c2 are between x and r and hence in the interval I. Using the
bounds for f 0 , f 00 and 1/f 0 above shows that (7) holds because
A
| NG(x) − r| ≤ B|c1 − r||x − r| ≤ M|x − r|2
C2
which is what we needed to prove.
For another illustration of how fast the Newton method converges when
NG0 (r) = 0 let us assume that the constant M in (7) is one or less. If we
make an initial guess x0 which is in I and approximates r to n decimal places,
i.e., |x0 − r| < 10−n then the first approximate x1 = NG(x0 ) satifies
Theorem 2. Suppose that f (r) = 0 and that f (x) = (x − r)m g(x) where
g(x) is differentiable, m is a positive integer and g(r) 6= 0. Then, NG(x) is
defined for all x 6= r which are sufficiently close to r and the iterates converge
to r.
Proof. Since f 0 (r) = 0 for m > 1 (check!) it is not clear that we can even
define NG(x) for x near r. But
and since g(r) 6= 0 it is easy to see that the bracketed expression above can
not be zero for all x near to r and hence the same is true of f 0 (x) provided
x 6= r.
78 CHAPTER 5. FINDING ROOTS USING COMPUTERS
Now using (8) to simplify NG(x) (do this using Derive) we get
When x = r, the right hand side above simplifies to (m − 1)/m which is less
the one. This means that we can find an interval I, with midpoint r, and a
number 0 < λ < 1 such that
| NG(x) − r| ≤ λ|x − r|
for all x in the interval I. As we argued earlier in this section, this means
that starting at x0 in I we get the sequence of approximates xn satisfying
|xn − r| ≤ λn |x0 − r| which shows that the xn ’s converge to r.
• Addition and Subtraction. You do the real and imaginary parts sepa-
rately:
(a + b i) + (c + d i) = (a + c) + (b + d) i
(a + b i) − (c + d i) = (a − c) + (b − d) i
∗
5.4. COMPLEX NUMBERS, FRACTALS AND CHAOS 79
The formulas may look a little complicated at first but nevertheless the
usual rules of algebra hold for complex numbers and we can use them just as
we would ordinary real numbers. Of course, it is easier to make computations
using Derive then it is to compute with the formulas above. You should
verify the above formulas using Derive.
You probably have already encountered complex numbers in Derive
when, for example, you try to Solve an equation such as x2 + 1 = 0 or
something more complicated, while trying to find extreme points. The result
is that Derive computes the two solutions x = ±i. Of course, in calculus we
usually ignore complex solutions since they are not relevant to max-min the-
ory or graphics. Nevertheless, they do play an important role in algebra since
they provide a complete theory for the solution to polynomial equations.
Complex numbers can be represented as points in the plane: (a, b) for
a + b i. When we are plotting complex numbers we usually refer to the
plane as the complex plane. Derive does not directly plot complex numbers
as points in the plane but we have a utility function DRAW COMPLEX which
converts a vector of complex numbers into a matrix of points in the plane,
using the RE and IM functions, which can then be plotted as usual.
and then simplify to get [i, −1, −i, 1, i]1 . Now we apply our utility function
to this vector and simplify to get the matrix
0 1
−1 0
0 −1
1 0
0 1
Notice that the row entries (x, y) corresponds to the complex numbers x + y i
from the above vector. Now plotting these five (x, y)-points as usual gives us
an interesting picture of how complex multiplication is quite different from
ordinary real numbers. Notice that the powers are rotating in the counter-
clockwise direction about the origin. Also notice that the last point is the
1
We could use vector(#i^n,n,1,5). to produce this vector
∗
5.4. COMPLEX NUMBERS, FRACTALS AND CHAOS 81
same as the first since i5 = (−1)2 i = i. Which real numbers x have the
property that xn = x for some integer n?
Using the same function f (x) = x2 − 2, let’s see what happens if we start
with a complex number for x0 like 3 + 2 i. The formula (1) 1 on page 70 for
Newton’s method involves arithmetic operations with the complex number
x0 = 3 + 2 i and hence is valid. This time we Author
and then we approximate this with precision set to 6 decimal places, to obtain
DRAW COMPLEX(NEWT(x^2-2,x,#i,25))
Notice that all the values are purely imaginary (they only have an i com-
ponent) and that they seem to bounce around randomly. Moreover, if you
author
DRAW COMPLEX(NEWT(x^2-2,x,1.01#i,25))
you’ll notice that the corresponding entries of the answers are approximately
the same for the first few terms but very quickly seem to have no relation to
82 CHAPTER 5. FINDING ROOTS USING COMPUTERS
each other. Here’s a nice way to do this: Define the first set of points to be
data1 and the second set to be data2. Author the vector [data1,data2]
and simplify. Then scroll though the matrix to compare entries.
In other words, even though the two starting points above; namely i
and 1.01i are quite close together their long-term behavior seem completely
different. The above phenomenon is what is known as chaos.
We can illustrate this last property graphically by looking at the √
equation
x √+ 2 = 0 rather than x − 2 = 0. The former equation has roots 2 i and
2 2
− 2 i. Just as before if we start with any point a+ib in the upper half of the
complex
√ plane (b > 0), the Newton iterates of the function x2 + 2 converge
√
to 2 i and any point in the lower half plane (b < 0) converges to − 2 i.
Test this by plotting
DRAW COMPLEX(NEWT(x^2+2,x,1+#i,5))
But we get chaos on the real axis. To see this chaos plot the function x2 + 2
and the output to DRAW NEWT(x^2+2,x,2,4) in connected mode, see Fig-
ure 5.4 on the next page.
∗
5.4. COMPLEX NUMBERS, FRACTALS AND CHAOS 83
http://www.math.hawaii.edu/lab/newton.html
This web page contains an interactive Java applet that shows the iterations in
Newton’s method. You click a complex starting point and the applet shows
the first five iterations. Try it!
84 CHAPTER 5. FINDING ROOTS USING COMPUTERS
Constructing the Julia set The set of points where Newton’s method
fails, that is, the set of points x0 where the sequence
fails to converge, is called the Julia set for NG(x). In the example f (x) =
x3 −1 these are the points on the edge or boundary of the basin of attraction.
As the picture in Figure 5.5 shows this set can be very complicated, it looks
a little like a necklace with infinitely many smaller and smaller loops coming
out in many different directions.
There are two basic methods for constructing this set. Since NG(x0 ) is
not even defined when f 0 (x0 ) = 0 this is a good place to start. If x is a
solution of
NG(x) = x0
then the third term of the sequence (9) is not defined and so x will be in
the Julia set. In the case when f (x) = x3 − 1 the equation above has
∗
5.4. COMPLEX NUMBERS, FRACTALS AND CHAOS 85
three solutions. For each of these there are three more obtained by solving a
similar equation or in other words finding the points where NG(NG(x)) = x0 .
Continuing in this way we get a close approximation to the Julia set. The
actual set is obtained by taking limits of these points. This method is called
the backward method and is done in the file F-JULIA-BACKWARD for the
polynomial x3 − 3x which has√critical points at ±1. This function has three
real roots at x = 0 and x = ± 3 and the Julia set somehow has to separate
the three basins of attraction corresponding to these roots. See Figure 5.6
for a picture of it’s Julia “necklace”.
Figure 5.6: Bad Newton starting points for x3 − 3x = 0 in the complex plane
The trouble with the backward method is that it uses the cubic formula
for solving 3rd order equations and this formula is pretty complicated. Even
worst is the fact that there is no analogous formula for degrees 5 or greater.
To get around this problem there is the “forward method” which involves
simply looking at the sequence:
and checking whether it gets closer and closer to root or else just wanders
around forever. Since you have to do this for each point or pixel in the graph
this can be a very lengthy computation. A number of shortcuts and tricks
are typically employed and you can study the file F-JULIA-FORWARD to
see how we did it. Or you can just check out the pictures; see Figure 5.7.
and apply the bisection method to it; otherwise we use [(a + b)/2, b]. At each
stage the root lies in an interval that is only half the size of the previous
stage. So we can eventually find the root to any number of decimal places.
We can automate this process by authoring two functions:
F(x) :=
BIS2(a,b) := IF(f(a)f((a+b)/2)<0, [a, (a+b)/2], [(a+b)/2, b])
BIS(v) := BIS2(v SUB 1, v SUB 2).
The main function is BIS(v) and BIS2(a,b) is a helper function. The ar-
gument v to BIS is a vector with two entries, e.g., [a, b]. The Derive
function SUB, which we discussed in the previous section, returns the parts
of a vector so that [a,b] SUB 1 = a and [a,b] SUB 2 = b. So BIS starts
with a vector like [a,b] and calls BIS2(a,b). This then uses the values f (a)
and f ((a + b)/2) to decide if there is a root in [a, (a + b)/2] or in [(a + b)/2, b].
In the discussion above we assumed that f (x) < 0 and that f (b) > 0. The
way we have defined BIS it will work also in the case f (x) > 0 and f (b) < 0.
To do this we test if the product f (a)f ((a + b)/2) is negative. If it is, then
one of f (a) and f ((a + b)/2) is negative and the other is positive. In this
case the points (a, f (a)) and ( a+b
2
, f ( a+b
2
)) lie on opposite sides of the x–axis
and so there must be a root in the interval [a, (a + b)/2]. In the other case,
f (a)f ((a + b)/2) is positive and so they have the same sign. In this case
f ((a + b)/2) and f (b) must have the opposite signs (why?) and so there is a
root in [(a + b)/2, b].
Let us try the equation ln x = 1 which has the (unique) solution x = e =
2.718 . . . . Of course we are finding the root of ln x − 1 so we author f(x)
:= ln(x) - 1 and apply BIS. Graphing f (x) shows that there is a root
between 2 and 3 so we author BIS([2,3]). This returns [2.5, 3], indicating
that 2.5 < e < 3.
Now we want to apply BIS to the answer [2.5, 3]. You can do this several
times by choosing author, typing BIS, and then inserting the highlighted
vector. Once again we have an iteration process and we can use the ITERATES
function that does this for you.
Using this technique, we author
ITERATES(BIS(v),v,[2,3],10)
and then approximate it to see how well this approximates e, see Figure 5.8
on the following page.
88 CHAPTER 5. FINDING ROOTS USING COMPUTERS
An easier way to see the bisection method in action is to use the function
BISECT(u,x,v,k) in the utility file ADD-UTIL. To get the above results
we would simply enter BISECT(ln x-1,x,[2,3],10) and the press the ap-
proximation button. It is interesting to compare the results of the
bisection method with the Newton method of the previous section. The bi-
section method is fairly fast at getting a good approximation but not nearly
as fast as the Newton method.
The bisection method will work for any f that is continuous on the interval
[a, b] and f (a) and f (b) have opposite signs. It is easy to see that after n
iterates the error is at most (b − a)/2n . (In fact this is the width of the
resulting interval. If we choose the midpoint as our estimate, the error will
be at most (b − a)/2n+1 .)
5.6. LABORATORY EXERCISES 89
x0 , x1 , x2 , x3 . . .
would look like? What are their signs? What do you think
limn→∞ |xn | is?
e. Choose any numbers a and b which satisfy
7. Use the definitions for adding and multiplying complex numbers given
in Section 5.4 on page 78 to answer these problems using pencil and
paper and then check your computations using Derive. Recall that
the definition of the complex conjugate is a + b i = a − b i and√that in
Derive you use the function CONJ. The modulus is |a + b i| = a2 + b2
and in Derive you use the function ABS, i.e. absolute value.
a. (1 + 2i) + (3 − i) b. (1 + 2i) + (3 − i)
c. (1 + 2i)(3 − i) d. (1 + 2i)(3 − i)
a. z + w = z̄ + w̄ b. z − w = z̄ − w̄
c. zw = wz d. zw = z̄ w̄
g. P (z) = P (z̄)
92 CHAPTER 5. FINDING ROOTS USING COMPUTERS
9. Use the discussion about complex numbers in Section 5.4 and the utility
function DRAW COMPLEX to answer these questions.
*10. Theorem 1 had the hypothesis that f (r) = 0 and f 0 (r) 6= 0. In this
problem we explore what happens to a function when f 0 (r) = 0. Let
f (x) = x(x2 − 2)2 .
11. The function f (x) = x1/3 has a root at x = 0. Find NG(x), NG0 (x), and
NG0 (0). Find 10 iterates of Newton’s method starting with x0 = 0.1.
(Note: Make sure that the Precision Mode is set to Exact or else there
may be problems with this exercise.)
holds for n = 0, 1, 2, . . . .
94 CHAPTER 5. FINDING ROOTS USING COMPUTERS
*14. The light area in Figure 5.5 on page 84 shows the basin of attraction
of the root 1 when using Newton’s method on x3 − 1. The origin of the
complex plane is in the middle of this figure. Note that most of the
negative real axis (the negative x–axis) is in the white area. This means
that starting with most negative real numbers, Newton’s method will
converge to 1. Try this for x0 = −1 and −2. If you look closely at
the figure you see that black pinches down on the negative real axis at
various points. Find the value of the first such point to the left of the
origin.
Chapter 6
Numerical Integration
Techniques
6.1 Introduction
This lab discusses numerical integration. Numerical integration is described
in most calculus books and is sometimes covered in second semester calculus.
You may want to look over this part of your calculus text.
A function is called elementary if it is made up of sums, products, powers,
and compositions of the trig functions and ln x and ex . Although the deriva-
tive of any elementary function is elementary, not all such functions have
elementary antiderivatives. For Rexample, there is no elementary function
whose derivative is sin(x2 ), i.e., sin(x2 ) dx is not an elementary function.
Consider the problem
Z 1
sin(x2 ) dx
−1
Even though sin(x2 ) has no elementary antiderivative, the area defined by the
integral certainly exists. So how do we find it? We use numerical integration.
Rb
Consider the integral a f (x) dx, and for simplicity assume f (x) ≥ 0 and
that a < b. The idea of numerical integration is to choose intermediate points
a = x0 < x1 < x2 < · · · < xn = b and estimate the area in the strip below
f (x) for xi ≤ x ≤ xi+1 and then add up these estimates; see Figure 6.2 on
page 100. Of course the width of this strip is xi+1 − xi . The height varies
with x. Some of the most common ways of estimating the area of the strip
95
96 CHAPTER 6. NUMERICAL INTEGRATION TECHNIQUES
are:
• Left endpoint: f (xi ) (xi+1 − xi )
• Right endpoint: f (xi+1 ) (xi+1 − xi )
xi+1 + xi
• Midpoint: f( ) (xi+1 − xi )
2
1
• Trapezoid: [f (xi+1 ) + f (xi )] (xi+1 − xi )
2
1 xi+1 + xi
• Simpson’s Rule: f (xi+1 ) + 4f ( ) + f (xi ) (xi+1 − xi )
6 2
The last one, Simpson’s Rule, is based on the best quadratic approximation
to f (x). This basic idea was derived in Exercise 4 on page 65 in Chapter 4.
Section 6.5 on page 104 has a detailed explanation.
Usually we choose the xi ’s equally spaced, so that
b−a
(1) xi = a + i
n
b−a
Of course, in this case, xi+1 − xi = . Thus, if we use the left endpoint
n
approximation, we get
Z b
b−aX
n−1
(2) f (x) dx ≈ f (xi )
a n i=0
b−a
Notice that we factor out the term and multiply by the sum rather
n
than multiplying every term.
6.2 An Example
Formula (2) suggests how we might do numerical integration with Derive.
Let u be the expression for f (x). We can define a Derive function for the
left endpoint method by
LEFT(u,x,n,a,b) :=
(b-a)/n * SUM(SUBST(u, x, a + k*(b-a)/n), k, 0, n-1)
6.2. AN EXAMPLE 97
SUBST(u, x, a+k(b-a)/n)
To use the left endpoint method with n = 10 intervals, we would just author
and then approximate
LEFT(1/x, x, 10, 1, 2)
Doing this gives the answer 0.718771. Similarly if we wanted to use the
trapezoid method we would author and approximate TRAP(1/x, x, 10, 1,
2) which gives 0.693771.
We want to compare the accuracy of these methods of approximation
and also see how much the accuracy is improved by increasing n. We will
try them for n = 10, 100, 1,000 and 10,000. A fancy way to see and compare
approximation values, using the left endpoint rule for a range of n is to start
by authoring the vector
[10^n,LN(2),LEFT(1/x,x,10^n,1,2),LEFT(1/x,x,10^n,1,2)-LN(2)].
method and the fourth column being the error. See Section 0.14 on page 18
for more discussion on the vector function.
Notice from Figure 6.1 that the accuracy in this method seems to be
roughly 1, 2, 3 and 4 digits respectively. This is in fact the case and it can
be proved that using 10n subdivisions yields an accuracy of n decimal places.
This is not very efficient since it requires a billion computations (109 ) to
achieve calculator accuracy of 9 digits. Try comparing computation times
for various powers of 10 to see how this rapidly becomes impractical. If we
try to obtain simple calculator accuracy of 8-12 decimal places, then this
can take hours on a PC which is impractical. It is for this reason that we
investigate the other methods for computational purposes.
By replacing the left endpoint method with the trapezoid method in the
computation in Figure 6.1 we see a remarkable difference. The accuracy now
appears to be approximately 2, 4, 6 and 8 digits respectively. Thus, the 4
decimal place accuracy achieved by the left endpoint method using 10,000
rectangles is equivalent to the trapezoid method using only 100 trapezoids.
We can summarize the theoretical error for these methods as follows. It
6.2. AN EXAMPLE 99
can be shown that error in using the left endpoint method is no greater than
(b − a)2 0 1
(3) max |f (x)| .
2 x∈[a,b] n
On the other hand, the error in using the trapezoid method is no greater
than
(b − a)3 00 1
(4) max |f (x)| 2 .
12 x∈[a,b] n
Notice the main difference between (4) and (5) is that we now have an
error which is roughly 1/n4 (the bracketed quantity in our example is 24/180).
Thus, with n = 10 we obtain the same accuracy as n = 100 in the trapezoid
method or n = 10, 000 in the left endpoint method. A table illustrating these
differences can be obtain by approximating
These functions are available by doing Load/Utility with the file ADD-UTIL.
Seeing the accuracy of SIMP(1/x, x, 104 , 1, 2) requires 16 digits of accuracy.
Recall from Section 0.6 how to increase the accuracy of a calculation.
To get a geometric feeling for why the trapezoid method is so much better
than the left endpoint method one need only draw a sketch comparing the two
methods. It’s possible to graphically represent these approximations using
Derive. Recall from Chapter 4 that one can plot a collection of points,
100 CHAPTER 6. NUMERICAL INTEGRATION TECHNIQUES
which is the desired result. This completes the proof of part (a).
The proof of part (b) is similar except it uses the Mean Value Theorem
three times. We estimate the error from approximating f (x) by the linear
function obtained from the endpoints values f (xk−1 ) and f (xk ). Thus, for
xk−1 ≤ x ≤ xk we have
f (x) − f (xk ) − f (xk−1 ) · (x − xk−1 ) + f (xk−1 )
xk − xk−1
f (xk ) − f (xk−1 )
= (f (x) − f (xk−1 )) − · (x − xk−1 )
xk − xk−1
= |f (c1 )(x − xk−1 ) − f (c2 )(x − xk−1 )| = |f 00 (c3 )||c1 − c2 ||x − xk−1 |
0 0
2
00 b−a
≤ max |f (x)|
x∈[a,b] n
102 CHAPTER 6. NUMERICAL INTEGRATION TECHNIQUES
and thus at each point the error is proportional to 1/n2 and so is the integral
over [a, b]. More precisely,
Z b
1
f (x) dx − TRAP(f (x), x, n, a, b) ≤ (b − a) max |f (x)| 2
3 00
x∈[a,b] n
a
We note that the error estimates above differ from (3) and (4) only in the
constant term and not the power of n. To obtain the better constant more
careful estimation needs to done in the above argument. On the other hand,
the constants obtained above suffice for most applications.
Use Derive to graph g(x) = 1/(1 + x3/2 ). Notice that the graph is pretty
tame; there are no wild oscillations and it would appear that the trapezoid
method could be used to obtain a good approximation of (6). In fact it does
give a good approximation.
In order to use (4) to estimate the error in using the trapezoid rule to
evaluate (6) we need to find g 00√ . Use Derive to do this. Note that g 00 (0)
is undefined; but that limx→0+ xg 00 (x) = − 34 . This means that g 00 (x) ≈
− 34 x−1/2 and hence is not bounded on [0, 1] so that (4) gives us no information
about the error.
We can work around this problem by noticing that for each n we can
apply (4) to the interval [ n1 , 1] instead and use a different technique for that
first interval. Thus, using |g 00 (1/n)| for the maximum on [1/n, 1] (check that
this is valid for all large n), we obtain from (4) that
Z 1
√ 1 c
g(x) dx − TRAP(g(x), x, n − 1, 1/n, 1) ≤ c n · ≈ .
(n − 1)2 n3/2
1/n
∗
6.4. MORE ON ERROR ESTIMATES 103
On the small interval we observe that g(x) is decreasing for x > 0 and that
g(0) − g(x) = x3/2 /(1 + x3/2 ) ≤ x3/2 . Thus, by comparing areas we see that
Z
1/n 1
1 1
g(x) dx − TRAP(g(x), x, 1, 0, 1/n) ≤ (g(0) − g( )) ≤ 5/2 .
0 n n n
Combining these estimates shows that the error obtained using the trapezoid
method is proportional to n−3/2 (which is the larger of the two errors). This is
a better result than 1/n but not as good as 1/n2 . Actually, one can improve
the 3/2-power a little by refining these estimates.
The next question is what can you do without explicit estimates like the
above but only using monotonicity or convexity of the graph. If f Ris increas-
b
ing on [a, b] notice that the left endpoint method of estimating a f (x) dx
always underestimates the integral while the right endpoint method overes-
timates it. Similarly, if f is decreasing the opposite inequalities hold. If we
let LEFT(f (x), x, n, a, b) and RIGHT(f (x), x, n, a, b) be the left and right
endpoint estimates then:
Z b
(7) LEFT(f (x), x, n, a, b) ≤ f (x) dx ≤ RIGHT(f (x), x, n, a, b)
a
if f 0 (x) ≥ 0 on [a, b]
and
Z b
(8) RIGHT(f (x), x, n, a, b) ≤ f (x) dx ≤ LEFT(f (x), x, n, a, b)
a
if f 0 (x) ≤ 0 on [a, b]
See Figure 8.2 on page 143 which makes these relations quite obvious.
A similar relation holds between the trapezoid and midpoint methods
but depends on the concavity, i.e., the second derivative of f rather than
the slope, i.e., the first derivative of f . If we let TRAP(f (x), x, n, a, b) and
MID(f (x), x, n, a, b) be the trapezoid and midpoint estimates then
Figure 6.3 shows why this is true. It has two graphs of the same function
which is concave up. the line in the left part shows the trapezoid used in the
trapezoid rule. Clearly it overestimates the integral. The midpoint rule is
illustrated in the right graph. The midpoint rule gives the area under the line
AB. The line CD is the tangent line through the midpoint. The area below
AB is the same as the area below CD (why?). So both are the midpoint
estimate. But clearly the area under CD is less than the area under the
curve.
A B
b. Plot the graph of f (x) = 4/(1 + x2 ) and verify that this function
is decreasing on the interval [0,1].
106 CHAPTER 6. NUMERICAL INTEGRATION TECHNIQUES
5. For the following integrals use the error estimate (4) described above
to find an n large enough so that the trapezoid method will give an
approximation of the integral with error at most 0.005. Give both the
approximate value of the integral and the smallest n which guarantees
(using formula (4)) that you will be within this error, and also give
M2 = max{|f 00 (x)| : a ≤ x ≤ b}.
Hints: Use Derive to find f 0 , f 00 , and f 000 . For the first integral below,
you can easily see that the maximum for |f 00 | occurs when x = 1. For
the second, solve f 000 (x) = 0; this tells you where the maximums of
|f 00 (x)| can occur, and, using this (and maybe some plotting), you can
find M2 . For the third integral don’t forget that M2 if the maximum of
the absolute value of f 00 (x) on [0, 2]. Once you have M2 , find n large
enough so that the error given in (4) is at most 0.005.
6.6. LABORATORY EXERCISES 107
Z e2 Z 2
1
a. ln x dx b. dx
1 1
2
1 + x2
Z 2
1
c. dx
0 1 + x2
6. Do the same as the last problem, but use Simpson’s Rule this time and
of course use formula (5) instead of (4).
1 2
SIMP(f (x), x, n, a, b) = TRAP(f (x), x, n, a, b) + MID(f (x), x, n, a, b)
3 3
108 CHAPTER 6. NUMERICAL INTEGRATION TECHNIQUES
SIMP(f(x),x,n,a,b)
(1/3) TRAP(f(x),x,n,a,b) + (2/3) MID(f(x),x,n,a,b)
Taylor Polynomials
One natural way to to do this is to require that f (0) = Pn (0), f 0 (0) = Pn0 (0),
f 00 (0) = Pn00 (0), etc., i.e., f (k) (0) = Pn (0) for k = 0, . . . , n. This gives n + 1
(k)
(1) P3 (x) = a0 + a1 x + a2 x2 + a3 x3
(2) P30 (x) = 1 · a1 + 2 · a2 x + 3 · a3 x2
(3) P300 (x) = 2 · 1 · a2 + 3 · 2 · a3 x
(4) P3000 (x) = 3 · 2 · 1 · a3
f (k) (0)
(5) ak = for 0≤k≤n
k!
109
110 CHAPTER 7. TAYLOR POLYNOMIALS
X
n
f (k) (0)
(6) Pn (x) = xk
k=0
k!
The coefficient of xk in Pn (x) is just f (k) (0)/k! (which is the same for all n
as long as n ≥ k). This quantity is called the k th –Taylor coefficient for f (x).
As our first application, notice that it follows from (5) that the graph of
y = P1 (x) is just the tangent line to the curve y = f (x) at the point (0, f (0)).
We studied this method of approximation extensively in Section 2.2. Thus,
since the tangent line yields the best degree-one approximation to the func-
tion, near the point x = 0, it is reasonable that guess that Pn (x) is the best
nth –degree approximation, near x = 0.
To have Derive compute a Taylor polynomial for a function first select
the Calculus/Taylor menu, then enter the function in the form, enter some
integer n, say n = 5, for the Degree and leave the Point1 value at its default
value of 0. This results in the expression TAYLOR(f(x),x,0,5). An alternate
approach after becoming familiar with its syntax is to simply author this
expression. See Figure 7.1 on the facing page for some of the basic examples
and a comparison of the graph of f (x) = 1/(1 − x) and it’s 5th degree Taylor
polynomial approximation. An interesting exercise is to load the file F-TAY0
which contains the expressions from Figure 7.1 and compare graphically the
various functions with their Taylor polynomials of different degrees.
7.2 Examples
As we see from Figure 7.1 on the next page, the formulas for the Taylor
polynomials for the function f (x) = 1/(1 − x) are pretty easy to guess after
looking at a few examples:
Pn (x) = 1 + x + x2 + x3 + · · · + xn .
1 Xn X∞
(7) = 1 + x + x2 + x3 + · · · = lim xk = xk .
1−x n→∞
k=0 k=0
The three dots in the above formula means that you are to keep adding more
and more terms forever. Since adding an infinite list of numbers together
looks like an impossible task, we define this in terms of limits. That is what
we mean by the last expression above. More generally, we denote an infinite
series by
X∞ Xn
ak = lim ak .
n→∞
k=0 k=0
that our conjecture in (7) above is in fact the most famous example of an
infinite series; namely, it is the geometric series.
We got the form of the Taylor polynomials above by the method of pattern
recognition, looking at lots of examples using Derive and then guessing at
the general result. To verify this result we must use (5) to compute the
Taylor coefficients. We will need to show f (k) (0) = k!. Using Derive we can
make a table of derivatives by authoring
VECTOR([k, DIF((1-x)^-1,x,k)], k, 0, 4)
and then simplifying. The answers seem to follow the pattern k!/(1 − x)k+1
which can be verified by having Derive check that:
1 + x + x2 /2! + · · · + xn /n!. Now, if we can also take the limit as in the case
of the geometric series, then
X
∞
xk x2 Xn
xk
(8) x
e = =1+x+ + · · · = lim .
k=0
k! 2! n→∞
k=0
k!
In fact, the series above does converge, for all values of x, to the exponential
function. Moreover, it is this series that forms the basis for numerical calcu-
lations of the exponential function on computers and calculators. Section 7.5
will give a complete explanation of this matter.
7.2. EXAMPLES 113
x3 x5 X ∞
x2n+1
(9) sin x = x − + + ··· = (−1)n
3! 5! n=0
(2n + 1)!
x2 x4 X ∞
x2n
(10) cos x = 1 − + +··· = (−1)n
2! 4! n=0
(2n)!
|x|n+1 |x|n+1
(11) |f (x) − Pn (x)| ≤ max |f (n+1) (t)| · =M·
0≤t≤x (n + 1)! (n + 1)!
Just put u = g(t) and v = (x − t)m+1 /(m + 1)! and apply the integration by
parts formula. Notice that the integrated terms, i.e., the uv|x0 terms, vanish
because g(0) = 0 at the left endpoint and (x − t)m+1 is zero when t = x.
116 CHAPTER 7. TAYLOR POLYNOMIALS
Proof. Put g(t) = f (t) − Pn (t) let M be the maximum of |f (n+1) | on the
interval [0, x]. By the definition of the Taylor polynomial, observe that
where the second fact follows since the (n + 1)st –derivative of any degree n
polynomial is zero (look back at (1) on page 109). Now we get to apply (12)
to g 0 , g 00 , . . . , g (n) with the result that
Z x Z x Z x
0 00 000 (x − t)2
g (t) dt = g (t)(x − t) dt = g (t) dt
2
0 0
Z x 0
(x − t)n
= ··· = g (n+1) (t) dt
0 n!
and hence that
Z x
f (x) − Pn (x) = g(x) − g(0) = g 0 (t) dt
Z x 0
Z x
(n+1) (x − t)n (x − t)n
= g (t) dt = f (n+1) (t) dt .
0 n! 0 n!
We take absolute values of the left and right-hand sides of the above to get
Z x
(x − t)n
|f (x) − Pn (x)| ≤ |f (n+1) (t)| dt
n!
Z x
0
(x − t)n xn+1
≤ M dt = M
0 n! (n + 1)!
the identity sin(π − x) = sin x to reduce the problem to the smaller interval
[0, π/2]. It’s even possible to reduce the interval to [0, π/4].
We can use formula (11) to estimate the error in using the Taylor poly-
nomial to estimate sin x. The computation of M = max |f (n+1) | might look
a little formidable at first but we observe that any derivative is equal to ei-
ther ± sin t or ± cos t and in either case M ≤ 1. Thus, we can take M = 1
and achieve 6 decimals of accuracy by determining the smallest integer n
satisfying
|x|n+1
(13) ≤ 10−6
(n + 1)!
For approximations on the interval [0, π/2] we could just take the worst case
by setting x = π/2 in the above.
We now have reduced the problem to solving (13) for the smallest possible
integer n. Unfortunately, the factorial expression means that we can’t use
simple algebra to solve this inequality. A simple numerical approach would
be to make a table with n in the first column and the above expression in the
second column. Examining the data will result in an answer provided n is
reasonably small. We did this earlier on page 138 when we studied the ratio
test. If this fails, as with the 1/k 2 series, you might try testing various powers
of 10. Both of these techniques are easy to do using the vector function. (In
the next section we present another way of finding n.) In Figure 7.4, see the
file F-TAY3, we analyze sin 100 by reducing the computation to a smaller
value of x (x = 0.530973), determining which n yields an error of less than
10−6 (n = 7) and then computing using P7 (0.530973). Observe that for the
sine function P2n+1 (x) = P2n+2 (x) and so for the error computation (13) we
use the higher power 2n + 3 instead 2n + 2 and hence
X
n
x2k+1 |x|2n+3
(14) | sin x − (−1)k |≤ for n = 0, 1, . . .
k=0
(2k + 1)! (2n + 3)!
Lastly, let us observe that the right hand side of (14) tends to zero for
any x. After all, for x fixed, the ratio of terms above is
|x|2n+5 (2n + 3)! |x|2 1
· = ≤
(2n + 5)! |x| 2n+3 (2n + 5)(2n + 4) 2
for all large n. Thus, |x|2n+3 /(2n + 3)! ≤ cx /2n (or for that matter |x|n /n! ≤
cx /2n ) for some constant cx and the sequence tends to zero because 1/2n → 0.
118 CHAPTER 7. TAYLOR POLYNOMIALS
By applying Theorem 1 we see that the Taylor series converges for all x and
we indeed have the representation
x3 x5 X∞
x2n+1
(15) sin x = x − + +··· = (−1)n
3! 5! n=0
(2n + 1)!
which is valid for all −∞ < x < ∞. In a similar manner we establish (10)
on page 114.
Thus,
x 2 n x n+1
(16) e − (1 + x + x + · · · + x ) < 3 x
2! n! (n + 1)!
and so we need only find n so that the right-hand side is sufficiently small.
We would like to define a function in Derive to determine the number of
terms n necessary to achieve 6 decimals places, rather than looking at tables.
First of all, recall from the previous section that
3x xn+1
lim =0
n→∞ (n + 1)!
for all values of x. Hence we are guaranteed that there is a first n for which
the above quantity is less than 10−6 . Moreover, this proves that the Taylor
series converges for all x, by Theorem 1, to ex . Thus, as stated earlier
X
∞
xk x2 Xn
xk
x
e = =1+x+ + · · · = lim
k=0
k! 2! n→∞
k=0
k!
X
∞
f (k) (c)
f (x) = (x − c)k
k=0
k!
x2 x3 X (−1)k+1
∞
(17) ln(1 + x) = x − + ··· = xk , −1 < x ≤ 1
2 3 k=1
k
although the proof that the series converges to ln x on the interval 0 < x ≤
1/2 is straightforward, see Exercise 2, is quite a bit harder than our earlier
examples to get the full interval 0 < x ≤ 1. You can load the file F-LOG
and try approximating ln x with higher degree Taylor polynomials to see if
you can confirm the above representation. The full convergence problem for
the logarithm function will be studied in Chapter 9.
122 CHAPTER 7. TAYLOR POLYNOMIALS
be an interval which is centered about the origin. Actually, there are four
possibilities for the interval of convergence: (−r, +r), (−r, +r], [−r, +r) or
[−r, +r] for some 0 ≤ r ≤ +∞.2 This number r is called the radius of
convergence.
Now consider the function 1/(1 + x2 ). We can obtain the Taylor series
for this function by substituting −x2 for x in (7):
1 X ∞
(18) = 1 − x2
+ x4
− x6
+ · · · = (−1)k x2k
1 + x2 k=0
As before we can use Derive to plot several of the Taylor polynomials for
this function; see Figure 7.7.
(i) Define f (x) to be the given function and plot its graph in the color
red.
(ii) Simplify the first vector function above to make a 7–vector which
has the degree n Taylor polynomial, expanded about x = 0, for
n = 4, . . . , 10 as its entries.
(iii) Graph this vector to plot each of these Taylor polynomial in suc-
cession. It is probably best to plot the first 2-3 polynomials indi-
vidually first and then the entire vector so that you can see how
7.8. LABORATORY EXERCISES 125
their graphs are getting closer and closer to the graph of f (x).
You might need to scale your pictures appropriately.
(iv) Simplify the second vector function to make a 7 × 2–table that
has the degrees n in the first column and Pn (x) in the second
column.
(v) Use your table to guess what the infinite Taylor series expansion
is.
*(vi) Prove that in each case, the Taylor series expansion converges to
the function and determine the interval of x’s for which it is valid.
In parts (b), (c) use the geometric series (2) on page 132 and the
fact that the series converges only for −1 < x < 1.
x8 x2
a. f (x) = +···+ +x
8 2
1
b. f (x) =
3−x
1
c. f (x) = 2 (Hint: For the pattern recognition you will need
x +2
to change the output mode to Rational. Use the Declare/Algebra
state menu to access the Output menu.)
2. Let f (x) = − ln(1 − x). We want to determine the Taylor series for
f (x) and prove that it converges to f (x) using Theorem 1 on page 115.
X
∞
1 xk x2 x3
(19) ln = =x+ + + ...
1−x k 2 3
k=1
c. Use Theorem 1 to show that the (1) holds for any −1 ≤ x ≤ 1/2.
(Hint: Carefully compute the right-hand side of (11) on page 115.
Then, show that the error estimate tends to zero as n → ∞ only
for −1 ≤ x ≤ 1/2.)
126 CHAPTER 7. TAYLOR POLYNOMIALS
1
a. b. tan−1 x
(x + 1)2
2
c. e−x
4. Load the file F-TAY3 and following the methods of Section 7.4 compute
sin 7. Which degree Taylor polynomial should you use to get an error
of less than 10−6 ?
6. The functions f (x) = sin x has only odd powers in its Taylor series ex-
pansion. This property can be explained by the fact that f (x) satisfies
the equation f (−x) = −f (x) as do all odd powers of x. It is because of
this that we call any such f (x) an odd function. Similarly, a function is
an even function if f (−x) = f (x) holds for all x, as do all even powers
of x.
*7. Recall from our discussion of complex numbers in Section 5.4 on page 78
that a complex number is of the form α = a + b i where a, b are real
numbers and i satisfies i2 = −1. Our goal in this exercise is define the
complex exponential eαx in such a way that the basic formula:
d αx
(20) e = αeαx
dx
is still valid. We will use Taylor polynomials as a guidepost to the
correct definition.
a. Find the 6th degree Taylor polynomial P6 (x) for the function eαx
(enter the α symbol using the symbol bar). Be sure to expand out
your answer into powers of x.
128 CHAPTER 7. TAYLOR POLYNOMIALS
In fact, this is the Euler Formula which is one of the most fun-
damental equations with complex numbers. Verify that using this
definition, (20) with α = i is valid.
f. Lastly, use the definition
(22) eαx = e(a+b i)x = eax eibx = eax (cos bx + i sin bx)
*8. In this problem we use the Taylor polynomials for the arc tangent
function tan−1 x to estimate π. Recall that tan−1 x is entered as atan x.
Series
8.1 Introduction
An infinite series is a sum with infinitely many terms:
X
∞
ai = a0 + a1 + a2 + · · ·
i=0
P∞
We define i=0 ai = s to mean that
X
n
lim ai = s,
n→∞
i=0
if this limit exists. If the limit does exist we say the series converges; oth-
erwise we say it diverges. There are two basic techniques for showing that
a series is convergent. One method is to show directly that the above limit
exists. There are not many examples when we can do this but a particularly
important one is geometric series which will be discussed in the next section.
The second method for showing convergence applies to series with non-
negative terms, i.e., the case that ai ≥ 0 for all i = 1, 2, . . . . In this case the
partial sums,
X
n
sn = a0 + a1 + · · · + an = ai , n = 0, 1, . . .
i=0
131
132 CHAPTER 8. SERIES
sn − xsn = (1 + x + x2 + · · · + xn )
− (x + x2 + · · · + xn + xn+1 )
= 1 − xn+1
X
n
1 − xn+1
(1) sn = xi = 1 + x + x2 + · · · + xn = , if x 6= 1
i=0
1−x
It’s instructive to verify this formula in DfW. You start by clicking the sum
button and enter x^k. Make sure the variable is k (not x) and set the
Start value to 0 and the End value to n. Click OK and edit the resulting
expression SUM( x^k, k, 0, n) by multiplying it by the factor (1 − x).
Lastly, use Simplify/Expand to get the desired 1 − xn+1 .
If |x| < 1, then limn→∞ xn+1 = 0. Thus, we get that limn→∞ sn exists
and so the series is convergent. In addition, the following simple formula for
evaluating geometric series holds:
X
∞
a
(2) axi = a + ax + ax2 + · · · = , if |x| < 1
i=0
1−x
If |x| ≥ 1 then the series diverges because limn→∞ sn does not exist.
We can verify this formula in Derive by entering, as we did above or
directly, the expression SUM( ax^k, k, 0, inf) which displays as the left
hand side of (2). Now we must declare that −1 < x < 1. We do this using
8.3. APPLICATIONS 133
the Declare/Variable Domain menu for the variable x and the open interval
(−1, 1). If you do this right the expression x :∈ Real (-1,1) will be the
result. If not you should be able to edit the expression until it is right.
Finally simplify the infinite series to get a/(1 − x) which is the desired result.
8.3 Applications
Geometric series are useful in several areas, for example, business and finance.
We will give some of the important examples.
s bal(k,p,r):=IF(k=0,p,(1+r)*s bal(k-1,p,r)).
Now suppose that the bank compounds your money quarterly instead of
annually. This means that they give you 6/4 = 1.5% interest four times a
year. So the amount of money in your account after n years is p(1.015)4n .
For a general rate r compounded k times a year, the amount of money after n
years is
r kn
(3) p 1+
k
This can be also be expressed as s bal(kn,p,r/k).
Now suppose that you deposit a dollars each year into a bank account
paying a rate r in interest, compounded annually. Suppose that you opened
the account with an amount p dollars. How much money will the account
have after n years? This is easy to do using a small modification in the s bal
function as follows:
s bal(k,p,r,a):=IF(k=0,p,(1+r)*s bal(k-1,p,r,a)+a).
In other words we need only account for the extra a dollars which are de-
posited each year. We can get a nice table of values by numerically approxi-
mating
VECTOR([k,s bal(k,1000,.06,100)], k, 0, 10)
to see how an initial balance of $1000 will grow over a ten year period, at
6% annual interest, if we add an extra $100 each year.
Now if you make the same table as above using the symbolic values for
p, r and a you get a sequence of expressions which don’t appear to follow
any clear pattern. On the other hand, if we substitute r1 = 1 + r everything
is much clearer. In DfW you would declare a variable r1 and use r1-1 as a
replacement for r, then the table obtained by entering
Here we used (1) with x = r1 = 1 + r. Thus, the geometric series arises nat-
urally in compound interest problems and provides us with a useful formula.
l bal(k,p,r,a):=IF(k=0,p,(1+r)*l bal(k-1,p,r,a)-a).
X
k−1
1 − (1 + r)k
(1 + r) p − a
k
(1 + r)i = (1 + r)k p − a
i=0
1 − (1 + r)
(1 + r)k − 1
(4) = (1 + r)k p − a
r
Using this formula we can easily get the general formula for a by solving
l bal(k,p,r,a) = 0 for a. Thus, the monthly payments a on a loan of p
dollars at a monthly interest rate r (divide the annual rate by 12) for a period
of n years (so k = 12n payments) is:
r(1 + r)k p
(5) a= where k = 12n.
(1 + r)k − 1
bal(48,20000,0.01,526.68) = 0.
d1 d2 d3 X dk ∞
(6) x = 0.d1 d2 · · · = + 2 + 3 +··· =
10 10 10 10k
k=1
Now clearly, the partial sums form an increasing sequence since the terms
are nonnegative numbers. However, maybe it is not completely obvious that
8.4. APPROXIMATING INFINITE SERIES 137
d1 d2 dn 9 9 9
+ 2 +···+ n ≤ + 2 +···+ n
10 10 10 10 10 10
9 9 1 9 1 9 1
= + + ···+
10 10 10 10 102 10 10n−1
9
9 9 1 9 1
≤ + + + · · · = 10
1 = 1.
10 10 10 10 102 1 − 10
Note
P∞ that thek key step above was recognizing that the geometric series
k=0 a(1/10) , where a = 9/10, sums to 1 by (2) on page 132.
Of course, we also notice that repeating decimals like 0.999 · · · = 1 and
0.333 · · · = 1/3 are all geometric series when represented as above. Try to
figure out the a, x in (2) in each case. This turns out to be true of any
repeating decimal and hence by the formula (2) these decimal numbers must
be fractions a/b where a, b are integers. In fact, the converse is also true,
namely, a decimal is a fraction if and only if it is eventually repeating.
1 X −3
∞
5 1 1 k
x= + 3 + 5 +··· = + 10 10−2
10 10 10 2
k=0
1 10−3 1 1 248
= + = + = .
2 1 − 10−2 2 990 495
We might notice that it is not possible to enter x in Derive as a decimal
but we can define it by means of the infinite series above. Then, simplifying
we get the above result.
(b) Suppose that 0 < x < 1 and that for some n, ai+1 /ai ≤ x for all i > n.
Then
X
∞ X
n
an+1
(7) 0≤ ai − ai ≤
i=1 i=1
1−x
Proof. If limi→∞ ai+1 /ai = λ < 1, then for any x satisfying λ < x < 1,
we know that the ratios will be less than x for all large i. Thus, given x
there is an n for which ai+1 ≤ xai whenever i > n. So an+2 ≤ xan+1 and
an+3 ≤ xan+2 ≤ x2 an+1 . In general, an+1+k ≤ xk an+1 and hence for all m > n
X
m X
n
0≤ ai − ai = an+1 + an+2 + · · · + am
i=1 i=1
an+1
≤ an+1 (1 + x + x2 + · · · ) = .
1−x
Thus, the partial sums {sm } are bounded and the series is convergent. More-
over, the inequality (7) follows by taking limits as m → ∞.
X
∞
2k 1 2 4
= + + + ...
k=0
k! 1 1 2
converges and estimate its value with an error of at most 10−6 . The first
step is to show that limk→∞ ak+1 /ai < 1. We think of the terms as a
8.4. APPROXIMATING INFINITE SERIES 139
for all k > 2. Thus, limk→∞ ak+1 /ak = limk→∞ 2/(k + 1) = 0 < 1 so the series
converges and furthermore we can take x = 1/2 and n = 2 in Theorem 1(b).
Now we must determine n so that
an+1
= 2an+1 < 10−6 .
1−x
We do this by authoring
VECTOR([n,2*term(n+1)],n,2,20) ,
approximating the result and then searching the entries (by scrolling) until
we find one smaller than 10−6 . It P turns out n = 13 works. The last step is
to compute the partial sum s13 = 13 k
k=0 2 /k! giving 7.38905.
We might observe how fortunate we were that k turned out to be so small.
Recall some of our computations using the trapezoid method or Simpson’s
rule where similar accuracy required thousands or even millions of computa-
tions (using the left endpoint method, for example). It is one of the funda-
mental properties of geometric series that they converge very rapidly. Think
about it, 6–decimal place accuracy with just 15 computations!
As it turns out this series is rather special since
X
∞
2k
= e2 = 7.38905 · · ·
k=0
k!
∗
Example . Now consider the harder problem of approximating
X∞
k!
k=1
kk
140 CHAPTER 8. SERIES
with error again of at most 10−6. Proceeding as before we author the formula
term(k):=k!/k^k. Now by first simplifying and then taking limits we see
that
ak+1 kk 1 1
lim = lim k
= <
k→∞ ak k→∞ (k + 1) e 2
since e > 2. Thus, the limit is less than 1 so the series converges. Further-
more, we can take x = 1/2 in Theorem 1(b). But now we need to find an
integer n so that ak+1 /ak = k k /(1 + k)k < 1/2 for all k > n.
This step is harder than before. If we graph f (x) = (x/(x+1))x it appears
to be decreasing for all x ≥ 0. See Figure 8.1.
after some rearrangement (see the file F-RATIO for the step-by-step pro-
cedure). Since xx /(x + 1)x is positive for x > 0, we need to show that
ln(x + 1) − ln x − 1/(x + 1) is positive. At this point Derive can’t help so
we need an idea from calculus. One quick way to solve this problem is to
use the Mean Value Theorem for the function g(x) = ln x. The Mean Value
Theorem says: g(b) − g(a) = g 0(c)(b − a) for some a < c < b. For a = x and
b = x + 1 this gives ln(x + 1) − ln x = 1/c for some x < c < x + 1. Thus
1 1 1
ln(x + 1) − ln x − = − >0
x+1 c x+1
where we obtain the desired inequality since c < x + 1. Thus, f 0 (x) < 0 for
all x > 0 and so f is decreasing.
Now that we have established that the ratios decrease we need to know
when they are less than 1/2. Since
ak+1 kk 1
= k
= f (k) ≤ f (1) = .
ak (k + 1) 2
for all k ≥ 1, it follows that we may apply the theorem for any n. Finally,
by (7), we must determine n so that
an+1 −6
1 = 2an+1 < 10 .
1− 2
As before we author
VECTOR([n,2*term(n+1)],n,2,20) ,
approximate the result and then search the entries until we find one smaller
than 10−6 . It turns out in this case that k = 16. Computing the partial sum
s15 gives 1.87985. P
Now suppose that your series ak satisfies lim ak+1 /ak = λ but that
λ ≥ 1. The case λ > 1 is pretty much like the case λ < 1 except that now
the series diverges. The idea is to pick 1 < x < λ and observe that
∞ = an + an x + an x2 · · · ≤ an + an+1 + an+2 + . . .
for some large n since now ak+1 /ak ≥ x for all k ≥ n. The case λ = 1 is
much harder since, as we see in the next section, there are examples in which
the series converges and examples where it diverges.
142 CHAPTER 8. SERIES
The Integral Test Suppose that f (x) is a decreasing positive P valued func-
∞
tion, for x ≥ 1. Let an = f (n). We want to approximate i=1 ai and
determine whether the series is convergent or divergent.
In Section 6.4 we saw that, for a decreasing function like f (x), the left
endpoint method of estimating a definite integral of f (x) always overesti-
mates the integral while the right endpoint method underestimates it. This
is quite obvious by looking at Figure 8.2 on the next page where we use
f (x) = 1/x as our function and apply the box drawing function from Chap-
ter 6; see the file F-SERINT for a demonstration. Now, we observe that since
the interval size is one, the area of the box with height f (n) is just an . From
this wePget that adding the area of boxes corresponds to partials sums of the
series ak . Thus, for any 1 ≤ n ≤ m
X
m+1 Z m X
m
(8) ai ≤ f (x) dx ≤ ai
i=n+1 n i=n
The sum on the left is the right endpoint estimate and the sum on the right
is the left endpoint estimate, when we use ∆x = 1 as the subinterval size.
From this inequality, we obtain the following theorem:
X
n X
n Z ∞ X
∞ X
n Z ∞
(9) ai ≤ ai + f (x) dx ≤ ai ≤ ai + f (x) dx
i=1 i=1 n+1 i=1 i=1 n
(c) The value of the series can be estimated using the following:
Z !
X
∞ X
n ∞
(10) 0≤ ai − ai + f (x) dx ≤ an
i=1 i=1 n+1
8.4. APPROXIMATING INFINITE SERIES 143
and hence the partial sums {sm } are bounded (the first n terms are irrelevant).
Thus, the series converges and the second inequality in (9) follows from
X∞ Xn X∞ Xn Z ∞
ai = ai + ai ≤ ai + f (x) dx .
i=1 i=1 i=n+1 i=1 n
A similar argument shows that the integral is convergent if the series is and
that the first inequality in (9) holds.
The first inequality in (10) is an immediate consequence of R(9) and simi-
n+1
larly it follows that the middle expression in (10) is bounded by n f (x) dx.
But since f (x) is decreasing this integral is less than or equal to f (n) = an
and the theorem is proved.
144 CHAPTER 8. SERIES
and the second one is the more refined estimate (10) which uses the quantity
X
n Z ∞
ai + f (x) dx
i=1 n+1
is convergent. First note that λ = 1 in the ratio test so we cannot use that
approach. Next, we take f (x) = 1/x2 and observe (say using Derive) that
Z n n
dx 1 1
2
= − = 1 − → 1 as n → ∞
1 x x 1 n
R∞
and hence the Pimproper integral 1 dx/x2 is convergent. Now, by the integral
test the series 1/i2 <P∞, that is the series is convergent. Actually, a similar
argument shows that 1/i < ∞ whenever p > 1.
p
You have to wonder how the π–term can possibly be involved in this compu-
tation. The proof of this fact is beyond the scope of this text but Derive
can help us believe this result. One way to do this is to have Derive simplify
the series and get π 2 /6 as the answer. It works, try it. A more independent
8.4. APPROXIMATING INFINITE SERIES 145
approach would be to compute the partial sums sn for several n and compare
with a decimal approximation to π 2 /6. In Figure 8.3 on the following page
we have Derive make these comparisons with n = 100, 1000, and 10,000.
We also observe that Derive knows about (11) and simplifies the series
accordingly.
Thus, to solve our problem we need only find n so that right-hand side of
(12) is less than 10−m , and then use
Xn
1 1 1 1 1 1
2
+ = 1+ + +···+ 2 +
i=1
i n+1 4 9 n n+1
for our estimate. For example, with m = 6 we need to take n2 > 10m or
n > 1000.
Something rather amazing occurred in this problem. In Figure 8.3 on the
next page we used the partial sums {sn } to approximate the sum s. This is
the natural thing to do since s = lim sn . But the accuracy in Figure 8.3 is
only 3 or 4 decimal places with n = 1000. This error is to expected since (8)
yields
X∞ Z ∞
1 dx 1
0 ≤ s − sn = ≤ =
i=n+1
i2 n x2 n
and the right-hand side is just less than 10−3 . On the other hand, adding
the term 1/(n + 1) (which is the integral in (10)) increases the accuracy to
1/n2 = 10−6. This accuracy is 1000 times better than the other estimate.
Put another way, suppose for example that both computations take about 3
seconds with n = 1000 on your PC, the amount of computation time needed
to produce 6 decimal place accuracy using the less efficient method is almost
an hour! See the file F-2-SER which contains a comparison of these methods.
This problem illustrates the potential value of a innovative approach to a
computation compared to the conventional solution.
146 CHAPTER 8. SERIES
P
Figure 8.3: Summing the series 1/i2
The Harmonic Series Let us apply the integral test to the harmonic
series, namely,
X ∞
1 1 1
= 1 + + + ...
i=1
i 2 3
so that any n > e100 ≈ 2.6 × 1043 is certainly large enough. On the other
hand, using (8) again we have
X
n X
n Z n−1
1 1 dt
100 ≤ =1+ ≤1+ < 1 + ln n
1
i 2
i 1 t
so that when combined with the above we see that the best n satisfies: e99 <
n < e100 .
2. Suppose you get a 30 year mortgage loan for $200, 000 which is to be
repaid in 30 · 12 equal monthly payments, based on an annual interest
rate of 7.5%.
b. How much do you still owe after your first payment? How much
of your first year’s payment went to interest and how much went
to paying off the principal?
c. Formula (4) gives the amount you still owe after k months. Re-
place the k with 12k in Formula (4) so that k now represents years,
approximate the resulting expression and then plot. (You’ll need
to adjust the range in such a way that the visible x-axis contains
the range 0 to 30 and the y–axis contains the range 0 to 200, 000.)
Notice that at the beginning the amount you owe changes slowly
but that near the end of the 30 years it changes quickly.
4. The bank says that it will give you a car loan of $6,000 provided you
make monthly payments of $135 for 5 years. What interest rate is the
bank charging? (Hint: You may need to be a little careful how you
compute this.)
X
∞
1
c.
k=1
1 + k2
9. Some of the following series converge and some diverge. Decide which
do which and state the required Theorem needed to prove your conclu-
sion.
X
∞ X
∞
1
a. 2k b.
k=0 k=2
k ln k
X∞
1 X
∞
2 · 4 · · · 2k
c. d.
k=0
ek k=1
(2k)!
P∞
10. Consider the series k=0 1/k!.
is known to hold for all t > 0. The formula is derived from an im-
portant technique in the theory of Fourier transforms called Poisson
summation. We will not attempt to prove this formula but instead try
to use it as a method of approximating π more efficiently than in an
earlier problem. It has a number of other useful applications too. We
will fix the value of t = 2 for the rest of this problem.
a. Using the ratio test, show that both infinite series in (13) are
convergent.
P
that ∞
2 −2k 2 π 2
b. Use Theorem 1 with x = e−6π and n = 1 to show √ k=1 e
is less than 10−8 . Thus, with an accuracy of 2 2π10−8 or roughly
7 decimal
√ places we can take the right hand side of equation (13)
to be 2π.
c. Using Theorem 1 again, show that
! !
X∞
2
X
6
2 /2
0< 1+2 e−k /2 − 1+2 e−k
k=1 k=1
X
∞
2 /2
=2 e−k < 10−10
k=7
and hence !2
1 X
6
−k 2 /2
π≈ 1+2 e .
2 k=1
9.1 Introduction
Z b Z b Z bX
n X
n
bk+1
f (x) dx ≈ Pn (x) dx = k
ak x dx = ak
0 0 0 k+1
k=0 k=0
153
154 CHAPTER 9. APPROXIMATING INTEGRALS
This technique
Rb works for integrals going from 0 to b. If you want to
approximate a f (x) dx, you can make the substitution u = x − a so the
R b−a
integral becomes 0 f (u) du.
Now, the idea is to approximate the integrand by its Taylor series but in
this case we recognize the connection with the geometric series; namely,
X
∞
1
tk = for −1 < x < 1.
k=0
1−t
We’ll actually use the following more refined estimate from Section 8.2:
1 X
n X∞
tn+1
(5) − tk = tk = tn+1 + tn+2 + · · · =
1 − t k=0 k=n+1
1−t
Taking absolute values of the above we need to evaluate the integral in (7).
Since this looks complicated, we instead try to obtain an upper bound. For
1
positive x, we uses the inequality 0 < 1−t ≤ 1−x
1
to obtain that
Z x k+1 Z x n+1
t t xn+2
(8)
dt ≤ dt = .
0 1−t 0 1−x (1 − x)(n + 2)
On the other hand, for negative x, we instead use 0 < 1−t 1
≤ 1 to get a similar
bound:
Z 0 k+1 Z 0
t |x|n+2
(9) dt ≤ |t| n+1
dt = .
x 1−t x (n + 2)
Hence, we have the desired approximation result because
xn+2
lim = 0 whenever 0≤x<1
n→∞ (1 − x)(n + 2)
and
|x|n+2
lim = 0 whenever −1 ≤ x ≤ 0
n→∞ (n + 2)
156 CHAPTER 9. APPROXIMATING INTEGRALS
One small point, the polynomial approximations in (6) look a little differ-
ent from the standard Taylor polynomials because the powers are expressed
with the index k + 1. This is just an artificial difference since
X
n+1 j
x x2 xn+1 X xk+1n
Pn+1 (x) = =x+ +···+ =
j=1
j 2 n + 1 k=0 k + 1
1 X xj ∞
(11) ln = for − 1 ≤ x < 1.
1−x j=1
j
x3 x5 X (−1)k
∞
sin x = x − + −··· = x2k+1
3! 5! k=0
(2k + 1)!
for all −∞ < x < +∞. Now for x 6= 0 we can divide both sides by x to get
sin x x2 x4
=1− + − ... .
x 3! 5!
9.4. AN INTEGRAL APPROXIMATION 157
Finally, since 1/(9 · 9!) ≈ 3.0619 × 10−7 we get 6 decimal place accuracy by
approximating the integral using 4 terms from the series.
In Figure 9.1 on the next page we have Derive approximate our integral
using 20 digit precision. This computation, which uses Simpson’s rule, is
actually quite slow, Load the file F-SININT and try this yourself. On the
other hand, we enter the partial sums of the series solution and make a table
comparing the first several sums with the answer from Derive. Notice that
the theoretical error that we calculated above is practically the same as the
actual error when n = 3.
158 CHAPTER 9. APPROXIMATING INTEGRALS
a. Using the error in Simpson’s rule, formula (5) on page 99, deter-
mine approximately how many subdivisions (and hence how many
computations) are needed to obtain 8 decimal place accuracy.
b. For completeness, also do part (a) using the left endpoint method
and the trapezoid method.
c. Show that the ratio of terms in the above series is less than 1/9.
d. Using formula (7) of Theorem 1 on page 138 with x = 1/9, deter-
mine how many terms are needed to approximate ln 2 to 8 decimal
places.
e. Now compare all four approximation techniques. Which method
is the most efficient?
a. Define and simplify P(x):= the 8th –degree Taylor polynomial for
ex .
√
b. Author P ( x) and integrate the result from 0 to 2. Simplify this
integral and then express the answer as a decimal.
c. Compute the√maximum value of the ninth derivative of ex on the
interval 0 to 2. Denote this maximum by M. (Note: This is the
M value associated with the Taylor polynomial √ p(x) in (11) on
page 115√ corresponding to the interval 0 ≤ x ≤ 2. The reason √
|ex − p(x)|
we use √ 2 and√not 2 is that if √ √ ≤ c for 0 ≤ x ≤ 2,
then |e x − p( x)| ≤ c for 0 ≤ x ≤ 2, i.e., for 0 ≤ x ≤ 2.)
d. In a manner similar to what was done in Section 9.4, find the error
in the approximation you obtained in part b.
R2 √
e. Have Derive evaluate 0 e x dx and then approximate it and
compare the answer with what you obtained in part b.
Instead of starting with the Taylor polynomial for ex , start with the
Taylor polynomial of e−x .
Chapter 10
10.1 Introduction
Graphs of the form y = f (x) or x = g(y) can be used to represent a wide
variety of curves in the plane, there are many important curves, such as
circles or ellipses, that cannot be represented by a single graph of this type.
More generally, imagine the curve traced out by an ant walking on a flat
surface. In this chapter we will introduce two techniques for plotting general
curves. One is the method of polar coordinates, which is a coordinate system
based on angles and distance from the origin. The other is the method of
parametric representation, which allows one to specify completely arbitrary
curves like the motion of a particles (or the ant).
161
162 CHAPTER 10. POLAR AND PARAMETRIC GRAPHS
[r, θ]
r
θ
other values of θ. Of course, Derive will plot the curve for us. We enter
either 1 + cos t or just highlight expression #1 in Figure 10.2, then click
the Plot button in the graphics window. So far, this is just like rectangular
plots except for the change in the coordinate mode. But now Derive will
prompt you for the parameter interval (interval of θ’s to use) and suggest the
default range of −π ≤ θ ≤ π. Since many of the standard examples of polar
curves involve the θ-variable only in the form of either cos θ or sin θ, it usually
suffices to only consider 0 ≤ θ ≤ 2π (or as Derive prefers −π ≤ θ ≤ π).
Of course, you can change it to whatever interval you want. For example, in
Figure 10.2 the range 0 ≤ θ ≤ π was used. You might want to plot the full
graph at this point by using the default range. The resulting graph heart
shaped curve is called a cardioid.
Tracing. It is important to actually see the curve being plotted but the
computer plots so quickly that it is nearly impossible to see it happen. De-
rive has an approach for “driving” around a curve called tracing. After
plotting the polar curves above select the Trace Mode option on the Options
164 CHAPTER 10. POLAR AND PARAMETRIC GRAPHS
menu (or just press F3 to toggle the Trace mode) and the cross will turn into
a box and it will be moved onto the last curve plotted. Now press and hold
down the right arrow key and watch the little car drive around the curve.
You can see the value of θ, which we can also interpret as time, as it increases,
as well as the r and θ coordinates, on the lower part of the screen. If you
have more than one graph you can switch between curves by using the up or
down arrow keys.
When plotting the cardioid a = 1 pay particular attention to the way the
plotting slows down as we approach the cusp. It turns out that the only way
for cusps or corners to occur in the graph, when r(θ) is differentiable, is for
the plotting to slow to a stop and then to start up again. This notion of
speed will be discussed in Section 10.5.
x2 y 2
− 2 = ±1.
a2 b
We need to discuss converting polar graphs to rectangular graphs and
vice versa. Figure 10.1 on page 162 makes it clear how to do this. The
algebraic relationship between the polar coordinates [r, θ] and the rectangular
coordinates (x, y) is given by the right triangle formed from the 3 points:
(0, 0), (x, y) and (x, 0). The equations are:
p
(1) x = r cos θ, y = r sin θ and r = x2 + y 2, θ = tan−1 (y/x)
but Derive does not simplify these by default. Instead we need to choose
Declare/Algebra State/Simplification and on the Trigonometry box we select
Expand. Simplifying these standard formulas will now yield the above.
Simplifying our rotated curve now yields: r 2 cos2 t−r 2 /2 = 1. Converting
back to rectangular coordinates we use (1) to replace r 2 cos2 t with x2 and
r 2 with x2 + y 2 . This yields the desired result; namely, rotating the graph
y = 1/x by 45◦ results in an equation x2 − y 2 = 2 which is a hyperbola.
cos2 (θ + π) = cos2 θ
and the angle θ comes from its polar coordinate representation. Note that
the modulus |z| is simply the r in the polar coordinate representation. Thus,
Of course, this makes the plotting problem harder since we probably wouldn’t
use the geometry of polar coordinates to plot points. The computer on the
other hand doesn’t use geometric consideration since it just plots lots of
points and connects them with line segments.
Let us consider the non-polar example
x(t) = 4 cos t, y(t) = sin t, where 0 ≤ t ≤ π.
We can plot n points in order by taking ti = ti−1 + ∆t where ∆t = (β − α)/n
and making a n × 2–matrix using the vector function. Enter [4*cos t,sin
t] and then use Calculus/Vector with Start: 0, End: π and Step: .2 (this
gives 16 points). Now we can plot this as usual in rectangular coordinates
(you will need to switch back to rectangular coordinates). To draw a curve,
select Options/Points again and set plotting mode to connected. Then, replot
the points. See Figure 10.4 and Load the file F-PARAM1.MTH for a
demonstration.
the parameter interval and then plot the curve. You might have thought that
Derive would plot the two functions 4 cos t and sin t since we know that this
happens for 3 or more functions in a vector. But when a vector contains only
two functions, it is treated as a parametric curve.
Looking at the picture you might have guessed that the curve in Fig-
ure 10.4 was an ellipse (even if you didn’t read the caption) because of it’s
oval shape. Of course, not all oval shaped curves are ellipses but indeed this
example is one since one easily checks that
x(t)2 y(t)2
+ 2 = sin2 t + cos2 t = 1 for all t
42 1
and hence the particle travels along the ellipse x2 /4 + y 2 = 1 centered at
the origin with semi-major axis 4 and semi-minor axis 1. Observe that this
information does not tell you how the particle travels around on this ellipse.
For instance, is it going clockwise or counterclockwise? Does it ever stop?
See Figure 10.4 again and try the slow down technique to understand
how the parametric curve can be thought of as a particle moving along a
curve (like a car traveling over a roadway). By using this technique it is
apparent that the motion is counterclockwise (as time t increases) and it
never stops. This interpretation will be extremely important in later courses
when Newton famous F = ma law is used to analyze the forces acting on a
moving particle.
equal time. Assuming this fact, then the particle needs to be faster near the
top and bottom since these points are closer to the origin and hence sweep
out less area. Whereas the left and right portions of the curve are further
from the origin and hence require less time to sweep out an equal amount of
area.
One can calculate the speed directly as follows: Over a small time interval
∆t the x-position changes by ∆x (= x(t + ∆t) − x(t)) and the y-position
changes by ∆y.p Thus, the distance traveled during that time interval is
approximately ∆x2 + ∆y 2 and hence the average speed is given by
p s 2 2
2
∆x + ∆y 2 ∆x ∆x p
= + ≈ x0 (t)2 + y 0 (t)2
∆t ∆t ∆t
Taking limits as ∆t → 0 leads to the formula:
p
(2) Speed at time t = x0 (t)2 + y 0(t)2 .
Derive has an alternate approach for curves described by [x(t),y(t)];
namely, one uses Calculus/Differentiate on the vector and then √ applies ABS
to the result. This works because ABS([a,b]) simplifies to a2 + b2 .
Another use of Derive’s tracing feature is for curves that retrace them-
selves and hence make motion on the curve difficult to see. Try the example,
x = sin t cos t and y = sin 2t for 0 ≤ t ≤ 2π. That is, plot the vector [sin
t cos t, sin(2t)]. Surprisingly, the picture is simply a line segment with
endpoints (−1/2, −1) and (1/2, 1). But how does the particle travel around
this curve? By pressing F3 and tracing the curve we see a back and forth
motion which reminds us of a swinging pendulum. In fact, by carefully ob-
serving the motion near the endpoints we see the particle slow down and
stop. Then, it turns around and goes back in the opposite direction gaining
speed as it approaches the center of the line segment and then slowing down
as it approaches the other endpoint. A point where the speed is zero is ac-
tually the only way a smoothly parametrized curve, i.e., one for which x(t)
and y(t) are continuously differentiable, can have cusps (like the cardioid) or
corners (as in this example) or otherwise exhibit nonsmooth behavior. Check
directly the speed at the endpoints.
As a last example, enter x = 2 cos2 t and y = 2 sin t cos t for 0 ≤ t ≤ π. In
this case we have another surprising picture of a circle, which we can verify
by showing
(x(t) − 1)2 + y(t)2 = 1 for all t.
10.5. PARAMETRIC CURVES 171
Two interesting features are that the complete circle is plotted with t in the
[0, π] (instead of requiring 0 ≤ t ≤ 2π) and also that a particle travels around
the curve with uniform speed. Observe this with the tracing technique and
then verify it directly using (2).
172 CHAPTER 10. POLAR AND PARAMETRIC GRAPHS
a. Plot the first two curves (e = .5 and e = .75) with −π < θ < π
and identify the conics.
b. Plot the curve with e = 1 with −3.10 < θ < 3.10. Can you
identify this conic?
10.6. LABORATORY EXERCISES 173
c. Plot the curve with e = 2 with −2.09 < θ < 2.09. Can you
identify this conic?
d. In the last plot, what is the significance of the number θ = 2.09?
What curves do you get when −π < θ < −2.10 or 2.10 < θ < π?
(Warning: If you try plotting with the default range [−π, π] it will
eventually graph the complete conic but it takes a very long time!
Press Esc if you can’t wait.)
The butterfly curve, Amer. Math. Monthly, vol. 96, May 1989, p. 442.)
It can be viewed on the Web as Figure 3 on our home page
http://www.math.hawaii.edu/lab/
10. Imagine a circle (or wheel) of radius one rolling along the x–axis at
unit speed. Now try to picture the path followed by a fixed point on
this circle as its rolls. This is the parametric curve in problem 9, it is
called a cycloid curve. It may seem a little surprising that the speed of
the point on the wheel is 0 once every time the wheel revolves even as
the center of the wheel travels at a constant speed.
a. Make a graph of the speed function (2) and determine how fast
the point on the wheel going when it is at its highest point? (Hint:
Plot the speed function and cycloid curve together on the same
graph.)
b. Load the file F-CYCL.MTH and plot expression #8 which contains
the parametric curves for 5 positions of the rolling wheel along
with a dot marking the particle’s position on the wheel. Then
plot the cycloid expression #3, see Figure 10.6 on the facing page.
11. Plot the parametric curve x = sin(π sin t) and x = cos(π sin t) for
−π ≤ t ≤ π.
a. What geometric object does this look like? Prove that your answer
is correct.
b. Using the trace feature to see how a particle following these para-
metric equations moves along this geometric object. Describe this
motion in words.
10.6. LABORATORY EXERCISES 175
c. Are there places where the point seems to have speed 0? Find a
formula for the speed of the particle at time t. At what times does
the particle have speed 0 and what is the position of the particle
at these times?
12. Two particles move in the plane. The motion of the first is described
by the parametric equations
Plot both of these curves. Find where the curves intersect. But just
because the curves cross does not mean the particles collide; they might
arrive at the intersection point at different times. Where do the parti-
cles collide?
176 CHAPTER 10. POLAR AND PARAMETRIC GRAPHS
Chapter 11
11.1 Introduction
Suppose that y is a function of x. A first order differential equation is an
equation which involves x, y and its derivative y 0 . An nth order differential
equation involves x, y, y 0, . . . , y (n) . For example, y 00 + xy 0 = x2 + 1 is a second
order differential equation.
Differential equations occur frequently in every field of science and engi-
neering, especially biology. Libraries have many volumes devoted to solving
differential equations (even for first order differential equations). In this
chapter we study first order differential equations and show some of the ap-
plications. One of the most important examples is population growth (of
humans, cells, radioactive substances, savings account balances, etc.)
We will show you how to get an exact solution to certain types of first
order differential equations and we will introduce slope fields and Euler’s
method for obtaining approximate solutions to more general first order dif-
ferential equations.
11.2 Examples
Population Growth. The standard model for population growth states
that the rate of change y 0(x), of the population size y(x), with respect to
time x is proportional to the population size at any given time. This means
177
178 CHAPTER 11. DIFFERENTIAL EQUATIONS
that y 0 (x) = ky(x) for some fixed constant k and all x. Now it is easy to
check that y(x) = y0 ekx satisfies the relations
dy
(1) = y 0 = ky y(0) = y0
dx
where we simplify our notation by dropping the explicit reference to the
variable x. Thus, the exponential function provides a model for population
growth. Recall from Problem 7 on page 66 that we compared a linear model
versus the exponential model above for the population of the US and found
a significant difference in the long run behavior with the exponential model
giving a much larger growth. This comparison was also made in Section 1.3
where it was shown that exponential growth eventually exceeds the growth
of any polynomial.
We now want to show that the above solution is the only solution to (1).
Suppose that y = u(x) is a solution then
Z x1 0 Z x1
u (x)
dx = k dx = kx1 .
0 u(x) 0
On the other hand, using the substitution y = u(x) in the integral above
yields that
Z x1 0 Z u(x1 )
u (x) dy u(x1 )
dx = = ln .
0 u(x) y0 y y0
Now, if we set the right hand sides of the two equations above equal to each
other we get
u(x1 )
y0 = kx1 .
Finally, we observe that since u(x1 ) > 0 and y0 > 0 we can remove the abso-
lute values in the above equation. Exponentiate both sides of the resulting
equation (remembering that eln a = a for a > 0) we get
u(x1 )
= ekx1
y0
which is equivalent to our solution above.
Less formally, this approach is quite simple. We rewrite
dy dy
= ky as =k
dx y
11.2. EXAMPLES 179
You should always check your answer to make sure it satisfies the differ-
ential equation and the initial conditions. In Derive you could differentiate
the answer and set that equal to xy 2 − 4x where you use your answer in
place of y. Now Simplify, you should get the value true which indicates that
Derive was able to verify that the equation is an identity. You might also
try this problem with pencil and paper ! The y-integral part involves some
work since this is a partial fractions integration problem.
dy
(2) = −k(y − ycool ) where y(0) = yhot > ycool
dt
and k > 0 is a constant which depends on the physical properties of the
pan, for example, copper cools faster than iron so the corresponding k-value
11.2. EXAMPLES 181
would be larger. Notice that the derivative dy/dt above is negative since the
temperature is decreasing.
We solve (2) by separating variables and integrating:
Z Z
dy
= −k dt + c .
y − ycool
Simplifying this equation in Derive, we then solve for the unknown constant
c by substituting in the initial conditions x = 0 and y = yhot . Note that we
assume the quantity y − ycool to be positive; hence, we can take its logarithm.
We now Substitute this value for c back into the relation above and solve for
y. This gives the solution
DF(r,x,x0,xm,m,y,y0,yn,n)
where the first argument r is f (x, y) and x0, xm, m represent the initial
and final x-values in a rectangular grid with m x-values plotted. Similarly,
y0, yn, n represent the initial and final y-values in a rectangular grid with
n y-values. Hence, the total number of line segments plotted will be m ·
n. In order that line segments are plotted, not just the endpoints, we put
the plotting window into connected mode by choosing Options/Points and
setting Connect to ‘Yes.’
As an example, we can take the cooling problem above, namely,
y0 + y = 1
df(1-y,x,0,4,8,y,0,4,8)
1
R
Notice that the simple differential equation y 0 = f (x) has solution y = f (x) dx so
that the class of differential equations contains all integration problems.
11.3. APPROXIMATION OF SOLUTIONS 183
to get the slope field. In the plot window select Option/Plot Color and set
it to ‘Off’ so that all slope lines will be in one color. Of course, if you like
colorful diagrams then you can skip that last step. Also choose Option/Points
to set the Connected Mode and to set Size to Small. Make sure to delete
all existing graphs and then plot the slope field. Try to picture the solution
though a given initial point (0, y0) by following the slope field. Finally, plot
some actual solutions that we obtained above using the DE function and see
how it conforms to the slope field. See Figure 11.1 where we have graphed
the solution y = 3e−x +1, which corresponds to the initial condition y(0) = 4.
Try several other initial conditions to see how the slope lines approximate
the solution.
∗
Another Population Model. In the model we used for population growth
we had
dP
= kP.
dt
184 CHAPTER 11. DIFFERENTIAL EQUATIONS
This works well for many populations. But the population cannot continue
to grow forever. When a country no longer has room for expansion the rate of
growth slows. For example, a bacteria culture in a petri dish will satisfy the
above differential equation for awhile, but as the dish fills the above equation
becomes invalid. Verhulst, a Belgian mathematician, proposed a model using
the differential equation
dP P
(5) = kP 1 − .
dt P1
so that the population starts out at a value of 1 and has a maximum value
of 5. We’ll guess as to how long of a time interval to look at, say 0 ≤ t ≤ 12.
Thus, we enter
df(p(1-p/5),t,0,12,12,p,0,6,6)
and Approximate to get a large direction field matrix. See Figure 11.2 on
the facing page for a picture of the direction field.
Looking at the picture of the direction field you should be able to see the
solution approximately by starting at the point (0, 1) and “following along”
in the direction of the tangent lines. The picture makes it appears that it
takes only about 6-7 time units before the population reaches its maximum
value of 5.
To solve (5) exactly we separate variables and integrate
Z Z
dP
= k dt + c .
P (1 − PP1 )
But, there are two problems with the above equation. We assume that
0 < P0 < P < P1 when t > 0 and so two of the logarithm have negative
arguments. As a result we modify the above to the equation
a power where the power is the above equation and then Simplifying. This
yields the equation
P P0 ekt
=
P1 − P P1 − P0
which Derive can solve for P . The result is
P0 P1 ekt
(6) P =
P0 ekt + (P1 − P0 )
Notice that P (0) = P0 , as we would expect, and that limt→∞ P (t) = P1 .
u(x) = y0 ekx
2
In Derive you can multiply both sides of an equation by an expression by right clicking
the equation, inserting with parentheses into the author box and then multiplying.
188 CHAPTER 11. DIFFERENTIAL EQUATIONS
as we claimed.
Equation (1) is a special case of the general equation
since (1) can be written as y 0 −ky = 0. Thus, in (8) the functions p(x) = −k,
q(x) = 0 and initial time x0 = 0. Any differential equation with the form
of (8) is called a linear first order differential equation. We will prove that
any such equation has a unique solution which is obtained in manner a similar
to the above. See Theorem 1 below for the formula for the solution.
The formula for the solution to (8) can be made into a Derive function
quite easily. This has been done in Derive’s utility file ODE1 with the
name LINEAR1. For convenience we have added this function to our utility
file ADD-HEAD but we use the shorter name DE. It has the form
where p and q are expressions in the variable x. The initial conditions are
y = y0 when x = x0.
y 0 = 2y where y(0) = 5
we rewrite the equation so that it has the form of the general first order linear
differential equation in (8):
and thus we can solve this equation with Derive by using the DE function.
We will use the variables yh and yc in place of yhot and ycool . So that these
11.5. LINEAR FIRST ORDER DIFFERENTIAL EQUATIONS 189
are treated as single variables (and not as y · h) we first Author the vector
[yh:=,yc:=]. Then we Author
de(k,k*yc,t,y,0,yh)
See Figure 11.4 for a demonstration of these functions and observe how
rapidly the temperature y(t) tends to the water temperature yc. Use Derive
to calculate limt→∞ y(t).
We now solve the general linear first order differential equation (8)
by proving the following theorem:
Theorem 1. Suppose that y(x) satisfies (8) where p(x) and q(x) are con-
tinuous functions of x. If y satisfies the initial condition y(x0 ) = y0 then
R
Z x Ru
− xx p(u) du p(v) dv
(12) y=e 0 · x
q(u)e 0 du + y0 .
x0
Rx
p(u) du
Proof. Let h(x) = e x0 . By hthe fundamental
i theorem of calculus,
d
Rx 0 d
Rx
p(u) du
x0
dx x0
p(u) du = p(x). So h (x) = dx e = p(x)h(x). Thus
If we multiply equation (8) by h(x) and use the above, we see that (h(x)y)0 =
R x sides of this from x0 to x and use the fact that
h(x)q(x). If we integrate both
h(x0 ) = 0, we get h(x)y = x0 h(u)q(u) du + C, or
Rx
Z x Ru
− p(u) du p(v) dv
y=e x0
· q(u)e x0
du + C .
x0
As we said, the solution (12) to the differential equation can be made into
a Derive function quite easily. You should look at the formula above and
see if you can write a Derive function that will produce the solution. Then
compare your answer with the following definition of the function DE.3
(13) de(p,q,x,y,x0,y0):= y =
e^(-int(p,x,x0,x)) * (int(q*^
^ e^int(p,x,x0,x),x,x0,x)+y0)
ma = mv 0 = mg − kv.
Solve this equation for v with v(0) = 0. Find limt→∞ v(t) (don’t in-
clude the v= part from above). Derive returns an expression contain-
ing SIGN(km) because it does not know that k and m are positive. Use
Declare/Variable Domain to declare that k is a positive real number,
and do the same for m. Now reevaluate the limit. Note that v never
exceeds this value, which is called the terminal velocity. No wind re-
sistance corresponds to k = 0. Find v in this case both by solving the
differential equation with k = 0, and by taking the limit of the general
solution for v found above as k → 0.
*6. Suppose the population growth of a small country satisfies (5) with
P1 = 10 and k = 0.05 (with population in millions). Plot the direction
field for this. (There are instructions for doing this in Section 11.3.)
Suppose P (0) = 2. Find P (20), P (50), and P (100). Graph P (t).
Adjust the scale of the graph so that you get a clear picture of the
nature of the population growth.
12.1 Introduction
In this chapter we consider second order differential equations of the form
It is easy to see that two initial conditions are needed because in the simply
case y 00 = f (t) one would simple integrate f (t) twice and that there would
be two constants of integration. Our applications involve moving physical
systems such as oscillating springs and swinging pendulums. These systems
can be analyzed using Newton’s Law of Motion: F = ma. Since the a in this
formula is acceleration which is the second derivative with respect to time
we see that Newton’s formula is actually a second order differential equation.
Free Fall The equation for the vertical height y of a free falling particle,
in the time variable t, is given by
where g is the force due to gravity (which we assume is constant). The initial
conditions in this case are the initial height y0 and the initial velocity v0 . By
integrating twice and solving for the initial conditions we get the solution
y = y0 + v0 t − 12 gt2 . Here are the steps:
195
196 CHAPTER 12. HARMONIC MOTION
Z
0
y = −g dt = −gt + c1 = −gt + v0
Z
1 1
y= −gt + v0 dt = − gt2 + v0 t + c2 = − gt2 + v0 t + y0
2 2
Notice that the constants are obtained by substituting t = 0, y = y0 and
y 0 = v0 .
In this chapter we will mainly be interested in second order linear differ-
ential equations of the form
(1) y 00 + by 0 + cy = 0
These equations are linear differential equations because of they only involve
y 00 , y 0 and y in a linear relation. It turns out that to solve these equations
we need to use the complex exponential function eαx which was discussed
earlier in Problem 6 on page 127 in Chapter 7. Here α = a + b i is in general
a complex number (see Section 5.4 on page 78 for a brief tutorial on complex
numbers).
12.2 Examples
Springs and Hooke’s Law. Suppose we have a mass m attached to the
end of a spring hanging from the ceiling. If we pull the mass down a little it
will bounce (oscillate) up and down. We image it moving along the y–axis
with y = 0 denoting the rest position. Newton’s law says F = ma where F
is the force on the mass and a = y 00 is the acceleration. A reasonably good
approximation of the force is given by Hooke’s Law which states
F = Fspring = −ky
where k is a positive constant. Since Fspring = F = ma = my 00 this gives y 00 =
−mk
y. After transposing the y term we get the following linear differential
equation:
k
(2) y 00 + y=0
m
As an example suppose we pull the particle down d units and let go.
Then, the initial condition would be: y0 = −d and v0 = 0
12.2. EXAMPLES 197
b 0 k
(3) y 00 + y + y = 0.
m m
The component of the force in the direction of the pendulum rod is can-
celled out by the rod. The force along the arc is −gm sin θ. Newton’s law,
F = ma says md2 s/dt2 = −gm sin θ. In terms of θ we get the differential
equation:
d2 θ g
(4) 2
= − sin θ
dt l
This is not a linear equation because of the sin θ. But the Taylor series is
sin θ = θ − θ3 /6 + · · · , so if θ is small we can approximate sin θ with θ. Using
this (4) becomes
d2 θ g
(5) + θ=0
dt2 l
If we start by pulling the pendulum back by an angle θ0 and letting go
we get the initial conditions:
dθ
θ = θ0 and = 0 when t = 0 .
dt
(6) r 2 + br + c = 0
Now, we take advantage of the fact that our differential equation is linear.
We claim that if y1 , y2 are solutions to (1) then so is any function of the form
y = C1 y1 + C2 y2 where C1 , C2 are arbitrary constants. This is easy to verify:
Finally, we take y1 = er1 t , y2 = er2 t above and then use the two arbitrary
constants C1 , C2 to solve for the two initial conditions:
r 2 − r − 2 = (r + 1)(r − 2) = 0 so r = −1, 2 .
Our general solution is therefore y = C1 e−t + C2 e2t and we need to solve for
the initial conditions. This leads to the pair of equations:
C1 + C2 = y(0) = 1
−C1 + 2C2 = y 0(0) = 0
r 2 + 1 = (r + i)(r − i) = 0 so r = ±i .
and in general,
e(α+β i)t = eαt (cos βt + i sin βt) .
It’s important to note (again from Problem 6 on page 127) that the differen-
tiation formula
d (α+β i)t
e = (α + β i)e(α+β i)t
dt
still holds, even for complex exponents, and hence
C1 + C2 = y(0) = 3
iC1 + (−i)C2 = y 0(0) = 0
C1 + C 2 = 3
C 1 − C2 = 0
mass-spring system (see (2)) and this couldn’t possibly have complex num-
bers in the solution. Let look at the computation:
3 3
y = eit + e−it
2 2
3 3
= (cos t + i sin t) + (cos t − i sin t)
2 2
= 3 cos t
and hence the i-term drops out of the solution. Although we were convinced
of this fact from a physical point of view it was far from obvious from an
algebraic point.
To see why this is always the case we make an important observation
about (1). Since we tacitly assume that the coefficients b, c are real numbers
observe that if y is a solution to (1) then so are the real part <y and the
imaginary part =y of y also solutions to this equation. Let’s verify this for
the real part (you may need to refer back to Section 5.4 and Problem 8 on
page 91). Recall that we have already shown that solutions can be combined
by multiplying them by constants and then adding them. Since <y = (y +
ȳ)/2 we need only show that given a solution y that it’s complex conjugate
ȳ is also a solution. Let verify this fact:
This means that instead of using complex exponentials for the two solutions
y1 , y2 and then solving for complex C1 , C2 so that C1 y1 + C2 y2 satisfies the
initial conditions we could instead solve the initial value equations
which would not involve complex numbers at all. Of course, since we ulti-
mately use Derive to do the algebra, i.e., DE2(b, c, t,t0, y0, v0) this
is mostly important for hand calculations.
Let us return to the general equation (1) with characteristic equation
r 2 + br + c = (r − r1 )(r − r2 ) = 0
Thus, as we showed above eαt cos(βt) and eαt sin(βt) are both solutions. Of
course, we can easily check this fact directly by substituting eαt cos(βt) for y
in (1) and show the left side simplifies to 0. Try this yourself.
When b2 − 4ac = 0, r1 = r2 we have a problem because we know longer
have two different solutions with which to solve for C1 , C2 . However, in this
case both er1 t and ter1 t are solutions; see Exercise 3 on page 213.
Summarizing, with r1 , r2 , α, and β as above, the general solution of (1)
are
Using this form of the general solution we solve for the two arbitrary
constants C1 , C2 using the initial conditions as above. Again, assuming you
have loaded the ADD-HEAD file, you can solve (1) subject to these initial
conditions by authoring
m are both positive numbers. You use the Declare/Variable Domain to tell
Derive that k is positive and do the same for m. Now simplifying the DE2
function gives the answer
√ !
k
−d cos √ t
m
Figure 12.2 shows the graph of this function when k = 2 and m = 1 and
d varies between −2 and 2 in increments of 0.5. Notice all of the graphs
cross the x–axis at the same place; that is, at the same time. So it doesn’t
matter how far the spring is pulled down it will takeq the same amount of
k
time to return to its original position. If we let ω = m then time to return
to the original position is 2π/ω independent of d. This is called the period
of oscillation. ω is called the angular frequency while the reciprocal of the
period, ω/2π is the frequency. In Exercise 4 on page 213 you investigate
what happens if we start with y0 = 0 but vary the velocity.
Damped oscillation. Recall that the frictional force due to the oil in the
vat is proportional to the velocity of the mass. This extra force changed the
differential equation (2) to
b 0 k
(14) y 00 + y + y=0
m m
where b, k, m are all positive constants.
Before we analyze the roots of the characteristic equation just imagine
the effect of going from a very small value of b, close to zero, to a large
value. Since the value of b represents how difficult it is to move in the oil
(this is essentially the viscosity of the oil) a small value would mean the
solution is similar the oscillating spring although it eventually slows down
and stops. A larger value will mean less oscillations and eventually there
will be no oscillations and the system will quickly come to a stop. Actually,
the mathematical model never actually stops, it just goes into exponential
decay!
Going back to the characteristic equation for this equation we have
b k
r2 + r+ =0
m m
Its roots are
√
−b ± b2 − 4km
2m
√
The solutions to (14) are given in (10)–(12). The sign of b2 − 4km
determines which equation applies. If b2 > 4km we say the spring is over
damped . In this case the solutions of (14) have the form of (10) with r1 and
r2 the solutions to the characteristic equation given above. Notice that since
b, k, and m are all positive
√
0 < b2 − 4km < b2 so b2 − 4km < b
and thus, both r1 and r2 are negative. So the general solution is the sum of
two decaying exponentials.
If b2 = 4km the spring is critically damped. The solutions are given
by (11). Again r1 is negative. If b2 ≤ 4km the spring is under damped. The
solutions are given by (12).
12.3. SOLVING LINEAR DIFFERENTIAL EQUATIONS 205
d2 θ g
(15) + sin θ = 0
dt2 l
Since this is not a linear equation we approximate sin θ with θ to get (5), i.e.,
d2 θ g
(16) + θ=0
dt2 l
In our illustration, we start by pulling the pendulum back by a small
angle θ0 and letting go we can solve the equation by authoring DE2(0, g/L,
t, 0, θ0 , 0). This gives the solution
r
g
θ0 cos t
l
Notice that the period depends on l but not on θ0 . This is why old clocks
often use pendulums. As the spring runs down the pendulum will continue to
swing with the same frequency (until it stops completely of course). You can
adjust the speed by making a small change in the length of the pendulum.
Later in this chapter we will show how to find approximate solutions to
exact pendulum equation (15).
y 0 = r(t, y, z)
(18)
z 0 = s(t, y, z)
These can be solved exactly if r(t, y, z) has the form ay + bz and s(t, y, z)
has a similar form. (When r and s have this form, the system of equations
(18) is called linear.) Since the examples we are interested are not linear, we
concentrate on finding approximate solution to (18).
In Chapter 11 we described Euler’s method for finding an approximate
solution of a single first order equation y 0 = f (t, y) subject to the initial
conditions y(t0 ) = y0 . We start with the point (t0 , y0). Since we know the
slope of y at this point is f (t0 , y0) we draw a short line segment from (t0 , y0 )
to (t1 , y1) = (t0 + h, y0 + hf (t0 , y0 )), where h is a small increment. The
(n + 1)st point is obtained the nth by
Figure 12.4 on the following page gives the direction field for the simple
differential equation y 0 = −4(t − 1). Of course we can find the solutions
by integration. If y(0) = 0 this gives y = 2t(2 − t), which we have also
graphed. If we use Euler’s method with h = 1/2 the first three points are
(0, 0), (1/2, 2), and (1, 3). As the graph indicates these points are not very
close to the true solution.
208 CHAPTER 12. HARMONIC MOTION
If instead of using the slope at (tn , yn ) we average this slope with the slope
at the next point (tn+1 , yn+1) we obtain a much more accurate approximation
of the solution. This is known as the second order Runge-Kutta method. The
precise formulae for tn+1 and yn+1 are
tn+1 = tn + h = t0 + (n + 1)h
(19) h
yn+1 = yn + f (tn , yn ) + f (tn + h, yn + hf (tn , yn ))
2
As we can see from Figure 12.4 this is much more accurate. If we also take into
account the slope at the midpoint of the two points we obtain the fourth order
Runge-Kutta method . This is usually just called the Runge-Kutta method .
where h is the step size and n is the total number of steps you want. When we
approXimate this we get a matrix of triples. To graph y(t) we use the function
extract 2 columns(m,1,2) (where m is the matrix we got). This gives a ma-
trix of pairs which we can plot. To plot z(t) we use extract 2 columns(m,1,3).
Returning to the predator-prey problem let’s look at the rabbits and foxes
problem with specific data for the constants in (17):2
We approximate this and then use extract 2 columns for columns 1 and 2
to see R(t). The result is graphed in the upper right window of Figure 12.5
on the following page. Extracting column 1 and 3 gives the fox population
graphed in the lower right. Notice both populations oscillate with the fox
population following the rabbit population. After the rabbits increase the
foxes will then increase but when the fox population gets large the rabbit
1
RK is the same as the one in Derive’s utility file ODE-APPR so the description of it
in Derive’s Help applies.
2
This example is taken from J. Callahan and K. Hoffman Calculus in Context,
W. H. Freeman, 1995.
210 CHAPTER 12. HARMONIC MOTION
population will decrease which in turn will cause the fox population to de-
crease and so on.
The window in the lower left of Figure 12.5 shows the results of extracting
columns 2 and 3. The point near the crosshair in that window is (2000, 10),
the initial rabbit and fox populations. At the beginning both the rabbit and
fox populations increase. When they reach the point furthest to the right
the rabbit population starts to decrease while the fox population continues
to increase. As we continue along the curve it spirals inward indicating that
the oscillation in the populations get smaller. In Exercise 6 on page 213 your
find the point to which the spiral approaches.
θ0 = w
(20) g
w 0 = − sin θ
l
As an example suppose g/l = 25. Then to get an approximate solution
of (15) we author and then approximate the following.
Figure 12.6 on the next page shows the resulting graphs for θ0 = π/8,
π/4, π/2, and 15π/16. The graph of the solution of (15), namely θ0 cos(5t), is
also shown on each graph. Looking at these graphs we can see several things.
First for θ0 = π/8 = 22.5◦ the curves are almost identical showing that using
(5) rather than (15) works well for small and even moderate angles. Even
for θ0 = π/4 = 45◦ , shown in the lower left, is fairly close to the true graph.
Since we are considering a pendulum without friction (undamped) we
expect that when we release it with an initial angle of θ0 it will swing to the
other side reaching the angle θ = −θ0 and then return back to the original
position with θ = θ0 . Then of course it will just repeat this. This means
that the solution of (15) will be periodic. The linear approximation (5) has
a shorter period than the true equation (15). This makes sense since the
magnitude of the force pushing the pendulum back towards its rest position
(θ = 0) is proportional to sin θ for the true equation and to θ in the linear
approximation and sin θ ≤ θ for θ > 0.
The lower right frame of Figure 12.6 on the following page gives the graphs
when θ0 = 15π/16. This corresponds to starting the pendulum almost at the
top. Notice that not only is the true period much greater than the linear
approximation but that the shape of curve is different.
a. y 00 − 2y 0 − 8y = 0 b. y 00 = −4y 0 − 4y
c. y 00 = −16y d. y 00 − 4y 0 + 5y = 0
3. Show that if the characteristic equation (6) has a double root, that is,
if r1 = r2 , then both y = er1 t and y = ter1 t satisfy (1) on page 196.
y 000 + by 00 + cy 0 + dy = 0
6. Suppose we want to find solutions to (17) such that both R(t) and
F (t) are constant. One (trivial) solution is R(t) = 0 = F (t) but we
would like something more interesting than this. If the populations are
constant then R0 (t) = 0 and F 0 (t) = 0.
a. Solve (17) for R and F when R0 (t) = F 0 (t) = 0. Hint: Use the
second equation of (17) to find R, substitute this into the first
equation and then solve for F .
214 CHAPTER 12. HARMONIC MOTION
b. The curve in the lower left window of Figure 12.5 on page 210
spirals inward. Use your answer from the previous part with the
constants on page 209 to guess where it is heading.
7. Suppose in our rabbits and foxes example instead of (17) we use the
simpler equations
R0 = kR − cRF
F 0 = dRF − eF
8. We saw in the lower right frame of Figure 12.6 on page 212 on page 212
that for θ0 = 15π/16 the solution to (15) and to the linear approxi-
mation (5) were quite different. In this exercise we will compare the
solution to (15) with a cos function of the same amplitude and period.
and w0 = 0 except that starting place for the graph will be different.
That is, the curves will be the same except one will be shifted to the
right. However if w0 is large enough, the pendulum will swing all the
way over the top.
a. Author
EXTRACT 2 COLUMNS(
RK([w,-25 SIN(θ)],[t,θ,w],[0,0,w0 ],.05,60),1,2)
EXTRACT 2 COLUMNS(
RK([w,-25 SIN(θ)],[t,θ,w],[0,0,10],.01,500),1,2)
Utility Files
This book defines 19 Derive functions. All of these functions are defined in
the utility file ADD-UTIL which is quietly loaded; which means the formulas
are not displayed, whenever you load the file ADD-HEAD. This file and its
companion, ADD-HEAD, can be downloaded from our web page
http://www.math.hawaii.edu/lab/
Listings of these files are given at the end of this appendix in the (hopefully
unlikely) event you have trouble downloading them. Additional information
on the use of the functions as well as examples are included on the web site
above.
How to use these files. When a student first starts to work on a lab
assignment he should:
1. Enter (Author) his name and the lab number as a comment. (Com-
ments are entered using the double quotes, ".)
2. Do File/Load/Math add-head.
3. Begin working on his assignment.
The file ADD-HEAD has only four lines. Two of these are comments
and one gives the variable syntax for the commands. The other line is
LOAD("add-util"). This automatically silently loads ADD-UTIL.1 . For
1
This assumes that this file is in the default directory. A default directory should be
established and all F-*.MTH files placed in that directory along with ADD-HEAD and
ADD-UTIL
217
218 APPENDIX A. UTILITY FILES
Home use of these files. If you have Derive for your own computer
you can install ADD-HEAD and ADD-UTIL. When Derive starts it has
a default directory it looks for MTH files. If you haven’t changed this it is
..\DfW\Math. That directory contains the utility files that come with the
system. You should add a new subdirectory, such as ..\DfW\Lab, using the
Windows file manager and then make that your default directory using the
File/Change Directory command in DfW. Next, you put both ADD-UTIL
and ADD-HEAD in this default directory and the above directions should
work fine. You should also add all the F-*.MTH files that are used in the
book’s figures to this directory. All of these files are available from our web
site.
the integral of the expression u, in the variable x over the interval [a, b], using
Simpson’s rule with n subdivisions. You Author SIMP(u,x,n,a,b) and press
.
Example 7. Suppose that you want to solve a Newton cooling type differ-
ential equation: y 0 = −(y − 2) with initial conditions y(0) = 4. You start by
manipulating the equation to the form y 0 + py = q where p = 1 and q = 2.
The function DE(p,q,x,y,x0,y0) solves this equation so we just substitute
in the right values which in this case means that we Author DE(1,2,x,y,0,4)
and press .
Example 8. Suppose that you want to look at the direction field for the
equation y 0 = r where r is an expression in x and y. You use the function
DF(r,x,x0,xm,m,y,y0,yn,n))
[SUBST(u,x,a):=,SECANT(u,x,a,h):=,TANGENT(u,x,a):=,CURVEFIT(x,data):=,SPLINE(~
x,data,m1):=,NEWT(u,x,x0,k):=,DRAW_NEWT(u,x,x0,k):=,DRAW_COMPLEX(v):=,BISECT(~
u,x,v0,k):=,LEFT(u,x,n,a,b):=,MID(u,x,n,a,b):=,RIGHT(u,x,n,a,b):=,TRAP(u,x,n,~
a,b):=,SIMP(u,x,n,a,b):=,DRAW_LEFT(u,x,n,a,b):=,DRAW_TRAP(u,x,n,a,b):=,DE(p,q~
,x,y,x0,y0):=,DF(r,x,x0,xm,m,y,y0,yn,n):=,EULER(r,x,y,x0,y0,xn,n):=]
LOAD("add-util")
SUBST(u,x,a):=LIM(u,x,a)
SECANT(u,x,a,h):=y=(SUBST(u,x,a+h)-SUBST(u,x,a))/h*(x-a)+SUBST(u,x,a)
TANGENT(u,x,a):=y=SUBST(u,x,a)+SUBST(DIF(u,x),x,a)*(x-a)
UNK(n):=RHS(VECTOR(SOLVE(upsilon=upsilon,upsilon),i,1,n,1)‘ SUB 1)
224 APPENDIX A. UTILITY FILES
ANS(x,data,ddata,a,n):=IF(DIMENSION(ans_:=SOLVE(EQNS(data,ddata,a,n),[a]))=0,~
[],POLY(x,(RHS(ans_)) SUB 1,n))
"Finds the polynomial of the right degree through the nx2-data matrix."
CURVEFIT1(x,data):=ANS(x,data,[],UNK(DIMENSION(data)),DIMENSION(data)-1)
CURVEFIT2(x,data,ddata):=ANS(x,data,ddata,UNK(DIMENSION(data)+DIMENSION(ddata~
)),DIMENSION(data)+DIMENSION(ddata)-1)
CURVEFIT(x,data,ddata):=IF(DIMENSION(ddata)>0,CURVEFIT2(x,data,ddata),CURVEFI~
T1(x,data),CURVEFIT1(x,data))
"Quadratic spline function interpolating data points with initial slope m1."
SLOPE(data,m1):=ITERATES([v SUB 1+1,2*(data SUB (v SUB 1+1) SUB 2-data SUB (v~
SUB 1) SUB 2)/(data SUB (v SUB 1+1) SUB 1-data SUB (v SUB 1) SUB 1)-v SUB 2]~
,v,[1,m1],DIMENSION(data)-1)
SPLINE(x,data,m1):=SPLINE_AUX(x,data,SLOPE(data,m1))
"Newton algorithm"
NEWT_ITERATES(u,x,x0,k):=ITERATES(x-u/DIF(u,x),x,x0,k)
NEWT(u,x,x0,k):=IF(k>0,NEWT_ITERATES(u,x,x0,k),?,SUBST(x-u/DIF(u,x),x,x0))
DRAW_NEWT(u,x,x0,k):=VECTOR([[v,0],[v,SUBST(u,x,v)],[NEWT(u,x,v),0]],v,ITERAT~
ES(NEWT(u,x,w),w,x0,k))
DRAW_COMPLEX(v):=VECTOR([RE(z),IM(z)],z,v)
BIS_AUX(u,x,a,b):=IF(SUBST(u,x,a)*SUBST(u,x,(a+b)/2)<0,[a,(a+b)/2],[(a+b)/2,b~
])
"Bisection method"
"Formula for the left-endpoint method for integrating u(x) over [a,b] with n ~
subdivisions."
LEFT(u,x,n,a,b):=(b-a)/n*SUM(SUBST(u,x,a+k*(b-a)/n),k,0,n-1)
"Formula for the midpoint method for integrating u(x) over [a,b] with n subdi~
A.2. LISTINGS OF THE UTILITY FILES 225
visions."
MID(u,x,n,a,b):=(b-a)/n*SUM(SUBST(u,x,a+(k+0.5)*(b-a)/n),k,0,n-1)
"Formula for the right-endpoint method for integrating u(x) over [a,b] with n~
subdivisions."
RIGHT(u,x,n,a,b):=(b-a)/n*SUM(SUBST(u,x,a+k*(b-a)/n),k,1,n)
"Formula for the trapezoid method for integrating u(x) over [a,b] with n subd~
ivisions."
TRAP(u,x,n,a,b):=(b-a)/(2*n)*(SUBST(u,x,a)+SUBST(u,x,b)+2*SUM(SUBST(u,x,a+k*(~
b-a)/n),k,1,n-1))
"Formula for Simpson’s rule for integrating u(x) over [a,b] with n subdivisio~
ns."
SIMP(u,x,n,a,b):=(b-a)/(6*n)*(SUBST(u,x,a)+SUBST(u,x,b)+2*SUM(SUBST(u,x,a+k*(~
b-a)/n),k,1,n-1)+4*SUM(SUBST(u,x,a+(k+1/2)*(b-a)/n),k,0,n-1))
"The box and trapezoid drawing functions used in the graphical demonstrations~
of numerical integration techniques."
D_BOX(x1,y1,x2,y2):=[[x1,y1],[x2,y1],[x2,y2],[x1,y2],[x1,y1]]
D_TRAP(x1,y1,x2,y2,x3,y3,x4,y4):=[[x1,y1],[x2,y2],[x3,y3],[x4,y4],[x1,y1]]
DRAW_LEFT(u,x,n,a,b):=VECTOR(D_BOX(t,0,t+(b-a)/n,SUBST(u,x,t)),t,a,b-(b-a)/n,~
(b-a)/n)
DRAW_RIGHT(u,x,n,a,b):=VECTOR(D_BOX(t,0,t+(b-a)/n,SUBST(u,x,t+(b-a)/n)),t,a,b~
-(b-a)/n,(b-a)/n)
DRAW_TRAP(u,x,n,a,b):=VECTOR(D_TRAP(t,0,t+(b-a)/n,0,t+(b-a)/n,SUBST(u,x,t+(b-~
a)/n),t,SUBST(u,x,t)),t,a,b-(b-a)/n,(b-a)/n)
DE(p,q,x,y,x0,y0):=y=(y0+INT(q*#e^INT(p,x,x0,x),x,x0,x))/#e^INT(p,x,x0,x)
"The direction field (DF) for the differential equation: y’=r(x,y) with a de~
termined by x0<x<xm with m subdivisions and y0<y<ym with n subdivisions."
DF(r,x,x0,xm,m,y,y0,yn,n):=VECTOR(VECTOR(SEG(LIM(r,[x,y],[x0+j*(xm-x0)/m,y0+k~
*(yn-y0)/n]),x0+j*(xm-x0)/m,y0+k*(yn-y0)/n,(xm-x0)/(4*m),(yn-y0)/(4*n)),j,0,m~
),k,0,n)
"The EULER function gives an approximate solution to: y’=r(x,y) with y(x0)=y0~
."
EULER(r,x,y,x0,y0,xn,n):=ITERATES(v+(xn-x0)/n*[1,LIM(r,[x,y],v)],v,[x0,y0],n)
226 APPENDIX A. UTILITY FILES
Appendix B
DfW Version 5
227
228 APPENDIX B. DFW VERSION 5
Next, notice that the comment text is now a specially formatted object
in the window and it is no longer necessary to use double quotes to add
remarks. Inserted text can be formatted in different fonts and colors so that
the algebra window can be given the appearance of a page in a textbook.
This is a powerful feature that makes it much easier to use the printout from
your Derive session in classroom notes. In addition, the plot window can be
embedded into the algebra window and hence be included with the printout of
the expressions and comments, making for a true workbook like appearance.
You can see the embedded graphics feature in Figure B.1. Notice that there
is no Plot window showing.
B.1. WHAT’S NEW IN DFW5 229
231
232 INDEX
series
convergent, 131
divergent, 131
geometric, 131, 132
harmonic, 146
Simpson’s rule
derivation, 104–105
error estimate, 99
slope field, 182
solving equations, 6–8
numerically, 30
spline functions, 62
subscripts, 20
super attractor, 77
tangent line, 40
trapezoid rule
error estimate, 99
vectors, 18
web page, xv