Pell Sequence

Math 242 Modern Computational Math, spring 2021

First we will define a module that generates a recursive sequence. The following module generateSequence takes three parameters and returns a list.

Specifically, generateSequence[(Subscript[a, 0]), (Subscript[a, 1]), n] returns a list of n≥ 2 terms. The first two terms are (Subscript[a, 0]) and (Subscript[a, 1]); all other terms are computed by the recurrence (Subscript[a, n])=2(Subscript[a, n-1])+(Subscript[a, n-2]).

generateSequence[a_,b_,n_]:=Module[
{seq={a, b}, val1=a, val2=b,next, m=2},
If[n ≤ 2, Return[seq]];
Do[
next=m*val2+val1;            (* compute the next value *)
seq=Append[seq,next];     (* append the next value to the sequence *)
val1=val2;                                (* shift val2 to val1 *)
val2=next;                                 (* shift next to val2 *)
,n-2];
seq      (* return the sequence *)
]

Note: This module doesn’t simply return the nth term of a sequence. It returns a list of n terms.

For example, we can generate the first ten terms of the Pell sequence:

generateSequence[0,1,20]

Or we can generate the first ten terms of the Pell-Lucas sequence:

generateSequence[2,2,20]

Compute a long list of Pell numbers and store them.

pellSeq=generateSequence[0,1,1000];

pellSeq[[1;;15]]

The first Pell number is usually denoted (Subscript[P, 0]). Let’s define an indexed variable to fix this:

p[n_]:=pellSeq[[n+1]]

p[0]

{p[0],p[1],p[2],p[3],p[4]}

Let' s find a pattern

{
p/@Range[0,18,3], (* every third Pell number *)
(p/@Range[0,6])^3, (* cube of every Pell number *)
p/@ Range[0,6]   (* every Pell numer *)
} //Transpose//TableForm

Finding Coefficients

Suppose we want to find a and b such that (Subscript[P, 3n])=((a((Subscript[P, n])))^3)+b (Subscript[P, n]). Let’s write equations for n=2 and n=4 and solve for a and b.

eq1=p[3*2]==a*(p[2])^3+b*p[2]  (* n=2 equation *)
eq2=p[3*4]==a*(p[4])^3+b*p[4]  (* n=4 equation *)

solution=Solve[{eq1,eq2},{a,b}]

Let' s now check that the identity (Subscript[P, 3n])=8((((Subscript[P, n])))^3)+3 (Subscript[P, n]) holds for lots of even values of n.

lhs=Table[p[3n],{n,0,300,2}];
rhs=Table[a(p[n])^3+b*p[n],{n,0,300,2}]/.solution[[1]];

lhs==rhs

Your turn: insert (((-1))^n) in the formula and check that it holds for odd n

Here is one solution:

lhs=Table[p[3n],{n,0,300}];
rhs=Table[a(p[n])^3+(-1)^n*b*p[n],{n,0,300}]/.solution[[1]];
lhs==rhs