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]

Let' s store a huge list of Pell numbers in memory:

pellSeq=generateSequence[0,1,10000];

pellSeq[[1;;20]]

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

p[n_]:=pellSeq[[n+1]]  (* return the Pell number (Subscript[P, n]) *)

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

Now generate three lists:

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

Finding Coefficients

Start with two equations (Subscript[P, 3n])=a ((Subscript[P, n])^3)+b (Subscript[P, n]) for n=2 and n=4:

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

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

Check that (Subscript[P, 3n])=8 ((Subscript[P, n])^3)+3 (Subscript[P, n]) holds for lots of even n besides 2 and 4.

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

lhs==rhs

We have confirmed the identity for 150 values of n.

Your turn: Now insert (((-1))^n) in the verification and check that it works 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