Math 166
 Computer Lab # 8
 Sequences and Series

New Mathematica  syntax introduced:  
  [ ]  for storing and accessing elements of a sequence.
  g[n_]:= g[n] = relationship for implicit sequences.
  Limit[ f[x], x->Infinity] to evaluate a limit.
  Sum[ f[n], {n,1,Infinity}] to evaluate an infinite series.
  Product [ f[i], {i,1,n}] to evaluate the product 
     f[1]*f[2]* ... *f[n].

For more information about how to use Mathematica  with sequences and series, refer to the discussion in the "Commands" file, available on the Mathematics Delta server in the Mac Lab.

Introduction

In this lab, we consider:
  a sequence, an unending list of numbers, and
  a series, the sum of an unending list of numbers.
 
We are interested in the behavior of sequences and series as we "approach" infinity.

We first look at three ways sequences can be defined, and how Mathematica  can help us.

Then we look at a typical series, and show several ways that Mathematica can help us apply the tests for convergence, and  determine the limiting value of the series.

Simple Sequences Defined by a Formula

A sequence is an unending list of numbers.  Each number in the sequence is called a term, and is identified by its order in the sequence.  Thus, in the sequence
  3, 5/2, 7/3, 9/4, 11/5, ...
the third term is 7/3.

The simplest way to describe a sequence is to write a formula for the n-th term.  Here, for example, the sequence above has n-th term (2 n+1)/n.  We use n instead of x as the argument, to remind ourselves that we are really only interested in the value of f at whole numbers:

Clear[f, n]

f[n_] := (2*n+1)/n

To find the value of the 10-th term of this sequence is easy:

N[f[10]]

We can use Mathematica commands like Table to print out the terms in the sequence:

Table[ {n, N[f[n]]}, {n,1,15}];
TableForm[%]

We can use ListPlot to visualize the terms in the sequence.

values=Table[ {n, f[n]}, {n,1,50}];
p2=ListPlot[values,PlotRange->{0,3}]

Because we have a formula, it's easy for Mathematica to check the limit of the sequence.

Limit[ f[n], n->Infinity]

Exercise:  Consider the sequence ((n+1)/n)^n. Make a table of the first 10 values, and a plot.  What is the limit of the sequence as n goes to Infinity?

Clear[a,n]
a[n_]:=( (n+1)/n)^n
Table[  {n,a[n]},{n,1,10}]//N;
TableForm[%]
values=Table[  {n,a[n]},{n,1,50}]//N;
ListPlot[values]
Limit[ a[n], n -> Infinity]

Implicit Sequences Defined Recursively

Sometimes a sequence is not defined by giving a formula for each term.  Instead, a sequence is defined implicitly, by giving several starting values, and a way to produce the next value from the previous values.  For instance:
 
  2, 5, 7, 12, 19, 31, ... 
    Term n is the sum of Term(n-1) and Term(n-2).
    
  5, 1, 3, 2, 2.5, 2.25, 2.375, ...
    Term n is the average of Term(n-1) and Term(n-2).
    
  1, 1, 2, 1+sqrt[2], ...
    Term n is the sum of the square roots of Term(n-1) and Term(n-2).

If we want to know the value of the n-th term of an implicit sequence like this, we actually have to work out the values of all the previous terms to do so.

For instance, to set up the last sequence, we must tell Mathematica  the first two values of the sequence, and the rule for getting the third value from the first two, the fourth value from the second and third, and so on:

Clear[a,n]

a[1] = 1;
a[2] = 1;

a[n_] :=  a[n-1]^(1/2) + a[n-2]^(1/2)

We can now generate the terms in the sequence.  For example, enter the commands defining a[n] above and the following command.

Table[{n,a[n]},{n,1,20}]//N;
TableForm[%]

You may notice that Mathematica  took some time in responding with the answer.  The reason for this is that when Mathematica  goes to calculate a[10] for example, it needs a[9] and a[8]; but to calculate a[9] it needs a[8] and a[7], etc.  This means that when it determines a[10] it recalculates a[1] through a[9].  There is really no need to do this since we have already calculated these values.

Mathematica  has a way that will allow us to define a relationship for such an implicit sequence and to store all values that it finds; thus we can  compute any term we want without recomputing the previous values.  To do so, we use a new notation.   Look at this notation carefully, and notice how it differs from the usual function definition!  

Instead of 
  a[n_] := formula
we have 
  a[n_] := a[n] = formula

(* These Mathematica commands set up a sequence *)
(* which needs two initial values to get started. *)

Clear[a,n]

a[1] = 1;
a[2] = 1;

a[n_] := a[n] = a[n-1]^(1/2) + a[n-2]^(1/2)

Table[{n,a[n]},{n,1,20}]//N;
TableForm[%]

Now enter the new definition for a[n] above and see how quickly  Mathematica  responds when you ask it to set up the table compared with defining a[n] in the usual manner.

Exercise: Consider the sequence defined as follows:
  a(1) = 1,
  a(n) = 1 + 1/a(n-1)
  
Define the nth term a[n] using the new notation defined in the box above.  Make a table of the first 20 terms of a, and a plot.

Clear[a,n]

a[1] = 1; 
a[n_] := a[n] = 1+1/a[n-1]

values=Table[{n,a[n]},{n,1,20}]//N;
TableForm[%]
ListPlot[values]

Newton's Method Produces a Sequence

Newton's method tries to find a solution of an equation f(x)=0.  For any starting point x0, Newton's method tells us how to produce a sequence x1, x2, x3, ...  

The behavior of Newton's method depends completely on the starting point x0 that we give it.  If x0 is close to a correct solution, the sequence will usually converge to a solution of the equation,
but a "bad" starting point will cause the sequence to diverge (increase or decrease without limit), or to wander aimlessly without ever converging.

Newton's method proceeds as follows:

  1) Choose an estimate x[0] for the root.
  2) The next estimate is 
        x[n+1] = x[n] - f[x[n]] / f'[x[n]] 
  3) If f[x[n+1]] is small enough, stop.
     If it's too big, stop, you're diverging.
     If you've taken too many steps, stop.
     Otherwise, try step 2 again.

Example:  Apply 5 steps of Newton's method to find a solution of the equation f(x)=0, where f(x)=cos(x)-x.

To get a good starting point, we first plot cos(x)-x:

Plot[ Cos[x]-x, {x,-Pi,2*Pi} ]

From the graph, it looks like a good starting estimate would be x=1.  The following Mathematica commands set up Newton's method, and will save any computed entries as a sequence.  Note that once again we have used the new notation for defining x[n] recursively.

(* Newton's method for f(x)=0. *)
(* Just define f, and a starting value x[0] *)

Clear[f,n,x]

f[x_]:=Cos[x]-x

x[0] = 1;

x[n_] := x[n] = x[n-1] - f[x[n-1]] / f'[x[n-1]]//N

This new notation makes tables easy.  Note that we have printed the table using more digits of accuracy.

results = N[ Table[ {n, x[n], f[x[n]]}, {n,0,5}], 20] ;

TableForm[ results, 
  TableHeadings->{ None, {"n", "x[n]", "f[x[n]]"} } ]

We'd expect a graph of the x values to show convergence to the solution:

xvals = Table[{n,x[n]}, {n,1,10}];
ListPlot[xvals,PlotRange->{.5,1}]

and a graph of the f values to converge to zero.

Exercise:  Use Newton's method to find a zero of f(x)=E^x-4x which is nearest 2.  Use  2 as your starting value.  Compare your results with the FindRoot command.

Clear[f,n,x]
f[x_]:=E^x-4*x 
x[0] = 2;
x[n_] :=   x[n] = x[n-1] - f[x[n-1]] / f'[x[n-1]] //N
xvals = N[ Table[{n,x[n]}, {n,1,5}], 20];
TableForm[%]
N[ FindRoot[f[x]==0,{x,2}],20 ]

Infinite Series

An infinite series is an expression that represents the sum of an infinite sequence of quantities.  There is a standard mathematical notation for infinite series.  For instance, to represent the infinite sum S of the quantities 1, 1/2, 1/4, 1/8, ..., we would write

É
S =  Ö  1 / (2^i)
    i=0

Using limits, it may be possible to assign a finite value to an infinite series.  We do this by considering the long-term behavior of partial sums.  For instance, let us write Sn to represent the sum of the first n terms of the series.  That is,

n
Sn =  Ö  1 / (2^i) = 1 + 1/2 + 1/4 + 1/8 + ... + 1/2^n
     i=0

Certainly, for any n, we can compute this quantity.  The first few values of Sn are 1, 1.5, 1.75, 1.875, ..., and we may guess (correctly) that the value of Sn never reaches 2, but can be made as close to 2 as we want by choosing n large enough.  In other words,

  limit  Sn = 2
  n->É

In some cases, the limit is infinite, and we say the series diverges  to infinity.  In other cases, there is no limit at all, and the series cannot be assigned any value.

Example: Use the idea of partial sums to find the limit of the series whose i-th term is 1/(2^i):

The sequence of partial sums is:

  S1 = 1                    = 1      = 2 - 1
  S2 = 1 + 1/2              = 3/2    = 2 - 1/2
  S3 = 1 + 1/2 + 1/4        = 7/4    = 2 - 1/4
  S4 = 1 + 1/2 + 1/4 + 1/8   = 15/8   = 2 - 1/8
  
and so, in general, we can see that Sn = 2 - 1 / 2^(n-1)

Therefore, we can ask Mathematica for the limit of these partial sums:

Limit[ 2-1/2^(n-1), n->Infinity]

And we conclude that the infinite series has a value of 2.

Instead of figuring out the partial sums, and then looking at the limit, we may be able to ask Mathematica to sum the series to infinity right away.

Example:  Try using the Sum command with an infinite upper limit to evaluate the series whose i-th term is 1/(2^i):

Sum[ 1/(2^i), {i,0,Infinity} ]

This seems to be getting us nowhere!  Actually, Mathematica has computed the answer, but we have to force it out with the N command!  (In cases where this trick doesn't work, you can sometimes still get a reasonable answer by using the NSum command).

N[Sum[ 1/(2^i), {i,0,Infinity} ]]

Therefore, we may conclude that the sum of the series is 2.

Exercise:  Ask Mathematica for the limit of the series whose n-th term is 3 / (4 n^2 - 1) and whose first term is 1.  (Be prepared to use the N[] command!)  Make a table of the sum of the first m terms for m=10,20,...,100.  Compare your answers.

Sum[ 3/(4*n^2-1), {n,1,Infinity} ]//N

Table[{m,Sum[ 3/(4*n^2-1), {n,1,m} ]},{m,10,100,10}]//N;
TableForm[%]

The Integral Test for Infinite Series

Again, if we are only interested in knowing whether a series has a limit, we can sometimes answer this question, and even get a ballpark figure for the value of the series, by using the integral test.  The integral test states that if f(x) is positive, continuous, and decreasing for x Ñ 1, then the infinite series

É
S =  Ö  f(i)
    i=1

and the improper integral

É
I = Ü    f(x)   dx
      1

must both diverge or both converge.

Example: Show that the infinite series ln(n)/n must have a infinite value.

We have f(x) = ln(x)/x.  Set up the improper integral

É
I = Ü    ln(x)/x   dx  = limit  Ü1b ln(x)/x   dx
      1                 b -> É

and note that for any finite upper limit b, we can integrate the problem by setting u=ln(x), du = 1/x dx:

b                    ln(b)                  ln(b)
 Ü    ln(x)/x   dx  = Ü   u du  = [ 1/2 u^2 ] 
   1                    0                      0

= 1/2 ( (ln(b))^2 - 0).

Now as b increases to infinity, so does ln(b), and so does 1/2(ln(b))^2, so the improper integral diverges.  The integral test tells us that the corresponding infinite series must also diverge.

We can ask Mathematica  to integrate and take the limit for us:

Limit[ Integrate[ Log[x]/x, {x,1,b} ],  b->Infinity ]

Exercise:  Perform the integral test on the series whose n-th term is 3 / (4 n^2 + 1).

Limit[ Integrate[3/(4*x^2+1), {x,1,b} ],  b->Infinity ] 
N[%]

Math 166
Computer Lab #8
Name:_ _ _ _ _ _ _ _ _ _
Instructor:_ _ _ _ _ _ _ _ _

1. On January 1, Leo buys a newborn gremlin at the pet store.  The owner tells Leo that after two months, the gremlin will give birth to a new gremlin baby every month for the rest of its life.  On February 1, he still has just 1 gremlin, but on March 1 his gremlin gives birth and he has 2.   In April, his first gremlin gives birth again and he has 3 gremlins.  In May, his first gremlin gives birth, but so does his second one (it's now two months old), so he has 5 gremlins.  Things begin to get out of hand.  In June, gremlins #1, #2, and #3 give birth, making 8.  Leo begins to worry, because he only has room for 10,000 gremlins.

a.  We want to write a sequence for the number of gremlins alive on the first of each month.  The sequence begins

  J  F  M  A  M  J  J  A  S  O  N  D  (the months)
  1, 1, 2, 3, 5, 8...                 (the gremlins)
  
What is the next term in the sequence, for July?  In June, there were 8 gremlins alive, so in July there will be 8 gremlins, plus any newborns.  How many newborns are there?  Every gremlin at least two months old gives birth.  There were 5 gremlins alive in May (two months ago) so exactly 5 gremlins will give birth on the first of July.  So in July, there are 13 = 8 (old) + 5 (baby) gremlins.

a. What is the gremlin population for August and September?

b. What is a simple rule for the monthly population, given the population for the two previous months?

c. The number of gremlins alive on December 1 (the 12-th term) is 233.  What does your rule predict for December 1?

d. Copy the Mathematica commands labeled 
  (* These Mathematica commands set up a sequence *)
  (* which needs two initial values to get started. *)
and modify them for the gremlin sequence.  When will Leo have 10,000 or more gremlins?

2. Let us consider using Newton's method to find a solution of the equation f(x) = x^3 - x.

a) Plot f(x).  Use the Solve command to find the exact values x where f(x) is zero.

b) Newton's method will fail whenever f'(x) is zero.  Why?  What are the exact points where this will occur?

c) Copy the lines of Mathematica commands labeled 
  (* Newton's method for f(x)=0. *).  (They're in a box.)
Modify the commands for our function f(x), and start the iteration at x[0]=2.  Make a table of n, x[n], and f[x[n]] for n from 0 to 5.

d)  Find a starting point which causes the iteration to converge to x=0.  You might start with a value very close to zero; experiment to see how much larger you can make your starting point before Newton's method stops converging to x=0.

Make a table of your final results.

3.  The numerical approximation of  º  is a problem that has concerned mathematicians from ancient times to the present.  Archimedes used inscribed and circumscribed polygons to estimate its value.  In later times many other mathematicians including Newton and Euler proposed methods yielding improved estimates.  The discovery of calculus in the 17th century resulted in a host of formulas from which approximate values of  º  can be calculated.  In this problem you will compute several approximations to º.  In another lab we will consider other approximations to  º.

A famous relation involving  º  was discovered by Euler in 1736, namely
    º2/6 = 1  + 1/4  +  1/9  +  1/16 +  ...  +  1/n2 + ...

a. Use Euler's formula for  º2/6  to obtain approximations to  º.  In particular choose n=100,200,400, 1000.   Make a table of n, your approximation of º2/6, the corresponding estimate for º, and your error in approximating º.  (Remember that Mathematica knows what Pi is).

b.  The area of a circle of radius one is  º.  Use this to give an expression for  º  in terms of a definite integral.  Approximate this integral using a left- and right- Riemann sum with n=100, 200, 400.  Find the average of the two sums. Make a table of your results and the error in each approximation.

c.  Another expression for  º  is Wallis' formula, discovered in 1655:
 
    º/2 = (2 / 1) (2 / 3) (4 / 3) (4 / 5)  (6 / 5) (6 / 7)  (8 / 7)...   
    
    To approximate  º  from this expression we need to carry out the indicated multiplication, terminating the process after a suitable number of terms.  The Mathematica  command Product can be used.  The simplest way is to realize that we can group the terms in pairs; in this case each pair of terms is of the form
      (2k)(2k) /  ((2k-1) (2k+1) ) .
      
Use Wallis' formula with 100, 200 and 400 terms to approximate  º.  Make a table of your results and the error in each case.

d.  Compare your three results.

4. Let us consider the series

É
S =  Ö  1  /  ( n (n+1) )
     n=1

a) Use the Integral Test  to show that the infinite series must converge to some finite value.  Note:  you may have to help Mathematica  evaluate the limit by using properties of logarithms!

b) Make a table of the partial sums of the series for n=50, 100 and 200.

c)  Use the Mathematica  Sum command to sum the infinite series.  If you don't get a numerical answer, use N[] to get one. Does your value agree with your result in b?   Does it make sense that the sum of your infinite series is different from the value  your improper integral in (a) converges to?

d) (A telescoping series): The n-th term of our sequence is 1/(n*(n+1)), but this can be rewritten as 1/n  - 1/(n+1).  (You can show this by using the method of Partial Fractions!). 

Therefore the terms of the sequence can be rewritten as
  1 - 1/2,  1/2 - 1/3, 1/3 - 1/4, 1/4 - 1/5, ...
Since S3 is the sum of the first three terms, we can see that
  S3 = 1-1/4. 
Think about what S4 would be, and S5.  Now use this idea to write the n-th partial sum Sn in a simple form.

Since the value of the infinite series is the limit of the partial sums, use your simple form of Sn to figure out the value of the infinite series.  Does it agree with your result in (b)?