Runge-Kutta-Nystrom Formula of order 10 for the HP-41
Overview
0°) 116 Constants
1°) 1 Second Order Conservative Equation
2°) System of 2 Second Order Conservative
Equations
3°) System of 3 Second Order Conservative
Equations
4°) System of n Second Order Conservative
Equations
-These routines apply a 13-stage 10th-order Runge-Kutta-Nystrom formula to solve ODEs of the form: y" = f(x,y)
-All s-stage explicit Runge-Kutta-Nystrom formulae without error estimate are determined by ( s2 + 5.s - 2 ) / 2 coefficients ai,j ; bi ; ci
-Here, s =13 , so 116 constants
k1 = h.f(x;y)
k2 = h.f(x+b2.h;y+b2.h.y'+h
a2,1.k1)
k3 = h.f(x+b3.h;y+b3.h.y'+h(a3,1.k1+a3,2.k2)
)
.................................................................................
k13 = h.f(x+b13.h;y+b13.h.y'+h (a13,1.k1+a13,2.k2+ .......... +a13,12.k12) )
then, y(x+h) = y(x) + h y'(x) + h [ c1.k1+
................ + c13.k13 ]
and y'(x+h) = y'(x) + [ c'1.k1+
................ + c'13.k13 ]
0°) 116 Constants
-All the programs listed in paragraphs 1-2-3 use 116 constants.
-"C116" stores these numbers into registers R65 thru R180
-Registers R07 to R64 are also used.
-Allways execute "C116" before executing the other programs... unless
you store your own 116 constants into R65 to R180 !
• R65 = a2,1
• R66 = b2
• R67 = a3,1
• R68 = a3,2
• R69 = b3
......................................................................... ...........
• R142= a13,1 ........................................ • R153 = a13,12 • R154 = b13
• R155 = c1 .........................................................................
• R167 = c13
• R168 = c'1 ........................................................................
• R180 = c'13
| 01 LBL "C116"
02 7.068 03 CLRGX 04 .0507388591 05 CHS 06 STO 07 07 .852583158 08 CHS 09 STO 09 10 .2925628199 11 STO 10 12 .4263130454 13 CHS 14 STO 11 15 .3084512679 16 STO 12 17 .8206806306 18 STO 13 19 .2635320056 20 STO 14 21 .0380029595 22 CHS 23 STO 15 24 .05083695 25 STO 16 26 .8584004334 27 STO 17 28 .855063419 29 CHS 30 STO 18 31 4.417296783 32 CHS 33 STO 20 |
34 1.462041474
35 STO 21 36 1.580060367 37 STO 22 38 2.150047309 39 CHS 40 STO 23 41 5.219288295 42 STO 24 43 .701222246 44 CHS 45 STO 25 46 .4067428472 47 STO 26 48 .109950164 49 CHS 50 STO 27 51 .0254989511 52 STO 28 53 .9592205307 54 STO 29 55 3.612420577 56 STO 30 57 19.61437642 58 STO 32 59 6.554095453 60 CHS 61 STO 33 62 5.363477518 63 CHS 64 STO 34 65 8.954920063 66 STO 35 |
67 22.11199958
68 CHS 69 STO 36 70 3.066641833 71 STO 37 72 1.297438034 73 CHS 74 STO 38 75 .5982268483 76 STO 39 77 .0241070789 78 CHS 79 STO 40 80 .0045319136 81 STO 41 82 1 83 STO 42 84 .0114450454 85 STO 43 86 .1518222881 87 STO 48 88 .0938333833 89 STO 49 90 .131388714 91 STO 50 92 .042452794 93 STO 51 94 .0466143591 95 STO 52 96 .0197393408 97 STO 53 98 .0027040753 99 STO 54 |
100 .0114450454
101 STO 56 102 .1896345577 103 STO 61 104 .0990150484 105 STO 62 106 .2237772082 107 STO 63 108 .1046681151 109 STO 64 110 .1535817253 111 STO 65 112 .1394025506 113 STO 66 114 .0663097234 115 STO 67 116 .0121660259 117 STO 68 118 7.119062 119 REGMOVE 120 7.06 121 CLRGX 122 .0006255035 123 STO 07 124 .0353695786 125 STO 08 126 .0008340047 127 STO 09 128 .0016680095 129 STO 10 130 .0707391571 131 STO 11 132 .0143775922 |
133 STO 12
134 .0255203896 135 CHS 136 STO 13 137 .0299554191 138 STO 14 139 .1939722747 140 STO 15 141 .005660601 142 STO 16 143 .0224381472 144 STO 18 145 .009619895 146 STO 19 147 .2746584906 148 STO 20 149 .0042497841 150 STO 21 151 .0137373896 152 STO 23 153 .0021090824 154 STO 24 155 .0002169817 156 CHS 157 STO 25 158 .1993954583 159 STO 26 160 .0008446992 161 STO 27 162 .0007628617 163 STO 29 164 .004527489 165 CHS |
166 STO 30
167 .0000947045 168 CHS 169 STO 31 170 .0043839568 171 STO 32 172 .0523320971 173 STO 33 174 .0402571887 175 STO 34 176 .2813255857 177 STO 36 178 .0990869733 179 CHS 180 STO 37 181 .0315586786 182 STO 38 183 .0762630093 184 STO 39 185 .2450911018 186 CHS 187 STO 40 188 .4128592672 189 STO 41 190 .5185220615 191 CHS 192 STO 42 193 4.022049675 194 CHS 195 STO 44 196 1.333954239 197 STO 45 198 .3628013762 |
199 CHS
200 STO 46 201 .446629326 202 CHS 203 STO 47 204 4.10933539 205 STO 48 206 .0833718601 207 STO 49 208 .5944056701 209 STO 50 210 .4552651504 211 STO 51 212 3.44102443 213 STO 53 214 1.142030463 215 CHS 216 STO 54 217 .3312250453 218 STO 55 219 .5251081451 220 STO 56 221 3.400720497 222 CHS 223 STO 57 224 .0149590255 225 STO 58 226 .017714834 227 STO 59 228 .6964849889 229 STO 60 230 7.065054 231 REGMOVE 232 END |
( 1302 bytes / SIZE 181 )
| STACK | INPUT | OUTPUT |
| X | / | / |
---Execution
time = 37s---
1°) 1 Differential Equation
-We want to solve y" = f(x,y)
with the initial conditions y(x0)
= y0 , y'(x0) = y'0
Data Registers: • R00: f name ( Registers R00 to R05 & R65 to R180 are to be initialized before executing "RKN10" )
• R01 = x0 • R04 =
h = stepsize
R06 & R09: temp
• R65 to • R180: the 116 coefficients stored by "C116"
• R02 = y0 • R05 =
N = number of steps R11 to R23: k1 to
k13
• R03 = y'0
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
| 01 LBL "RKN10"
02 RCL 05 03 STO 07 04 11.001 05 STO 06 06 GTO 01 07 LBL 00 08 RCL IND 08 09 RCL IND 09 10 * 11 + 12 ISG 09 13 CLX 14 ISG 08 15 GTO 00 16 RCL 03 |
17 +
18 RTN 19 LBL 01 20 .01 21 STO 08 22 65.154 23 STO 09 24 RCL 02 25 RCL 01 26 XEQ IND 00 27 RCL 04 28 * 29 STO 11 30 LBL 02 31 RCL 08 32 FRC |
33 RCL 06
34 + 35 STO 08 36 CLX 37 LBL 03 38 RCL IND 08 39 RCL IND 09 40 * 41 + 42 ISG 09 43 ISG 08 44 GTO 03 45 RCL 03 46 RCL IND 09 47 * 48 + |
49 RCL 04
50 * 51 RCL 02 52 + 53 RCL 04 54 RCL IND 09 55 * 56 RCL 01 57 + 58 XEQ IND 00 59 RCL 04 60 * 61 STO IND 08 62 ISG 09 63 GTO 02 64 11.999 65 ST- 08 |
66 CLX
67 XEQ 00 68 RCL 04 69 ST+ 01 70 * 71 ST+ 02 72 13 73 ST- 08 74 CLX 75 XEQ 00 76 STO 03 77 DSE 07 78 GTO 01 79 RCL 03 80 RCL 02 81 RCL 01 82 END |
( 138 bytes / SIZE 181 )
| STACK | INPUTS | OUTPUTS |
| Z | / | y'(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 ,
y'(0) = 0 , y" = - y ( x2 + y2
) 1/2
>>> Calculate y(1) & y'(1)
-Load this short routine:
| 01 LBL "T"
02 X^2 03 X<>Y 04 STO Z 05 X^2 06 + 07 SQRT 08 * 09 CHS 10 RTN |
"T" ASTO 00
0 STO 01 STO 03 1 STO 02
-With h = 0.1 and N = 10 , 0.1 STO 04 10 STO 05
XEQ "RKS1" >>>> x
= 1
---Execution time = 8m30s---
RDN y(1) = 0.536630617
RDN y'(1) = -0.860171927
Note:
-Simply press R/S to continue with the same h & N or modify registers
R04 & R05 if need be.
2°) System of 2 Differential Equations
-"2RKN10" solves the system:
y" = f(x,y,z)
y(x0) = y0 y'(x0)
= y'0
z" = g(x,y,z)
z(x0) = z0 z'(x0)
= z'0
Data Registers: • R00: f name ( Registers R00 to R07 & R65 to R180 are to be initialized before executing "2RKN10" )
• R01 = x0 • R04 =
h = stepsize
• R06 = y'0
R11 to R23: k1 to k13
R08 to R10 & R37: temp
• R02 = y0 • R05 =
N = number of steps • R07 = z'0
R24 to R36: l1 to l13
• R03 = z0
• R65 to • R180: the 116 coefficients stored by "C116"
Flags: /
Subroutine: a program which calculates 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 "2RKN10"
02 RCL 05 03 STO 37 04 11.001 05 STO 38 06 GTO 01 07 LBL 00 08 X<>Y 09 RCL IND 10 10 RCL IND 09 11 * 12 + 13 X<>Y 14 RCL IND 08 15 RCL IND 09 16 * 17 + 18 ISG 08 19 ISG 09 20 TEXT0 21 ISG 10 22 GTO 00 |
23 RTN
24 LBL 01 25 .01 26 STO 10 27 65.154 28 STO 09 29 RCL 03 30 RCL 02 31 RCL 01 32 XEQ IND 00 33 RCL 04 34 ST* Z 35 * 36 STO 24 37 X<>Y 38 STO 11 39 LBL 02 40 RCL 10 41 FRC 42 RCL 38 43 + 44 STO 10 |
45 24.9
46 STO 08 47 CLST 48 XEQ 00 49 RCL 07 50 RCL IND 09 51 * 52 + 53 RCL 04 54 * 55 RCL 03 56 + 57 X<>Y 58 RCL 06 59 RCL IND 09 60 * 61 + 62 RCL 04 63 * 64 RCL 02 65 + 66 RCL 04 |
67 RCL IND 09
68 * 69 RCL 01 70 + 71 XEQ IND 00 72 RCL 04 73 ST* Z 74 * 75 STO IND 08 76 X<>Y 77 STO IND 10 78 ISG 09 79 GTO 02 80 11.999 81 ST- 10 82 ST- 08 83 CLST 84 XEQ 00 85 RCL 07 86 + 87 RCL 04 88 * |
89 ST+ 03
90 CLX 91 RCL 06 92 + 93 RCL 04 94 ST+ 01 95 * 96 ST+ 02 97 13 98 ST- 10 99 ST- 08 100 CLST 101 XEQ 00 102 ST+ 07 103 X<>Y 104 ST+ 06 105 DSE 37 106 GTO 01 107 RCL 03 108 RCL 02 109 RCL 01 110 END |
( 180 bytes / SIZE 181 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z
y(0) = 2 y'(0) = 1
Evaluate y(1) , z(1) , y'(1) , z'(1) with
h = 0.1 & N = 10
d2z/dx2 = x.(
y + z ) z(0) = 1
z'(0) = 1
| 01 LBL "T"
02 RDN 03 STO Z 04 X<>Y 05 CHS 06 ST* Z 07 - 08 R^ 09 * 10 RTN |
"T" ASTO 00
0 STO 01
1 STO 06
2 STO 02
1 STO 07
1 STO 03
0.1 STO 04
10 STO 05
XEQ "2RKN10" >>>> x
= 1
---Execution time = 13m58s---
RDN y(1) = 1.531356645 and
R06 = y'(1) = -2.312840138
RDN z(1) = 2.620254282
R07 = z'(1) = 2.941748401
Note:
-Simply press R/S to continue with the same h & N or modify registers
R04 & R05 if need be.
3°) System of 3 Differential Equations
-"3RKN10" solves the system:
y" = f(x,y,z,t)
y(x0) = y0 y'(x0)
= y'0
z" = g(x,y,z,t)
z(x0) = z0 z'(x0)
= z'0
t" = h(x,y,z,t)
t(x0) = t0 t'(x0)
= t'0
Data Registers: • R00: f name ( Registers R00 to R09 & R65 to R180 are to be initialized before executing "3RKN10" )
• R01 = x0 • R05 =
h = stepsize
• R07 = y'0
R16 to R28: k1 to k13
R10 to R15 & R55: temp
• R02 = y0 • R06 =
N = number of steps • R08 = z'0
R29 to R41: l1 to l13
• R03 = z0
• R09 = t'0
R42 to R54: m1 to m13
• R04 = t0
• R65 to • R180: the 116 coefficients stored by "C116"
Flags: /
Subroutine: a program which calculates f(x;y;z;t)
in Z-register , g(x;y;z;t) in Y-register and h(x;y;z;t) in X-register
assuming x is in X-register , y is in Y-register , z is in Z-register and
t is in T-register upon entry.
| 01 LBL "3RKN10"
02 RCL 06 03 STO 55 04 16.001 05 STO 10 06 GTO 01 07 LBL 00 08 RCL IND 11 09 RCL IND 14 10 * 11 ST+ 15 12 X<> Z 13 RCL IND 12 14 RCL IND 14 15 * 16 + 17 X<>Y 18 RCL IND 13 19 RCL IND 14 20 * 21 + 22 ISG 14 23 TEXT0 24 ISG 13 25 ISG 12 26 ISG 11 27 GTO 00 28 RTN 29 LBL 01 30 .015 31 STO 11 |
32 65.154
33 STO 14 34 RCL 04 35 RCL 03 36 RCL 02 37 RCL 01 38 XEQ IND 00 39 RCL 05 40 ST* T 41 ST* Z 42 * 43 STO 42 44 RDN 45 STO 29 46 X<>Y 47 STO 16 48 LBL 02 49 RCL 11 50 FRC 51 RCL 10 52 + 53 STO 11 54 13.9 55 + 56 STO 12 57 LASTX 58 INT 59 + 60 STO 13 61 CLST 62 STO 15 |
63 XEQ 00
64 RCL 09 65 RCL IND 14 66 * 67 + 68 RCL 05 69 * 70 RCL 04 71 + 72 X<>Y 73 RCL 08 74 RCL IND 14 75 * 76 + 77 RCL 05 78 * 79 RCL 03 80 + 81 RCL 07 82 RCL IND 14 83 * 84 RCL 15 85 + 86 RCL 05 87 * 88 RCL 02 89 + 90 STO 15 91 CLX 92 RCL 05 93 RCL IND 14 |
94 *
95 RCL 01 96 + 97 RCL 15 98 X<>Y 99 XEQ IND 00 100 RCL 05 101 ST* T 102 ST* Z 103 * 104 STO IND 13 105 RDN 106 STO IND 12 107 X<>Y 108 STO IND 11 109 ISG 14 110 GTO 02 111 11.999 112 ST- 11 113 ST- 12 114 ST- 13 115 CLST 116 STO 15 117 XEQ 00 118 RCL 09 119 + 120 RCL 05 121 * 122 ST+ 04 123 CLX 124 RCL 08 |
125 +
126 RCL 05 127 * 128 ST+ 03 129 RCL 15 130 RCL 07 131 + 132 RCL 05 133 ST+ 01 134 * 135 ST+ 02 136 13 137 ST- 11 138 ST- 12 139 ST- 13 140 CLST 141 STO 15 142 XEQ 00 143 ST+ 09 144 X<>Y 145 ST+ 08 146 RCL 15 147 ST+ 07 148 DSE 55 149 GTO 01 150 RCL 04 151 RCL 03 152 RCL 02 153 RCL 01 154 END |
( 242 bytes / SIZE 181 )
| STACK | INPUT | OUTPUTS |
| T | / | t(x0+N.h) |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z.t
y(0) = 1 y'(0) = 1
Evaluate y(1) , z(1) , t(1) , y'(1) , z'(1) , t'(1)
with h = 0.1 & N = 10
d2z/dx2 = x.( y + z
- t ) z(0) = 1
z'(0) = 1
d2u/dx2 = x.y - z.t
t(0) = 2 t'(0) = 1
| 01 LBL "T"
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 |
"T" ASTO 00
0 STO 01
1 STO 02 STO 03 STO 07
STO 08 STO 09
2 STO 04
0.1 STO 05
10 STO 06 XEQ "3RKN10"
>>>> x = 1
---Execution time = 20m16s---
RDN y(1) = 0.439524100
y'(1) = -2.101122880 = R07
RDN z(1) = 2.070940654
z'(1) = 1.269596951 = R08
RDN t(1) = 1.744524962
t'(1) = -1.704234756 = R09
Notes:
-The errors are about E-9 ( no more than roundoff-errors )
-You can use registers R56 to R64 to write the subroutine.
-Simply press R/S to continue with the same h & N or modify registers
R05 & R06 if need be.
4°) System of n Differential Equations
-We have to solve the system:
y"1 = f1(x,y1,
... ,yn)
y"2 = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
& y'i(x0)
i = 1 , ... , n
y"n = fn(x,y1,
... ,yn)
Data Registers: • R00 = subroutine name
( Registers R00 thru R03 , R20 thru R20+2n and R21+17n thru R136+17n are to be initialized before executing "NRKN10" )
• R01 = n = number of equations = number of functions
R04 thru R11 & R19: temp R12 to R18: unused
• R02 = h = stepsize
• R03 = N = number of steps
• R20 = x0
• R21 = y1(x0)
• R21+n = y'1(x0)
R21+2n thru R20+17n: temp
• R22 = y2(x0)
• R22+n = y'2(x0)
............... ( initial
values ) .......................
• R20+n = yn(x0)
• R20+2n = y'n(x0)
• R21+17n = a2,1
• R22+17n = b2
• R23+17n = a3,1
• R24+17n = a3,2
• R25+17n = b3
......................................................................... ...........
• R98+17n = a13,1 ................................. • R109+17n = a13,12 • R110+17n = b17
• R111+17n = c1 ........................................................................
• R123+17n = c17
• R124+17n = c'1 .......................................................................
• R136+17n = c'17
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+2n, ... , R20+3n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
-Line 162 is a three-byte GTO 01
| 01 LBL "NRKN10"
02 RCL 03 03 STO 11 04 RCL 01 05 .1 06 % 07 + 08 STO 05 09 FRC 10 21.02 11 + 12 STO 04 13 RCL 05 14 + 15 STO 10 16 LASTX 17 ST+ X 18 + 19 STO 06 20 GTO 01 21 LBL 02 22 XEQ IND 00 23 RCL 04 24 RCL 05 25 ST+ X 26 + 27 RCL 07 28 LBL 00 |
29 RCL IND Y
30 RCL 02 31 * 32 STO IND Y 33 ISG Y 34 RDN 35 ISG Y 36 GTO 00 37 RTN 38 LBL 03 39 RCL IND Y 40 STO IND Y 41 ISG Y 42 RDN 43 ISG Y 44 GTO 03 45 RTN 46 LBL 07 47 RCL IND 07 48 RCL IND 08 49 * 50 RCL 02 51 * 52 ST+ IND 04 53 ISG 04 54 CLX 55 ISG 07 56 GTO 07 |
57 RCL 01
58 ST- 04 59 RCL 05 60 FRC 61 ST+ 07 62 ISG 08 63 GTO 07 64 RTN 65 LBL 01 66 RCL 05 67 RCL 06 68 + 69 STO 07 70 XEQ 02 71 RCL 04 72 RCL 06 73 XEQ 03 74 RCL 01 75 17 76 * 77 21 78 + 79 .1 80 % 81 + 82 STO 08 83 3.014 84 STO 09 |
85 RCL 20
86 STO 19 87 LBL 04 88 XEQ 07 89 RCL 02 90 RCL IND 08 91 * 92 ST+ 20 93 SIGN 94 LBL 06 95 RCL IND 10 96 LASTX 97 * 98 ST+ IND 04 99 ISG 04 100 CLX 101 ISG 10 102 GTO 06 103 RCL 01 104 ST- 04 105 ST- 10 106 XEQ 02 107 RCL 06 108 RCL 04 109 XEQ 03 110 RCL 19 111 STO 20 112 RCL 09 |
113 INT
114 E3 115 / 116 ISG X 117 ST+ 08 118 RCL 05 119 RCL 06 120 + 121 STO 07 122 ISG 09 123 GTO 04 124 RCL 02 125 ST+ 20 126 LBL 05 127 RCL IND 10 128 RCL 02 129 * 130 ST+ IND 04 131 ISG 10 132 CLX 133 ISG 04 134 GTO 05 135 RCL 01 136 ST- 04 137 ST- 10 138 XEQ 07 139 .013 140 ST+ 08 |
141 RCL 05
142 RCL 06 143 + 144 STO 07 145 LBL 08 146 RCL IND 07 147 RCL IND 08 148 * 149 ST+ IND 10 150 ISG 10 151 CLX 152 ISG 07 153 GTO 08 154 RCL 01 155 ST- 10 156 RCL 05 157 FRC 158 ST+ 07 159 ISG 08 160 GTO 08 161 DSE 11 162 GTO 01 163 RCL 23 164 RCL 22 165 RCL 21 166 RCL 20 167 END |
( 273 bytes / SIZE 137+17.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:
d2y1/dx2
= - y1y2y3
d2y2/dx2
= x ( y1+y2-y3 )
d2y3/dx2
= xy1- y2y3
with the initial values: y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2 & y'1(0) = 1 ; y'2(0) = 1 ; y'3(0) = 1
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 R27 ; R28 ; R29
-For instance, the following subroutine:
| 01 LBL "T"
02 RCL 21 03 RCL 22 04 RCL 23 05 * 06 * 07 CHS 08 STO 27 09 RCL 21 10 RCL 22 11 + 12 RCL 23 13 - 14 RCL 20 15 * 16 STO 28 17 RCL 20 18 RCL 21 19 * 20 RCL 22 21 RCL 23 22 * 23 - 24 STO29 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 116 constants into R21+17n to R136+17n
-Execute for instance the following routine:
| 01 LBL "INIT"
02 XEQ "C116" 03 RCL 01 04 17 05 * 06 21.116 07 + 08 E3 09 / 10 65 11 + 12 REGMOVE 13 END |
-Then the initial values:
x = 0 STO 20 and
y1 = 1 STO 21
y'1 = 1 STO 24
y2 = 1 STO 22
y'2 = 1 STO 25
y3 = 2 STO 23
y'3 = 1 STO 26
-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 "NRKN10" >>>> x = 1 ---Execution time = 4s--- with V41 in turbo mode
RDN y1(1) = 0.439524097 = R21
y'1(1) = -2.101122877 = R24
RDN y2(1) = 2.070940658 = R22
and y'2(1) =
1.269596947 = R25
RDN y3(1) = 1.744524962 = R23
y'3(1) = -1.704234764 = R26
Notes:
-With a real HP41, execution time is about 40 minutes, but with free42,
the results are returned almost instantaneously !
-Since SIZE = 137+17n, "NRKN10" can solve at most a system of n = 10
differential equations with SIZE 307.
-But with free42, n may reach 50 with SIZE 990 ( beyond register 999,
the control numbers would not work )
Reference:
[1] Philip Sharp - http://www.math.auckland.ac.nz/~sharp/rkn1.dat
-The programs listed in this page employ formulae without error-estimate,
so you'll have to calculate again
the solutions with a smaller stepsize - say h/2 - to get
an idea of the precision.
-But in reference [1], you will find the coefficients of
pairs
of embedded Runge-Kutta-Nystrom formulas
of various orders to control the stepsize more directly.