# hp41programs

Implicit-Runge-Kutta Implicit Runge-Kutta methods for the HP-41

Overview

1°) An Implicit Runge-Kutta Method with order 8
2°) A General Implicit-Runge-Kutta Program  ( with an example of a 12th-order method )
3°) A Slighly less general Implicit-Runge-Kutta Program  ( with an example of a 13th-order method )

-These methods are applied to solve:

a) A first order differential equation:       dy/dx = f(x,y)                                          with the initial value:    y(x0) = y0
b) A system of 2 differential equations:   dy/dx = f(x,y,z)  ;   dz/dx = g(x,y,z)       with the initial values:   y(x0) = y0  ;   z(x0) = z0

1°) An implicit Runge-Kutta Method with order 8

a) First order differential equations:  dy/dx = f (x;y)

-The formulae are called explicit Runge-Kutta methods if the successive k1 , k2 , .... , kn  may be directly computed.
-In the following program , we have:

k1 = h.f ( x + b1h ; y + a1,1 k1 +  a1,2 k2 +  ..... + a1,5 k5 )
.....................................................................................          where ai,j # 0  for j > i-1
k5 = h.f ( x + b5h ; y + a5,1 k1 +  a5,2 k2 +  ..... + a5,5 k5 )

and   ci = a5,i

-The advantage is that only 5 stages ( instead of 11 ) are required to obtain an 8th-order method: the program is shorter and less data registers are needed.
-But we have to solve a 5x5 non-linear system within each step ! Speed is not a characteristic of implicit methods, especially with an HP-41...
-"IRK8" uses an iterative algorithm based on the "fixed-point theorem"
-Some implicit Runge-Kutta formulae are also satisfactory for stiff problems, but the elementary iterative method used by "IRK8" will seldom lead to convergence
in such a case, unless h is very small: another algorithm ( like the Newton's method ) would be better to solve the 5x5 non-linear system, if you're in the mood...

Data Registers:     •   R00:  f name             ( Registers R00 to R04 and R16 to R45 are to be initialized before executing "IRK8" )

•   R01 = x0      •   R03 = h = stepsize                         R05 to R10: temp   ;    R11 = k1 ......   R15 = k5
•   R02 = y0      •   R04 = N = number of steps       •   R16 thru  R45 = the 30 following numbers:

R16 = 1                             R26 = 29/180                         R36 = 29/180
R17 = 1/20                        R27 = -3/140                         R37 = -0.06901154103                      the decimal numbers are actually
R18 = 49/180                    R28 = 1/2                               R38 = 0.05200216599                       rational functions of sqrt(21)
R19 = 16/45                      R29 = 1/20                             R39 = -3/140                                     ( all decimals are correct )
R20 = 49/180                    R30 = 0.2813091833             R40 = 0
R21 = 1/20                        R31 = 73/360                         R41 = 1/20
R22 = 0.8273268354        R32 = -0.05283696110         R42 = -7/60
R23 = 1/20                        R33 = 3/160                           R43 = 2/15
R24 = 0.2702200562        R34 = 0.1726731646             R44 = -7/60
R25 = 0.3674242394        R35 = 1/20                             R45 = 1/20

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

 01  LBL "IRK8"    02  RCL 04 03  STO 09 04  CLX 05  STO 11 06  STO 12 07  STO 13 08  STO 14 09  STO 15 10  LBL 01 11  11.015 12  STO 05 13  45 14  STO 07 15  CLX 16  STO 08 17  LBL 02 18  15.01 19  STO 06 20  CLX 21  LBL 03 22  RCL IND 06 23  RCL IND 07 24  * 25  + 26  DSE 07 27  DSE 06 28  GTO 03 29  RCL 02 30  + 31  STO 10 32  RCL IND 07 33  RCL 03 34  * 35  RCL 01 36  + 37  XEQ IND 00 38  RCL 03 39  * 40  ENTER^ 41  X<> IND 05 42  - 43  ABS 44  ST+ 08 45  DSE 07 46  ISG 05 47  GTO 02 48  RCL 08 49  VIEW X 50   E-9      51  X

( 99 bytes / SIZE 046 )

 STACK INPUTS OUTPUTS Y / y(x0+N.h) X / x0+N.h

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(0.5)

-Initialize registers R16 thru R45.

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.1
0.1  STO 03
5    STO 04    XEQ "IRK8"   yields          x = 0.5                       ( in 5mn57 )
X<>Y                  y(0.5) = 0.7788007830     ( exact result is  0.7788007831 )

-If you replace lines 50-51 with  X#0? , you'll get the exact value 0.7788007831 but this severe test might sometimes lead to an infinite loop!

-The HP-41 displays  | k1 - k'1 | + ............ +  | k5 - k'5 |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

-Now, we have to solve a 10x10 non-linear system within each step.

Data Registers:      •   R00:  f name           ( Registers R00 to R05 and R25 to R54 are to be initialized before executing "IRK8B" )

•   R01 = x0      •   R04 = h = stepsize                         R06 to R14: temp   ;    R15 , ..... , R24 = ki
•   R02 = y0      •   R05 = N = number of steps       •   R25 thru  R54 = the same 30 numbers used in the previous program.
•   R03 = z0
Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-If you don't have an HP-41CX,
replace lines 04-05 with  CLX  STO 15  STO 16  STO 17  STO 18  STO 19  STO 20  STO 21  STO 22  STO 23  STO 24
-Line 80 is a three-byte GTO 01

 01  LBL "IRK8B" 02  RCL 05 03  STO 12 04  15.024 05  CLRGX  06  LBL 01  07  15.019 08  STO 06 09  20 10  STO 07 11  54 12  STO 10 13  CLX 14  STO 11 15  LBL 02 16  19.014 17  STO 08 18  24 19  STO 09 20  CLST 21  LBL 03 22  X<>Y 23  RCL IND 08 24  RCL IND 10 25  * 26  + 27  X<>Y 28  RCL IND 09 29  RCL IND 10 30  * 31  + 32  DSE 10 33  DSE 09 34  DSE 08 35  GTO 03 36  RCL 03 37  + 38  STO 14 39  X<>Y 40  RCL 02 41  + 42  STO 13 43  RCL IND 10 44  RCL 04 45  * 46  RCL 01 47  + 48  XEQ IND 00 49  RCL 04 50  ST* Z 51  * 52  ENTER^ 53  X<> IND 07 54  - 55  ABS 56  X<>Y 57  ENTER^ 58  X<> IND 06 59  - 60  ABS 61  + 62  ST+ 11 63  SIGN 64  ST+ 07 65  ST- 10 66  ISG 06 67  GTO 02 68  RCL 11 69  VIEW X 70   E-9  71  X

( 138 bytes / SIZE 055 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(0.5) and z(0.5)

-Initialize registers R25 thru R54.

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03   and if we choose  h = 0.1
0.1  STO 04
5    STO 05    XEQ "IRK8B"   yields          x = 0.5                          execution time = 12mn
RDN                      y(0.5) =  0.7788007830      ( exact result is  0.7788007831 )
RDN                      z(0.5) = -0.7788007830      ( exact result is -0.7788007831 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | k10 - k'10 |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

2°)  A General implicit-Runge-Kutta Program

a) First order differential equations:  dy/dx = f (x;y)

-A n-stage implicit Runge-Kutta method is defined by   n(n+2)  coefficients ai,j ; bi ; ci

k1 = h.f ( x + b1h ; y + a1,1 k1 +  a1,2 k2 +  ..... + a1,n kn )
.....................................................................................            (  actually,     ai,1 + ai,2 + .......... + ai,n = bi  )
kn = h.f ( x + bnh ; y + an,1 k1 +  an,2 k2 +  ..... + an,n kn )

then:     y(x+h) = y(x) + c1.k1+ ................ + cn.kn

-The following program will work for any implicit method.

Data Registers:         •   R00:  f name            (  Registers R00 to R05 and    Rn+11 to Rn2+3n+10     are to be initialized before executing "IRK" )

•   R01 = x0                                      R06 to R10:  temp             •  Rn+11 to Rn2+3n+10 =  the ( n2 + 2n ) coefficients ai,j ; bi ; ci
•   R02 = y0                                      R11 = k1                                                                     which are to be stored as shown below:
•   R03 = h = stepsize                        R12 = k2
•   R04 = N = number of steps          ............
•   R05 = n  = number of stages        Rn+10 = kn

•  Rn+11   = a1,1        •  Rn+12   = a1,2   ...........................    •  R2n+10 = a1,n     •  R2n+11 = b1
•  R2n+12 = a2,1        •  R2n+13 = a2,2   ...........................    •  R3n+11 = a2,n     •  R3n+12 = b2

..................................................................................................................................................

•  Rn2+n+10 = an,1     •  Rn2+n+11 = an,2   .......................   •  Rn2+2n+9 = an,n   •  Rn2+2n+10 = bn

•  Rn2+2n+11 = c1      •  Rn2+2n+12 = c2   .......................   •  Rn2+3n+10 = cn

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

-If you don't have an HP-41CX , replace line 12 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00

 01  LBL "IRK"     02  RCL 04 03  STO 06 04  RCL 05 05  10 06  + 07  .1 08  % 09  11 10  + 11  STO 07 12  CLRGX 13  LBL 01 14  RCL 05 15  STO 08 16  RCL 07 17  FRC 18  11.9 19  ST+ 08 20  INT 21  + 22  STO 07 23  CLX 24  STO 10 25  LBL 02 26  RCL 07 27  FRC 28  11 29  + 30  STO 09 31  CLX 32  LBL 03 33  RCL IND 08 34  RCL IND 09 35  * 36  + 37  ISG 08 38  ISG 09 39  GTO 03 40  RCL 02 41  + 42  RCL 03 43  RCL IND 08 44  * 45  RCL 01 46  + 47  XEQ IND 00 48  RCL 03 49  * 50  ENTER^ 51  X<> IND 07 52  - 53  ABS 54  ST+ 10 55  ISG 08 56  ISG 07 57  GTO 02 58  RCL 10 59  VIEW X 60   E-9   61  X

( 126 bytes / SIZE n2+3n+11 )

 STACK INPUTS OUTPUTS Y / y(x0+N.h) X / x0+N.h

-Implicit methods are very accurate:  n-stage 2n-order methods do exist which would be impossible with explicit methods,
but, on the other hand, we have to solve a nxn non-linear system within each step!
-For instance, here are the 48 coefficients of a 6-stage 12-order method, based on Gauss-Legendre quadrature.
-These coefficients are given with 20 decimals in case you want to use them with a calculator working with more than 10 digits:

R17 = a1,1 =  0.04283 11230 94792 58626          R24 = a2,1 =  0.09267 34914 30378 86319          R31 = a3,1 =  0.08224 79226 12843 87381
R18 = a1,2 = -0.01476 37259 97197 41248          R25 = a2,2 =  0.09019 03932 62034 65189          R32 = a3,2 =  0.19603 21623 33245 00606
R19 = a1,3 =  0.00932 50507 06477 75119          R26 = a2,3 = -0.02030 01022 93239 58595          R33 = a3,3 =  0.11697 84836 43172 76185
R20 = a1,4 = -0.00566 88580 49483 51191          R27 = a2,4 =  0.01036 31562 40246 42374          R34 = a3,4 = -0.02048 25277 45656 09763
R21 = a1,5 =  0.00285 44333 15099 33514          R28 = a2,5 = -0.00488 71929 28037 67146          R35 = a3,5 =  0.00798 99918 99662 33580
R22 = a1,6 = -0.00081 27801 71264 76211          R29 = a2,6 =  0.00135 55610 55485 06178          R36 = a3,6 = -0.00207 56257 84866 33419
R23 =  b1  =  0.03376 52428 98423 98609          R30 =  b2  =  0.16939 53067 66867 74317          R37 =   b3 =   0.38069 04069 58401 54568

R38 = a4,1 =  0.08773 78719 74451 50671          R45 = a5,1 =  0.08430 66851 34100 11074          R52 = a6,1 =  0.08647 50263 60849 93463
R39 = a4,2 =  0.17239 07946 24406 96799          R46 = a5,2 =  0.18526 79794 52106 97525          R53 = a6,2 =  0.17752 63532 08969 96865
R40 = a4,3 =  0.25443 94950 32001 62132          R47 = a5,3 =  0.22359 38110 46099 09996          R54 = a6,3 =  0.23962 58253 35829 03560
R41 = a4,4 =  0.11697 84836 43172 76185          R48 = a5,4 =  0.25425 70695 79585 10965          R55 = a6,4 =  0.22463 19165 79867 77250
R42 = a4,5 = -0.01565 13758 09175 70227         R49 = a5,5 =  0.09019 03932 62034 65189           R56 = a6,5 =  0.19514 45125 21266 71626
R43 = a4,6 =  0.00341 43235 76741 29871          R50 = a5,6 = -0.00701 12452 40793 69067          R57 = a6,6 =  0.04283 11230 94792 58626
R44 =  b4  =  0.61930 95930 41598 45432          R51 =  b5  =  0.83060 46932 33132 25683          R58 =  b6  =  0.96623 47571 01576 01391

R59 =  c1  =  0.08566 22461 89585 17252          R61 =  c3  =  0.23395 69672 86345 52369           R63 =  c5  =  0.18038 07865 24069 30378
R60 =  c2  =  0.18038 07865 24069 30378          R62 =  c4  =  0.23395 69672 86345 52369           R64 =  c6  =  0.08566 22461 89585 17252

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

-Initialize registers R17 thru R64

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.5
0.5  STO 03
2    STO 04    ( the number of steps )
6    STO 05    ( we are using a 6-stage method )

XEQ "IRK"   yields          x = 1                         ( in 5mn36s )
X<>Y                  y(1) = 0.3678794411     ( error = -10-10 )

( With  h = N = 1 the error is only  -8.4 10-9 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | kn - k'n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

-The same formulae may be generalized to a system of 2 differential equations, and we have to solve a 2nx2n non-linear system within each step.

Data Registers:      •   R00:  f name        ( Registers R00 to R06 and    R2n+14 to Rn2+4n+13     are to be initialized before executing "IRKB" )

•   R01 = x0                                      R07 to R13:  temp             •  R2n+14 to Rn2+4n+13 =  the ( n2 + 2n ) coefficients ai,j ; bi ; ci
•   R02 = y0                                      R14 to R2n+13 = ki                                                      which are to be stored as shown below:
•   R03 = z0
•   R04 = h   = stepsize
•   R05 = N  = number of steps
•   R06 = n   = number of stages

•  R2n+14 = a1,1        •  R2n+15 = a1,2   ...........................    •  R3n+13 = a1,n         •  R3n+14 = b1
•  R3n+15 = a2,1        •  R3n+16 = a2,2   ...........................    •  R4n+14 = a2,n         •  R4n+15 = b2

.....................................................................................................................................................

•  Rn2+2n+13 = an,1   •  Rn2+2n+14 = an,2   ......................   •  Rn2+3n+12 = an,n   •  Rn2+3n+13 = bn

•  Rn2+3n+14 = c1      •  Rn2+3n+15 = c2   ........................   •  Rn2+4n+13 = cn

Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-If you don't have an HP-41CX , replace line 13 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00
-Line 113 is a three-byte GTO 01.

 01  LBL "IRKB"   02  RCL 05   03  STO 07   04  RCL 06   05  ST+ X   06  13   07  +   08  .1   09  %   10  14   11  +   12  STO 09   13  CLRGX    14  LBL 01   15  RCL 09   16  FRC   17  14.9   18  STO 10   19  INT   20  +   21  STO 08   22  RCL 06   23  ST+ 10   24  ST+ 10 25  +   26  STO 09   27  CLX   28  STO 13   29  LBL 02   30  RCL 09   31  FRC   32  14   33  +   34  STO 11   35  RCL 06   36  +   37  STO 12   38  CLST   39  LBL 03   40  X<>Y   41  RCL IND 10   42  RCL IND 11   43  *   44  +   45  X<>Y   46  RCL IND 10   47  RCL IND 12   48  * 49  +   50  ISG 10   51  ISG 11   52  ISG 12   53  GTO 03   54  RCL 03   55  +   56  X<>Y   57  RCL 02   58  +   59  RCL 04   60  RCL IND 10   61  *   62  RCL 01   63  +   64  XEQ IND 00   65  RCL 04   66  ST* Z   67  *   68  ENTER^   69  X<> IND 09   70  -   71  ABS   72  X<>Y 73  ENTER^   74  X<> IND 08   75  -   76  ABS   77  +   78  ST+ 13   79  ISG 10   80  ISG 08   81  ISG 09   82  GTO 02   83  RCL 13   84  VIEW X   85   E-9      86  XY   94  RCL IND 10   95  RCL IND 11   96  * 97  +   98  X<>Y   99  RCL IND 10 100  RCL IND 12 101  * 102  + 103  ISG 10 104  ISG 11 105  ISG 12 106  GTO 04 107  ST+ 03 108  X<>Y 109  ST+ 02 110  RCL 04 111  ST+ 01 112  DSE 07 113  GTO 01 114  RCL 03 115  RCL 02 116  RCL 01 117  CLD 118  END

( 177 bytes / SIZE n2+4n+14 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

-If, for instance, you are using the 6-stage 12-order method described above, the same 48 coefficients are to be stored into registers R26 thru R73

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

-Initialize registers R26 thru R73.

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03   and if we choose  h = 0.5
0.5  STO 04
2    STO 05    ( the number of steps )
6    STO 06    ( the number of stages )

XEQ "IRKB"   yields          x  =  1                           execution time = 10mn43s
RDN                         y(1) =  0.3678794411      ( error = -10-10 )
RDN                      z(0.5) = -0.7357588823      ( error = 0 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | k2n - k'2n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

3°)  A slighly less general Implicit-Runge-Kutta Program

a) First order differential equations:  dy/dx = f (x;y)

-Implicit Runge-Kutta methods are usually more "stable" than any explicit method, but n-stage (2n-1) or (2n-2) order methods are even more stable
than n-stage 2n-order methods. ( For a detailed description of stability properties of Runge-Kutta methods, cf Butcher's work )
-That's precisely what happens with the Radau IIA methods ( n-stage (2n-1)-order methods ) and Lobatto IIIC methods ( n-stage (2n-2)-order methods )
which also have another feature in common:  ci = an,i  ( i = 1 , 2 , .... , n ) , thus, less data registers are needed.
-The 2 following programs take advantage of this fact but the 2 programs presented above would work too.

Note:   The programs listed in paragraph 1 use a 5-stage 8th-order Lobatto IIIC method.

Data Registers:     •   R00:  f name          (  Registers R00 to R05   and    Rn+12 to Rn2+2n+11     are to be initialized before executing "I2RK" )

•   R01 = x0                                      R06 to R11:  temp             •  Rn+12 to Rn2+2n+11 =  the ( n2 + n ) coefficients ai,j ; bi
•   R02 = y0                                      R12 = k1                                                                     which are to be stored as shown below:
•   R03 = h = stepsize                        R13 = k2
•   R04 = N = number of steps          ............
•   R05 = n  = number of stages        Rn+11 = kn

•  Rn+12   = a1,1        •  Rn+13   = a1,2   ...........................    •  R2n+11 = a1,n       •  R2n+12 = b1
•  R2n+13 = a2,1        •  R2n+14 = a2,2   ...........................    •  R3n+12 = a2,n       •  R3n+13 = b2

..................................................................................................................................................

•  Rn2+n+11 = an,1     •  Rn2+n+12 = an,2   ........................   •  Rn2+2n+10 = an,n  •  Rn2+2n+11 = bn

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

-If you don't have an HP-41CX , replace line 12 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00

 01  LBL "I2RK"   02  RCL 04 03  STO 06 04  RCL 05 05  11 06  + 07  .1 08  % 09  12 10  + 11  STO 07 12  CLRGX  13  LBL 01 14  RCL 05 15  STO 08 16  RCL 07 17  FRC 18  12 19  ST+ 08         20  + 21  STO 07 22  CLX 23  ST0 10 24  LBL 02 25  RCL 07 26  FRC 27  12 28  + 29  STO 09 30  CLX 31  LBL 03 32  RCL IND 08 33  RCL IND 09 34  * 35  + 36  ISG 08 37  CLX 38  ISG 09 39  GTO 03 40  RCL 02 41  + 42  STO 11 43  RCL 03 44  RCL IND 08 45  * 46  RCL 01 47  + 48  XEQ IND 00 49  RCL 03 50  * 51  ENTER^ 52  X<> IND 07 53  - 54  ABS 55  ST+ 10 56  ISG 08 57  CLX 58  ISG 07 59  GTO 02 60  RCL 10 61  VIEW X 62   E-9    63  X

( 109 bytes / SIZE n2+2n+12 )

 STACK INPUTS OUTPUTS Y / y(x0+N.h) X / x0+N.h

-For instance, here are the 56 coefficients of a 7-stage 13-order method ( a Radau IIA method )
-These coefficients are given with 20 decimals in case you want to use them with a more accurate calculator:

R19 = a1,1 =  0.03754 62649 93921 33133          R27 = a2,1 =  0.08014 75965 15618 96780          R35 = a3,1 =  0.07206 38469 41881 90211
R20 = a1,2 = -0.01403 93345 56460 40154          R28 = a2,2 =  0.08106 20639 85891 53668          R36 = a3,2 =  0.17106 83549 83886 61942
R21 = a1,3 =  0.01035 27896 00742 30094          R29 = a2,3 = -0.02123 79921 20711 03494          R37 = a3,3 =  0.10961 45640 40072 10923
R22 = a1,4 = -0.00815 83225 40275 01191          R30 = a2,4 =  0.01400 02912 38817 11898          R38 = a3,4 = -0.02461 98717 28984 05386
R23 = a1,5 =  0.00638 84138 79534 68494          R31 = a2,5 = -0.01023 41857 30090 16383          R39 = a3,5 =  0.01476 03770 43950 81707
R24 = a1,6 = -0.00460 23267 79148 65550          R32 = a2,6 =  0.00715 34651 51364 59050          R40 = a3,6 = -0.00957 52593 96791 40056
R25 = a1,7 =  0.00182 89425 61470 64370          R33 = a2,7 = -0.00281 26393 72406 72334          R41 = a3,7 =  0.00367 26783 97138 30567
R26 =  b1  =  0.02931 64271 59784 89197          R34 =  b2  =  0.14807 85996 68484 29185          R42 =   b3 =   0.33698 46902 81154 29910

R43 = a4,1 =  0.07570 51258 19824 42042          R51 = a5,1 =  0.07391 23421 63191 84654          R59 = a6,1 =  0.07470 55620 59796 23017
R44 = a4,2 =  0.15409 01551 42171 14465          R52 = a5,2 =  0.16135 56076 15942 43219          R60 = a6,2 =  0.15830 72238 72468 70066
R45 = a4,3 =  0.22710 77366 73202 38641          R53 = a5,3 =  0.20686 72415 52104 19782          R61 = a6,3 =  0.21415 34232 67200 03111
R46 = a4,4 =  0.11747 81870 37024 78199          R54 = a5,4 =  0.23700 71153 42694 23476          R62 = a6,4 =  0.21987 78470 31860 03999
R47 = a4,5 = -0.02381 08271 53044 17358         R55 = a5,5 =  0.10308 67935 33813 44662           R63 = a6,5 =  0.19875 21216 80635 26980
R48 = a4,6 =  0.01270 99855 33661 20563          R56 = a5,6 = -0.01885 41391 52580 44884          R64 = a6,6 =  0.06926 55016 05509 13323
R49 = a4,7 = -0.00460 88442 81289 63344         R57 = a5,7 =  0.00585 89009 74888 79182           R65 = a6,7 = -0.00811 60081 97728 29011
R50 =  b4  =  0.55867 15187 71550 13208          R58 =  b5  =  0.76923 38620 30054 50092           R66 =  b6  =  0.92694 56713 19741 11485

R67 = a7,1 =  0.07449 42355 56010 31793
R68 = a7,2 =  0.15910 21157 33650 74087
R69 = a7,3 =  0.21235 18895 02977 80420
R70 = a7,4 =  0.22355 49145 07283 23475                         (   and    ci = a7,i  )
R71 = a7,5 =  0.19047 49368 22115 57690
R72 = a7,6 =  0.11961 37446 12656 20289
R73 = a7,7 =  0.02040 81632 65306 12245 = 1/49
R74 =  b7  =  1

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

-Initialize registers R19 thru R74

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.5
0.5  STO 03
2    STO 04    ( the number of steps )
7    STO 05    ( we are using a 7-stage method )

XEQ "I2RK"   yields          x = 1                         ( in 7mn13s )
X<>Y                      y(1) = 0.3678794410     ( error = -2 10-10 )

( With  h = N = 1 the error is only  -1.5 10-9 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | kn - k'n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:      •   R00:  f name          ( Registers  R00 to R06  and   R2n+14 to Rn2+4n+13  are to be initialized before executing "I2RKB" )

•   R01 = x0                                      R07 to R15:  temp             •  R2n+16 to Rn2+3n+15 =  the ( n2 + n ) coefficients ai,j ; bi
•   R02 = y0                                      R16 to R2n+15 = ki                                                      which are to be stored as shown below:
•   R03 = z0
•   R04 = h   = stepsize
•   R05 = N  = number of steps
•   R06 = n   = number of stages

•  R2n+16 = a1,1        •  R2n+17 = a1,2   ...........................    •  R3n+15 = a1,n         •  R3n+16 = b1
•  R3n+17 = a2,1        •  R3n+18 = a2,2   ...........................    •  R4n+16 = a2,n         •  R4n+17 = b2

.....................................................................................................................................................

•  Rn2+2n+15 = an,1   •  Rn2+2n+16 = an,2   ......................   •  Rn2+3n+14 = an,n   •  Rn2+3n+15 = bn

Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-If you don't have an HP-41CX , replace line 13 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00
-Line 98 is a three-byte GTO 01

 01  LBL "I2RKB"   02  RCL 05   03  STO 07   04  RCL 06   05  ST+ X   06  15   07  +   08  .1   09  %   10  16   11  +   12  STO 09   13  CLRGX     14  LBL 01   15  RCL 09   16  FRC   17  16   18  STO 10   19  +   20  STO 08   21  RCL 06 22  ST+ 10   23  ST+ 10   24  +   25  STO 09   26  CLX   27  STO 13   28  LBL 02   29  RCL 09   30  FRC   31  16   32  +   33  STO 11   34  RCL 06   35  +   36  STO 12   37  CLST   38  LBL 03   39  X<>Y   40  RCL IND 10   41  RCL IND 11   42  * 43  +   44  X<>Y   45  RCL IND 10   46  RCL IND 12   47  *   48  +   49  ISG 10   50  CLX   51  ISG 11   52  ISG 12   53  GTO 03   54  RCL 03   55  +   56  STO 15   57  X<>Y   58  RCL 02   59  +   60  STO 14   61  RCL 04   62  RCL IND 10   63  * 64  RCL 01   65  +   66  XEQ IND 00   67  RCL 04   68  ST* Z   69  *   70  ENTER^   71  X<> IND 09   72  -   73  ABS   74  X<>Y   75  ENTER^   76  X<> IND 08   77  -   78  ABS   79  +   80  ST+ 13   81  ISG 10   82  CLX   83  ISG 08   84  ISG 09 85  GTO 02   86  RCL 13   87  VIEW X          88   E-9       89  X

( 146 bytes / SIZE  n2+3n+16 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

-If, for instance, you are using the 7-stage 13-order method described above, the same 56 coefficients are to be stored into registers R30 thru R85.

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

-Initialize registers R30 thru R85.

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03   and if we choose  h = 0.5
0.5  STO 04
2    STO 05    ( the number of steps )
7    STO 06    ( the number of stages )

XEQ "I2RKB"   yields        x  =  1                           execution time = 13mn54s
RDN                         y(1) =  0.3678794411      ( error = -10-10 )
RDN                      z(0.5) = -0.7357588822      ( error =  10-10 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | k2n - k'2n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

Reference:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5