The Museum of HP Calculators

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
 

-The purpose of the following programs is to solve differential equations with various Runge-Kutta methods:

    1°) RK4  is the classical 4th order Runge-Kutta formula
    2°) RKF is a Runge-Kutta-Fehlberg method of order 4 ( embedded within 5th order )
    3°) RK6 uses a 6th order Runge-Kutta -----
    4°) RK8 is an 8th order method
    5°) ERK is a general Runge-Kutta program suitable to all explicit formulae
    6°) IRK8 uses an implicit Runge-Kutta method of order 8.

-All 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)          -------------------    y(x0) = y0  ;   z(x0) = z0
 

1°) The classical 4th order Runge-Kutta method
 

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

RK4 uses the following formulae:       k1 = h.f(x;y)
                                                         k2 = h.f(x+h/2;y+k1/2)
                                                         k3 = h.f(x+h/2;y+k2/2)
                                                         k4 = h.f(x+h;y+k3)

  and then,   y(x+h) = y(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6      This formula duplicates the Taylor serie through the term in h4
 

Data Registers:     Registers R00 to R04 are to be initialized before executing "RK4"

    R00:  f name      R01 = x0         R03 = h/2 = half of the stepsize          R05 & R06: temp            ( storing h/2 instead of h saves several bytes )
                               R02 = y0         R04 = N = number of steps

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 "RK4"
02  RCL 04
03  STO 06
04  LBL 01
05  RCL 02
06  RCL 01
07  XEQ IND 00
08  RCL 03
09  ST+ 01
10  *
11  STO 05
12  RCL 02
13  +
14  RCL 01
15  XEQ IND 00
16  RCL 03
17  *
18  ST+ 05
19  ST+ 05
20  RCL 02
21  +
22  RCL 01
23  XEQ IND 00
24  RCL 03
25  ST+ 01
26  ST+ X
27  *
28  ST+ 05
29  RCL 02
30  +
31  RCL 01
32  XEQ IND 00
33  RCL 03
34  *
35  RCL 05
36  +
37  3
38  /
39  ST+ 02
40  DSE 06
41  GTO 01
42  RCL 02
43  RCL 01
44  END

( 65 bytes / SIZE 007 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)              ( The solution is y = exp ( -x2 ) )
 

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.05  STO 03
10     STO 04    XEQ "RK4"   yields          x = 1                                 ( in 20s )
                                X<>Y                    y(1) = 0.367881                   ( the exact result is 1/e =  0.367879441 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 24 = 16 )
 

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

-The preceding formulae are now generalized to a system of 2 differential equations ( see also "systems of differential equations for the HP-41" ):

                         k1 = h.f(x;y;z)                                             l1 = h.g(x;y;z)
                         k2 = h.f(x+h/2;y+k1/2;z+l1/2)                      l2 = h.g(x+h/2;y+k1/2;z+l1/2)
                         k3 = h.f(x+h/2;y+k2/2;z+l2/2)                      l3 = h.g(x+h/2;y+k2/2;z+l2/2)
                         k4 = h.f(x+h;y+k3;z+l3)                               l4 = h.g(x+h;y+k3;z+l3)

  then,   y(x+h) = y(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6   and   z(x+h) = z(x) + ( l1 + 2.l2 + 2.l3 + l4 ) / 6     duplicate the Taylor series through the terms in h4
 

Data Registers:     Registers R00 to R05  are to be initialized before executing RK4B

    R00:  f name      R01 = x0         R04 = h/2 = half of the stepsize              R06 to R08: temp            ( storing h/2 instead of h saves several bytes )
                               R02 = y0         R05 = N = number of steps
                               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.

                                                                                                                  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 "RK4B"
02  RCL 05
03  STO 08
04  LBL 01
05  RCL 03
06  RCL 02
07  RCL 01
08  XEQ IND 00
09  RCL 04
10  ST+ 01
11  ST* Z
12  *
13  STO 07
14  RCL 03
15  +
16  X<>Y
17  STO 06
18  RCL 02
19  +
20  RCL 01
21  XEQ IND 00
22  RCL 04
23  ST* Z
24  *
25  ST+ 07
26  ST+ 07
27  RCL 03
28  +
29  X<>Y
30  ST+ 06
31  ST+ 06
32  RCL 02
33  +
34  RCL 01
35  XEQ IND 00
36  RCL 04
37  ST+ 01
38  ST+ X
39  ST* Z
40  *
41  ST+ 07
42  RCL 03
43  +
44  X<>Y
45  ST+ 06
46  RCL 02
47  +
48  RCL 01
49  XEQ IND 00
50  RCL 04
51  ST* Z
52  *
53  RCL 07
54  +
55  3
56  /
57  ST+ 03
58  X<>Y
59  RCL 06
60  +
61  3
62  /
63  ST+ 02
64  DSE 08
65  GTO 01
66  RCL 03
67  RCL 02
68  RCL 01
69  END

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

Example:     y(0) = 1 ; y'(0) = 0 ;  y'' + 2.x.y' + y = 0   Evaluate y(1) and y'(1)            (  The solution is y = exp ( -x2 ) ;  y' = -2.x exp ( -x2 )  )

-This second order equation is equivalent to the system:             y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    where z = y'
 

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.05  STO 04
10     STO 05    XEQ "RK4B"   yields       x =  1                                        ( in 32s )
                              RDN                      y(1) =  0.367881                            the exact result is  y(1) = 1/e =  0.367879441
                              RDN                      z(1) = -0.735762                            -----------------  z(1) = -2/e = -0.735758883
 

-Press R/S to continue with the same h and N.
 

N.B:  -To obtain an error estimate, start over again with a smaller stepsize:
           when h is divided by 2 , errors are approximately divided by 16

     -In the following program, a fourth order Runge-Kutta method is imbedded within a fifth order Runge-Kutta method:
 the difference of the results provides an error estimate for the 4th order method.
 
 

2°) A 4th order Runge-Kutta-Fehlberg method
 

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

Data Registers:     Registers R00 to R05 are to be initialized before executing "RKF"

    R00:  f name      R01 = x0         R03 = h = stepsize                         R05 = the error estimate  ( clear this register before the first execution )
                               R02 = y0         R04 = N = number of steps             R06 thru R11: temp

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

001  LBL "RKF"
002  RCL 04
003  STO 06
004  LBL 01
005  RCL 02
006  RCL 01
007  XEQ IND 00
008  RCL 03
009  *
010  STO 07
011  ST+ X
012  9
013  /
014  RCL 02
015  +
016  RCL 03
017  ST+ X
018  9
019  /
020  RCL 01
021  +
022  XEQ IND 00
023  RCL 03
024  *
025  4
026  /
027  STO 08
028  RCL 07
029  12
030  /
031  +
032  RCL 02
033  +
034  RCL 03
035  3
036  /
037  RCL 01
038  +
039  XEQ IND 00
040  RCL 03
041  *
042  STO 09
043  270
044  *
045  RCL 08
046  972
047  *
048  -
049  RCL 07
050  69
051  *
052  +
053  128
054  /
055  RCL 02
056  +
057  RCL 03
058  .75
059  *
060  RCL 01
061  +
062  XEQ IND 00
063  RCL 03
064  ST+ 01
065  *
066  64
067  *
068  STO 10
069  RCL 07
070  85
071  *
072  -
073  60
074  /
075  RCL 08
076  RCL 09
077  5
078  /
079  -
080  27
081  *
082  +
083  RCL 02
084  +
085  RCL 01
086  XEQ IND 00
087  RCL 03
088  *
089  15
090  *
091  RCL 10
092  +
093  STO 11
094  RCL 09
095  351
096  *
097  +
098  RCL 08
099  540
100  *
101  -
102  RCL 07
103  65
104  *
105  +
106  432
107  /
108  RCL 02
109  +
110  RCL 01
111  RCL 03
112  6
113  /
114  -
115  XEQ IND 00
116  RCL 03
117  *
118  72
119  *
120  RCL 11
121  -
122  RCL 09
123  9
124  *
125  +
126  RCL 07
127  ST+ X
128  -
129  3
130  1/X
131  %
132  ST- 05               if you replace this line by   ABS   ST+ 05   you'll get a less optimistic error estimate ( in absolute value )
133  RCL 11
134  RCL 09
135  81
136  *
137  +
138  PI
139  R-D
140  /
141  RCL 07
142  9
143  /
144  +
145  ST+ 02
146  DSE 06
147  GTO 01            a three-byte GTO
148  RCL 02
149  RCL 01
150  END
 

( 197 bytes / SIZE 012 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 

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

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

 >>> Don't forget to clear register R05:       CLX  STO 05

"T" ASTO 00
 0    STO 01
 1    STO 02    and if we choose  h = 0.1
0.1  STO 03
10   STO 04    XEQ "RKF"   yields          x = 1
                              X<>Y                    y(1) = 0.367879263
and the error estimate    RCL 05      >>>       -9.7 10-8           ( or  5.4 10-7  with   ABS  ST+ 05  at line 132  )
-The error is actually  -1.8 10-7
 

-Press R/S to continue with the same h and N.
 
 

      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 R07  are to be initialized before executing RKFB  ( clear registers R06 & R07 before the first execution )

    R00:  f name      R01 = x0         R04 = h = stepsize                         R06 = y-error estimate           R08 to R18: temp
                               R02 = y0         R05 = N = number of steps           R07 = z-error estimate
                               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.
 
 

001  LBL "RKFB"
002  RCL 05
003  STO 18
004  GTO 01
005  LBL 00
006  XEQ IND 00
007  RCL 04
008  ST* Z
009  *
010  RTN
011  LBL 01
012  RCL 03
013  RCL 02
014  RCL 01
015  XEQ 00
016  STO 13
017  ST+ X
018  9
019  /
020  RCL 03
021  +
022  X<>Y
023  STO 08
024  ST+ X
025  9
026  /
027  RCL 02
028  +
029  RCL 04
030  ST+ X
031  9
032  /
033  RCL 01
034  +
035  XEQ 00
036  STO 14
037  4
038  /
039  RCL 13
040  12
041  /
042  +
043  RCL 03
044  +
045  X<>Y
046  STO 09
047  4
048  /
049  RCL 08
050  12
051  /
052  +
053  RCL 02
054  +
055  RCL 04
056  3
057  /
058  RCL 01
059  +
060  XEQ 00
061  STO 15
062  270
063  *
064  RCL 14
065  243
066  *
067  -
068  RCL 13
069  69
070  *
071  +
072  128
073  /
074  RCL 03
075  +
076  X<>Y
077  STO 10
078  270
079  *
080  RCL 09
081  243
082  *
083  -
084  RCL 08
085  69
086  *
087  +
088  128
089  /
090  RCL 02
091  +
092  RCL 04
093  .75
094  *
095  RCL 01
096  +
097  XEQ 00
098  64
099  ST* Z
100  *
101  STO 16
102  RCL 15
103  324
104  *
105  -
106  RCL 14
107  405
108  *
109  +
110  RCL 13
111  85
112  *
113  -
114  60
115  /
116  RCL 03
117  +
118  X<>Y
119  STO 11
120  RCL 10
121  324
122  *
123  -
124  RCL 09
125  405
126  *
127  +
128  RCL 08
129  85
130  *
131  -
132  60
133  /
134  RCL 02
135  +
136  RCL 04
137  RCL 01
138  +
139  STO 01
140  XEQ 00
141  15
142  ST* Z
143  *
144  RCL 16
145  +
146  STO 17
147  RCL 15
148  351
149  *
150  +
151  RCL 14
152  135
153  ST* 09
154  *
155  -
156  RCL 13
157  65
158  *
159  +
160  432
161  /
162  RCL 03
163  +
164  X<>Y
165  RCL 11
166  +
167  STO 12
168  RCL 10
169  351
170  *
171  +
172  RCL 09
173  -
174  RCL 08
175  65
176  *
177  +
178  432
179  /
180  RCL 02
181  +
182  RCL 04
183  6
184  /
185  RCL 01
186  X<>Y
187  -
188  XEQ 00
189  72
190  ST* Z
191  *
192  RCL 17
193  -
194  RCL 15
195  9
196  *
197  +
198  RCL 13
199  ST+ X
200  -
201  3
202  1/X
203  %
204  ST- 07               if you replace this line by   ABS   ST+ 07   you'll get a less optimistic error estimate ( in absolute value )
205  R^
206  RCL 12
207  -
208  RCL 10
209  9
210  *
211  +
212  RCL 08
213 ST+ X
214  -
215  3
216  1/X
217  %
218  ST- 06               if you replace this line by   ABS   ST+ 06   you'll get a less optimistic error estimate ( in absolute value )
219  RCL 17
220  RCL 15
221  81
222  ST* 10
223  *
224  +
225  PI
226  R-D
227  /
228  RCL 13
229  9
230  ST/ 08
231  /
232  +
233  ST+ 03
234  RCL 12
235  RCL 10
236  +
237  PI
238  R-D
239  /
240  RCL 08
241  +
242  ST+ 02
243  DSE 18
244  GTO 01            a three-byte GTO
245  RCL 03
246  RCL 02
247  RCL 01
248  END
 

( 342 bytes / SIZE 019)
 
 
 
      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(1) and z(1)
 

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

 >>> Don't forget to clear registers R06 & R07:       CLX  STO 06  STO 07

"U" ASTO 00
 0    STO 01
 1    STO 02
 0    STO 03   and if we choose  h = 0.1
0.1  STO 04
10   STO 05    XEQ "RKFB"   yields       x  =  1                      ( in 2mn16s )
                              RDN                      y(1) =  0.367879517
                              RDN                      z(1) = -0.735759034

and the error estimates    RCL 06     >>>    -8.7 10-8    ( or  6.5 10-7  with   ABS  ST+ 05 at line 218 )      actually   y-error =    7.6 10-8
                                      RCL 07     >>>    -2.1 10-7    ( or    8  10-7  ----    ABS  ST+ 06  at line 204 )    actually    z-error =  -1.5 10-7
 

-Press R/S to continue with the same h and N.
 

3°) A 6th order method
 

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

-The formula duplicates the Taylor serie through the term in h6
-Methods of this order require 7 stages ( 7 evaluations of the function f within each step)
 

Data Registers:     Registers R00 to R04 are to be initialized before executing "RK6"

    R00:  f name      R01 = x0         R03 =  h = the stepsize                   R05: counter
                               R02 = y0         R04 = N = the number of steps      R06 , ...... , R11 = k1 , ..... , k6

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

001  LBL "RK6"
002  RCL 04
003  STO 05
004  LBL 01
005  RCL 02
006  RCL 01
007  XEQ IND 00
008  RCL 03
009  *
010  STO 06
011  3
012  /
013  RCL 02
014  +
015  RCL 03
016  3
017  /
018  RCL 01
019  +
020  XEQ IND 00
021  RCL 03
022  *
023  STO 07
024  1.5
025  /
026  RCL 02
027  +
028  RCL 03
029  1.5
030  /
031  RCL 01
032  +
033  XEQ IND 00
034  RCL 03
035  *
036  STO 08
037  CHS
038  RCL 07
039  4
040  *
041  +
042  RCL 06
043  +
044  12
045  /
046  RCL 02
047  +
048  RCL 03
049  3
050  /
051  RCL 01
052  +
053  XEQ IND 00
054  RCL 03
055  *
056  STO 09
057  90
058  *
059  RCL 08
060  35
061  *
062  +
063  RCL 07
064  110
065  *
066  -
067  RCL 06
068  25
069  *
070  +
071  48
072  /
073  RCL 02
074  +
075  RCL 03
076  1.2
077  /
078  RCL 01
079  +
080  XEQ IND 00
081  RCL 03
082  *
083  STO 10
084  10
085  /
086  RCL 09
087  2
088  /
089  +
090  RCL 08
091  8
092  /
093  -
094  RCL 07
095  55
096  *
097  RCL 06
098  18
099  *
100  -
101  120
102  /
103  -
104  RCL 02
105  +
106  RCL 03
107  6
108  /
109  RCL 01
110  +
111  XEQ IND 00
112  RCL 03
113  ST+ 01
114  *
115  STO 11
116  80
117  *
118  RCL 09
119  118
120  *
121  -
122  RCL 07
123  99
124  *
125  +
126  39
127  /
128  RCL 10
129  128
130  *
131  RCL 08
132  215
133  *
134  +
135  RCL 06
136  783
137  *
138  -
139  780
140  /
141  +
142  RCL 02
143  +
144  RCL 01
145  XEQ IND 00
146  RCL 03
147  *
148  RCL 06
149  +
150  13
151  *
152  RCL 11
153  RCL 10
154  +
155  32
156  *
157  +
158  RCL 09
159  RCL 08
160  +
161  55
162  *
163  +
164  .5
165  %
166  ST+ 02
167  DSE 05
168  GTO 01              a three-byte GTO
169  RCL 02
170  RCL 01
171  END
 

( 226 bytes / SIZE 012 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 

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

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
10     STO 04    XEQ "RK6"   yields          x = 1                       ( in 89s )
                                X<>Y                    y(1) = 0.367879436    ( error =  -5 10-9 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 26 = 64 )
 
 

      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 R05  are to be initialized before executing RK6B

    R00:  f name      R01 = x0         R04 = h = stepsize                         R06 to R18: temp
                               R02 = y0         R05 = N = number of steps
                               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.
 
 

001  LBL "RK6B"
002  RCL 05
003  STO 18
004  GTO 01
005  LBL 00
006  XEQ IND 00
007  RCL 04
008  ST* Z
009  *
010  RTN
011  LBL 01
012  RCL 03
013  RCL 02
014  RCL 01
015  XEQ 00
016  STO 12
017  3
018  /
019  RCL 03
020  +
021  X<>Y
022  STO 06
023  3
024  /
025  RCL 02
026  +
027  RCL 04
028  3
029  /
030  RCL 01
031  +
032  XEQ 00
033  STO 13
034  1.5
035  /
036  RCL 03
037  +
038  X<>Y
039  STO 07
040  1.5
041  /
042  RCL 02
043  +
044  RCL 04
045  1.5
046  /
047  RCL 01
048  +
049  XEQ 00
050  STO 14
051  CHS
052  RCL 13
053  4
054  *
055  +
056  RCL 12
057  +
058  12
059  /
060  RCL 03
061  +
062  X<>Y
063  STO 08
064  CHS
065  RCL 07
066  4
067  *
068  +
069  RCL 06
070  +
071  12
072  /
073  RCL 02
074  +
075  RCL 04
076  3
077  /
078  RCL 01
079  +
080  XEQ 00
081  STO 15
082  90
083  *
084  RCL 14
085  35
086  *
087  +
088  RCL 13
089  110
090  *
091  -
092  RCL 12
093  25
094  *
095  +
096  48
097  /
098  RCL 03
099  +
100  X<>Y
101  STO 09
102  90
103  *
104  RCL 08
105  35
106  *
107  +
108  RCL 07
109  110
110  *
111  -
112  RCL 06
113  25
114  *
115  +
116  48
117  /
118  RCL 02
119  +
120  RCL 04
121  1.2
122  /
123  RCL 01
124  +
125  XEQ 00
126  STO 16
127  10
128  /
129  RCL 15
130  2
131  /
132  +
133  RCL 14
134  8
135  /
136  -
137  RCL 13
138  11
139  *
140  24
141  /
142  -
143  RCL 12
144  .15
145  *
146  +
147  RCL 03
148  +
149  X<>Y
150  STO 10
151  10
152  /
153  RCL 09
154  2
155  /
156  +
157  RCL 08
158  8
159  /
160  -
161  RCL 07
162  11
163  *
164  24
165  /
166  -
167  RCL 06
168  .15
169  *
170  +
171  RCL 02
172  +
173  RCL 04
174  6
175  /
176  RCL 01
177  +
178  XEQ 00
179  STO 17
180  80
181  *
182  RCL 15
183  118
184  *
185  -
186  RCL 13
187  99
188  *
189  +
190  20
191  *
192  RCL 16
193  128
194  *
195  +
196  RCL 14
197  215
198  *
199  +
200  RCL 12
201  783
202  *
203  -
204  780
205  /
206  RCL 03
207  +
208  X<>Y
209  STO 11
210  80
211  *
212  RCL 09
213  118
214  *
215  -
216  RCL 07
217  99
218  *
219  +
220  20
221  *
222  RCL 10
223  128
224  *
225  +
226  RCL 08
227  215
228  *
229  +
230  RCL 06
231  783
232  *
233  -
234  780
235  /
236  RCL 02
237  +
238  RCL 04
239  RCL 01
240  +
241  STO 01
242  XEQ 00
243  RCL 12
244  +
245  13
246  *
247  RCL 17
248  RCL 16
249  +
250  32
251  *
252  +
253  RCL 15
254  RCL 14
255  +
256  55
257  *
258  +
259  .5
260  %
261  ST+ 03
262  R^
263  RCL 06
264  +
265  13
266  *
267  RCL 11
268  RCL 10
269  +
270  32
271  *
272  +
273  RCL 09
274  RCL 08
275  +
276  55
277  *
278  +
279  .5
280  %
281  ST+ 02
282  DSE 18
283  GTO 01         a three-byte GTO
284  RCL 03
285  RCL 02
286  RCL 01
287  END

( 390 bytes / SIZE 019 )
 
 
      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(1) and z(1)
 

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
10     STO 05    XEQ "RK6B"   yields      x =  1                                ( in 2mn44s )
                              RDN                      y(1) =  0.367879433             ( y-error = -9 10-9 )
                              RDN                      z(1) = -0.735758865             ( z-error = 17 10-9 )
 

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 26 = 64 )
 

4°) A method with order 8
 

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

-Methods of this order require 11 stages. The following one was discovered in 1972 by Cooper and Verner.
-The formula duplicates the Taylor serie through the term in h8
-The following program uses 41 coefficients which are of the form ( a + b 211/2 )/c  where a , b , c are integers.
-They are to be stored in registers R11 thru R51  ( all the digits are correct )
 

Data Registers:     Registers R00 to R04  and R11 to R51 are to be initialized before executing "RK8"

    R00:  f name      R01 = x0         R03 =  h = the stepsize                      R05  to R10 and R52: temp
                               R02 = y0         R04 = N = the number of steps        R11 to R51:  the 41 coefficients listed below

   R11 = 0.8273268354                   R21 = 0.5766714727                   R31 = 0.1289862930                           R41 = 2.031083139
   R12 = 0.5                                     R22 = 0.1855068535                   R32 = -0.01186868389                       R42 = -0.6379313502
   R13 = 0.1726731646                   R23 = 0.3865246267                   R33 = 0.002002165993                      R43 = 0.9401094520
   R14 = 9                                        R24 = -0.4634553896                 R34 = 0.9576053440                           R44 = -2.229158210
   R15= 14                                       R25 = 0.3772937693                   R35 = -0.6325461607                         R45 = 7.553840442
   R16 = 0.25                                   R26 = 0.1996369936                   R36 = 0.1527777778                          R46 = -7.164951553
   R17 = 0.8961811934                   R27 = 0.09790045616                 R37 = -0.009086961101                     R47 = 2.451380432
   R18 = -0.2117115009                  R28 = 0.3285172131                  R38 = 0.03125                                     R48 = -0.5512205631
   R19 = 7                                        R29 = -0.3497052863                 R39 = 1.062498447                             R49 = 0.05
   R20 = 0.06514850915                 R30 = -0.03302551131               R40 = -1.810863083                           R50 = 0.2722222222
                                                                                                                                                                        R51 = 2.8125

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

001  LBL "RK8"
002  RCL 04
003  STO 52
004  GTO 01
005  LBL 00
006  XEQ IND 00
007  RCL 03
008  *
009  RTN
010  LBL 01
011  RCL 02
012  RCL 01
013  XEQ 00
014  STO 05
015  RCL 12
016  *
017  RCL 02
018  +
019  RCL 03
020  RCL 12
021  *
022  RCL 01
023  +
024  XEQ 00
025  STO 06
026  RCL 05
027  +
028  RCL 16
029  *
030  RCL 02
031  +
032  RCL 03
033  RCL 12
034  *
035  RCL 01
036  +
037  XEQ 00
038  STO 07
039  RCL 17
040  *
041  RCL 06
042  RCL 18
043  *
044  +
045  RCL 05
046  RCL 19
047  /
048  +
049  RCL 02
050  +
051  RCL 03
052  RCL 11
053  *
054  RCL 01
055  +
056  XEQ 00
057  STO 08
058  RCL 20
059  *
060  RCL 07
061  RCL 21
062  *
063  +
064  RCL 05
065  RCL 22
066  *
067  +
068  RCL 02
069  +
070  RCL 03
071  RCL 11
072  *
073  RCL 01
074  +
075  XEQ 00
076  STO 09
077  RCL 23
078  *
079  RCL 08
080  RCL 24
081  *
082  +
083  RCL 07
084  RCL 25
085  *
086  +
087  RCL 05
088  RCL 26
089  *
090  +
091  RCL 02
092  +
093  RCL 03
094  RCL 12
095  *
096  RCL 01
097  +
098  XEQ 00
099  STO 10
100  RCL 27
101  *
102  RCL 09
103  RCL 28
104  *
105  +
106  RCL 08
107  RCL 29
108  *
109  +
110  RCL 07
111  RCL 30
112  *
113  +
114  RCL 05
115  RCL 31
116  *
117  +
118  RCL 02
119  +
120  RCL 03
121  RCL 13
122  *
123  RCL 01
124  +
125  XEQ 00
126  STO 06
127  RCL 14
128  /
129  RCL 10
130  RCL 32
131  *
132  +
133  RCL 09
134  RCL 33
135  *
136  +
137  RCL 05
138  RCL 15
139  /
140  +
141  RCL 02
142  +
143  RCL 03
144  RCL 13
145  *
146  RCL 01
147  +
148  XEQ 00
149  STO 07
150  RCL 34
151  *
152  RCL 06
153  RCL 35
154  *
155  +
156  RCL 10
157  RCL 36
158  *
159  +
160  RCL 09
161  RCL 37
162  *
163  +
164  RCL 05
165  RCL 38
166  *
167  +
168  RCL 02
169  +
170  RCL 03
171  RCL 12
172  *
173  RCL 01
174  +
175  XEQ 00
176  STO 08
177  RCL 39
178  *
179  RCL 07
180  RCL 40
181  *
182  +
183  RCL 06
184  RCL 41
185  *
186  +
187  RCL 10
188  RCL 42
189  *
190  +
191  RCL 09
192  RCL 14
193  /
194  +
195  RCL 05
196  RCL 15
197  /
198  +
199  RCL 02
200  +
201  RCL 01
202  RCL 03
203  ST+ 01
204  RCL 11
205  *
206  +
207  XEQ 00
208  X<> 09
209  RCL 48
210  *
211  RCL 08
212  RCL 44
213  *
214  +
215  RCL 07
216  RCL 45
217  *
218  +
219  RCL 06
220  RCL 46
221  *
222  +
223  RCL 10
224  RCL 47
225  *
226  +
227  RCL 09
228  RCL 43
229  *
230  +
231  RCL 02
232  +
233  RCL 01
234  XEQ 00
235  RCL 05
236  +
237  RCL 49
238  *
239  RCL 09
240  RCL 07
241  +
242  RCL 50
243  *
244  +
245  RCL 08
246  RCL 51
247  /
248  +
249  ST+ 02
250  DSE 52
251  GTO 01            a three-byte GTO
252  RCL 02
253  RCL 01
254  END
 

( 329 bytes / SIZE 053 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 

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

 -Initialize registers R11 thru R51.

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
10   STO 04    XEQ "RK8"   yields            x = 1                        ( in 2mn )
                              X<>Y                     y(1) = 0.3678794412   correct to 10 D!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 28 = 256 )
 
 

      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 R05 and R18 to R58 are to be initialized before executing "RK8B"

    R00:  f name      R01 = x0         R04 = h = stepsize                         R06 to R17: temp    R59: counter
                               R02 = y0         R05 = N = number of steps          R18 thru  R58 = the same 41 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.
 

001  LBL "RK8B"
002  RCL 05
003  STO 59
004  GTO 01
005  LBL 00
006  XEQ IND 00
007  RCL 04
008  ST* Z
009  *
010  RTN
011  LBL 01
012  RCL 03
013  RCL 02
014  RCL 01
015  XEQ 00
016  STO 12
017  RCL 19
018  *
019  RCL 03
020  +
021  X<>Y
022  STO 06
023  RCL 19
024  *
025  RCL 02
026  +
027  RCL 04
028  RCL 19
029  *
030  RCL 01
031  +
032  XEQ 00
033  STO 13
034  RCL 12
035  +
036  RCL 23
037  *
038  RCL 03
039  +
040  X<>Y
041  STO 07
042  RCL 06
043  +
044  RCL 23
045  *
046  RCL 02
047  +
048  RCL 04
049  RCL 19
050  *
051  RCL 01
052  +
053  XEQ 00
054  STO 14
055  RCL 24
056  *
057  RCL 13
058  RCL 25
059  *
060  +
061  RCL 12
062  RCL 26
063  /
064  +
065  RCL 03
066  +
067  X<>Y
068  STO 08
069  RCL 24
070  *
071  RCL 07
072  RCL 25
073  *
074  +
075  RCL 06
076  RCL 26
077  /
078  +
079  RCL 02
080  +
081  RCL 04
082  RCL 18
083  *
084  RCL 01
085  +
086  XEQ 00
087  STO 15
088  RCL 27
089  *
090  RCL 14
091  RCL 28
092  *
093  +
094  RCL 12
095  RCL 29
096  *
097  +
098  RCL 03
099  +
100  X<>Y
101  STO 09
102  RCL 27
103  *
104  RCL 08
105  RCL 28
106  *
107  +
108  RCL 06
109  RCL 29
110  *
111  +
112  RCL 02
113  +
114  RCL 04
115  RCL 18
116  *
117  RCL 01
118  +
119  XEQ 00
120  STO 16
121  RCL 30
122  *
123  RCL 15
124  RCL 31
125  *
126  +
127  RCL 14
128  RCL 32
129  *
130  +
131  RCL 12
132  RCL 33
133  *
134  +
135  RCL 03
136  +
137  X<>Y
138  STO 10
139  RCL 30
140  *
141  RCL 09
142  RCL 31
143  *
144  +
145  RCL 08
146  RCL 32
147  *
148  +
149  RCL 06
150  RCL 33
151  *
152  +
153  RCL 02
154  +
155  RCL 04
156  RCL 19
157  *
158  RCL 01
159  +
160  XEQ 00
161  STO 17
162  RCL 34
163  *
164  RCL 16
165  RCL 35
166  *
167  +
168  RCL 15
169  RCL 36
170  *
171  +
172  RCL 14
173  RCL 37
174  *
175  +
176  RCL 12
177  RCL 38
178  *
179  +
180  RCL 03
181  +
182  X<>Y
183  STO 11
184  RCL 34
185  *
186  RCL 10
187  RCL 35
188  *
189  +
190  RCL 09
191  RCL 36
192  *
193  +
194  RCL 08
195  RCL 37
196  *
197  +
198  RCL 06
199  RCL 38
200  *
201  +
202  RCL 02
203  +
204  RCL 04
205  RCL 20
206  *
207  RCL 01
208  +
209  XEQ 00
210  STO 13
211  RCL 21
212  /
213  RCL 17
214  RCL 39
215  *
216  +
217  RCL 16
218  RCL 40
219  *
220  +
221  RCL 12
222  RCL 22
223  /
224  +
225  RCL 03
226  +
227  X<>Y
228  STO 07
229  RCL 21
230  /
231  RCL 11
232  RCL 39
233  *
234  +
235  RCL 10
236  RCL 40
237  *
238  +
239  RCL 06
240  RCL 22
241  /
242  +
243  RCL 02
244  +
245  RCL 04
246  RCL 20
247  *
248  RCL 01
249  +
250  XEQ 00
251  STO 14
252  RCL 41
253  *
254  RCL 13
255  RCL 42
256  *
257  +
258  RCL 17
259  RCL 43
260  *
261  +
262  RCL 16
263  RCL 44
264  *
265  +
266  RCL 12
267  RCL 45
268  *
269  +
270  RCL 03
271  +
272  X<>Y
273  STO 08
274  RCL 41
275  *
276  RCL 07
277  RCL 42
278  *
279  +
280  RCL 11
281  RCL 43
282  *
283  +
284  RCL 10
285  RCL 44
286  *
287  +
288  RCL 06
289  RCL 45
290  *
291  +
292  RCL 02
293  +
294  RCL 04
295  RCL 19
296  *
297  RCL 01
298  +
299  XEQ 00
300  STO 15
301  RCL 46
302  *
303  RCL 14
304  RCL 47
305  *
306  +
307  RCL 13
308  RCL 48
309  *
310  +
311  RCL 17
312  RCL 49
313  *
314  +
315  RCL 16
316  RCL 21
317  /
318  +
319  RCL 12
320  RCL 22
321  /
322  +
323  RCL 03
324  +
325  X<>Y
326  STO 09
327  RCL 46
328  *
329  RCL 08
330  RCL 47
331  *
332  +
333  RCL 07
334  RCL 48
335  *
336  +
337  RCL 11
338  RCL 49
339  *
340  +
341  RCL 10
342  RCL 21
343  /
344  +
345  RCL 06
346  RCL 22
347  /
348  +
349  RCL 02
350  +
351  RCL 04
352  RCL 18
353  *
354  RCL 01
355  +
356  XEQ 00
357  X<> 16
358  RCL 55
359  *
360  RCL 15
361  RCL 51
362  *
363  +
364  RCL 14
365  RCL 52
366  *
367  +
368  RCL 13
369  RCL 53
370  *
371  +
372  RCL 17
373  RCL 54
374  *
375  +
376  RCL 16
377  RCL 50
378  *
379  +
380  RCL 03
381  +
382  X<>Y
383  X<> 10
384  RCL 55
385  *
386  RCL 09
387  RCL 51
388  *
389  +
390  RCL 08
391  RCL 52
392  *
393  +
394  RCL 07
395  RCL 53
396  *
397  +
398  RCL 11
399  RCL 54
400  *
401  +
402  RCL 10
403  RCL 50
404  *
405  +
406  RCL 02
407  +
408  RCL 04
409  RCL 01
410  +
411  STO 01
412  XEQ 00
413  RCL 12
414  +
415  RCL 56
416  *
417  RCL 16
418  RCL 14
419  +
420  RCL 57
421  *
422  +
423  RCL 15
424  RCL 58
425  ST/ 09
426  /
427  +
428  ST+ 03
429  X<>Y
430  RCL 06
431  +
432  RCL 56
433  *
434  RCL 10
435  RCL 08
436  +
437  RCL 57
438  *
439  +
440  RCL 09
441  +
442  ST+ 02
443  DSE 59
444  GTO 01             a three-byte GTO
445  RCL 03
446  RCL 02
447  RCL 01
448  END
 

( 593 bytes / SIZE 060 )
 
 
 
      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(1) and z(1)

 -Initialize registers R18 thru R58.

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
10   STO 05    XEQ "RK8B"   yields        x = 1                         in 3mn10s
                              RDN                      y(1) =  0.3678794412   correct to 10D.
                              RDN                      z(1) = -0.7357588824   z-error = -10-10
 

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 28 = 256 )
 

5°) A general Explicit-Runge-Kutta program

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

-All n stage explicit Runge-Kutta formulae are determined by ( n2 + 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci

    k1 = h.f(x;y)
    k2 = h.f(x+b2.h;y+a2,1.k1)
    k3 = h.f(x+b3.h;y+a3,1.k1+a3,2.k2)
    .............................................................                                         (  actually,     ai,1 + ai,2 + .......... + ai,i-1 = bi  )

    kn = h.f(x+bn.h;y+an,1.k1+an,2.k2+ .......... +an,n-1.kn-1)

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

-Therefore, a single program may be used for all Runge-Kutta methods, provided the coefficients are stored in appropriate registers.
 

Data Registers:     Registers R00 to R05 and    Rn+10 to R(n2+5.n+16)/2     are to be initialized before executing "ERK"

    R00:  f name      R01 = x0                                      R06 to R09:  temp               Rn+10 to R(n2+5.n+16)/2 =  the ( n2 + 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci
                               R02 = y0                                      R10 = k1                                                                            which are to be stored as shown below:
                               R03 = h = stepsize                        R11 = k2
                               R04 = N = number of steps          ............
                               R05 = n  = number of stages        Rn+9 = kn
 

                                                    Rn+10 = a2,1                                                                                  Rn+11 = b2
                                                    Rn+12 = a3,1          Rn+13 = a3,2                                                  Rn+14 = b3

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

                                                    R(n2+n+18)/2 = an,1  ................     R(n2+3n+14)/2 = an,n-1         R(n2+3n+16)/2 = bn

                                                    R(n2+3n+18)/2 = c1 ............................................................         R(n2+5n+16)/2 = cn
 

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

01  LBL "ERK"
02  RCL 04
03  STO 06
04  LBL 01
05  RCL 05
06  10.9
07  STO 07
08  +
09  STO 08
10  RCL 05
11  9
12  +
13  .1
14  %
15  11
16  +
17  STO 09
18  RCL 02
19  RCL 01
20  XEQ IND 00
21  RCL 03
22  *
23  STO 10
24  LBL 02
25  RCL 07
26  INT
27  .1
28  %
29  LASTX
30  1/X
31  +
32  STO 07
33  CLX
34  LBL 03
35  RCL IND 07
36  RCL IND 08
37  *
38  +
39  ISG 08
40  ISG 07
41  GTO 03
42  RCL 02
43  +
44  RCL 03
45  RCL IND 08
46  *
47  RCL 01
48  +
49  XEQ IND 00
50  RCL 03
51  *
52  STO IND 09
53  ISG 08
54  ISG 09
55  GTO 02
56  RCL 05
57  ST- 09
58  CLX
59  LBL 04
60  RCL IND 08
61  RCL IND 09
62  *
63  +
64  ISG 08
65  ISG 09
66  GTO 04
67  ST+ 02
68  RCL 03
69  ST+ 01
70  DSE 06
71  GTO 01
72  RCL 02
73  RCL 01
74  END

( 110 bytes / SIZE ( n2+5.n+18 )/2 )
 
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 

-For example, if you are using the classical four stage method:      4   STO 05   and the 13 coefficients:

                                                           1/2                 1/2                                             R14                   R15
                                                            0    1/2          1/2       are to be stored into        R16  R17          R18        respectively
                                                            0     0     1      1                                               R19  R20  R21  R22
                                                           1/6  1/3  1/3  1/6                                              R23  R24  R25  R26

-With the 7 stage method used in "RK6":                                    7  STO 05    and the 34 coefficients:

                                1/3                                                                        1/3                                                R17                                            R18
                                 0           2/3                                                          2/3                                                R19  R20                                   R21
                                1/12       1/3    -1/12                                              1/3                                               R22  R23  R24                           R25
                              25/48   -55/24  35/48       15/8                                5/6      are to be stored into          R26  R27  R28  R29                   R30
                                3/20   -11/24   -1/8          1/2       1/10                   1/6                                              R31  R32  R33  R34  R35          R36
                           -261/260  33/13  43/156  -118/39  32/195   80/39      1                                                R37  R38  R39  R40  R41  R42  R43
                              13/200     0       11/40      11/40     4/25      4/25    13/200                                          R44  R45  R46  R47  R48  R49  R50

                                                      ................ and so on ..........

-For the other instructions, cf "RK6"
-With the same differential equation, execution time = 2mn50s for a 7 stage method ( instead of 89s )
 ( the number of bytes is decreased, but execution time is increased ... )
 

       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+13 to R(n2+7.n+22)/2     are to be initialized before executing "ERKB"

    R00:  f name      R01 = x0                                      R07 to R12:  temp               R2n+13 to R(n2+7n+22)/2 =  the ( n2 + 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci
                               R02 = y0                                      R13 to Rn+12 = kiy                                                            which are to be stored as shown below:
                               R03 = z0                                      Rn+13 to R2n+12 = kiz
                               R04 = h = stepsize
                               R05 = N  = number of steps
                               R06 = n  = number of stages
 

                                                    R2n+13 = a2,1                                                                                  R2n+14 = b2
                                                    R2n+15 = a3,1          R2n+16 = a3,2                                                R2n+17 = b3

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

                                                    R(n2+3n+24)/2 = an,1  ................     R(n2+5n+20)/2 = an,n-1         R(n2+5n+22)/2 = bn

                                                    R(n2+5n+24)/2 = c1 ............................................................           R(n2+7n+22)/2 = 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.
 
 
 

001  LBL "ERKB"
002  RCL 05
003  STO 07
004  LBL 01
005  13.9
006  STO 08
007  RCL 06
008  +
009  STO 09
010  LASTX
011  +
012  STO 10
013  LASTX
014  12
015  +
016  .1
017  %
018  14
019  +
020  STO 11
021  RCL 06
022  +
023  STO 12
024  RCL 03
025  RCL 02
026  RCL 01
027  XEQ IND 00
028  RCL 04
029  ST* Z
030  *
031  STO IND 09
032  X<>Y
033  STO 13
034  LBL 02
035  RCL 08
036  INT
037  .1
038  %
039  13
040  +
041  STO 08
042  RCL 06
043  +
044  STO 09
045  CLST
046  LBL 03
047  X<>Y
048  RCL IND 08
049  RCL IND 10
050  *
051  +
052  X<>Y
053  RCL IND 09
054  RCL IND 10
055  *
056  +
057  ISG 10
058  ISG 09
059  CLX
060  ISG 08
061  GTO 03
062  RCL 03
063  +
064  X<>Y
065  RCL 02
066  +
067  RCL 04
068  RCL IND 10
069  *
070  RCL 01
071  +
072  XEQ IND 00
073  RCL 04
074  ST* Z
075  *
076  STO IND 12
077  X<>Y
078  STO IND 11
079  ISG 12
080  CLX
081  ISG 10
082  ISG 11
083  GTO 02
084  RCL 06
085  ST- 11
086  ST- 12
087  CLST
088  LBL 04
089  X<>Y
090  RCL IND 10
091  RCL IND 11
092  *
093  +
094  X<>Y
095  RCL IND 10
096  RCL IND 12
097  *
098  +
099  ISG 12
100  CLX
101  ISG 10
102  ISG 11
103  GTO 04
104  ST+ 03
105  X<>Y
106  ST+ 02
107  RCL 04
108  ST+ 01
109  DSE 07
110  GTO 01          a three-byte GTO
111  RCL 03
112  RCL 02
113  RCL 01
114  END
 

( 167 bytes / SIZE ( n2+ 7.n + 24 )/2 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 

-For example, if you are using the classical four stage method:      4   STO 06   and the 13 coefficients:

                                                           1/2                 1/2                                             R21                   R22
                                                            0    1/2          1/2       are to be stored into        R23  R24          R25        respectively
                                                            0     0     1      1                                               R26  R27  R28  R29
                                                           1/6  1/3  1/3  1/6                                              R30  R31  R32  R33

-With the 7 stage method used in "RK6B":                                    7  STO 06    and the 34 coefficients:

                                1/3                                                                        1/3                                                R27                                            R28
                                 0           2/3                                                          2/3                                                R39  R30                                   R31
                                1/12       1/3    -1/12                                              1/3                                               R32  R33  R34                           R35
                              25/48   -55/24  35/48       15/8                                5/6      are to be stored into          R36  R37  R38  R39                   R40
                                3/20   -11/24   -1/8          1/2       1/10                   1/6                                              R41  R42  R43  R44  R45          R46
                           -261/260  33/13  43/156  -118/39  32/195   80/39      1                                                R47  R48  R49  R50  R51  R52  R53
                              13/200     0       11/40      11/40     4/25      4/25    13/200                                          R54  R55  R56  R57  R58  R59  R60

                                                      ................ and so on ..........

-For the other instructions, cf "RK6B"
-With the same differential system, execution time = 4mn36s for a 7 stage method ( instead of 2mn44s )
 

6°) An implicit Runge-Kutta method of order 8

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

-The formulae used by "RK4" to "RK8" are explicit Runge-Kutta methods in that 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 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 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 function 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 ( when h is divided by 2 , errors are approximately divided by 28 = 256 )
 

       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 ( when h is divided by 2 , errors are approximately divided by 28 = 256 )
 
 

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