"Mathematics with Maple"

S.Duzhin

Lesson 2: Number Sequences and Limits

Mathematics: sequence, limit, rate of growth.

Maple: seq, limit, proc.

An important remark: if you remember, in lesson 1 we had a command

plot3d(cos(x)+cos(y)+cos(x+y),x=0..2*Pi,y=0..2*Pi);

that didn't work well, because the letter y had been assigned some value

(we had y:=x^x somewhere earlier in the same worksheet).

This is a typical error in Maple. To correct it, use the following command of unassignment:

y:='y';

After this command, the variable y 'forgets' its previous value and becomes simply y. In this lesson, we will use this anassignment command quite often, because after every loop sentence (like "for i from 1 to 10 do ...") the loop variable (i) preserves some constant value. If you see that some command does not work correctly, try to see the current values of the variables that it contains.

Example:

> n:=7;

> sum(n,n=1..7);

Now let us begin!

In this lesson, we will review the notion of limit of a numeric sequence and study how to work with sequences and find limits with Maple.

Sequences in mathematics and in computer

First of all, let me explain the difference between the mathematical notion of a sequence and the notion of a sequence in Maple.

Sequences in mathematics are, by definition, infinite: x_1, x_2, x_3,...,x_n,... -- where n is assumed to go to infinity.

Of course, we cannot insert an infinite object into a computer.

In Maple, a sequence is a finite collection of expressions (e.g., numbers) separated by commas. The Maple function seq is used to create sequences:

> s:=seq(sin(1/n),n=1..15);

To understand the behaviour of a sequence a[n], it is useful to plot it as the set of points (n,a[n]) in the plane.

Before giving the plot command to Maple, let us prepare the sequence that consists of pairs:

> sp:=seq([n,sin(1/n)],n=1..15);

Now we can plot it:

> plot([sp],style=POINT);

You see that points in the picture get closer and closer to the line y=0, so it seems that

> n:='n':

> Limit(sin(1/n),n=infinity)=0;

However, this is not a proof. To prove this fact, you have to use the epsilon-delta argument studied in the first-year course of calculus. Maple also knows some standard methods of evaluating limits, so you can use the following command to find the limit:

> s:=limit(sin(1/n),n=infinity);

Now you may be 99% sure that the sequence really has limit 0. (99%, not 100% because every software package has bugs).

In the previous command, we have used the word 'infinity'.

Let us check whether Maple understands it correctly:

> infinity;

Yes, that's OK. It does understande the infinity.

But can it make an infinite sequence?

> n:='n':

> seq(n,n=1..infinity);

No, it cannot. This can't be helped.

Rate of growth

The limit of a numeric sequence can be not only an ordinary number, it also can be infinity.

Different sequences that grow to infinity may grow with different rates. For example, the sequence n^2 grows faster than log(n) but slower than 2^n, because

> limit(n^2/log(n),n=infinity);

> limit(n^2/2^n,n=infinity);

The next example is not so evident, because in the beginning (for n<1000) the values of 1000^n are much bigger than the values of n! But the second sequence has bigger rate of growth and finally it overcomes:

> limit(1000^n/n!,n=infinity);

However, in some books on celestial mechanics the series 1000^n/n! is considered to be divergent, and the series n!/1000^n convergent (see H.Poincare, "La Science et l'Hypothese")!

Two sequences are said to have the same rate of growth if their ratio tends to a constant different from 0 and from infinity. Example:

> n:='n':

> limit((2^n+3^n)/3^n,n=infinity);

The sequence 2^n+3^n has the same growth as 3^n, because the summand 2^n is negligibly small in comparison with 3^n.

The notion of the rate of growth of sequences is very important in computer science. Typically, the sequence considered is the number of operations N(p) in the algorithm as a function of the size of input data p.

In any book or article about complexity of algorithms you can see such phrases:

"This algorithm is polynomial" (meaning that N(p) < C*p^k for some k).

or:

"This algorithm has exponential complexity" (meaning that N(p) < C*exp(p)).

"Double exponent" means something like 2^(2^p).

Polynomials, logarithms and exponents in various combinations form a set of standard rates of growth, so if you have to describe the rate of growth of a sequence, you have to find a function from this standard set which has the same rate of growth as the given function.

Fibonacci numbers

Let us find the rate of growth of the famous sequence of Fibonacci numbers defined by the conditions that its first and second numbers are equal to 1 and every next number is the sum of two preceding numbers.

Let us write the definition of Fibonacci numbers as a Maple procedure (this is our first example of Maple programming):

> f:= proc(n) if n=1 or n=2 then 1 else f(n-1)+f(n-2) fi end;

Let us check it:

> seq(f(n),n=1..15);

Looks correct. The numbers are growing rather fast. How fast? To get an idea, let us plot the ratios f[n+1]/f[n] (do not forget to unassign n first, i.e. execute the second line, then the first line!):

> plot([seq([f(n),f(n+1)],n=1..10)],style=POINT,scaling=constrained);

Aha. The points are on one line -- so it seems.

Let us compute the ratios f[n+1]/f[n] numerically:

> n:='n':

> seq(evalf(f(n+1)/f(n)),n=1..20);

It certainly looks very suspicious: the values approach the number 1.618033...

Actually, this is true indeed, and the exact value of this limit is the famous "golden ratio"

> (1+5^(1/2))/2;

> evalf(%);

There is a very nice proof of this fact based on the notion of linear dependence in the space of numeric sequences and we will study it some other day, but now I will show you how to rigorously prove this fact using Maple.

In Maple, there is a procedure "rsolve" that solves recurrence relations in much the same way as the procedure "solve" solves usual equations.

If you write the defining relations for Fibonacci numbers as the first argument of "rsolve" and the name of the unknown function as the second argument, you will get the solution:

> rs:=rsolve({F(n)=F(n-1)+F(n-2),F(0)=0,F(1)=1},F(n));

Oh, we have forgotten to unassign n again! Let's do it, and then return to "rsolve" once more:

The formula looks a bit complicated, so let's try to simplify it:

> normal(rs);

Now it only remains to note that the second summand tends to 0, because

2(1+5^(1/2)) < 1, and 2/(-1+5^(1/2)) is exactly the golden ration. Let us see if Maple can do it, too.

The quotient of two successive elements of the sequence:

> q:=subs(n=n+1,rs)/rs;

> q:=simplify(q);

> limit(q,n=infinity);

> evalf(%);

Sequences of complex numbers

One of the most beautiful formulas of elementary mathematics is the following formula due to L.Euler:

> E:=exp(1);

> E^(Pi*I);

> simplify(%);

It relates three most famous and completely independent numbers: E (the base of natural logarithms), Pi (the length of circumference of radius 1) and square root of -1.

We will try to visualize this formula, using the definition of the exponent as a power series expansion:

> E^x=Sum(x^n/n!,n=0..infinity);

Let s(n) be the partial sum of this series from 0 to n:

> s:= n -> sum((Pi*I)^k/k!,k=0..n);

Here are the first 15 numeric values of s(n):

> for n from 1 to 15 do evalf(s(n)) od;

This sequence certainly looks like tending to -1+0*I. To make this observation more transparent, let us plot the complex numbers in the plane. Here is the list of pairs (x,y) consisting of the real (Re) and imaginary (Im) parts of our numbers:

> pts:= [seq([Re(evalf(s(n))),Im(evalf(s(n)))],n=1..15)];

Now:

> plot(pts);

The first point is in the upper right-hand corner of the picture. Every point is connected to the next point by a straight line.

To see the behaviour of this sequence in more detail near the limit point, let us eliminate the first point from the list (the picture will show the remaining points in a "zoom"). The Maple construction used below, op(2..15,pts), means that we pick elements from 2 to 15 from the list.

> plot([op(2..15,pts)]);

To make bigger zoom, move the cursor back into the line above and change 2 by 3, then by 4 etc. and repeat the command.

Finally, let us see whether Maple can find the limit honestly:

> limit(s(m),m=infinity);

Oops! That didn't work. Let us try to evaluate that limit numerically at least:

> evalf(%);

>

That's better than nothing!

Exercises

1. Square root extraction.

Let a be a positive real number.

Set x[1]:=1, x[2]:=(x[1]+a/x[1])/2, ..., x[n+1]:=(x[n]+a/x[n])/2,...

Prove that limit(x[k],k=infinity) = a^(1/2);

2. Let a[n]:=(1+1/2)*(1+1/4)*(1+1/16)*...*(1+1/2^(2^n)).

Find limit(a[n]).

Hint. It is easier to solve the problem without the computer:

consider (1-1/2)*a[n].

3. Prove that the sequence a[1]:=0; a[n+1]:=(a[n]+6)^(1/2)

has limit and find it. (Try to guess the answer with the help of Maple, then prove it mathematically).

4. Let f_n(x) be the n-th iteration of the function sin(x):

f_n(x):=sin(sin(...(sin(x)...))) (the function sin is applied n times). By experimenting with Maple check that the sequence f_n(x) for any x between 0 and 1 has growth rate n^(-1/2):

f_n(x)=C*n^(-1/2).

Find a numeric value for C.