The Museum of HP Calculators

Implicit Runge-Kutta methods for the HP-41

This program is by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

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)          -------------------s:   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:                           ( Registers R00 to R04 and R16 to R45 are to be initialized before executing "IRK8" )

    R00:  f name      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                       or another "small" number:
51  X<Y?
52  GTO 01
53  RCL 03
54  ST+ 01
55  RCL 10
56  STO 02
57  DSE 09
58  GTO 01
59  RCL 01
60  CLD
61  END

( 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:                          ( Registers R00 to R05 and R25 to R54 are to be initialized before executing "IRK8B" )

    R00:  f name      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.
 
 

01  LBL "IRK8B"
02  RCL 05
03  STO 12
04  15.024
05  CLRGX             if you don't have an HP-41CX , replace lines 04-05 with  CLX  STO 15  STO 16  STO 17  STO 18  STO 19
06  LBL 01                                                                                                               STO 20  STO 21  STO 22  STO 23  STO 24
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                or another "small" number
71  X<Y?
72  GTO 01
73  RCL 04
74  ST+ 01
75  RCL 14
76  STO 03
77  RCL 13
78  STO 02
79  DSE 12
80  GTO 01         a three-byte GTO
81  RCL 01
82  CLD
83  END

( 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:                                              (  Registers R00 to R05 and    Rn+11 to Rn2+3n+10     are to be initialized before executing "IRK" )

    R00:  f name      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.

     >>>>>       In other words, this subroutine has to change the stack from   Y = y
                                                                                                                    X = x        into       X' = f(x;y)

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               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
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                     or another "small" number
61  X<Y?
62  GTO 01
63  RCL 05
64  ST- 07
65  CLX
66  LBL 04
67  RCL IND 07
68  RCL IND 08
69  *
70  +
71  ISG 08
72  ISG 07
73  GTO 04
74  ST+ 02
75  RCL 03
76  ST+ 01
77  DSE 06
78  GTO 01
79  RCL 02
80  RCL 01
81  CLD
82  END

( 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:                                           ( Registers R00 to R06 and    R2n+14 to Rn2+4n+13     are to be initialized before executing "IRKB" )

    R00:  f name      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.

                                                                                                                  Z = z
      >>>>>    In other words, this subroutine has to change the stack from   Y = y       into     Y' = f(x;y;z)
                                                                                                                  X = x                  X' = g(x;y;z)

  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               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
  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                     or another "small" number
  86  X<Y?
  87  GTO 01
  88  RCL 06
  89  ST- 11
  90  ST- 12
  91  CLST
  92  LBL 04
  93  X<>Y
  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           a three-byte GTO
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:                                        (  Registers R00 to R05   and    Rn+12 to Rn2+2n+11     are to be initialized before executing "I2RK" )

    R00:  f name      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.
 

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               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
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                     or another "small" number
63  X<Y?
64  GTO 01
65  RCL 03
66  ST+ 01
67  RCL 11
68  STO 02
69  DSE 06
70  GTO 01
71  RCL 01
72  CLD
73  END

( 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:                                         ( Registers  R00 to R06  and   R2n+14 to Rn2+4n+13  are to be initialized before executing "I2RKB" )

    R00:  f name      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.
 
 

  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               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
  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                     or another "small" number
  89  X<Y?
  90  GTO 01
  91  RCL 04
  92  ST+ 01
  93  RCL 15
  94  STO 03
  95  RCL 14
  96  STO 02
  97  DSE 07
  98  GTO 01               a three-byte GTO
  99  RCL 01
100  CLD
101  END

( 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:     J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
 
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall