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:

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