Optimal Runge-Kutta Method of Order 4 for the HP-41

Overview
 

-This program solves differential equations:   dy/dx = f(x,y)

  systems of 2 differential equations:   dy/dx = f(x,y,z)
                                                        dz/dx = g(x,y,z)

  and systems of 3 differential equations:  dy/dx = f(x,y,z,t)
                                                              dz/dx = g(x,y,z,t)
                                                              dt/dx = h(x,y,z,t)

  by an optimal Runge-Kutta method of order 4 given in reference [1] page 280.

-It is a 4-step method:

    k₁ = h.f(x;y)
    k₂ = h.f(x+b₂.h;y+a_2,1.k₁)
    k₃ = h.f(x+b₃.h;y+a_3,1.k₁+a_3,2.k₂)
    k₄ = h.f(x+b₄.h;y+a_4,1.k₁+a_4,2.k₂+a_4,3.k₃)

then,     y(x+h) = y(x) + c₁.k₁+ ........... + c₄.k₄

-And the coefficients are:

          0.3716151060     0.3716151060
          0.6                     -0.1180444797        0.7180444797
          1                          0.5173871366       -0.5608902997         1.043503163

          0.1474734369     0.3125088197        0.3903768538         0.1496408895
 

Program Listing
 

-The coefficients - except 1 - are stored in registers R16 to R26 ( line 02 to 25 )
-Then the differential equation is solved.

-Clear flags F02 & F03 for a simple differential equation    CF 02  CF 03
-Set F02 & clear F03 for a system of 2 ODEs                   SF 02  CF 03
-Set F03 for a system of 3 differential equations                 SF 03
 

Data Registers:           •  R00 = function name               ( Registers R00 thru R04 or R05 or R06 are to be initialized before executing "RK4OP" )

    1 equation                 •  R01 = x          •  R03 = h                            R05 to R08: temp    R09 to R15: unused   R16 to R26: coeff   R27-R28: unused
                                      •  R02 = y          •  R04 = N
                                      •  R03 = z

    2 equations               •  R01 = x          •  R04 = h                            R06 to R12: temp    R13 to R15: unsed   R16 to R26: coeff    R27-R28: unused
                                      •  R02 = y          •  R05 = N
                                      •  R03 = z

    3 equations               •  R01 = x          •  R05 = h                             R07 to R28: temp   ( R16 to R26: coeff )
                                      •  R02 = y          •  R06 = N
                                      •  R03 = z
                                      •  R04 = t

Flags:   CF 02  CF 03  <->  1 differential equation
              SF 02  CF 03  <->  2 differential equations
                   SF 03         <->  3 differential equations

Subroutine:   CF 02  CF 03  A program that takes x , y in registers  X , Y  and returns  f(x,y)  in X-register
                        SF 02  CF 03  ---------------------  x , y , z  in  X , Y , Z    -----------  f(x,y,z) in Y-register and g(x,y,z) in X-registe
                             SF 03         ---------------------  x , y , z , t  in X , Y , Z , T  -------  f(x,y,z,t) , g(x,y,z,t) , h(x,y,z,t)  in registers  Z , Y , X  respectively.
 
 

-Lines 27-234-239-428-434 are three-byte GTOs
 
 

 
     ( 669 bytes / SIZE 027 or 029 )
 
 

 
Example1:    y' = -2 x.y   y(0) = 1  y'(0) = 0      Evaluate  y(1)
 
 

 
  CF 02  CF 03

 "T"  ASTO 00

  0     STO 01
  1     STO 02    and if we choose  h = 0.1
 0.1   STO 03
 10    STO 04    XEQ "RK4OP"   yields          x = 1                                       ---Execution time = 35s---
                                X<>Y                        y(1) = 0.367879270

-The exact result is 1/e =  0.367879441...

-Press R/S to continue with the same h and N or modify R03-R04
 

Example2:     y(0) = 1 ; y'(0) = 0 ;  y'' + 2.x.y' + 2.y = 0   Evaluate y(1) and y'(1)

-This second order equation is equivalent to the system:             y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    where z = y'
 
 

 
  SF 02  CF 03

 "U"  ASTO 00

  0      STO 01
  1      STO 02
  0      STO 03   and if we choose  h = 0.1
 0.1    STO 04
 10     STO 05    XEQ "RK4OP"   yields       x =  1                                               ---Execution time = 55s---
                                  RDN                      y(1) =  0.367879816
                                  RDN                      z(1) = -0.735759631

-The exact results are  y(1) = 1/e =  0.367879441... & z(1) = -2/e = -0.735758883...

-Press R/S to continue with the same h and N or modify  R04-R05

Example3:

     dy/dx = -y.z.t                    y(0) = 1         Evaluate  y(1) , z(1) , t(1)  with  h = 0.1   ( whence N = 10 )
     dz/dx = x.( y + z - t )         z(0) = 1
     dt/dx = x.y - z.t                  t(0) = 2
 
 

 
  SF 03

  0     STO 01
  1     STO 02  STO 03
  2     STO 04
 0.1   STO 05
 10    STO 06    XEQ "RK4OP"   >>>>      x   =  1                                                       ---Execution time = 88s---
                                                     RDN     y(1) =  0.258209755
                                                     RDN     z(1) =  1.157620123
                                                     RDN     t(1) =  0.842178165

-The exact values are:

      y(1) = 0.258207906                      ( rounded to 9 decimals )
      z(1) = 1.157623981
      t(1) = 0.842178312

-Press R/S to continue with the same h and N or modify   R05-R06
 

Notes:

-When you press R/S after the 1st execution, registers R16 to R26 are unchanged,
  so the constants are still stored in these registers and lines 02 to 29 are no more executed.

-The results are often - but not always - more accurate than with the "classical" 4th-order Runge-Kutta method !
 
 

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  Anthony Ralston - "Runge-Kutta Methods with Minimum Error Bounds"

-In reference [2], you can find coefficients for another optimal method