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:
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 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 = x0
• R21 = y1(x0)
R21+n thru R20+3n: temp
• R22 = y2(x0)
............... ( initial
values )
• R20+n = yn(x0)
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 "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 | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example: Solve 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
| 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 y1(1) = 0.258210908 = R21
y1(1) = 0.258207906
RDN y2(1) = 1.157620520 = R22
the exact results, rounded to 9D are
y2(1) = 1.157623981
RDN y3(1) = 0.842179307 = R23
y3(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:
y1(1) = 0.258208088
y2(1) = 1.157623789
y3(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