Runge-Kutta-Nystrom10

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 ( s² + 5.s - 2 ) / 2 coefficients a_i,j ; b_i ; c_i

-Here, s =13 , so 116 constants

    k₁ = h.f(x;y)
    k₂ = h.f(x+b₂.h;y+b₂.h.y'+h a_2,1.k₁)
    k₃ = h.f(x+b₃.h;y+b₃.h.y'+h(a_3,1.k₁+a_3,2.k₂) )
    .................................................................................

k₁₃ = h.f(x+b₁₃.h;y+b₁₃.h.y'+h (a_13,1.k₁+a_13,2.k₂+ .......... +a_13,12.k₁₂) )

then, y(x+h) = y(x) + h y'(x) + h [ c₁.k₁+ ................ + c₁₃.k₁₃ ]
and y'(x+h) = y'(x) + [ c'₁.k₁+ ................ + c'₁₃.k₁₃ ]

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 = a_2,1 • R66 = b₂
• R67 = a_3,1 • R68 = a_3,2 • R69 = b₃

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

• R142= a_13,1 ........................................ • R153 = a_13,12 • R154 = b₁₃

• R155 = c₁ ......................................................................... • R167 = c₁₃
• R168 = c'₁ ........................................................................ • R180 = c'₁₃

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(x₀) = y₀ , y'(x₀) = y'₀

Data Registers: • R00: f name ( Registers R00 to R05 & R65 to R180 are to be initialized before executing "RKN10" )

                                     • R01 = x₀      • R04 = h = stepsize                 R06 & R09: temp            • R65 to • R180: the 116 coefficients stored by "C116"
                                     • R02 = y₀      • R05 = N = number of steps    R11 to R23: k₁ to k₁₃
                                     • R03 = y'₀

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'(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example: y(0) = 1 , y'(0) = 0 , y" = - y ( x² + y² )^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(x₀) = y₀ y'(x₀) = y'₀
z" = g(x,y,z) z(x₀) = z₀ z'(x₀) = z'₀

Data Registers: • R00: f name ( Registers R00 to R07 & R65 to R180 are to be initialized before executing "2RKN10" )

                                     • R01 = x₀      • R04 = h = stepsize                 • R06 = y'₀             R11 to R23: k₁ to k₁₃     R08 to R10 & R37: temp
                                     • R02 = y₀      • R05 = N = number of steps    • R07 = z'₀             R24 to R36: l₁ to l₁₃
                                     • R03 = z₀                      • 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(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example:

d²y/dx² = -y.z y(0) = 2 y'(0) = 1 Evaluate y(1) , z(1) , y'(1) , z'(1) with h = 0.1 & N = 10
d²z/dx² = 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(x₀) = y₀     y'(x₀) = y'₀
     z" = g(x,y,z,t)         z(x₀) = z₀     z'(x₀) = z'₀
     t" = h(x,y,z,t)          t(x₀) = t₀      t'(x₀) = t'₀

Data Registers: • R00: f name ( Registers R00 to R09 & R65 to R180 are to be initialized before executing "3RKN10" )

                                     • R01 = x₀      • R05 = h = stepsize                 • R07 = y'₀             R16 to R28: k₁ to k₁₃     R10 to R15 & R55: temp
                                     • R02 = y₀      • R06 = N = number of steps    • R08 = z'₀             R29 to R41: l₁ to l₁₃
                                     • R03 = z₀                                                        • R09 = t'₀             R42 to R54: m₁ to m₁₃
                                     • R04 = t₀                                     • 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(x₀+N.h)
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example:

     d²y/dx² = -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
     d²z/dx² = x.( y + z - t )         z(0) = 1        z'(0) = 1
     d²u/dx² = 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"₁ = f₁(x,y₁, ... ,y_n)
     y"₂ = f₂(x,y₁, ... ,y_n)
             ..................                                        with the initial values:     y_i(x₀) & y'_i(x₀)      i = 1 , ... , n

y"_n = f_n(x,y₁, ... ,y_n)

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 = x₀
                                      • R21 = y₁(x₀)                                   • R21+n = y'₁(x₀)                                   R21+2n thru R20+17n: temp
                                      • R22 = y₂(x₀)                                   • R22+n = y'₂(x₀)
                                          ...............         ( initial values )           .......................

• R20+n = y_n(x₀) • R20+2n = y'_n(x₀)

• R21+17n = a_2,1 • R22+17n = b₂
• R23+17n = a_3,1 • R24+17n = a_3,2 • R25+17n = b₃

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

• R98+17n = a_13,1 ................................. • R109+17n = a_13,12 • R110+17n = b₁₇

• R111+17n = c₁ ........................................................................ • R123+17n = c₁₇
• R124+17n = c'₁ ....................................................................... • R136+17n = c'₁₇

Flags: /
Subroutine: A program which calculates and stores f₁ ( x, y₁ , ... , y_n) , ... , f_n ( x, y₁ , ... , y_n) in registers R21+2n, ... , R20+3n respectively
with x, y₁, ... , y_n 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	/	y₃(x₀+N.h)
Z	/	y₂(x₀+N.h)
Y	/	y₁(x₀+N.h)
X	/	x₀+N.h

Example: Let's solve again the system:

      d²y₁/dx² = - y₁y₂y₃
      d²y₂/dx² = x ( y₁+y₂-y₃ )
      d²y₃/dx² = xy₁- y₂y₃

with the initial values: y₁(0) = 1 ; y₂(0) = 1 ; y₃(0) = 2 & y'₁(0) = 1 ; y'₂(0) = 1 ; y'₃(0) = 1

x ; y₁ ; y₂ ; y₃ 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

     y₁ = 1   STO 21     y'₁ = 1 STO 24
     y₂ = 1   STO 22     y'₂ = 1 STO 25
     y₃ = 2   STO 23     y'₃ = 1 STO 26

-Suppose we want to know y₁(1) , y₂(1) , y₃(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     y₁(1) = 0.439524097 = R21                     y'₁(1) = -2.101122877 = R24
                                                 RDN     y₂(1) = 2.070940658 = R22       and         y'₂(1) = 1.269596947 = R25
                                                 RDN     y₃(1) = 1.744524962 = R23                      y'₃(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.

hp41programs

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