Differential Equations & Gravitational n-body Problem

 Overview
 

-'DIFF'  is a directory ( binary file / #55319d / 9713.5 bytes ) that works with an HP48S/SX/G/GX
 ( In most of the routines, the <<  >>  have been deleted to save bytes ).

-However, the programs in the 'RkF' subdirectory doesn't work with an HP48S/SX

-'DIFF50' is an ascii file ( #47108d / 9848 bytes ) with the same programs that should be transfered to an HP50G too.

-These functions may solve any differential equation or system of differential equations ( provided it's not too large ).
-They are also applied to solve the gravitational n-body problem in an inertial frame of reference.
-However, the effects of general relativity are neglected.
 
 
 

Functions

Description

RK4  RK8  RK10  GM      rk4      rk8      rk10      N      h      GM      G      M1, M2 ...  RKV8E9  RkF      X      Y      F      TOL      rkf      GM2         T         Y         TOL         G         GM         rkf         M1, M2       InFo  IRK12  IRK13

Runge-Kutta "classical" 4th order formula  Runge-Kutta 8th order formula  Runge-Kutta 10th order formula  Directory  (  gravitational n-body problem )  Runge-Kutta 4th order formula  Runge-Kutta 8th order formula  Runge-Kutta 10th order formula  Number of steps  Stepsize expressed in days  Main subroutine  ( Gaussian gravit. contant )^2 = ( 0.01720209895 )^2  Mass of the 1st celestial body, 2nd body , ...  Runge-Kutta-Verner 8-9 method ( with error estimate )  Directory ( doesn't work with an HP48S/SX )  Independant Variable  Y = [ Y1 .... Yn ] = n-vector  A program that computes dY/dX = [ Y1' .... Yn' ]  Tolerance = small positive number  Uses the built-in Runge-Kutta-Fehlberg 4-5 function  Subdirectory  (  gravitational n-body problem by RKF )  Time  Vector [ X1 Y1 Z1 ... Xn Yn Zn X'1 Y'1 Z'1 ... X'n Y'n Z'n ]  Tolerance  ( Gaussian gravit. Contant )^2 = ( 0.01720209895 )^2  Main subroutine  Uses the built-in Runge-Kutta-Fehlberg 4-5 function  Mass of the 1st celestial body, 2nd body , ...  Info  Implicit Runge-Kutta method of order 12  Implicit Runge-Kutta method of order 13

 

Runge-Kutta "classical" 4th order formula
 

  RK4 uses the following formulae:

         k₁ = h.f(x;y)
         k₂ = h.f(x+h/2;y+k₁/2)
         k₃ = h.f(x+h/2;y+k₂/2)
         k₄ = h.f(x+h;y+k₃)

  and then,   y(x+h) = y(x) + ( k₁ + 2.k₂ + 2.k₃ + k₄ ) / 6

-This formula duplicates the Taylor series up to the term in h⁴

-In order to save bytes, the equations are written  dY/dX = f(Y)  with  Y = [ X  Y₁ , ......... , Y_n ]
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f >>	<< f >>
Level 3	h	h
Level 2	N	N
Level 1	[ x y1 ... yn ]	[ x y1 ... yn ]*

  where

    << f >>  is a program ( or its name ) that computes the derivative  of the (n+1)-vector  [ x y₁ ... y_n ]   Since dx/dx = 1 the first component is always 1
         h = stepsize
        N = number of steps
  [ x y₁ ... y_n ] = initial value  and  [ x y₁ ... y_n ]* is the result with  x* = x + N.h

Example1:    y(0) = 1 and  dy/dx = 2 x.y     Evaluate  y(1)              ( The exact solution is y = exp ( x² ) )

-With  h = 0.1  &  N = 10

  <<  OBJ->  DROP  *  2  *  1  SWAP  2  ->ARRY   >>   ENTER
                                          .1                                              ENTER
                                         10                                              ENTER
                                     [ 0   1 ]                                            RK4             gives        [ 1  2.71827017536 ]

-So,  y(1) = 2.71827017536
-The exact value is exp(1) = 2.71828182846   rounded to 12D

-Press  RK4 again to get  y(2)  ... and so on ...

Example2:

         dy/dx = -yzu         y(0)=1
         dz/dx = x(y+z-u)   z(0)=1
         du/dx = xy-zu       u(0)=2

-Evaluate y(1) , z(1) , u(1)  with h = 0.1 , N = 10:

 << OBJ-> DROP -> X Y Z U
  << 1 Y Z U * * NEG
     Y Z + U - X *
     X Y * Z U * - 4 ->ARRY
  >>
 >>              ENTER
 0.1             ENTER
 10              ENTER
 [ 0 1 1 2 ]   RK4               yields           [ 1   0.258209385512   1.15761955337   0.842178650981 ]

-Press   RK4 again to get y(2)  z(2)  u(2)  .....

Notes:

-Since there is no error estimate, the calculations must be performed with a smaller stepsize, for instance  h = 0.05
-To solve, say a 2nd order differential equation like y" = f( x , y , y' ), replace it by the system

  dy/dx = z
  dz/dx = f( x , y , z )
 

Runge-Kutta 8th order formula
 

-RK8  employs an 8th-order formula that is given in reference [4]
-It duplicates the Taylor series up to the term in h⁸
-11 evaluations of the function are required at each step.

-RK8 & RK4 are used in the same way.
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f >>	<< f >>
Level 3	h	h
Level 2	N	N
Level 1	[ x y1 ... yn ]	[ x y1 ... yn ]*

  where

    << f >>  is a program ( or its name ) that computes the derivative  of the (n+1)-vector  [ x y₁ ... y_n ]   Since dx/dx = 1 the first component is always 1
         h = stepsize
        N = number of steps
  [ x y₁ ... y_n ] = initial value  and  [ x y₁ ... y_n ]* is the result with  x* = x + N.h

Example1:    y(0) = 1 and  dy/dx = 2 x.y     Evaluate  y(1)              ( The exact solution is y = exp ( x² ) )

-With  h = 0.1  &  N = 10

  <<  OBJ->  DROP  *  2  *  1  SWAP  2  ->ARRY   >>   ENTER
                                          .1                                              ENTER
                                         10                                              ENTER
                                     [ 0   1 ]                                            RK8             gives        [ 1  2.7182818285 ]

-Press  RK8 again to get  y(2)  ... and so on ...

Example2:

         dy/dx = -yzu         y(0)=1
         dz/dx = x(y+z-u)   z(0)=1
         du/dx = xy-zu       u(0)=2

-Evaluate y(1) , z(1) , u(1)  with h = 0.1 , N = 10:

 << OBJ-> DROP -> X Y Z U
  << 1 Y Z U * * NEG
     Y Z + U - X *
     X Y * Z U * - 4 ->ARRY
  >>
 >>              ENTER
 0.1             ENTER
 10              ENTER
 [ 0 1 1 2 ]   RK8               yields           [ 1   0.258207906459   1.1576239808   0.842178311686 ]

-Press   RK8 again to get y(2)  z(2)  u(2)  .....

Note:

-Since there is no error estimate, the calculations must be performed with a smaller stepsize, for instance  h = 0.05
 

Runge-Kutta 10th order formula
 

-RK10 duplicates the Taylor series up to the term in h¹⁰
-17 evaluations of the function are required at each step.

-RK10 , RK8 & RK4 are used in the same way.
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f >>	<< f >>
Level 3	h	h
Level 2	N	N
Level 1	[ x y1 ... yn ]	[ x y1 ... yn ]*

  where

    << f >>  is a program ( or its name ) that computes the derivative  of the (n+1)-vector  [ x y₁ ... y_n ]   Since dx/dx = 1 the first component is always 1
         h = stepsize
        N = number of steps
  [ x y₁ ... y_n ] = initial value  and  [ x y₁ ... y_n ]* is the result with  x* = x + N.h

Example1:    y(0) = 1 and  dy/dx = 2 x.y     Evaluate  y(1)              ( The exact solution is y = exp ( x² ) )

-With  h = 0.1  &  N = 10

  <<  OBJ->  DROP  *  2  *  1  SWAP  2  ->ARRY   >>   ENTER
                                          .1                                              ENTER
                                         10                                              ENTER
                                     [ 0   1 ]                                           RK10             gives        [ 1  2.71828182846 ]

-Press  RK10 again to get  y(2)  ... and so on ...

Example2:

         dy/dx = -yzu         y(0)=1
         dz/dx = x(y+z-u)   z(0)=1
         du/dx = xy-zu       u(0)=2

-Evaluate y(1) , z(1) , u(1)  with h = 0.1 , N = 10:

 << OBJ-> DROP -> X Y Z U
  << 1 Y Z U * * NEG
     Y Z + U - X *
     X Y * Z U * - 4 ->ARRY
  >>
 >>              ENTER
 0.1             ENTER
 10              ENTER
 [ 0 1 1 2 ]   RK10               yields           [ 1   0.258207906453   1.15762398081   0.842178311706 ]

-Press   RK10 again to get y(2)  z(2)  u(2)  .....

Notes:

-Since there is no error estimate, the calculations must be performed with a smaller stepsize, for instance  h = 0.05
-Of course,the execution time increases like the number of functions in the system.
 

Runge-Kutta-Verner 8th order formula enbedded within 9th order
 

-RKV8E9  employs an 8th order formula and a 9th order formula.
-The difference between these gives an error estimate of the 8th order result
-This difference is stored in a variable 'ERROR'
-Always PURGE the variable 'ERROR' before the 1st execution.

-The successive "errors" are simply added to 'ERROR', so the estimation of the error may be underestimated...
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f >>	<< f >>
Level 3	h	h
Level 2	N	N
Level 1	[ x y1 ... yn ]	[ x y1 ... yn ]*

  where

    << f >>  is a program ( or its name ) that computes the derivative  of the (n+1)-vector  [ x y₁ ... y_n ]   Since dx/dx = 1 the first component is always 1
         h = stepsize
        N = number of steps
  [ x y₁ ... y_n ] = initial value  and  [ x y₁ ... y_n ]* is the result with  x* = x + N.h

Example:

         dy/dx = -yzu         y(0)=1
         dz/dx = x(y+z-u)   z(0)=1
         du/dx = xy-zu       u(0)=2

 with h = 0.1 , N = 10:

*** PURGE 'ERROR' before the first execution.

 << OBJ-> DROP -> X Y Z U
  << 1 Y Z U * * NEG
     Y Z + U - X *
     X Y * Z U * - 4 ->ARRY
  >>
 >>   ENTER
 0.1  ENTER
 10   ENTER
 [ 0 1 1 2 ]  and press RK8E9

-It returns in level 1          [ 1 0.258207906449 1.15762398084 0.842178311722 ]

 and 'ERROR' contains [ 0 2.7631E-11 3.27392E-11 7.636E-11 ]

-Press   RK8E9   again to continue ( after changing h and/or N if need be )
 

Implicit Runge-Kutta formula of order 12
 

-IRK12 & IRK13 use implicit Runge-Kutta formulae:

-A n-stage implicit Runge-Kutta method is defined by   n(n+2)  coefficients a_i,j ; b_i ; c_i

                             k₁ = h.f ( x + b₁h ; y + a_1,1 k₁ +  a_1,2 k₂ +  ..... + a_1,n k_n )
                              .....................................................................................            (  actually,     a_i,1 + a_i,2 + .......... + a_i,n = b_i  )
                             k_n = h.f ( x + b_nh ; y + a_n,1 k₁ +  a_n,2 k₂ +  ..... + a_n,n k_n )

                                 then:     y(x+h) = y(x) + c₁.k₁+ ................ + c_n.k_n

-So, the calculator must solve a non-linear system at each step !
-It is solved by a simple successive approximation method which may ... not converge !
-Even if it does converge, the execution time makes IRK12 & IRK13 very slow.

-These 2 programs are given for completeness but I'm not sure they will be very useful...
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f >>	<< f >>
Level 3	h	h
Level 2	N	N
Level 1	[ x y1 ... yn ]	[ x y1 ... yn ]*

  where

    << f >>  is a program ( or its name ) that computes the derivative  of the (n+1)-vector  [ x y₁ ... y_n ]   Since dx/dx = 1 the first component is always 1
         h = stepsize
        N = number of steps
  [ x y₁ ... y_n ] = initial value  and  [ x y₁ ... y_n ]* is the result with  x* = x + N.h

Example:    y(0) = 1 and  dy/dx = 2 x.y     Evaluate  y(1)              ( The exact solution is y = exp ( x² ) )

-With  h = 0.1  &  N = 10

  <<  OBJ->  DROP  *  2  *  1  SWAP  2  ->ARRY   >>   ENTER
                                          .1                                              ENTER
                                         10                                              ENTER
                                     [ 0   1 ]                                           IRK12             gives        [ 1  2.71828182846 ]

-Press  IRK12 again to get  y(2)  ... and so on ...

Note:

-The HP48 displays successive values so that you can check if the process converges or not.
 

Implicit Runge-Kutta formula of order 13
 

-Same remarks.
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f >>	<< f >>
Level 3	h	h
Level 2	N	N
Level 1	[ x y1 ... yn ]	[ x y1 ... yn ]*

  where

    << f >>  is a program ( or its name ) that computes the derivative  of the (n+1)-vector  [ x y₁ ... y_n ]   Since dx/dx = 1 the first component is always 1
         h = stepsize
        N = number of steps
  [ x y₁ ... y_n ] = initial value  and  [ x y₁ ... y_n ]* is the result with  x* = x + N.h

Example:    y(0) = 1 and  dy/dx = 2 x.y     Evaluate  y(1)              ( The exact solution is y = exp ( x² ) )

-With  h = 0.1  &  N = 10

  <<  OBJ->  DROP  *  2  *  1  SWAP  2  ->ARRY   >>   ENTER
                                          .1                                              ENTER
                                         10                                              ENTER
                                     [ 0   1 ]                                           IRK13             gives        [ 1  2.71828182846 ]

-Press  IRK13 again to get  y(2)  ... and so on ...

Note:

-The HP48 displays successive values so that you can check if the process converges or not.
 

Gravitational n-body Problem
 

-You can choose between  'rk4'  'rk8'  'rk10'  ( which call  'RK4'  'RK8'  'RK10' and 'GM' )  to solve the gravitational n-body problem.

-Let   M₁( x₁, y₁,z₁) ; ....... ; M_n( x_n, y_n,z_n)  be n celestial bodies  with initial velocities   V₁( x'₁, y'₁,z'₁) ; ....... ; V_n( x'_n, y'_n,z'_n)  at an instant t
 and    m₁ , ...... , m_n   the n masses of these bodies.

-We want to know their future positions and velocities at an instant  t + N.h    ( h = step size in days ; N = number of steps )
-So, we have to solve the system of (6n+1) differential equations:

              d²x_i /dt² = SUM_j#i   Gm_j ( x_j - x_i ) / [ ( x_i - x_j )²+ ( y_i - y_j )² + ( z_i - z_j )² ]^3/2
              d²y_i /dt² = SUM_j#i   Gm_j ( y_j - y_i ) / [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^3/2          i , j = 1 , ..... , n
              d²z_i /dt² = SUM_j#i   Gm_j ( z_j - z_i ) / [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^3/2

 where G is the gravitational constant:  G = k²  where  k = 0.01720209895 is Gaussian gravitational constant.

-The distances are to be expressed in Astronomical Units, the times in days and the masses in Solar mass.

1°)  Store the n masses into  'M1'  'M2'  .....
2°)  Store the stepsize in 'h'
3°)  Store the number of steps in 'N'

4°)  Place in level 1 the (6n+1)-vector  [ t  x₁, y₁,z₁ ............... x_n, y_n,z_n , x'₁, y'₁,z'₁ ............... x'_n, y'_n,z'_n ]

      and press 'rk4'  or  'rk8'  or  'rk10' ,  you'll get the new vector in level 1

5°)  Press again  'rk4'  or  'rk8'  or  'rk10'  to get the position & velocities at the instant  t + 2N h
 
 
 

STACK	INPUT	OUTPUT
Level 1	[ t  x1, y1,z1 ............... xn, yn,zn , x'1, y'1,z'1 ............... x'n, y'n,z'n ]	[ t  x1, y1,z1 ............... xn, yn,zn , x'1, y'1,z'1 ............... x'n, y'n,z'n ]*

  where  [ t  x₁, y₁,z₁ ............... x_n, y_n,z_n , x'₁, y'₁,z'₁ ............... x'_n, y'_n,z'_n ]   is the initial vector
    and    [ t  x₁, y₁,z₁ ............... x_n, y_n,z_n , x'₁, y'₁,z'₁ ............... x'_n, y'_n,z'_n ]*  is the position-velocity vector at  t + N.h

Example:

-Let's imagine a system of three stars      M₁ ( 2 ; 0 ;0 )       M₂ ( 0 ; 4 ; 0 )      M₃ ( 0 ; 0 ; 1 )                     unit = 1 AU
              with  initial velocities                 V₁ ( 0 ; 0.03 ; 0 )  V₂ ( 0 ; 0 ; 0.01 )  V₃ ( -0.02 ; 0 ; 0 )               unit = 1 AU per day
                     and masses                        m₁ = 2  ;  m₂ = 1  ;  m₃ = 3                                                           unit = 1 Solar mass

-Calculate the positions and velocities 10 days later.

-Masses:   2   'M1'  STO   1  'M2'  STO   3  'M3'  STO

-With  N = 1  &  h = 10 days    1  'N'  STO    10   'h'   STO   and if you want to use  'rk4'

   [ 0  2  0  0  0  4  0  0  0  1  0  0.03  0  0  0  0.01  -0.02  0  0 ]    rk4   yields  ( rounded to 9 decimals )

   [ 10   1.992077551   0.300333856   0.003673779   0.000661665   3.996080593   0.100603408   -0.194938922   0.001083898   0.997349678
            -0.001550089  0.030038159   0.000706688   0.000132598   -0.000790383   0.010117548   -0.019010806   0.000238022   -0.000510308  ]
 

-With N = 2 & h = 5 days, the results are given in the last column hereunder
 
 
 

-	coord	input	h = 10 ; N=1	h = 5 ; N=2
p   o   s.	x1      y1      z1 --------      x2      y2      z2 --------      x3      y3      z3	2    0    0 -----    0    4    0 -----    0    0    1	1.992077551 0.300333856 0.003673779  -------------- 0.000661665 3.996080593 0.100603408  -------------- -0.194938922  0.001083898  0.997349678	1.992077584     0.300333570     0.003673683      --------------     0.000661669     3.996080575     0.100603412     --------------    -0.194938946     0.001084095     0.997349741
v   e   l   .	x'1      y'1      z'1 --------      x'2      y'2      z'2 --------      x'3      y'3      z'3	0  0.03    0 -----    0    0  0.01 ----- -0.02    0    0	-0.001550089  0.030038159  0.000706688  --------------  0.000132598 -0.000790383  0.010117548  -------------- -0.019010806  0.000238022 -0.000510308	-0.001550083     0.030038158     0.000706684     --------------     0.000132598    -0.000790384     0.010117549     --------------    -0.019010811     0.000238023    -0.000510306

 

Notes:

-Of course, 'rk8' & 'rk10' will produce more accurate results but the execution time will be increased too.
-The directory contains 10 variables 'M1'  'M2'  ..........  'M10'
-Add other 'Mk'  if there are more than 10 celestial bodies.
 

Runge-Kutta-Fehlberg 4-5 formula
 

-Here, we use the built-in RKF function on the HP48G/GX,  HP49G,  HP50G
-The usage is different from that of RK4 , RK8 ....

-Store the initial x-value in 'X' , the initial y-value in 'Y' , and a program that computes  dy/dx in 'F'
-Store a tolerance value in 'TOL'

-Place the new x-value in level 1 and press the 'rkf'  key
-The y-result will be in level 1 ( and in 'Y' ) and the new x-value will be in 'X'
 
 
 

STACK	INPUT	OUTPUT
Level 1	x	[ y1 ... yn ]

  where           x         =  new x-value
    and     [ y₁ ... y_n ]  =  new y-value

Example1:    y(0) = 1 and  dy/dx = 2 x.y     Evaluate  y(1)  with a tolerance  10^(-10)

  -10  ALOG  'TOL'   STO              or  press the left shit key and TOL
    0   'X'  STO                                 ............................
    1   'Y'  STO

  <<  X  Y  *  2  *   >>   ENTER   'F'    STO

    1    rkf    gives    2.71828182955

-Press  2    rkf   to get  y(2)  ... and so on ...

Example2:       ( the required data are already stored in the DIFF file )

         dy/dx = -yzu         y(0)=1
         dz/dx = x(y+z-u)   z(0)=1
         du/dx = xy-zu       u(0)=2

-Evaluate y(1) , z(1) , u(1)  with  the same tolerance

   -10  ALOG  'TOL'   STO
      0          'X'  STO
[ 1  1  2 ]   'Y'   STO

 << Y  OBJ-> DROP ->  Y Z U
  << Y Z U * * NEG
     Y Z + U - X *
     X Y * Z U * -  3  ->ARRY
  >>  >>           ENTER            'F'  STO

       1    rkf   yields           [ 1   0.258207906452   1.15762398091   0.842178311783 ]    in level1  and in  'Y'  too.

-Press    2    rkf    to get y(2)  z(2)  u(2)  .....
 

Gravitational n-body Problem Solved by the built-in function RKF
 

-Here, we use the Runge-Kutta-Fehlberg 4-5 method on the HP48G/GX,  HP49G,  HP50G

-Store the initial t-value in 'T' , the initial Y-value = [  x₁, y₁,z₁ ............... x_n, y_n,z_n , x'₁, y'₁,z'₁ ............... x'_n, y'_n,z'_n ]  in 'Y'
-Store a tolerance value in 'TOL'  and the n masses into  'M1'  'M2'  .......

-Place the new T-value in level 1 and press the 'rkf'  key
-The y-result will be in level 1 ( and in 'Y' ) and the new T-value will be in 'T'
 
 
 

STACK	INPUT	OUTPUT
Level 1	T	[  x1, y1,z1 ............... xn, yn,zn , x'1, y'1,z'1 ............... x'n, y'n,z'n ]*

  where  [  x₁, y₁,z₁ ............... x_n, y_n,z_n , x'₁, y'₁,z'₁ ............... x'_n, y'_n,z'_n ]*   is the position-velocity vector at  the instant  T

Example:    the same as above:

             A system of three stars              M₁ ( 2 ; 0 ;0 )       M₂ ( 0 ; 4 ; 0 )      M₃ ( 0 ; 0 ; 1 )                     unit = 1 AU
              with  initial velocities                 V₁ ( 0 ; 0.03 ; 0 )  V₂ ( 0 ; 0 ; 0.01 )  V₃ ( -0.02 ; 0 ; 0 )               unit = 1 AU per day
                     and masses                        m₁ = 2  ;  m₂ = 1  ;  m₃ = 3                                                           unit = 1 Solar mass

-Calculate the positions and velocities 10 days later.

-Masses:   2   'M1'  STO   1  'M2'  STO   3  'M3'  STO
-Store 0 into 'T'
-With a tolerance of 10^(-7) , store this value into  'TOL'

-Store   [ 2  0  0  0  4  0  0  0  1  0  0.03  0  0  0  0.01  -0.02  0  0 ]   into  'Y'          ( all these variables are already initialized in the 'DIFF' file )

-Then,   10   rkf   yields  ( rounded to 9 decimals )

          [  1.992077586   0.300333553   0.003673677   0.000661669   3.996080574   0.100603412   -0.194938947   0.001084107   0.997349745
            -0.001550083  0.030038158   0.000706684   0.000132598   -0.000790385   0.010117549   -0.019010811   0.000238023   -0.000510306  ]

-Press   20   rkf   to get the positions & velocities 10 days later....

Notes:

-The directory contains 10 variables 'M1'  'M2'  ..........  'M10'
-Add other 'Mk'  if there are more than 10 celestial bodies.
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4
[3]  F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[4]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[5]  T. Feagin - "A tenth-order Runge-Kutta Method"