Sequences 215 problem 02.mws

1. Recursively defined sequences.

A sequence is recursively defined if it the value of a term in the sequence is determined by a rule which determines the value of a term in terms of the values of earlier terms, with some of the first few terms of the sequence possibly being given in advance. For example, if a[1] = 1 , and a[n+1] = a[n]+1 defines the sequence whose value at any positive integer n is simply a[n] = n . To define this sequence recursively in Maple, we use a procedure with a remember option. The remember option makes the procedure run more efficiently. You might want to look up remember in the help to understand what is going on.

> a:=proc(n)

> option remember;

> a(n-1)+1;

> end;

a := proc (n) option remember; a(n-1)+1 end proc

However, this recursion will not work unless we also tell Maple the initial value of the sequence.

> a(1):=1;

a(1) := 1

Let us verify that Maple will compute the correct value for n = 10 .

> a(10);

10

A very famous recursively defined sequence is the Fibonacci sequence, which is determined by a[n+1] = a[n-1]+a[n] , with initial values a[1] = 1 , a[2] = 1 . Note that to define the sequenc in Maple, you want to convert the definition to the form a[n] = a[n-2]+a[n-1] . Then all we do is enter

> a:=proc(n)

> option remember;

> a(n-2)+a(n-1);

> end;

a := proc (n) option remember; a(n-2)+a(n-1) end pr...

> a(1):=1;a(2):=1;

a(1) := 1

a(2) := 1

Let us calculate a[5] .

> a(5);

5

How about graphing the Fibonacci sequence.

> plot([seq([n,a(n)],n=1..100)],style=point);

[Maple Plot]

It looks like the Fibonacci sequence is increasing and diverges. However, because this sequence is defined recursively, Maple is unable to compute the limit.

> limit(a(n),n=infinity);

Error, (in limit) invalid limiting point

Submission:

For the sequences below, do the following.

a) Define a procedure which gives the sequence.

b) Plot the first 100 terms in the sequence.

c) Decide whether the sequence appears to converge, and if so, estimate the value to which it converges.

1) a[n+1] = a[n]/(n+1) , a[1] = 1 .

2) a[n+1] = n*a[n]/(n+1), a[1] = -2 .

Submission worksheet:

>

2. Sequences and recursion.

In a problem analysis it is often convenient to express a value at one point in time in terms of values at earlier times - this is called recursion. The following program generates a sequence, { p[n] } where p[n+1] = k*p[n](1-p[n]) which gives the size of a population at generation n which is of interest to an ecologist, p[n] is the fraction of the maximum size of the population, so it is always between 0 and 1.

> p := proc(k,p0,n)

> a := [p0]:

> for i from 2 to n do

> a:= [op(a),k*a[i-1]*(1-a[i-1])]:

> od;

> RETURN(a);

> end:

Warning, `a` is implicitly declared local to procedure `p`

Warning, `i` is implicitly declared local to procedure `p`

> popmodel := p(1.1,0.5,30);

popmodel := [.5, .275, .2193125, .1883359801, .1681...
popmodel := [.5, .275, .2193125, .1883359801, .1681...
popmodel := [.5, .275, .2193125, .1883359801, .1681...

If we want to look at a graph of this we would just type

> plot([seq([n,popmodel[n]],n=1..30)],style=point);

[Maple Plot]

Submission:

(a) Calculate 30 terms of the sequence { p[n] } for p[0] = 1/2 and for for two values of k such that 1< k < 3 . Graph the sequences. Do they appear to converge? Repeat for a different value of p[0] between 0 and 1. Does the limit depend on the choice of p[0] ? Does it depend on the choece of k ?

(b) Calculate terms of the sequence for a value of k between 3 and 3.4 and plot them. What do you notice about the behavior of the terms?

(c) Experimaent with values of k between 3.4 and 3.5. what happens to the terms?

(d) For values of k between 3.6 and 4, compute and plot at least 100 terms and comment on the behavior of the sequence. What hapens if you change p[0] by 0.001? this type of behavior is called chaotic and is exhibited by insect populations under certain conditions. Say some words about why you think this behavior is called chaotic .

Submission worksheet:

>

3. Picard's method for finding roots.

Picard's method of finding roots of the equation f(x) = 0 is to consider the equation function g(x) = f(x)+x , and solve the equation g(x) = x instead. The solutions to both equations are the same. At first, this seems to be no benefit, but Picard came up with an ingenious method for finding the roots of the second equation by means of a recursively defined sequence. Let x[0] be any number in the domain of the function, and define the sequence x[n] by x[n+1] = g(x[n]) . In many cases, this sequence will converge to some number x, and it can be shown that when it does converge to a number x , then it happens that g(x) = x , so that f(x) = 0 .

Let us consider a simple example. Suppose that f(x) = -3/4*x+3 . Then g(x) = x/4+3 . It is easy to determine that g(x) = x precisely when x = 4 . In fact, we can solve the equation very easily, because it is linear. Let us construct the Picard sequence for this function, starting with x[0] = 5 .

> g:=x->x/4+3;

g := proc (x) options operator, arrow; 1/4*x+3 end ...

> x:=proc(n)

> option remember;

> g(x(n-1));

> end;

x := proc (n) option remember; g(x(n-1)) end proc

> x(1):=5;

x(1) := 5

Let us plot the first 100 terms of this sequence to see what it seems to converge to.

> plot([seq([n,x(n)],n=1..100)],style=point);

[Maple Plot]

It does appear that the sequence is converging to x = 4 , as we would hope.

Let us solve the equation cos(x) = x . In this case, we have f(x) = cos(x)-x , or g(x) = cos(x) .

> g:=cos;

g := cos

> x:=proc(n)

> option remember;

> g(x(n-1));

> end;

x := proc (n) option remember; g(x(n-1)) end proc

Let us use x[1] = 0 .

> x(1):=0;

x(1) := 0

> plot([seq([n,x(n)],n=1..100)],style=point);

[Maple Plot]

This time the sequence jumps around a bit before approaching a value which we can estimate by clicking on the graph. x = .74 is an approximate solution. We could get a better approximation by evaluating the value of x[100] .

> evalf(x(100));

.7390851332

It is interesting to note that Maple cannot solve this problem exactly.

> x:='x';

x := 'x'

> solve(cos(x)=x);

RootOf(_Z-cos(_Z))

Nevertheless, Maple can find an estimate of the solution.

> fsolve(cos(x)=x);

.7390851332

The answer concurs with our estimate from the Picard sequence. I wonder what method Maple might have used to come up with this estimate?

Submission:

Use Picard's method to solve the following equations. First define the Picard sequence for the problem, then graph it, and use the graph to estimate the solution. Then determine the solution with Maple.

1) sqrt(x) = x .

2) cos(x) = x+1

Submission worksheet:

>

>

>