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.