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"