Pell Identities

First we will define a Module that generates recursive sequences.

generateSequences[n_, val0_, val1_]:=Module[{seq={val0, val1}, acc0=val0, acc1=val1,next, m=2},
If[n ≤ 1, Return[seq]];
Do[
next=m*acc1+acc0;            (* compute the next value *)
seq=Append[seq,next];     (* append to the sequence *)
acc0=acc1;     (* shift *)
acc1=next;      (* shift *)
,{i,2,n}];
seq                           (* return the sequence *)
]

Generate the Pell sequence to term (Subscript[P, 10]):

generateSequences[10,0,1]

Generate the Lucas sequence to term (Subscript[Q, 10]):

generateSequences[10,2,2]

Write some code to conveniently produce Pell numbers

Generate a long list of Pell numbers:

M=1000
pellSeq=generateSequences[M,0,1];

Make a function to easily access the Pell number (Subscript[P, n]). Note that we adjust the indexing so that (Subscript[P, 0])=0.

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

Try it out:

P /@ Range[0,8]

Looks good.

The 3n Pell identity

Produce some subsets of Pell numbers and their cubes:

P/@Range[0,15,3]     (* every third Pell number *)
(P/@Range[0,5])^3   (* cubes of Pell numbers *)
P/@Range[0,5]             (* Pell numbers *)

Observe:

0==8*0+3*0
5==8*1-3*1
70==8*8+3*2
985==8*125-3*5
13860==8*1728+3*12

It looks like (Subscript[P, 3n])=8(Subscript[P, n]^3)+3(((-1))^n)(Subscript[P, n]).

Now let’s pretend we didn’t know the coefficient values of 8 and ±3, and instead discover them from scratch.

Start with even n. We want to find A and B such that (Subscript[P, 3n])=((Subscript[AP, n])^3)+(Subscript[BP, n]). We will set up two equations to solve: