Runge-Kutta-Nystrom

Runge-Kutta-Nystrom Methods for the HP-41

Overview

1°) A 4th-Order Method

   a) Second-order differential equations
   b) Systems of 2 second-order differential equations
   c) Systems of 3 second-order differential equations

2°) A Runge-Kutta-Nystrom Method with order 6

a) Second order differential equations
b) Systems of 2 second-order differential equations

3°) S-Stage Explicit Runge-Kutta-Nystrom formulae

a) Second order differential equations
b) Systems of 2 second-order differential equations

-A differential equation of the type y" = f(x,y) may always be replaced by a system of 2 first order equations:

y' = z
z' = f(x,y)

-However, the same order of accuracy may be obtained with less evaluations of the function if we use formulae specially devised for these problems.

    3 evaluations instead of 4 for a 4th-order method
    5 evaluations instead of 7 for a 6th-order method
   13 evaluations instead of 17 for a 10th-order method ... etc ...

Notes:

-The programs of paragraph 1 are already listed elsewhere in this website
( cf "Explicit Runge-Kutta Methods" & "Systems of Differential Equations" )
-They are just included here for the sake of consistency.

1°) A 4th-Order Method

a) Second-order differential Equations

-Here, we have to solve: y" = f(x,y) with the initial values: y(x₀) = y₀ y'(x₀) = y'₀

-This program uses the 4th-order Runge-Kutta formula:

y(x+h) = y(x) + h ( y'(x) + k₁/6 + 2k₂/3 )
y'(x+h) = y'(x) + k₁/6 + 2k₂/3 + k₃/6 where k₁ = h.f (x,y) ; k₂ = h.f(x+h/2,y+h.y'/2+h.k₁/8) ; k₃ = h.f(x+h,y+h.y'+h.k₂/2)

-Note that only 3 evaluations of the function are needed ( for each step ).

Data Registers: • R00: subroutine name ( Registers R00 thru R05 are to be initialized before executing "RKS1" )

                               •   R01 = x₀      •   R04 = h/2 = half of the stepsize        R06 to R08: temp
                               •   R02 = y₀      •   R05 = N = number of steps
                               •   R03 = y'₀
Flags: /
Subroutine: a program which calculates f(x,y) in X-register, assuming x and y are in registers X and Y upon entry.

01 LBL "RKS1"
02 RCL 05
03 STO 08
04 LBL 01
05 RCL 02
06 RCL 01
07 XEQ IND 00
08 RCL 04
09 ST+ 01
10 *
11 STO 06
12 2

13 /
14 RCL 03
15 +
16 RCL 04
17 *
18 RCL 02
19 +
20 RCL 01
21 XEQ IND 00
22 RCL 04
23 ST+ 01
24 *

25 STO 07
26 RCL 03
27 +
28 RCL 04
29 ST+ X
30 *
31 RCL 02
32 +
33 RCL 01
34 XEQ IND 00
35 RCL 04
36 *

37 RCL 06
38 +
39 RCL 07
40 ST+ X
41 ST+ 06
42 ST+ X
43 +
44 3
45 ST/ 06
46 /
47 X<> 03
48 ST+ 03

49 RCL 06
50 +
51 RCL 04
52 ST+ X
53 *
54 ST+ 02
55 DSE 08
56 GTO 01
57 RCL 03
58 RCL 02
59 RCL 01
60 END

( 85 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Z	/	y'(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example: y(0) = 1 , y'(0) = 1 , 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.05 STO 04 10 STO 05

XEQ "RKS1" >>>>      x   = 1                                            ---Execution time = 29s---
                          RDN   y(1) = 0.536630911
                          RDN   y'(1) = -0.860172085

-With h = 0.02 & N = 50 , it yields

y(1) = 0.536630617
y'(1) = -0.860171928

Note:

-Simply press R/S to continue with the same h &N

b) Systems of 2 Second-order differential Equations

-"RKS2" 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: subroutine name ( Registers R00 thru R07 are to be initialized before executing "RKS2" )

                               •   R01 = x₀      •   R04 = h/2 = half of the stepsize       •   R06 = y'₀          R08 to R12: temp
                               •   R02 = y₀      •   R05 = N = number of steps            •   R07 = z'₀
                               •   R03 = z₀
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.

-Line 96 is a three-byte GTO 01

01 LBL "RKS2"
02 RCL 05
03 STO 12
04 LBL 01
05 RCL 03
06 RCL 02
07 RCL 01
08 XEQ IND 00
09 RCL 04
10 ST+ 01
11 ST* Z
12 *
13 STO 09
14 2
15 /
16 RCL 07
17 +
18 RCL 04
19 *
20 RCL 03

21 +
22 X<>Y
23 STO 08
24 2
25 /
26 RCL 06
27 +
28 RCL 04
29 *
30 RCL 02
31 +
32 RCL 01
33 XEQ IND 00
34 RCL 04
35 ST+ 01
36 ST* Z
37 *
38 STO 11
39 RCL 07
40 +

41 X<>Y
42 STO 10
43 RCL 06
44 +
45 RCL 04
46 ST+ X
47 ST* Z
48 *
49 RCL 03
50 ST+ Z
51 CLX
52 RCL 02
53 +
54 RCL 01
55 XEQ IND 00
56 RCL 04
57 ST* Z
58 *
59 RCL 09
60 +

61 RCL 11
62 ST+ X
63 ST+ 09
64 ST+ X
65 +
66 3
67 ST/ 09
68 /
69 X<> 07
70 ST+ 07
71 RCL 09
72 +
73 X<>Y
74 RCL 08
75 +
76 RCL 10
77 ST+ X
78 ST+ 08
79 ST+ X
80 +

81 3
82 ST/ 08
83 /
84 X<> 06
85 ST+ 06
86 RCL 08
87 +
88 RCL 04
89 ST+ X
90 ST* Z
91 *
92 ST+ 02
93 X<>Y
94 ST+ 03
95 DSE 12
96 GTO 01
97 RCL 03
98 RCL 02
99 RCL 01
100 END

( 139 bytes / SIZE 013 )

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.05 STO 04
   10    STO 05

XEQ "RKS2" >>>>     x   = 1                                                                                                  ( in 41 seconds )
                         RDN   y(1) = 1.531358015     and   R06 = y'(1) = -2.312838895
                         RDN   z(1) = 2.620254480              R07 = z'(1) = 2.941751649

-With h = 0.05 it yields y(1) = 1.531356736 y'(1) = -2.312840085
z(1) = 2.620254295 z'(1) = 2.941748608

-Actually, the exact results - rounded to 9D - are

y(1) = 1.531356646 y'(1) = -2.312840137
z(1) = 2.620254281 z'(1) = 2.941748399

-To continue the calculations, simply press R/S ( after changing h/2 & N in registers R04 & R05 if needed )

c) Systems of 3 Second-order differential Equations

-"RKS3" solves the system:

     y" = f(x,y,z,u)         y(x₀) = y₀     y'(x₀) = y'₀
     z" = g(x,y,z,u)         z(x₀) = z₀     z'(x₀) = z'₀
     u" = h(x,y,z,u)        u(x₀) = u₀     u'(x₀) = u'₀

Data Registers: • R00: subroutine name ( Registers R00 thru R09 are to be initialized before executing "RKS3" )

                               •   R01 = x₀      •   R05 = h/2 = half of the stepsize       •   R07 = y'₀
                               •   R02 = y₀      •   R06 = N = number of steps            •   R08 = z'₀
                               •   R03 = z₀                   R10 to R16: temp                  •   R09 = u'₀
                               •   R04 = u₀
Flags: /
Subroutine: a program which calculates f(x;y;z;u) in Z-register , g(x;y;z;u) in Y-register and h(x;y;z;u) in X-register
                       assuming x is in X-register , y is in Y-register , z is in Z-register and u is in T-register upon entry.

-Line 136 is a three-byte GTO 01

01 LBL "RKS3"
02 RCL 06
03 STO 16
04 LBL 01
05 RCL 04
06 RCL 03
07 RCL 02
08 RCL 01
09 XEQ IND 00
10 RCL 05
11 ST+ 01
12 ST* Z
13 ST* T
14 *
15 STO 12
16 2
17 /
18 RCL 09
19 +
20 RCL 05
21 *
22 RCL 04
23 +
24 X<>Y
25 STO 11
26 2
27 /
28 RCL 08
29 +

30 RCL 05
31 *
32 RCL 03
33 +
34 R^
35 STO 10
36 2
37 /
38 RCL 07
39 +
40 RCL 05
41 *
42 RCL 02
43 +
44 RCL 01
45 XEQ IND 00
46 RCL 05
47 ST+ 01
48 ST* Z
49 ST* T
50 *
51 STO 15
52 RCL 09
53 +
54 X<>Y
55 STO 14
56 RCL 08
57 +
58 R^

59 STO 13
60 RCL 07
61 +
62 RCL 05
63 ST+ X
64 ST* Z
65 ST* T
66 *
67 RCL 04
68 ST+ T
69 CLX
70 RCL 03
71 ST+ Z
72 CLX
73 RCL 02
74 +
75 RCL 01
76 XEQ IND 00
77 RCL 05
78 ST* Z
79 ST* T
80 *
81 RCL 12
82 +
83 RCL 15
84 ST+ X
85 ST+ 12
86 ST+ X
87 +

88 3
89 ST/ 12
90 /
91 X<> 09
92 ST+ 09
93 RCL 12
94 +
95 X<>Y
96 RCL 11
97 +
98 RCL 14
99 ST+ X
100 ST+ 11
101 ST+ X
102 +
103 3
104 ST/ 11
105 /
106 X<> 08
107 ST+ 08
108 RCL 11
109 +
110 R^
111 RCL 10
112 +
113 RCL 13
114 ST+ X
115 ST+ 10
116 ST+ X

117 +
118 3
119 ST/ 10
120 /
121 X<> 07
122 ST+ 07
123 RCL 10
124 +
125 RCL 05
126 ST+ X
127 ST* Z
128 ST* T
129 *
130 ST+ 02
131 X<>Y
132 ST+ 03
133 R^
134 ST+ 04
135 DSE 16
136 GTO 01
137 RCL 04
138 RCL 03
139 RCL 02
140 RCL 01
141 END

( 194 bytes / SIZE 017 )

STACK	INPUTS	OUTPUTS
T	/	u(x₀+N.h)
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example:

     d²y/dx² = -y.z.u                    y(0) = 1        y'(0) = 1         Evaluate y(1) , z(1) , u(1) ,   y'(1) , z'(1) , u'(1)      with h = 0.1 & N = 10
     d²z/dx² = x.( y + z - u )         z(0) = 1        z'(0) = 1
     d²u/dx² = x.y - z.u                 u(0) = 2       u'(0) = 1

LBL "T"       R^               RCL Z         X<> Z         LASTX            "T" ASTO 00
R^               ST* Y         R^                ST+ Y         *
CHS            ST+ L         ST+ L          ST* Z          X<>Y
STO L         CLX           ST* Y          X<> T         RTN

0     STO 01
1     STO 02 STO 03 STO 07 STO 08 STO 09
2     STO 04
0.05 STO 05
10    STO 06    XEQ "RKS3" >>>>     x   = 1                                                                                             ( in 65 seconds )
                                                 RDN   y(1) = 0.439528419           y'(1) = -2.101120400 = R07
                                                 RDN   z(1) = 2.070938499           z'(1) =   1.269599239 = R08
                                                 RDN   u(1) = 1.744522976           u'(1) = -1.704232092 = R09

-With h = 0.05 it yields

      y(1) = 0.439524393          y'(1) = -2.101122784
      z(1) = 2.070940521          z'(1) =   1.269597110
      u(1) = 1.744524843          u'(1) = -1.704234567

-Actually, the exact results rounded to 9D are:

      y(1) = 0.439524100          y'(1) = -2.101122880
      z(1) = 2.070940654          z'(1) =   1.269596950
      u(1) = 1.744524964          u'(1) = -1.704234756

-Press R/S to continue the calculations ( after changing h & N in registers R05 & R06 if needed )

2°) A Runge-Kutta-Nystrom Method with order 6

a) Second-order differential Equations

-Five stages are sufficient to get a 6-th order formula.
-The following one is due to Albrecht:

   f₀ = f(x₀,y₀)
   f₁ = f(x₀+ h/4 , y₀+ (h/4) y'₀+ (h²/32) f₀)
   f₂ = f(x₀+ h/2 , y₀+ (h/2) y'₀+ h² (-f₀/24 + f₁/6 )
   f₃ = f(x₀+ 3h/4 , y₀+ (3h/4) y'₀+ h² (3f₀/32 + f₁/8 + f₂/16 )
   f₄ = f(x₀+ h , y₀+ h y'₀+ h² (3f₁/7 - f₂/14 + f₃/7 )

-Then, y = y₀ + h y'₀+ h² ( 7 f₀/90 + 4 f₁/15 + f₂/15 + 4 f₃/45 ) + O(h⁷)
and y' = y'₀ + h ( 7 f₀/90 + 16 f₁/45 + 2 f₂/15 + 16 f₃/45 + 7 f₄/90 ) + O(h⁷)

Data Registers: • R00: subroutine name ( Registers R00 thru R05 are to be initialized before executing "RKN6" )

                               •   R01 = x₀      •   R04 = h = stepsize        R06 to R10: temp
                               •   R02 = y₀      •   R05 = N = number of steps
                               •   R03 = y'₀
Flags: /
Subroutine: a program which calculates f(x,y) in X-register, assuming x and y are in registers X and Y upon entry.

-Line 144 is a three-byte GTO 01

01 LBL "RKN6"
02 RCL 05
03 STO 10
04 LBL 01
05 RCL 02
06 RCL 01
07 XEQ IND 00
08 RCL 04
09 *
10 STO 06
11 8
12 /
13 RCL 03
14 +
15 RCL 04
16 *
17 4
18 /
19 RCL 02
20 +
21 RCL 04
22 4
23 /
24 RCL 01
25 +
26 XEQ IND 00
27 RCL 04
28 *
29 STO 07
30 6

31 /
32 RCL 06
33 24
34 /
35 -
36 RCL 03
37 2
38 /
39 +
40 RCL 04
41 *
42 RCL 02
43 +
44 RCL 04
45 2
46 /
47 RCL 01
48 +
49 XEQ IND 00
50 RCL 04
51 *
52 STO 08
53 ST+ X
54 RCL 07
55 4
56 *
57 +
58 RCL 06
59 3
60 *

61 +
62 8
63 /
64 RCL 03
65 3
66 *
67 +
68 RCL 04
69 *
70 4
71 /
72 RCL 02
73 +
74 RCL 04
75 .75
76 *
77 RCL 01
78 +
79 XEQ IND 00
80 RCL 04
81 *
82 STO 09
83 ST+ X
84 RCL 08
85 -
86 RCL 07
87 6
88 *
89 +
90 14

91 /
92 RCL 03
93 +
94 RCL 04
95 ST+ 01
96 *
97 RCL 02
98 +
99 RCL 01
100 XEQ IND 00
101 RCL 04
102 *
103 RCL 06
104 +
105 7
106 *
107 RCL 07
108 RCL 09
109 +
110 32
111 *
112 +
113 RCL 08
114 12
115 *
116 +
117 RCL 09
118 8
119 *
120 RCL 08

121 6
122 *
123 +
124 RCL 07
125 24
126 *
127 +
128 RCL 06
129 7
130 *
131 +
132 90
133 ST/ Z
134 /
135 RCL 03
136 +
137 RCL 04
138 *
139 ST+ 02
140 X<>Y
141 ST+ 03
142 DSE 10
143 GTO 01
144 RCL 03
145 RCL 02
146 RCL 01
147 END

( 178 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
Z	/	y'(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example: y(0) = 1 , y'(0) = 1 , 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 "RKN6"   >>>>     x   = 1                                            ---Execution time = 76s---
                           RDN   y(1) = 0.536630617
                           RDN   y'(1) = -0.860171927

-The exact results, rounded to 10D are y(1) = 0.5366306164 & y'(1) = -0.8601719268

Note:

-Press R/S to continue the calculations ( after changing h & N in registers R04 & R05 if need be )

b) Systems of 2 Second-order differential Equations

-"RKN6B" 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: subroutine name ( Registers R00 thru R07 are to be initialized before executing "RKN6B" )

                               •   R01 = x₀      •   R04 = h = stepsize                        •   R06 = y'₀          R08 to R16: temp
                               •   R02 = y₀      •   R05 = N = number of steps            •   R07 = z'₀
                               •   R03 = z₀
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.

-Line 255 is a three-byte GTO 01 ( or replace it by GTO 16 and line 04 by LNL 16 )

01 LBL "RKN6B"
02 RCL 05
03 STO 16
04 LBL 01
05 RCL 03
06 RCL 02
07 RCL 01
08 XEQ IND 00
09 RCL 04
10 ST* Z
11 *
12 STO 12
13 8
14 /
15 RCL 07
16 +
17 4
18 /
19 RCL 04
20 *
21 RCL 03
22 +
23 X<>Y
24 STO 08
25 8
26 /
27 RCL 06
28 +
29 4
30 /
31 RCL 04
32 *
33 RCL 02
34 +
35 RCL 04
36 4
37 /
38 RCL 01
39 +
40 XEQ IND 00
41 RCL 04
42 ST* Z
43 *
44 STO 13

45 6
46 /
47 RCL 12
48 24
49 /
50 -
51 RCL 07
52 2
53 /
54 +
55 RCL 04
56 *
57 RCL 03
58 +
59 X<>Y
60 STO 09
61 6
62 /
63 RCL 08
64 24
65 /
66 -
67 RCL 06
68 2
69 /
70 +
71 RCL 04
72 *
73 RCL 02
74 +
75 RCL 04
76 2
77 /
78 RCL 01
79 +
80 XEQ IND 00
81 RCL 04
82 ST* Z
83 *
84 STO 14
85 ST+ X
86 RCL 13
87 4
88 *

89 +
90 RCL 12
91 3
92 *
93 +
94 32
95 /
96 RCL 07
97 .75
98 STO 15
99 *
100 +
101 RCL 04
102 *
103 RCL 03
104 +
105 X<>Y
106 STO 10
107 ST+ X
108 RCL 09
109 4
110 *
111 +
112 RCL 08
113 3
114 *
115 +
116 32
117 /
118 RCL 06
119 RCL 15
120 *
121 +
122 RCL 04
123 *
124 RCL 02
125 +
126 RCL 04
127 RCL 15
128 *
129 RCL 01
130 +
131 XEQ IND 00
132 RCL 04

133 ST* Z
134 *
135 STO 15
136 ST+ X
137 RCL 14
138 -
139 RCL 13
140 6
141 *
142 +
143 14
144 /
145 RCL 07
146 +
147 RCL 04
148 *
149 RCL 03
150 +
151 X<>Y
152 STO 11
153 ST+ X
154 RCL 10
155 -
156 RCL 09
157 6
158 *
159 +
160 14
161 /
162 RCL 06
163 +
164 RCL 04
165 ST+ 01
166 *
167 RCL 02
168 +
169 RCL 01
170 XEQ IND 00
171 RCL 04
172 ST* Z
173 *
174 RCL 12
175 +
176 7

177 *
178 RCL 13
179 RCL 15
180 +
181 32
182 *
183 +
184 RCL 14
185 12
186 *
187 +
188 X<>Y
189 RCL 08
190 +
191 7
192 *
193 RCL 09
194 RCL 11
195 +
196 32
197 *
198 +
199 RCL 10
200 12
201 *
202 +
203 X<> 11
204 8
205 *
206 RCL 10
207 6
208 *
209 +
210 RCL 09
211 24
212 *
213 +
214 RCL 08
215 7
216 *
217 +
218 90
219 ST/ Z
220 ST/ 11

221 /
222 RCL 06
223 +
224 RCL 04
225 *
226 ST+ 02
227 X<> 11
228 ST+ 06
229 X<>Y
230 X<> 15
231 8
232 *
233 RCL 14
234 6
235 *
236 +
237 RCL 13
238 24
239 *
240 +
241 RCL 12
242 7
243 *
244 +
245 90
246 /
247 RCL 07
248 +
249 RCL 04
250 *
251 ST+ 03
252 X<> 15
253 ST+ 07
254 DSE 16
255 GTO 01
256 RCL 03
257 RCL 02
258 RCL 01
259 END

( 314 bytes / SIZE 017 )

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 "RKN6B" >>>>     x   = 1                                                                                            ---Execution time = 116s---
                            RDN   y(1) = 1.531356647     and   R06 = y'(1) = -2.312840139
                            RDN   z(1) = 2.620254282              R07 = z'(1) = 2.941748401

-The precision is almost perfect, since the exact results - rounded to 9D - are

y(1) = 1.531356646 y'(1) = -2.312840137
z(1) = 2.620254281 z'(1) = 2.941748399

-To continue the calculations, simply press R/S ( after changing h & N in registers R04 & R05 if need be )

3°) S-Stage Explicit Runge-Kutta-Nystrom formulae

a) Second-order differential Equations

-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

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

k_s = h.f(x+c_s.h;y+c_s.h.y'+h (a_s,1.k₁+a_s,2.k₂+ .......... +a_s,s-1.k_s-1) )

then, y(x+h) = y(x) + h y'(x) + h [ b₁.k₁+ ................ + b_s.k_s ] and y'(x+h) = y'(x) + [ b'₁.k₁+ ................ + b'_s.k_s ]

-Therefore, a single program may be used for all explicit Runge-Kutta-Nystrom methods,
provided the coefficients are stored in appropriate registers.

Data Registers: • R00: f name ( Registers R00 to R06 and Rs+10 to R(s²+7.s+16)/2 are to be initialized before executing "ERKN" )

               •   R01 = x₀                                      R07 to R09: temp          • Rs+10 to R(s²+7.s+16)/2 = the ( s² + 5.s - 2 ) / 2 coefficients a_i,j ; b_i ; c_i
               •   R02 = y₀                                      R10 = k₁                                                                            which are to be stored as shown below:
               •   R03 = y'₀                                     R11 = k₂
               •   R04= h = stepsize                        ............
               •   R05= N = number of steps
               •   R06 = s = number of stages        Rs+9 = k_s

• Rs+10 = a_2,1 • Rs+11 = c₂
• Rs+12 = a_3,1 • Rs+13 = a_3,2 • Rs+14 = c₃

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

• R(s²+s+18)/2 = a_s,1 ................ • R(s²+3s+14)/2 = a_s,s-1 • R(s²+3s+16)/2 = c_s

• R(s²+3s+18)/2 = b₁ ............................................................ • R(s²+5s+16)/2 = b_s
• R(s²+5s+18)/2 = b'₁ ........................................................... • R(s²+7s+16)/2 = b'_s

Flags: /
Subroutine: a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

-Line 97 is a three-byte GTO 01

01 LBL "ERKN"
02 RCL 05
03 STO 07
04 LBL 01
05 10
06 STO 08
07 3
08 RCL 06
09 ST+ Y
10 *
11 2
12 /
13 8
14 +
15 E3
16 /
17 +
18 RCL 06
19 +
20 STO 09
21 RCL 02

22 RCL 01
23 XEQ IND 00
24 RCL 04
25 *
26 STO 10
27 LBL 02
28 RCL 08
29 INT
30 E3
31 /
32 10
33 +
34 STO 08
35 CLX
36 LBL 03
37 RCL IND 08
38 RCL IND 09
39 *
40 +
41 ISG 09
42 ISG 08

43 GTO 03
44 RCL 03
45 RCL IND 09
46 *
47 +
48 RCL 04
49 *
50 RCL 02
51 +
52 RCL 04
53 RCL IND 09
54 *
55 RCL 01
56 +
57 XEQ IND 00
58 RCL 04
59 *
60 STO IND 08
61 ISG 09
62 GTO 02
63 RCL 06

64 1.001
65 -
66 ST- 08
67 CLX
68 LBL 04
69 RCL IND 08
70 RCL IND 09
71 *
72 +
73 ISG 09
74 CLX
75 ISG 08
76 GTO 04
77 RCL 03
78 +
79 RCL 04
80 ST+ 01
81 *
82 ST+ 02
83 RCL 06
84 ST- 08

85 CLX
86 LBL 05
87 RCL IND 08
88 RCL IND 09
89 *
90 +
91 ISG 09
92 CLX
93 ISG 08
94 GTO 05
95 ST+ 03
96 DSE 07
97 GTO 01
98 RCL 03
99 RCL 02
100 RCL 01
101 END

( 149 bytes / SIZE ( s^2+7.s+18 )/2 )

STACK	INPUTS	OUTPUTS
Z	/	y'(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example: y(0) = 1 , y'(0) = 1 , y" = - y ( x² + y² )^1/2

>>> Calculate y(1) & y'(1)

-Load this 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

-If we use the same ( Albrecht ) 5-stage method, store the following coefficients into R15 thru R38

     1/32                              1/4                             R15                            R16
    -1/24   1/6                     1/2                             R17   R18                   R19
     3/32    1/8   1/16           3/4            =              R20 R21 R22            R23
       0       3/7 -1/14   1/7    1                              R24 R25 R26 R27   R28
     7/90 4/15   1/15 4/45    0                             R29 R30 R31 R32   R33
     7/90 16/45 2/15 16/45 7/90                          R34 R35 R36 R37   R38

-Since it is a 5-stage method, 5 STO 06 and

XEQ "ERKN"   >>>>     x   = 1                                            ---Execution time = 147s---
                           RDN   y(1) = 0.536630617
                           RDN   y'(1) = -0.860171927

Notes:

-Press R/S to continue the calculations ( after changing h & N in registers R04 & R05 if need be )
-LBL 04 & LBL 05 are identical, so they could form a subroutine:
-The program would be shorter but slower.

Control Number of the Array of Coefficients for an S-stage method:

S	bbb.eee
5 6 7 8 9 10 11 12 13 14 15 16 17	15.038 16.047 17.057 18.068 19.080 20.093 21.107 22.122 23.138 24.155 25.173 26.192 27.212

-For instance, we can use the coefficients below - taken from reference [3] - to get a 13-stage 10th-order method.
-These 116 coefficients are to be stored into R23 thru R138 in the following order:

a( 2, 1) = 6.2550354381669436926370260206784041201048E-04 = R23
c( 2) = 3.5369578561715839852580662156128003290245E-02

a( 3, 1) =     8.3400472508892582568493680275712054934730E-04
a( 3, 2) =     1.6680094501778516513698736055142410986946E-03
c( 3)    =     7.0739157123431679705161324312256006580491E-02

a( 4, 1) =     1.4377592173651130455265029764307572975032E-02
a( 4, 2) =    -2.5520389586889438033660194426507730127969E-02
a( 4, 3) =     2.9955419091121746433503478890580675956400E-02
c( 4)    =     1.9397227470895648005755022968171580307125E-01

a( 5, 1) =     5.6606010103291956408355747332339448084576E-03
a( 5, 2) =     0.0000000000000000000000000000000000000000E+00
a( 5, 3) =     2.2438147205568485440414320313187753497836E-02
a( 5, 4) =     9.6198950134107937593128496677833016937067E-03
c( 5)    =     2.7465849059990290000000000000000000000000E-01

a( 6, 1) =     4.2497841145700712790685880544080818730288E-03
a( 6, 2) =     0.0000000000000000000000000000000000000000E+00
a( 6, 3) =     1.3737389558120657344878988260209037684744E-02
a( 6, 4) =     2.1090823734166889521084117969028056034572E-03
a( 6, 5) =    -2.1698166033104208350129690333992516122959E-04
c( 6)    =     1.9939545825206940000000000000000000000000E-01

a( 7, 1) =     8.4469918165805880106087924734363680346449E-04
a( 7, 2) =     0.0000000000000000000000000000000000000000E+00
a( 7, 3) =     7.6286170426770036487508428629504122228445E-04
a( 7, 4) =    -4.5274889794326362037916138257410810321144E-03
a( 7, 5) =    -9.4704498758374659029865654611073300453711E-05
a( 7, 6) =     4.3839567862160003018135640613696763068192E-03
c( 7)    =     5.2332097109723180000000000000000000000000E-02

a( 8, 1) =     4.0257188676144292921381323534074832686999E-02
a( 8, 2) =     0.0000000000000000000000000000000000000000E+00
a( 8, 3) =     2.8132558570671430000000000000000000000000E-01
a( 8, 4) =    -9.9086973311055342850364525233745309401480E-02
a( 8, 5) =     3.1558678646031991227831775524256371039974E-02
a( 8, 6) =     7.6263009309052972315471957912384344318273E-02
a( 8, 7) =    -2.4509110177537917817170493451316523864377E-01
c( 8)    =     4.1285926718800690000000000000000000000000E-01

a( 9, 1) =    -5.1852206149158582964016598801731374371401E-01
a( 9, 2) =     0.0000000000000000000000000000000000000000E+00
a( 9, 3) =    -4.0220496748396490000000000000000000000000E+00
a( 9, 4) =     1.3339542387088720000000000000000000000000E+00
a( 9, 5) =    -3.6280137620198678475540384353108923505414E-01
a( 9, 6) =    -4.4662932597538443824830112336373381632573E-01
a( 9, 7) =     4.1093353899979717196794462089910810944472E+00
a( 9, 8) =     8.3371860107714029551447509208700700646694E-02
c( 9)    =     5.9440567007045230000000000000000000000000E-01

a(10, 1) =     4.5526515039304022037674417630889739137541E-01
a(10, 2) =     0.0000000000000000000000000000000000000000E+00
a(10, 3) =     3.4410244296399106992261498105244921193431E+00
a(10, 4) =    -1.1420304634021581044575027059094514599562E+00
a(10, 5) =     3.3122504529344460041163434350588640458224E-01
a(10, 6) =     5.2510814510725835145952044841619177157062E-01
a(10, 7) =    -3.4007204966989148959577640759197726181960E+00
a(10, 8) =     1.4959025547025800806021814497275304454960E-02
a(10, 9) =     1.7714834024190792778113334621606086825758E-02
c(10)    =     6.9648498893199050000000000000000000000000E-01

a(11, 1) =    -5.0738859071291316865190644559724141522647E-02
a(11, 2) =     0.0000000000000000000000000000000000000000E+00
a(11, 3) =    -8.5258315802766612821674340523993103859370E-01
a(11, 4) =     2.9256281989303017562028747765972648095241E-01
a(11, 5) =    -4.2631304538837970000000000000000000000000E-01
a(11, 6) =     3.0845126792446679274386181736916000785280E-01
a(11, 7) =     8.2068063064681728553592877499937152955118E-01
a(11, 8) =     2.6353200561785124318114669620069412596405E-01
a(11, 9) =    -3.8002959536498107924589290701514760515073E-02
a(11,10) =     5.0836949957876193531014830503462796310983E-02
c(11)    =     8.5840043338316930000000000000000000000000E-01

a(12, 1) =    -8.5506341904459305448936612697062868735670E-01
a(12, 2) =     0.0000000000000000000000000000000000000000E+00
a(12, 3) =    -4.4172967833111218167905547111397401957735E+00
a(12, 4) =     1.4620414738250776987869127049027576880060E+00
a(12, 5) =     1.5800603670627463865955782231377536485336E+00
a(12, 6) =    -2.1500473087668603894015427990085825826253E+00
a(12, 7) =     5.2192882951843358632420497699580187954903E+00
a(12, 8) =    -7.0122224598841812446074882425859149461800E-01
a(12, 9) =     4.0674284722476557322266460477065027849355E-01
a(12,10) =    -1.0995016400771572122089410865399735394434E-01
a(12,11) =     2.5498951087297871672803054502111066906939E-02
c(12)    =     9.5922053070763063933434705195204984252787E-01

a(13, 1) =     3.6124205767243278950403294321213352485132E+00
a(13, 2) =     0.0000000000000000000000000000000000000000E+00
a(13, 3) =     1.9614376424164347318861614225759179537283E+01
a(13, 4) =    -6.5540954527470471068838634201654381264077E+00
a(13, 5) =    -5.3634775178894355119823336908885299369771E+00
a(13, 6) =     8.9549200633414078297191926517465420210063E+00
a(13, 7) =    -2.2111999575685298066957085805778892539418E+01
a(13, 8) =     3.0666418328538943745146953495973014350614E+00
a(13, 9) =    -1.2974380337600327021222956536808230517716E+00
a(13,10) =     5.9822684827188431105752038905995447885348E-01
a(13,11) =    -2.4107078892562108246647516780668823095046E-02
a(13,12) =     4.5319136185137669988740390100397569527517E-03
c(13)    =     1.0000000000000000000000000000000000000000E+00

b( 1)    =     1.1445045431083081161076675575316257764110E-02
b( 2)    =     0.0000000000000000000000000000000000000000E+00
b( 3)    =     0.0000000000000000000000000000000000000000E+00
b( 4)    =     0.0000000000000000000000000000000000000000E+00
b( 5)    =     0.0000000000000000000000000000000000000000E+00
b( 6)    =     1.5182228814165001267592100569234113490800E-01
b( 7)    =     9.3833383282371058262139365435233240765713E-02
b( 8)    =     1.3138871401731356660181720771852165705253E-01
b( 9)    =     4.2452793993460570894743314052986770767710E-02
b(10)    =     4.6614359052634087726314462069090004222681E-02
b(11)    =     1.9739340751760337610363200447178885100858E-02
b(12)    =     2.7040753297272850676247690093320494184034E-03
b(13)    =     0.0000000000000000000000000000000000000000E+00

bp( 1)   =     1.1445045431083081161076675575316257764110E-02
bp( 2)   =     0.0000000000000000000000000000000000000000E+00
bp( 3)   =     0.0000000000000000000000000000000000000000E+00
bp( 4)   =     0.0000000000000000000000000000000000000000E+00
bp( 5)   =     0.0000000000000000000000000000000000000000E+00
bp( 6)   =     1.8963455766836141912201425856209518713321E-01
bp( 7)   =     9.9015048411147152930074891966312987547087E-02
bp( 8)   =     2.2377720821386995442668671882695678061627E-01
bp( 9)   =     1.0466811506175315699681616849938105031123E-01
bp(10)   =     1.5358172529459859690396485461784888933400E-01
bp(11)   =     1.3940255061073114260364669651015505753606E-01
bp(12)   =     6.6309723413523447970758037743751771753948E-02
bp(13)   =     1.2166025894932047884961697698182018004080E-02 = R138

-With the same differential equation, h = 0.1 & N = 10 ( don't forget to store 13 into R06 ), it yields:

XEQ "ERKN"   >>>>     x   = 1                                            ---Execution time = 9m22s---
                           RDN   y(1) = 0.5366306165
                           RDN   y'(1) = -0.8601719269

-The exact results, rounded to 10D are y(1) = 0.5366306164 & y'(1) = -0.8601719268
-So, the error is only 1 unit in the last decimal.

-Of course, it's not very fast with a real HP-41, but with a good emulator or with a 41CL...

b) Systems of 2 Second-order differential Equations

-A similar program may be written to solve the system y" = f(x,y,z) z" = g(x,y,z)

Data Registers: • R00: fg name ( Registers R00 to R08 and R13+2s to R(s²+9.s+22)/2 are to be initialized before executing "ERKN2" )

               •   R01 = x₀                                      R09 to R12: temp          • R13+2s to R(s²+9.s+22)/2 = the ( s² + 5.s - 2 ) / 2 coefficients a_i,j ; b_i ; c_i
               •   R02 = y₀                                      R13 to R12+2s = k^y_i & k^z_i                                                which are to be stored as shown below:
               •   R03 = z₀
               •   R04= h = stepsize
               •   R05= N = number of steps
               •   R06 = y'₀
               •   R07 = z'₀
               •   R08 = s = number of stages

• R13+2s = a_2,1 • R14+2s = c₂
• R15+2s = a_3,1 • R16+2s = a_3,2 • R17+2s = c₃

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

• R(s²+3s+24)/2 = a_s,1 ................ • R(s²+5s+20)/2 = a_s,s-1 • R(s²+5s+22)/2 = c_s

• R(s²+5s+24)/2 = b₁ ............................................................ • R(s²+7s+22)/2 = b_s
• R(s²+7s+24)/2 = b'₁ ........................................................... • R(s²+9s+22)/2 = b'_s

Flags: /
Subroutine: a program that 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.

-Line 153 is a three-byte GTO 01

01 LBL "ERKN2"
02 RCL 05
03 STO 09
04 LBL 01
05 13
06 STO 10
07 RCL 08
08 +
09 STO 11
10 5
11 LASTX
12 ST+ Y
13 *
14 2
15 /
16 11
17 +
18 E3
19 /
20 +
21 RCL 08
22 +
23 STO 12
24 RCL 03
25 RCL 02
26 RCL 01
27 XEQ IND 00
28 RCL 04
29 ST* Z
30 *
31 STO IND 11
32 X<>Y

33 STO 13
34 LBL 02
35 RCL 10
36 INT
37 E3
38 /
39 13
40 +
41 STO 10
42 RCL 08
43 +
44 STO 11
45 CLST
46 LBL 03
47 X<>Y
48 RCL IND 10
49 RCL IND 12
50 *
51 +
52 X<>Y
53 RCL IND 11
54 RCL IND 12
55 *
56 +
57 ISG 12
58 ISG 11
59 CLX
60 ISG 10
61 GTO 03
62 RCL 07
63 RCL IND 12
64 *

65 +
66 RCL 04
67 *
68 RCL 03
69 +
70 X<>Y
71 RCL 06
72 RCL IND 12
73 *
74 +
75 RCL 04
76 *
77 RCL 02
78 +
79 RCL 04
80 RCL IND 12
81 *
82 RCL 01
83 +
84 XEQ IND 00
85 RCL 04
86 ST* Z
87 *
88 STO IND 11
89 X<>Y
90 STO IND 10
91 ISG 12
92 GTO 02
93 RCL 08
94 1.001
95 -
96 ST- 10

97 ST- 11
98 CLST
99 LBL 04
100 X<>Y
101 RCL IND 10
102 RCL IND 12
103 *
104 +
105 X<>Y
106 RCL IND 11
107 RCL IND 12
108 *
109 +
110 ISG 12
111 CLX
112 ISG 11
113 CLX
114 ISG 10
115 GTO 04
116 RCL 07
117 +
118 RCL 04
119 *
120 ST+ 03
121 CLX
122 RCL 06
123 +
124 RCL 04
125 ST+ 01
126 *
127 ST+ 02
128 RCL 08

129 ST- 10
130 ST- 11
131 CLST
132 LBL 05
133 X<>Y
134 RCL IND 10
135 RCL IND 12
136 *
137 +
138 X<>Y
139 RCL IND 11
140 RCL IND 12
141 *
142 +
143 ISG 12
144 CLX
145 ISG 11
146 CLX
147 ISG 10
148 GTO 05
149 ST+ 07
150 X<>Y
151 ST+ 06
152 DSE 09
153 GTO 01
154 RCL 03
155 RCL 02
156 RCL 01
157 END

( 225 bytes / SIZE ( s^2+9.s+24 )/2 )

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

-If we again use the same ( Albrecht ) 5-stage method, store the following coefficients into R23 thru R46

     1/32                              1/4                             R23                            R24
    -1/24   1/6                     1/2                             R25 R26                    R27
     3/32    1/8   1/16           3/4            =              R28 R29 R30            R31
       0       3/7 -1/14   1/7    1                              R32 R33 R34 R35   R36
     7/90 4/15   1/15 4/45    0                             R37 R38 R39 R40   R41
     7/90 16/45 2/15 16/45 7/90                          R42 R43 R44 R45   R46

-Since it is a 5-stage method, 5 STO 08 and

XEQ "ERKN2" >>>>     x   = 1                                                                                            ---Execution time = 3m41s---
                            RDN   y(1) = 1.531356646     and   R06 = y'(1) = -2.312840140
                            RDN   z(1) = 2.620254282              R07 = z'(1) = 2.941748401

Note:

-To continue the calculations, simply press R/S ( after changing h & N in registers R04 & R05 if need be )

Control Number of the Array of Coefficients for an S-stage method:

S	bbb.eee
5 6 7 8 9 10 11 12 13 14 15 16 17	23.046 25.056 27.067 29.079 31.092 33.106 35.121 37.137 39.154 41.172 43.191 45.211 47.232

-To use the 13-stage 10th-order method mentionned in the paragraph above,
store the same 116 coefficients into R39 to R154.

-With the same system of 2 differential equations, h = 0.1 & N = 10 ( don't forget to store 13 into R08 ), it yields:

XEQ "ERKN2"   >>>>     x   = 1                                                                                            ---Execution time = 14m58s---
                             RDN   y(1) = 1.531356645     and   R06 = y'(1) = -2.312840138
                             RDN   z(1) = 2.620254282              R07 = z'(1) = 2.941748401

-Albrecht's method already produced an almost perfect accuracy with this example,
so the superiority of the 10th-order method is not very apparent !
-In most cases, however, it's really worthwhile to use high-order formulae...

Notes:

-These programs employs 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 references [2] & [3], you will find the coefficients of pairs of embedded Runge-Kutta-Nystrom formulas
of various orders to control the stepsize more directly.

References:

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] Erwin Fehlberg - "Classical eighth- and lower-order Runge-Kutta-Nystrom formulas with stepsize control
for special second-order differential equations" - NASA TR R-381 & TR R-410
[3] Philip Sharp - http://www.math.auckland.ac.nz/~sharp/rkn1.dat

hp41programs

Runge-Kutta-Nystrom Methods for the HP-41