Runge-Kutta10B

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

Overview

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

-Thanks to Valentin Albillo, I've found on the web a 16-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 = 16 so there are 151 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  )

k16 = h.f(x+b16.h;y+a16,1.k1+a16,2.k2+ .......... +a16,15.k15)

then,     y(x+h) = y(x) + c1.k1+ ................ + c16.k16

0°)  151 Constants

-All the programs listed in paragraphs 1-2-3 use 151 constants.
-"C151" stores these numbers into registers R65 thru R215
-Registers R01 to R99 are also used.

-Allways execute "C151" before executing the other programs... unless you store your own 151 constants into R65 to R215

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

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

•  R184= a16,1  ........................................    •  R198 = a16,15       •  R199 = b16

•  R200 = c1 .........................................................................       •  R215 = c16

 01 LBL "C151"  02 .0688809661  03 STO 01  04 STO 02  05 STO 09  06 -.8381052035  07 STO 03  08 1.253690049  09 STO 04  10 .4155848452  11 STO 05  12 -.0049094868  13 STO 06  14 .0823217303  15 STO 07  16 -.0085312774  17 STO 08  18 1.048931954  19 STO 10  20 -.7538281773  21 STO 11  22 .8052281597  23 STO 12  24 -.3844632207  25 STO 13  26 .7158687154  27 STO 14  28 -.2399238343  29 STO 15  30 -.0636426116  31 STO 16  32 .1996754314  33 STO 17  34 .6103537951  35 STO 18  36 .3785994723  37 STO 19  38 .8850622527  39 STO 20  40 .0177883395  41 STO 21  42 -.0110502191  43 STO 22 44 -.0043934286  45 STO 23  46 .1059752729  47 STO 24  48 .0040508907  49 STO 25  50 -.0016792971  51 STO 26  52 .1106915584  53 STO 27  54 .2356604623  55 STO 28  56 .0889336269  57 STO 29  58 .0413883471  59 STO 30  60 -.852903036  61 STO 31  62 -.0230808755  63 STO 32  64 .0096859219  65 STO 33  66 .8014497793  67 STO 34  68 .301134226  69 STO 35  70 .0904869376  71 STO 36  72 .032044272  73 STO 37  74 .1243376822  75 STO 38  76 -.3073152176  77 STO 39  78 .1644784368  79 STO 40  80 -.0410068117  81 STO 41  82 .4066198993  83 STO 42  84 .1866998712  85 STO 43  86 .6563450699 87 STO 44  88 -.1281584914  89 STO 45  90 -.0949429224  91 STO 46  92 -.1445034451  93 STO 47  94 .921349407  95 STO 48  96 -.1326530105  97 STO 49  98 .0166148587  99 STO 50 100 -.6446177324 101 STO 51 102 .4714846368 103 STO 52 104 .1510115445 105 STO 53 106 .4155848452 107 STO 54 108 STO 65 109 -.3539426289 110 STO 55 111 -.1299925006 112 STO 56 113 .0729671747 114 STO 57 115 1.234092882 116 STO 58 117 .2432547931 118 STO 59 119 -.0488641534 120 STO 60 121 -.5144237141 122 STO 61 123 .0708850045 124 STO 62 125 -.2055533399 126 STO 63 127 .0471613277 128 STO 64 129 -1.288835964 130 STO 66 131 -.2802086976 132 STO 67 133 2.955025311 134 STO 68 135 2.687294528 136 STO 69 137 4.94062343 138 STO 70 139 -1.10669748 140 STO 71 141 -1.014991843 142 STO 72 143 2.537509826 144 STO 73 145 -3.050854355 146 STO 74 147 -4.48248369 148 STO 75 149 -1.212578432 150 STO 76 151 .6838026329 152 STO 77 153 -.2644034429 154 STO 78 155 -.0668111881 156 STO 79 157 .2298317881 158 STO 80 159 .6407414965 160 STO 81 161 .4294291108 162 STO 82 163 -.0096707888 164 STO 83 165 .0161366667 166 STO 84 167 -.0359631432 168 STO 85 169 -.0490267009 170 STO 86 171 .0138853767 172 STO 87 173 -.0232947511 174 STO 88 175 .0041107147 176 STO 89 177 .8849651386 178 STO 90 179 1.289959279 180 STO 91 181 .135032652 182 STO 92 183 -1.575337353 184 STO 93 185 -1.136141441 186 STO 94 187 -2.606299763 188 STO 95 189 -.1422551486 190 STO 96 191 -1.300502794 192 STO 97 193 2.707971631 194 STO 98 195 2.435542969 196 STO 99 197 1.065099 198 REGMOVE 199 -.2645818252 200 STO 01 201 .2843174728 202 STO 02 203 -.015673365 204 STO 03 205 .6035525316 206 STO 04 207 .4155848452 208 STO 05 209 STO 20 210 -.8217014524 211 STO 06 212 1.255912216 213 STO 07 214 -.0161528412 215 STO 08 216 -.0229703899 217 STO 09 218 -.0279603438 219 STO 10 220 -.7509429873 221 STO 11 222 -.0083102909 223 STO 12 224 .0278552386 225 STO 13 226 .0402142154 227 STO 14 228 1.617136464 229 STO 15 230 -1.396283777 231 STO 16 232 -.0142443916 233 STO 17 234 .7570905434 235 STO 18 236 -.2240573576 237 STO 19 238 .2772640185 239 STO 21 240 .0945381837 241 STO 22 242 .841727543 243 STO 23 244 -.9066525983 245 STO 24 246 -.0933485803 247 STO 25 248 4.088877589 249 STO 26 250 .795399865 251 STO 27 252 -.0485917515 253 STO 28 254 .1481985146 255 STO 29 256 -1.770211481 257 STO 30 258 1.838210845 259 STO 31 260 .0490934446 261 STO 32 262 -3.803788231 263 STO 33 264 .3315723987 265 STO 34 266 -.8422897608 267 STO 35 268 1 269 STO 36 270 .0318192746 271 STO 37 272 CLX 273 STO 38 274 STO 40 275 STO 41 276 .0468136929 277 STO 39 278 1.375535362 279 STO 42 280 .1750656714 281 STO 43 282 .1492479865 283 STO 44 284 .2718099231 285 STO 45 286 -.1707794051 287 STO 46 288 .3003551977 289 STO 47 290 -.010142544 291 STO 48 292 -1.191778158 293 STO 49 294 .0363913935 295 STO 50 296 -.0468331609 297 STO 51 298 .0324947663 299 STO 52 300 1.164052 301 REGMOVE 302 END

( 1929 bytes / SIZE 216 )

 STACK INPUT OUTPUT X / /

-Execution time = 55 seconds.

Note:

-There are of course less constants than with a 17-stage Runge-Kutta method,
but there are much more non-zero coefficients: so "C151" is longer & slower than "C169" !

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 R215 are to be initialized before executing "RK10" )

•  R01 = x0      •  R03 = h  = stepsize                 R05 & R10: temp            •  R65 to •  R215:  the 151 coefficients stored by "C151"
•  R02 = y0      •  R04 = N = number of steps    R11 to R26:  k1 to k16

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.026  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 16  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

( 115 bytes / SIZE 216 )

 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 "C151"

"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 = 10m01s---
X<>Y                       y(1) = 0.3678794410

-The exact result is 1/e =  0.3678794412...
-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 R215 are to be initialized before executing "2RK10" )

•  R01 = x0      •  R04 = h  = stepsize                 R06 & R10: temp            •  R65 to •  R215:  the 151 coefficients stored by "C151"
•  R02 = y0      •  R05 = N = number of steps    R11 to R26:  k1 to k16
•  R03 = z0                                                        R27 to R42:  l1 to l16
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 43  06 14.999  07 STO 44  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 15  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 43  51 +  52 STO 08  53 16.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 44  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 216 )

 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 "C151"

"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 = 17m01s---
RDN                      y(1) =  0.3678794411
RDN                      z(1) = -0.7357588821

-The exact results are  y(1) = 1/e =  0.3678794412... & z(1) = -2/e = -0.7357588823...
-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 R215 are to be initialized before executing "3RK10" )

•  R01 = x0      •  R05 = h  = stepsize                 R07 & R10: temp            •  R65 to •  R215:  the 151 coefficients stored by "C151"
•  R02 = y0      •  R06 = N = number of steps    R11 to R26:  k1 to k16         R59 to R61: temp
•  R03 = z0                                                        R27 to R42:  l1 to l16
•  R04 = t0                                                        R43 to R58:  m1 to m16
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 59  04 GTO 01  05 LBL 00  06 RCL IND 07  07 RCL IND 10  08 *  09 ST+ 61  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 61  31 RCL 02  32 +  33 RTN  34 LBL 01  35 15  36 STO 60  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 43  51 RDN  52 STO 27  53 X<>Y  54 STO 11  55 LBL 02  56 RCL 07  57 FRC  58 11.001  59 +  60 STO 07  61 16.9  62 +  63 STO 08   64 LASTX  65 INT  66 +  67 STO 09          68 CLST  69 STO 64 70 XEQ 00  71 STO 61  72 CLX  73 RCL 05  74 RCL IND 10  75 *  76 RCL 01  77 +  78 RCL 61  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 60 92 GTO 02  93 14.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 59 108 GTO 01 109 RCL 04  110 RCL 03        111 RCL 02 112 RCL 01 113 END

( 189 bytes / SIZE 216 )

 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 "C151"

"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 = 25m11s---
RDN     y(1) =  0.258207906
RDN     z(1) =  1.157623979
RDN     t(1) =  0.842178313

-The exact values are:

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

-Press R/S to continue with the same h and N or modify   R05-R06
-Synthetic registers M , N , O , P , Q  and registers R61-R62-R63-R64 may also be used in the subroutine to calculate f , g , h
( though register R61 is also used in "3RK10" )

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

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

•  R140+20n = a16,1  ..............................    •  R154+20n = a16,15       •  R155+20n = b16

•  R156+20n = c1 ........................................................................       •  R171+20n = c16

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.017  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 172+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

-Then we have to store the 151 constants into R21+20n to R171+20n

SIZE 232   XEQ "C151"  and since  n = 3

65.081151   REGMOVE

-We store the name of the subroutine in R00    "T"  ASTO 00
-The numbers of equations in R01                     3     STO 01
-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 = 40m12s---
RDN     y(1) =  0.258207906
RDN     z(1) =  1.157623981
RDN     t(1) =  0.842178312

Notes:

-With a real HP41, this program is very slow, but with free42, the results are returned almost instantaneously !
-Since SIZE = 172+20n, "NRK10" can solve at most a system of n = 7 differential equations with SIZE 312.
-But with free42, n may reach 41 with SIZE 992 ( 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]  David K. Zhang - "Discovering New Runge-Kutta Methods Using Unstructured Numerical Search"