Gill-Runge-Kutta

Gill Runge-Kutta Method for the HP-41

Overview

-The Gill's method is a 4-stage 4th order Runge-Kutta method to solve differential equations,
but the coefficients are choosen so that the formulae may be written
to use only 3 registers per variable instead of 4 registers in the "classical" RK4.
-So, we can solve larger systems of ODE.

Program Listing

-We solve the system:

     dy₁/dx = f₁(x,y₁, ... ,y_n)
     dy₂/dx = f₂(x,y₁, ... ,y_n)
             ..................                         with the initial values:     y_i(x₀)         i = 1 , ... , n

dy_n/dx = f_n(x,y₁, ... ,y_n)

Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "GRK" )

                                      • R01 = n = number of equations = number of functions              R04 thru R17: temp        R18 & R19 are unused
                                      • R02 = h = stepsize
                                      • R03 = N = number of steps

                                      • R20 = x₀
                                      • R21 = y₁(x₀)                                                                           R21+n thru R20+3n: temp
                                      • R22 = y₂(x₀)
                                          ...............         ( initial values )

                                      • R20+n = y_n(x₀)
Flags: /
Subroutine:   A program which calculates and stores f₁ ( x, y₁ , ... , y_n) , ... , f_n ( x, y₁ , ... , y_n) in registers R21+n, ... , R20+2n respectively
                        with x, y₁, ... , y_n in regiters R20, R21, ... , R20+n

01 LBL "GRK"
02 20
03 RCL 01
04 +
05 STO 14
06 LASTX
07 +
08 STO 15
09 LASTX
10 +
11 STO 16
12 RCL 14
13 E3
14 /
15 ST+ 15
16 SIGN
17 2
18 1/X
19 STO 06
20 SQRT
21 -
22 STO 07
23 ST+ X
24 STO 09
25 1/X
26 STO 10
27 18
28 SQRT

29 4
30 +
31 CHS
32 STO 11
33 1/X
34 CHS
35 STO 08
36 8
37 SQRT
38 4
39 +
40 STO 12
41 6
42 1/X
43 STO 13
44 LBL 00
45 RCL 03
46 STO 17
47 LBL 01
48 XEQ IND 00
49 RCL 14
50 RCL 15
51 RCL 02
52 RCL 06
53 *
54 ST+ 20
55 LBL 02
56 CLX

57 RCL 02
58 ST* IND Y
59 CLX
60 RCL IND Y
61 STO IND 16
62 LASTX
63 *
64 ST+ IND Z
65 DSE 16
66 DSE Z
67 DSE Y
68 GTO 02
69 RCL 01
70 ST+ 16
71 XEQ IND 00
72 RCL 07
73 STO 04
74 XEQ 03
75 RCL 08
76 STO 04
77 RCL 09
78 STO 05
79 XEQ 04
80 XEQ IND 00
81 RCL 10
82 STO 04
83 XEQ 03
84 RCL 11

85 STO 04
86 RCL 12
87 STO 05
88 XEQ 04
89 RCL 02
90 RCL 06
91 *
92 ST+ 20
93 XEQ IND 00
94 RCL 13
95 STO 04
96 XEQ 03
97 DSE 17
98 GTO 01
99 RCL 23
100 RCL 22
101 RCL 21
102 RCL 20
103 RTN
104 GTO 00
105 LBL 03
106 RCL 02
107 ST* IND 15
108 RCL IND 15
109 RCL IND 16
110 -
111 RCL 04
112 *

113 ST+ IND 14
114 DSE 14
115 DSE 16
116 DSE 15
117 GTO 03
118 RCL 01
119 ST+ 14
120 ST+ 15
121 ST+ 16
122 RTN
123 LBL 04
124 RCL IND 16
125 RCL 04
126 *
127 RCL IND 15
128 RCL 05
129 *
130 +
131 STO IND 16
132 DSE 16
133 DSE 15
134 GTO 04
135 RCL 01
136 ST+ 15
137 ST+ 16
138 RTN
139 END

( 205 bytes / SIZE 3n+21 )

STACK	INPUTS	OUTPUTS
T	/	y₃(x₀+N.h)
Z	/	y₂(x₀+N.h)
Y	/	y₁(x₀+N.h)
X	/	x₀+N.h

Example: Solve the system:

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

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

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 LASTX
18 RCL 21
19 *
20 RCL 22
21 RCL 23
22 *
23 -
24 STO 26
25 RTN

"T" ASTO 00

Initial values:

           0   STO 20
           1   STO 21    STO 22
           2   STO 23

n = 3     STO 01    ( 3 functions )
h = 0.1   STO 02
N = 10   STO 03    XEQ "GRK" >>>>      x   = 1                     = R20                      ( in 2mn32s )
                                                     RDN    y₁(1) = 0.258210908 = R21                                                                         y₁(1) = 0.258207906
                                                     RDN    y₂(1) = 1.157620520 = R22         the exact results, rounded to 9D are          y₂(1) = 1.157623981
                                                     RDN    y₃(1) = 0.842179307 = R23                                                                         y₃(1) = 0.842178312

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

-There is no built-in error estimate, so we must do the calculations again with a smaller stepsize h.
-With h = 0.05 & N = 20, it yields:

    y₁(1) = 0.258208088
    y₂(1) = 1.157623789
    y₃(1) = 0.842178380

Notes:

-Since the SIZE = 3n+21 instead of 4n+11 with RK4, we save registers only for n > 10
-Theoretically, "GRK" could solve a system of 99 differential equations.
-The method seems slightly more accurate than RK4 to solve the gravitational n-body problem...

References:

[1] F. Scheid - "Numerical Analysis" - Mc Graw Hill ISBN 07-055197-9
[2] J.C. Butcher - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons - ISBN 0-471-91046-5

hp41programs

Gill Runge-Kutta Method for the HP-41