Difference

Difference Equations for the HP-41

Overview

1°) Fibonacci Numbers
2°) A Similar Sequence
3°) Linear Equations
4°) General Case

1°) Fibonacci Numbers

-This famous sequence is defined by

u₀ = 0 , u₁ = 1 & u_n+1 = u_n + u_n-1

Data Registers: /
Flags: /
Subroutines: /

01 LBL "FIB"
02 ENTER^
03 SIGN
04 ENTER^
05 CHS
06 X<>Y
07 ST+ Z
08 LBL 01
09 ST+ Y
10 X<>Y
11 DSE Z
12 GTO 01
13 END

( 25 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	/	u_n-1
X	n	u_n
L	/	n

where n is a non-negative integer

Examples:

49 XEQ "FIB" >>>> u₄₉ = 7778742049 X<>Y u₄₈ = 4807526976
480 R/S >>>> u₄₈₀ = 9.216845715 E99 X<>Y u₄₇₉ = 5.696323921 E99

Notes:

-This routine uses a dummy u_-1 = 1
-The difference equation u_n+1 = u_n + u_n-1 which may be re-written u_n+1 - u_n - u_n-1 = 0 has a characteristic equation

r² - r - 1 = 0 which has 2 real roots r₁ = ( 1 + sqrt(5) ) / 2 and r₂ = ( 1 - sqrt(5) ) / 2

-So, the solution may be written u_n = c₁ r₁ⁿ + c₂ r₂ⁿ u₀ = 0 & u₁ = 1 lead to c₁ = 1/sqrt(5) & c₂ = -1/sqrt(5)
-And we could write the program

Data Registers: R00 = n , R01 = sqrt(5)
Flags: /
Subroutines: /

01 LBL "FIB"
02 STO 00
03 5
04 SQRT
05 STO 01
06 1
07 +
08 2
09 /
10 X<>Y
11 Y^X
12 1
13 RCL 01
14 -
15 2
16 /
17 RCL 00
18 Y^X
19 -
20 RCL 01
21 /
22 END

( 30 bytes / SIZE 002 )

STACK	INPUTS	OUTPUTS
X	n	u_n

where n is a non-negative integer

Examples:

49 XEQ "FIB" >>>> u₄₉ = 7778742108
480 R/S >>>> OUT OF RANGE

Remarks:

-Though the second version is faster, there are also larger roundoff-errors.
-Moreover, an OUT OF RANGE occurs before u_n overflows !
-Using M-code routines that don't check for overflow would give the result 9.216846403 E99 which is also much less accurate than the first version !

2°) A Similar Sequence

-If the formula that defines a sequence is not too much complicated, we can write SIZE 000 programs.

-For instance, let (u_n) defined by u₀ = 0 , u₁ = 1 , u₂ = 2 and u_n = 2 u_n-1 - 3 u_n-2 + u_n-3

Data Registers: /
Flags: /
Subroutines: /

01 LBL "SEQ"
02 SIGN
03 ST+ L
04 3
05 X<>Y
06 0
07 LBL 01
08 ENTER^
09 ST+ X
10 ST+ T
11 CLX
12 RCL Z
13 ST- T
14 ST+ X
15 ST- T
16 X<> T
17 DSE L
18 GTO 01
19 END

( 37 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Z	/	u_n-2
Y	/	u_n-1
X	n	u_n

where n is a non-negative integer

Example:

41 XEQ "SEQ" >>>> u₄₁ = -25048924
RDN u₄₀ = -1541375
RDN u₃₉ = +9734175

-Likewise, 149 R/S >>>> u₁₄₉ = 1.145414665 E27 ... and so on ...

Notes:

-This program employs: u_-3 = 3 , u_-2 = 1 , u_-1 = 0 to get correctly u₀ after 1 loop, ...
-If you want to save n in L-register, replace line 17 by DSE M , line 03 by ST+ M and add STO M after line 01.

-Here again, the difference equation u_n = 2 u_n-1 - 3 u_n-2 + u_n-3 may be re-written: u_n - 2 u_n-1 + 3 u_n-2 - u_n-3 = 0

-It leads to the characteristic equation: r³ - 2 r² + 3 r - 1 = 0 which has 1 real root and 2 complex conjugate roots

     r₁ = 0.4301597090
     r₂ = 0.7849201455 + 1.307141279 i
     r₃ = 0.7849201455 - 1.307141279 i

-The "initial values" u₀ = 0 , u₁ = 1 , u₂ = 2 lead to

     c₁ = 0.234486766
     c₂ = -0.117243393 - 0.414334183 i          and we have   u_n = c₁ r₁ⁿ + c₂ r₂ⁿ + c₃ r₃ⁿ
     c₃ = -0.117243393 + 0.414334183 i

-This last formula gives for instance u₄₁ = -25048924.04 & u₁₄₉ = 1.145414799 E27
-We can also avoid complex operations and it yields:

u_n = 0.234486766 ( 0.430159709 )ⁿ + 0.861205727 ( 1.524702580 )ⁿ cos ( 59°01576958 n - 105°7998242 )

3°) Linear Equations

-Let ( u_n ) defined by u_n = a₁ u_n-p + a₂ u_n-p+1 + ................ + a_p u_n-1 + b

and by p initial values, for instance u₁ , u₂ , ......... , u_p or u₀ , u₁ , ............ , u_p-1
or , more generally: u_m-p+1 , ............... , u_m

-We want to compute u_n for a given n > m

Data Registers: • R00 = p ( Registers R00 thru R2p+1 are to be initialized before executing "LDE" )

• R01 = u_m-p+1 , • R02 = u_m-p+2 , ............... , • Rpp = u_m
• Rp+1 = a₁ , • Rp+2 = a₂ , ............... , • R2p = a_p , R2p+1 = b

>>>> When the program stops, Rpp = u_n = X-register

Flags: /
Subroutines: /

01 LBL "LDE"
02 X<=Y?
03 SF 41
04 STO M
05 X<>Y
06 -
07 STO N

08 LBL 01
09 RCL 00
10 ST+ 00
11 ISG 00
12 CLX
13 RCL IND 00
14 ENTER^

15 DSE 00
16 LBL 02
17 X<> IND Z
18 RCL IND 00
19 X<>Y
20 *
21 ST+ Y

22 X<> L
23 DSE 00
24 DSE Z
25 GTO 02
26 X<>Y
27 STO IND 00
28 DSE N

29 GTO 01
30 0
31 X<> M
32 STO L
33 X<> Y
34 END

( 61 bytes / SIZE ppp+2 )

STACK	INPUTS	OUTPUTS
Y	m	n
X	n	u_n
L	/	n

with n > m

Example: Calculate u₁₀ & u₁₄ where (u_n) is defined by u₁ = 1 , u₂ = -3 , u₃ = 2 , u₄ = 5 and u_n = 2 u_n-4 - 4 u_n-3 + u_n-2 + 7 u_n-1 - 6

p = 4 STO 00

   1   STO 01          2 STO 05          -6   STO 09
-3 STO 02         -4 STO 06
   2 STO 03           1 STO 07
   5 STO 04           7 STO 08

m = 4 ENTER^
n = 10 XEQ "LDE" >>>> u₁₀ = 748401 = R04

and R03 = u₉ = 105902 , R02 = u₈ = 14985 , R01 = u₇ = 2123

10 ENTER^ ( or RDN if the stack is unchanged )
14 R/S >>>> u₁₄ = 1866782181 = R04

and R03 = u₁₃ = 264151943 , R02 = u₁₂ = 37377820 , R01 = u₁₁ = 5289009

Notes:

-Registers Rp+1 to R2p+1 are undisturbed
-Register R00 is temporarily modified, but its original content = p is restored at the end.

-Lines 02-03 only produce an error message if n <= m
-Otherwise, they could be deleted.
-Line 32 is almost superfluous.

4°) General Case

-Now we assume that a sequence is defined by a recurrence relation:

u_n = f ( u_n-p , u_n-p+1 , .......... , u_n-1 ; n ) where p is a positive integer

and by p initial values, for instance u₁ , u₂ , ......... , u_p or u₀ , u₁ , ............ , u_p-1
or , more generally: u_m-p+1 , ............... , u_m

-"DFE" will calculate u_n for a given n > m

Data Registers: • R00 = function name ( Registers R00 thru Rpp are to be initialized before executing "DFE" )

• R01 = u_m-p+1 , • R02 = u_m-p+2 , ............... , • Rpp = u_m ; Rp+1 = 1+m

>>>> When the program stops, Rpp = u_n = X-register

Flags: /
Subroutine: A function that computes f ( u_n-p , .......... , u_n-1 ; n ) assuming R01 = u_n-p , ............. , Rpp = u_n-1 ; Rp+1 = n

01 LBL "DFE"
02 X<=Y?
03 SF 41
04 STO M
05 X<>Y
06 -
07 STO N

08 X<>Y
09 STO O
10 ISG X
11 CLX
12 LASTX
13 STO IND Y
14 LBL 01

15 RCL O
16 1
17 +
18 ISG IND X
19 CLX
20 XEQ IND 00
21 RCL O

22 X<>Y
23 LBL 02
24 X<> IND Y
25 DSE Y
26 GTO 02
27 DSE N
28 GTO 01

29 RCL O
30 RCL M
31 STO L
32 RCL IND Y
33 CLA
34 END

( 61 bytes / SIZE ppp+2 )

STACK	INPUTS	OUTPUTS
Z	p	p
Y	m	n
X	n	u_n
L	/	n

with n > m

Example1: (u_n) is defined by u₀ = 0 , u₁ = 1 , u₂ = 2 and u_n = 2 u_n-1 - 3 u_n-2 + u_n-3 ; Calculate u₄₁ & u₁₄₉

-Load the following routine, any global LBL ( maximum 6 characters ):

01 LBL "T"
02 RCL 01
03 RCL 02
04 3
05 *
06 -
07 RCL 03
08 ST+ X
09 +
10 END

T ASTO 00 0 STO 01 1 STO 02 2 STO 03

p = 3    ENTER^
m = 2    ENTER^
n = 41 XEQ "DFE" >>>>   u₄₁ = -25048924 = R03   and   R02 = u₄₀ = -1541375 , R01 = u₃₉ = +9734175

-And to get u₁₄₉

   3     ENTER^           ( or RDN 149 R/S if the stack is unchanged )
41    ENTER^
149      R/S          >>>>    u₁₄₉ = 1.145414662 E27 = R03 and   R02 = u₁₄₈ = 1.090817476 E27   , R01 = u₁₄₇ = 2.438982157 E26

Example2: (u_n) is defined by u₁ = 1 , u₂ = 2 , u₃ = 1 and u_n = sqrt ( u_n-1 u_n-2 ) + 2 u_n-3 - Ln (n) ; Calculate u₁₀ & u₄₉

01 LBL "W"
02 RCL 02
03 RCL 03
04 *
05 SQRT
06 RCL 01
07 ST+ X
08 +
09 RCL 04
10 LN
11 -
12 END

"W" ASTO 00 1 STO 01 STO 03 2 STO 02

3 ENTER^ ENTER^
10 XEQ "DFE" >>>> u₁₀ = 18.96582848 .....

RDN 49 R/S >>>> u₄₉ = 1189084024 = R03

Example3: (u_n) is defined by u₀ = u₁ = 0 and u_n = (1/4) (n-2)² - u_n-1 - (1/4) u_n-2 Calculate u₁₀

01 LBL "U"
02 RCL 03
03 2
04 -
05 X^2
06 RCL 01
07 -
08 4
09 /
10 RCL 02
11 -
12 END

"U" ASTO 00 CLX STO 01 STO 02

   2   ENTER^
   1   ENTER^
10 XEQ "DFE" >>>>   u₁₀ = 8.296875 = R02

-The exact solution is u_n = (1/27) [ (-1/2)ⁿ ( 2n-4 ) + 3 n² - 8 n + 4 ]

Notes:

-Lines 02-03 only produce an error message if n <= m
-Otherwise, they could be deleted.
-Line 31 is almost superfluous.

Reference:

[1] F. Scheid - "Numerical Analysis" - Mc Graw Hill ISBN 07-055197-9

hp41programs

Difference Equations for the HP-41