Runge-Kutta10

# Runge-Kutta Formula of order 10 for the HP-41

Overview

0°)  169 Constants
1°)  1 Differential Equation
2°)  System of 2 Differential Equations
3°)  System of 3 Differential Equations
4°)  System of n Differential Equations

-These routines apply a 17-stage 10th-order Runge-Kutta formula to solve ODEs

-All s-stage explicit Runge-Kutta formulae without error estimate are determined by ( s2 + 3.s - 2 ) / 2  coefficients ai,j ; bi ; ci
-Here, s = 17 so  169 constants

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  )

k17 = h.f(x+b17.h;y+a17,1.k1+a17,2.k2+ .......... +a17,16.k16)

then,     y(x+h) = y(x) + c1.k1+ ................ + c17.k17

0°)  169 Constants

-All the programs listed in paragraphs 1-2-3 use 169 constants.
-"C169" stores these numbers into registers R65 thru R233
-Registers R07 to R64 are also used.

-Allways execute "C169" before executing the other programs... unless you store your own 169 constants into R65 to R233 !

•  R65 = a2,1                                                                                •  R66 = b2
•  R67 = a3,1        •  R68 = a3,2                                                    •  R69 = b3

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

•  R200= a17,1  ........................................    •  R215 = a17,16       •  R216 = b17

•  R217 = c1 .........................................................................       •  R233 = c17

 01 LBL "C169"  02 8.085  03 CLRGX  04 .8950802958  05 STO 07  06 .1979668312  07 STO 09  08 .0729547847  09 CHS  10 STO 10  11 .8512362396  12 CHS  13 STO 12  14 .3983201123  15 STO 13  16 3.639372632  17 STO 14  18 1.54822877  19 STO 15  20 2.122217147  21 CHS  22 STO 16  23 1.583503985  24 CHS  25 STO 17  26 1.715616082  27 CHS  28 STO 18  29 .0244036406  30 CHS  31 STO 19  32 .3090367612  33 STO 20  34 .9151765614  35 CHS  36 STO 21  37 1.454534402  38 STO 22  39 .7773336436 40 STO 34  41 CHS  42 STO 25  43 .0910895662  44 STO 33  45 CHS  46 STO 27  47 .5393578408  48 STO 35  49 .1  50 STO 36  51 STO 51  52 STO 61  53 .1571786658  54 STO 50  55 CHS  56 STO 38  57 .1817813007  58 STO 52  59 .675  60 STO 53  61 CHS  62 STO 67  63 .3427581598  64 STO 54  65 CHS  66 STO 66  67 .2591112146  68 STO 56  69 CHS  70 STO 65  71 .3582789667  72 CHS  73 STO 57  74 1.045948959  75 CHS  76 STO 58  77 .9303278454  78 STO 59 79 1.779509594  80 STO 60  81 .2825475695  82 CHS  83 STO 62  84 .1593273501  85 CHS  86 STO 63  87 .1455158946  88 CHS  89 STO 64  90 1  91 STO 68  92 30  93 1/X  94 STO 69  95 STO 71  96 STO 85  97 CHS  98 STO 83  99 40 100 1/X 101 STO 70 102 CHS 103 STO 84 104 .05 105 STO 73 106 CHS 107 STO 82 108 .04 109 STO 75 110 CHS 111 STO 81 112 .1892374781 113 STO 77 114 STO 80 115 .2774291885 116 STO 78 117 STO 79 118 7.155079 119 REGMOVE 120 9.096 121 CLRGX 122 .1 123 STO 07 124 STO 08 125 STO 83 126 .9151765614 127 CHS 128 STO 09 129 1.454534402 130 STO 10 131 .5393578408 132 STO 11 133 .2022591903 134 ENTER 135 STO 12 136 3 137 * 138 STO 14 139 + 140 STO 15 141 .5 142 - 143 STO 20 144 .1840247147 145 STO 16 146 .1979668312 147 STO 18 148 .0729547847 149 CHS 150 STO 19 151 .087900734 152 STO 21 153 .4104597025 154 STO 24 155 .4827137537 156 STO 25 157 .9810741902 158 STO 26 159 .0859700505 160 STO 27 161 .330885963 162 STO 30 163 .4896629573 164 STO 31 165 .0731856375 166 CHS 167 STO 32 168 1.2 169 1/X 170 STO 33 171 STO 96 172 .1209304491 173 STO 34 174 .2601246758 175 STO 38 176 .0325402622 177 STO 39 178 .0595780212 179 CHS 180 STO 40 181 .3540173659 182 STO 41 183 .1108543796 184 STO 42 185 .0605761488 186 CHS 187 STO 47 188 .3217637056 189 STO 48 190 .5104857256 191 STO 49 192 .882527662 193 STO 50 194 .1120544148 195 STO 51 196 .1449427759 197 CHS 198 STO 56 199 .3332697191 200 CHS 201 STO 57 202 .4992692296 203 STO 58 204 .5095046089 205 STO 59 206 .6426157582 207 STO 60 208 .113976784 209 STO 61 210 .0768813364 211 CHS 212 STO 66 213 .2395273603 214 STO 67 215 .3977746624 216 STO 68 217 .0107558957 218 STO 69 219 .3277691242 220 CHS 221 STO 70 222 .3573842418 223 STO 71 224 .0798314528 225 STO 72 226 .0520329687 227 CHS 228 STO 77 229 .0576954146 230 CHS 231 STO 78 232 .1947819157 233 STO 79 234 .1453849232 235 STO 80 236 .078294271 237 CHS 238 STO 81 239 .1145032994 240 CHS 241 STO 82 242 .117472338 243 STO 83 244 .9851156101 245 STO 84 246 .330885963 247 STO 87 248 .4896629573 249 STO 88 250 1.378964866 251 CHS 252 STO 89 253 .861164195 254 CHS 255 STO 90 256 5.784288136 257 STO 91 258 3.28807762 259 STO 92 260 2.386339051 261 CHS 262 STO 93 263 3.254793425 264 CHS 265 STO 94 266 2.163435417 267 CHS 268 STO 95 269 7.06509 270 REGMOVE 271 END

( 1271 bytes / SIZE 234 )

 STACK INPUT OUTPUT X / /

-Execution time = 35 seconds.

1°)  1 Differential Equation

-We want to solve  y' = f(x,y)  with the initial conditions    y(x0) = y0

Data Registers:          •  R00:  f name                                       ( Registers R00 to R04 & R65 to R233 are to be initialized before executing "RK10" )

•  R01 = x0      •  R03 = h  = stepsize                 R05 & R10: temp            •  R65 to •  R233:  the 169 coefficients stored by "C169"
•  R02 = y0      •  R04 = N = number of steps    R11 to R27:  k1 to k17

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

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

( 117 bytes / SIZE 234 )

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

Example:    y' = -2 x.y   y(0) = 1  y'(0) = 0      Evaluate  y(1)

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

XEQ "C169"

"T"  ASTO 00

0     STO 01
1     STO 02    and if we choose  h = 0.1
0.1   STO 03
10    STO 04    XEQ "RK10"   yields          x = 1                                       ---Execution time = 11m12s---
X<>Y                       y(1) = 0.3678794412

-The exact result is 1/e =  0.3678794412... so there is no error at all !
-Press R/S to continue with the same h and N or modify R03-R04

2°)  System of 2 Differential Equations

-"2RK10" solves the system

y' = f(x,y,z)      with the initial values     y(x0) = y0
z' = g(x,y,z)                                         z(x0) = z0

Data Registers:          •  R00:  f name                                     ( Registers R00 to R05 & R65 to R233 are to be initialized before executing "2RK10" )

•  R01 = x0      •  R04 = h  = stepsize                 R06 & R10: temp            •  R65 to •  R233:  the 169 coefficients stored by "C169"
•  R02 = y0      •  R05 = N = number of steps    R11 to R27:  k1 to k17
•  R03 = z0                                                        R28 to R44:  l1 to l17
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 "2RK10"  02 RCL 05  03 STO 07  04 11.001  05 STO 45  06 15.999  07 STO 46  08 GTO 01  09 LBL 00  10 X<>Y  11 RCL IND 08  12 RCL IND 10  13 *  14 +  15 X<>Y  16 RCL IND 09  17 RCL IND 10  18 * 19 +  20 ISG 10  21 ISG 09  22 ISG 08  23 GTO 00  24 RCL 03   25 +  26 X<>Y  27 RCL 02  28 +  29 RTN  30 LBL 01   31 .01  32 STO 08          33 16  34 STO 06  35 65.9  36 STO 10 37 RCL 03  38 RCL 02  39 RCL 01  40 XEQ IND 00  41 RCL 04  42 ST* Z  43 *  44 STO 28  45 X<>Y  46 STO 11  47 LBL 02  48 RCL 08  49 FRC  50 RCL 45  51 +  52 STO 08  53 17.9  54 + 55 STO 09  56 CLST  57 XEQ 00  58 RCL 04  59 RCL IND 10  60 *  61 RCL 01  62 +  63 XEQ IND 00  64 RCL 04  65 ST* Z  66 *  67 STO IND 09  68 X<>Y  69 STO IND 08  70 ISG 10  71 DSE 06 72 GTO 02  73 RCL 46  74 ST- 08  75 ST- 09  76 CLST  77 XEQ 00  78 STO 02  79 X<>Y  80 STO 03  81 RCL 04  82 ST+ 01  83 DSE 07  84 GTO 01  85 RCL 03   86 RCL 02           87 RCL 01  88 END

( 150 bytes / SIZE 234 )

 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' + 2.y = 0   Evaluate y(1) and y'(1)

-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

XEQ "C169"

"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 "2RK10"   yields       x =  1                                               ---Execution time = 19m00s---
RDN                      y(1) =  0.3678794412
RDN                      z(1) = -0.7357588823

-The exact results are  y(1) = 1/e =  0.3678794412... & z(1) = -2/e = -0.7357588823...  again no error at all !
-Press R/S to continue with the same h and N or modify  R04-R05

3°)  System of 3 Differential Equations

-"3RK10" solves the system of 3 ODEs:

dy/dx = f(x,y,z,t)                                          y(x0) = y0
dz/dx = g(x,y,z,t)   with the initial values:       z(x0) = z0
dt/dx = h(x,y,z,t)                                          t(x0) = t0

Data Registers:          •  R00:  f name                                     ( Registers R00 to R06 & R65 to R233 are to be initialized before executing "3RK10" )

•  R01 = x0      •  R05 = h  = stepsize                 R07 & R10: temp            •  R65 to •  R233:  the 169 coefficients stored by "C169"
•  R02 = y0      •  R06 = N = number of steps    R11 to R27:  k1 to k17         R62 to R64: temp
•  R03 = z0                                                        R28 to R44:  l1 to l17
•  R04 = t0                                                        R45 to R61:  m1 to m17
Flags: /
Subroutine:  A program that takes  x , y , z , t  in X , Y , Z , T and returns  f(x,y,z,t) , g(x,y,z,t) , h(x,y,z,t)  in registers  Z , Y , X  respectively.

-Line 108 is a three-byte GTO 01

 01 LBL "3RK10"  02 RCL 06  03 STO 62  04 GTO 01  05 LBL 00  06 RCL IND 07  07 RCL IND 10  08 *  09 ST+ 64  10 X<> Z  11 RCL IND 08  12 RCL IND 10  13 *  14 +  15 X<>Y  16 RCL IND 09  17 RCL IND 10  18 *  19 +  20 ISG 10  21 ISG 09  22 ISG 08  23 ISG 07 24 GTO 00  25 RCL 04  26 +  27 X<>Y  28 RCL 03  29 +  30 RCL 64  31 RCL 02  32 +  33 RTN  34 LBL 01  35 16  36 STO 63  37 .01  38 STO 07  39 65.9  40 STO 10  41 RCL 04  42 RCL 03  43 RCL 02  44 RCL 01  45 XEQ IND 00  46 RCL 05 47 ST* T  48 ST* Z  49 *  50 STO 45  51 RDN  52 STO 28  53 X<>Y  54 STO 11  55 LBL 02  56 RCL 07  57 FRC  58 11.001  59 +  60 STO 07  61 17.9  62 +  63 STO 08          64 LASTX  65 INT  66 +  67 STO 09   68 CLST  69 STO 64 70 XEQ 00  71 STO 64  72 CLX  73 RCL 05  74 RCL IND 10  75 *  76 RCL 01  77 +  78 RCL 64  79 X<>Y  80 XEQ IND 00  81 RCL 05  82 ST* T  83 ST* Z  84 *  85 STO IND 09  86 RDN  87 STO IND 08  88 X<>Y  89 STO IND 07  90 ISG 10  91 DSE 63 92 GTO 02  93 15.999  94 ST- 07  95 ST- 08  96 ST- 09  97 CLST  98 STO 64  99 XEQ 00 100 STO 02 101 RDN 102 STO 03 103 X<>Y 104 STO 04 105 RCL 05 106 ST+ 01 107 DSE 62 108 GTO 01 109 RCL 04        110 RCL 03  111 RCL 02 112 RCL 01 113 END

( 189 bytes / SIZE 234 )

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

Example:

dy/dx = -y.z.t                    y(0) = 1         Evaluate  y(1) , z(1) , t(1)  with  h = 0.1   ( whence N = 10 )
dz/dx = x.( y + z - t )         z(0) = 1
dt/dx = x.y - z.t                  t(0) = 2

 01 LBL "V"  02 R^  03 CHS  04 STO L  05 R^  06 ST* Y  07 ST+ L  08 CLX  09 RCL Z  10 R^  11 ST+ L  12 ST* Y  13 X<> Z  14 ST+ Y  15 ST* Z  16 X<> T  17 LASTX  18 *  19 X<>Y  20 RTN  21 END

XEQ "C169"

"V"   ASTO 00

0     STO 01
1     STO 02  STO 03
2     STO 04
0.1   STO 05
10    STO 06    XEQ "3RK10"   >>>>      x   =  1                                                       ---Execution time = 27m55s---
RDN     y(1) =  0.258207907
RDN     z(1) =  1.157623979
RDN     t(1) =  0.842178312

-The exact values are:

y(1) = 0.258207906                      ( rounded to 9 decimals )
z(1) = 1.157623981
t(1) = 0.842178312

-Here, we get a ( small ) error.
-Press R/S to continue with the same h and N or modify   R05-R06

-Synthetic registers M , N , O , P , Q  and register R64 may also be used in the subroutine to calculate f , g , h
-If it's not enough, store the constants in - for example - R71 to R239 and replace line 39 ( 65.9 ) by 71.9

Notes:

-These routines are slow with a real HP41
-But with V41 or 41CL, the results are obtained in a few seconds.

4°)  System of n Differential Equations

-Here, we want to solve the system:

dy1/dx = f1(x,y1, ... ,yn)
dy2/dx = f2(x,y1, ... ,yn)
..................                                        with the initial values:     yi(x0)         i = 1 , ... , n

dyn/dx = fn(x,y1, ... ,yn)

Data Registers:           •  R00 = subroutine name

( Registers R00 thru R03 , R20 thru R20+n and R21+20n thru R189+20n are to be initialized before executing "NRK10" )

•  R01 = n  = number of equations = number of functions         R04 thru R10 & R19: temp    R11 to R18: unused
•  R02 = h  = stepsize
•  R03 = N = number of steps

•  R20 = x0
•  R21 = y1(x0)                                                                           R21+n thru R20+20.n: temp
•  R22 =  y2(x0)
...............         ( initial values )

•  R20+n = yn(x0)

•  R21+20n = a2,1                                                                                •  R22+20n = b2
•  R23+20n = a3,1        •  R24+20n = a3,2                                            •  R25+20n = b3

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

•  R156+20n = a17,1  ..............................    •  R171+20n = a17,16       •  R172+20n = b17

•  R173+20n = c1 ........................................................................       •  R189+20n = c17

Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+n, ... , R20+2n  respectively
with  x, y1, ... , yn in regiters R20, R21, ... , R20+n

 01 LBL "NRK10"  02 RCL 03  03 STO 10  04 RCL 01  05 .1  06 %  07 +  08 STO 05  09 FRC  10 21.02  11 +  12 STO 04  13 RCL 05  14 ST+ X  15 +  16 STO 06  17 GTO 01  18 LBL 02  19 XEQ IND 00  20 RCL 04  21 RCL 05  22 +  23 RCL 07 24 LBL 00  25 RCL IND Y  26 RCL 02  27 *  28 STO IND Y  29 ISG Y  30 RDN  31 ISG Y  32 GTO 00  33 RTN  34 LBL 03  35 RCL IND Y  36 STO IND Y  37 ISG Y  38 RDN  39 ISG Y  40 GTO 03  41 RTN  42 LBL 07  43 RCL IND 07   44 RCL IND 08  45 *  46 ST+ IND 04 47 ISG 04  48 TEXT0  49 ISG 07  50 GTO 07  51 RCL 01  52 ST- 04  53 RCL 05  54 FRC  55 ST+ 07  56 ISG 08  57 GTO 07  58 RTN  59 LBL 01  60 RCL 05  61 RCL 06  62 +  63 STO 07  64 XEQ 02  65 RCL 04           66 RCL 06   67 XEQ 03  68 RCL 01  69 20 70 *  71 21  72 +  73 .1  74 %  75 +  76 STO 08  77 3.018  78 STO 09  79 RCL 20  80 STO 19  81 LBL 04  82 XEQ 07  83 RCL 02  84 RCL IND 08   85 *  86 ST+ 20  87 XEQ 02  88 RCL 06  89 RCL 04  90 XEQ 03  91 RCL 19  92 STO 20 93 RCL 09  94 INT  95  E3  96 /  97 ISG X  98 ST+ 08  99 RCL 04 100 RCL 05 101 3 102 * 103 + 104 STO 07 105 ISG 09 106 GTO 04 107 RCL 02 108 ST+ 20 109 XEQ 07 110 DSE 10 111 GTO 01 112 RCL 23  113 RCL 22          114 RCL 21 115 RCL 20 116 END

( 191 bytes / SIZE 190+20.n )

 STACK INPUTS OUTPUTS T / y3(x0+N.h) Z / y2(x0+N.h) Y / y1(x0+N.h) X / x0+N.h

Example:     Let's solve again the system:

dy1/dx = - y1y2y3
dy2/dx =  x ( y1+y2-y3 )
dy3/dx =  xy1- y2y3

with the initial values:   y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2

x ; y1 ; y2 ; y3  will be in R20 ; R21 ; R22 ; R23 respectively, and we have to write a program to compute the 3 functions
and store the results in R24 ; R25 ; R26   ( just after the "initial-value registers" ).

-For instance, the following subroutine:

 01  LBL "T"  02  RCL 21  03  RCL 22  04  RCL 23  05  *  06  *  07  CHS  08  STO 24  09  RCL 21  10  RCL 22  11  +  12  RCL 23  13  -  14  RCL 20  15  *  16  STO 25  17  RCL 20  18  RCL 21  19  *  20  RCL 22  21  RCL 23  22  *  23  -  24  STO26  25  RTN

We store the name of the subroutine in R00    "T"  ASTO 00
the numbers of equations in R01                      3    STO 01

-Then we have to store the 169 constants into R21+20n to R189+20n
-Execute for instance the following routine:

 01 LBL "INIT"  02 XEQ "C169"  03 RCL 01  04 20  05 *  06 21.169  07 +  08  E3  09 /  10 65  11 +  12 REGMOVE  13 END

-Then the initial values:

x     in R20                      0   STO 20
y1   in R21                      1   STO 21
y2   in R22                      1   STO 22
y3   in R23                      2   STO 23

-Suppose we want to know   y1(1) , y2(1) ,  y3(1)  with a step size h = 0.1 and therefore N = 10

h is stored in R02                 0.1   STO 02
N is stored in R03                 10    STO 03

and  XEQ "NRK10"  >>>>        x   =  1                                                ---Execution time = 5s---            with V41 in turbo mode
RDN     y(1) =  0.258207907
RDN     z(1) =  1.157623977
RDN     t(1) =  0.842178313

Notes:

-With a real HP41, execution time is about 50 minutes, but with free42, the results are returned almost instantaneously !
-Since SIZE = 190+20n, "NRK10" can solve at most a system of n = 6 differential equations with SIZE 310.
-But with free42, n may reach 40 with SIZE 990 ( beyond register 999, the control numbers would not work )

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  T. Feagin - "A Tenth-Order Runge-Kutta Method with Error Estimate"

-These programs employ the coefficients given in reference [2] without taking into account the error estimate.