Handling Sequences

This is based on Chapter 1.3.

We will illustrate rounding errors, and at the same time, look at ways to get Mathematica to compute sequences.

We will focus on the following sequence:

  p[n] = 10/3 p[n-1] - p[n-2],
  p[0] = 1
  p[1] = 1/3.
  
Let us first figure out the theoretical approach.  One guesses a solution:

  p[n] = A r^n

and we solve for r.  We end up with the equation:

Solve[(r^2)==((10 r)/3)-1,r]

which has roots r = 1/3 and r = 3.  Thus the general solution is

  p[n] = A (1/3)^n + B 3^n.

We use the conditions p[0] = 1, p[1] = 1/3 :

p[n_]=A ((((1/3)))^n)+B (3^n);
Solve[{p[0]==1,p[1]==(1/3)}]

to conclude that

  p[n] = 1/3^n.

Now, after having done this theoretically, let us have Mathematica check it out for us.

There are two ways to compute these sequences.  The first way is the recursive method.  This method seems mathematically slicker, but is actually fraught with troubles.  First, if one forgets to put := instead of =, the computer locks up.  Secondly, it is computationally very inefficient.  In fact, the computation time is exponential in n.

Clear[p]
p[n_]:=(10/3) p[n-1]-p[n-2]
p[0]=1;
p[1]=(1/3);
p[30]

The second way is the iterative method.  This is much more effective.  We will use this method from now on.

Clear[p]
p[0]=1;
p[1]=(1/3);
Do[p[n]=(10/3) p[n-1]-p[n-2],{n,2,100}]
p[30]
p[100]

This is all well and good, but most computers cannot handle the infinite precision we have here.  Let us look at what happens if we do the same computation using floating point arithmetic.  Let us compare the answer with the theoretical answer.

As you can see, the answers are very inaccurate for large enough n.  Why is this?

Clear[p]
p[0]=1.;
p[1]=N[(1/3)];
Do[p[n]=N[(10/3) p[n-1]-p[n-2]],{n,2,100}]
TableForm[Table[{n,N[(1/(3^n))],p[n]},{n,0,100}]]