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

u0 = 0  ,   u1 = 1  &   un+1 = un + un-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 / un-1 X n un L / n

where n is a non-negative integer

Examples:

49  XEQ "FIB"  >>>>   u49  = 7778742049           X<>Y   u48  = 4807526976
480      R/S        >>>>   u480 = 9.216845715 E99   X<>Y   u479 = 5.696323921 E99

Notes:

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

r2 - r - 1 = 0   which has 2 real roots  r1 = ( 1 + sqrt(5) ) / 2  and  r2 = ( 1 - sqrt(5) ) / 2

-So, the solution may be written   un = c1 r1n + c2 r2n       u0 = 0  &  u1 = 1 lead to   c1 = 1/sqrt(5)  &  c2 = -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 un

where n is a non-negative integer

Examples:

49  XEQ "FIB"  >>>>   u49 = 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 un 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  (un)  defined by  u0 = 0 , u1 = 1 , u2 = 2  and  un = 2 un-1 - 3 un-2 + un-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 / un-2 Y / un-1 X n un

where n is a non-negative integer

Example:

41  XEQ "SEQ"  >>>>  u41 = -25048924
RDN   u40 =  -1541375
RDN   u39 =  +9734175

-Likewise,   149  R/S  >>>>  u149 = 1.145414665 E27   ... and so on ...

Notes:

-This program employs:    u -3 = 3 , u -2 = 1 , u -1 = 0   to get correctly  u 0  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   un = 2 un-1 - 3 un-2 + un-3  may be re-written:    un - 2 un-1 + 3 un-2 - un-3 = 0

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

r1 = 0.4301597090
r2 = 0.7849201455 + 1.307141279 i
r3 = 0.7849201455 - 1.307141279 i

-The "initial values"   u0 = 0 , u1 = 1 , u2 = 2  lead to

c1 = 0.234486766
c2 = -0.117243393 - 0.414334183 i          and  we have   un = c1 r1n + c2 r2n + c3 r3n
c3 = -0.117243393 + 0.414334183 i

-This last formula gives for instance  u41 = -25048924.04 &  u149 = 1.145414799 E27
-We can also avoid complex operations and it yields:

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

3°)  Linear Equations

-Let ( un )  defined by  un = a1 un-p + a2 un-p+1 + ................ + ap un-1 + b

and by p initial values, for instance  u1 , u2 , ......... , up or  u0 , u1 , ............ , up-1
or , more generally:  um-p+1 , ............... , um

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

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

•  R01 = um-p+1 ,  •  R02 = um-p+2 ,  ...............  ,  •  Rpp = um
•  Rp+1 = a1     ,  •  Rp+2 = a2      ,  ...............  ,  •  R2p = ap   ,   R2p+1 = b

>>>>  When the program stops,  Rpp =  un = 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 un L / n

with n > m

Example:   Calculate u10 & u14  where (un)  is defined by   u1 = 1 , u2 = -3 ,  u3 = 2 , u4 = 5 and  un = 2 un-4 - 4 un-3 + un-2 + 7 un-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"   >>>>   u10 = 748401 = R04

and  R03 = u9 = 105902 , R02 = u8 = 14985 , R01 = u7 = 2123

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

and   R03 = u13 = 264151943 , R02 = u12 = 37377820 , R01 = u11 = 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:

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

and by p initial values, for instance  u1 , u2 , ......... , up or  u0 , u1 , ............ , up-1
or , more generally:  um-p+1 , ............... , um

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

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

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

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

Flags: /
Subroutine:   A function that computes f ( un-p , .......... , un-1 ; n )  assuming  R01 = un-p , ............. , Rpp = un-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 un L / n

with n > m

Example1:         (un)  is defined by  u0 = 0 , u1 = 1 , u2 = 2  and  un = 2 un-1 - 3 un-2 + un-3  ;   Calculate  u41 & u149

-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"  >>>>   u41 = -25048924 = R03   and   R02 =  u40 =  -1541375  ,  R01 =  u39 =  +9734175

-And to get  u149

3     ENTER^           ( or  RDN  149  R/S  if the stack is unchanged )
41    ENTER^
149      R/S          >>>>    u149 = 1.145414662 E27 = R03  and   R02 =  u148 = 1.090817476 E27   ,  R01 =  u147 = 2.438982157 E26

Example2:         (un)  is defined by  u1 = 1 , u2 = 2 , u3 = 1  and  un = sqrt ( un-1 un-2 ) + 2 un-3 - Ln (n) ;   Calculate  u10 & u49

 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"   >>>>    u10 = 18.96582848  .....

RDN   49   R/S   >>>>   u49 = 1189084024 = R03

Example3:         (un)  is defined by  u0 = u1 = 0  and  un = (1/4) (n-2)2 - un-1 - (1/4) un-2      Calculate  u10

 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"  >>>>   u10 = 8.296875 = R02

-The exact solution is  un = (1/27) [ (-1/2)n ( 2n-4 ) + 3 n2 - 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