The Bulirsch-Stoer Method for the HP-41

Overview
 

 1°) First-order Differential Equation                         dy/dx  =  f(x,y)
 2°) System of 2 first-order Differential Equations       dy/dx  =  f(x,y,z)  ;  dz/dx = g(x,y,z)
 3°) Second-order Conservative Equation                  d²y/dx² =  f(x,y)
 
 

1°) First-order Differential Equation
 

-The following program solves the differential equation   dy/dx = f(x,y)   with the initial value  y(x₀) = y₀
-It uses an extrapolation to the limit and the modified midpoint formula to compute y(x₀+H)
-Let

                       z₀ = y₀
                       z₁ = z₀ + h.f(x₀,z₀)
                       z₂ = z₀ + 2h.f(x₀+h,z₁)
                       z₃ = z₁ + 2h.f(x₀+2h,z₂)                                     where  h = H/n
                       ..................................
                       z_n = z_n-2 + 2h.f(x₀+(n-1).h,z_n-1)

    y(x₀+H) ~  y_n = (1/2).[ z_n + z_n-1 + h.f(x₀+H,z_n)

-The numbers y_n are computed for n = 2 , 4 , 6 , ..... , 16  ( at most )  and the Lagrange extrapolation polynomial is used
  until the estimated error is smaller than a specified tolerance ( tol )

-If the estimated error is still greater than tol with  n = 16 ,  H is halved ( test line 50 and LBL 00 line 08 )
-On the other hand, if the desired accuracy is satisfied with n < 8 , H is doubled ( lines 21-22 )
-Other ( and much better ) methods for stepsize control may be used ( cf "Numerical Recipes" )

-The Bulirsch-Stoer technique is quite similar to the Romberg integration:
  the modified midpoint formula is second-order only,
  but, for instance, with n = 16 and the deferred approach to the limit, we actually use a 16th-order method!
 

Data Registers:           •  R00 = f name                                                   ( Registers R00 thru R03 are to be initialized before executing "BST" )

                                      •  R01 = x₀      R04 = x₁     R07 = h = H/n             R13 = next to last H
                                      •  R02 = y₀      R05 = H     R08 to R12: temp        R14, R15, ........ , R21 = y-approximations
                                      •  R03 = tol     R06 = n = 2 , 4 , .... , 16

                          >>>>   where tol is a small positive number that defines the desired accuracy

Flags:  F09 - F10
Subroutine:  A program that computes  f(x,y)  assuming x is in X-register and y is in Y-register upon entry.

-Lines 130-136-144-148 are three-byte GTOs
 
 

 
   ( 203 bytes / SIZE 022 )
 
 

 
Example:     dy/dx =   x.(y/2)²           y(0) = 1         Evaluate  y(2)  and  y(2.5)
 
 

 "T"  ASTO 00

  0  STO 01
  1  STO 02   and if we choose  tol = 10 ^-7   STO 03

  2  XEQ "BST"  >>>>        x   =  2                              ( in 1mn49s )                             y(1)  has been computed with  H = 1 , n = 10
                           RDN      y(2) =  2.000000018      the exact result is 2                         y(2)  ------------------------  H = 1 , n = 14

  2.5  R/S    >>>>       x    = 2.5                                   ( in 3mn17s )
                   RDN   y(2.5) = 4.571428682     the last 3 decimals should be  571             H = 1 has been replaced by  0.5  but the tolerance 10 ^-7  cannot be
                                                                      the error is 1.11 10 ^-7                               satisfied with this stepsize. So, H finally became 0.25 ( and n = 14 )

-In fact, the solution is  y = 1/(1-x²/8)   so we'd need to increase tol in R03 and smaller and smaller stepsizes as x get closer to sqrt(8)
 

Notes:

-Do not choose a too small tolerance , it could produce  an infinite loop.
-The initial H-value is always +1 ( or -1 if  x₁ <  x₀ ) - lines 05-06 and 37-39
-Alternatively, place your own H in Y-register  after replacing lines 03-06 by  X<>Y   STO 05   9   STO 06

-Of course, if  x₁-x₀ is smaller than H , H is replaced by  x₁-x₀  ( lines 25-39 )

-The next to last H-value is stored in R13 to avoid losing the benefits of the previous calculations:
  For instance, if you have computed y(1.001)  with  H_initial = 1 , the last H-value will be 0.001
  but if you want to know y(x) for another x-value, the next step will be  H = 1  instead of  0.001
 
 

2°) System of 2 first-order Differential Equations
 

-The same method is employed to solve the system

     dy/dx = f(x,y,z)
     dz/dx = g(x,y,z)       with initial values    y(x₀) = y₀  ,    z(x₀) = z₀
 

Data Registers:           •  R00 = f name                                         ( Registers R00 thru R04 are to be initialized before executing "BST2" )

                                      •  R01 = x₀      R05 = H                              R09 to R17: temp
                                      •  R02 = y₀      R06 = x₁                              R18 = next to last H-value
                                      •  R03 = z₀      R07 = n = 2 , 4 , .... , 16      R19, R20, ........ , R26 = y-approximations
                                      •  R04 = tol     R08 = h = H/n                      R27, R28, ........ , R34 = z-approximations

                             >>>>   where tol is a small positive number that defines the desired accuracy

Flags:  F09 - F10
Subroutine:  A program that computes  f(x,y,z) in Y-register
                                                       and  g(x,y,z) in Z-register   assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Lines 172-180-188-192 are three-byte GTOs
 
 

 
   ( 267 bytes / SIZE 035 )
 
 

 
Example:    d²y/dx² = -2y - 2x.dy/dx    y(0) = 1    y'(0) = 0       Calculate  y(1) and y'(1)        [ the exact solution is  y = exp ( -x² ) ]

  This second-order equation is equivalent to the system

    dy/dx = z                      y(0) = 1
    dz/dx = -2y - 2x.z         z(0) = 0
 
 

 "T"  ASTO 00

  0  STO 01
  1  STO 02
  0  STO 03    and with  tol = 10 ^-7   STO 04

  1  XEQ "BST2"  >>>>        x   =  1                                     ( in 2mn08s )
                RDN                   y(1) =  0.367879446         the last decimal should be a 1
                RDN       z(1) =  y'(1) = -0.735758909        the last 3 decimals should be 882      z-error ~ 27 E-9  <  E-7
 
 

3°) Second-order Conservative Equation
 

-The following program solves the differential equation     d²y/dx² = f(x,y)     with initial values    y(x₀) = y₀  ,   y'(x₀) = y'₀
  where the first derivative does not appear explicitly.

-This kind of equations may be re-written as a system which could be solved by "BST2" but the following method is faster.

-Like "BST", "STM" also uses the deferred approach to the limit but - instead of the modified midpoint method - the Stoermer's rule is employed:
-Let

    d₀ = h.[ y'₀ + (h/2).f(x₀,y₀) ]
    y₁ = y₀ + d₀
    d_k = d_k-1 + h²f(x₀+k.h,y_k)                        d_k = y_k+1 - y_k          k = 1 , 2 , ..... , n-1         h = H/n
    y'_n = d_n-1/h + (h/2).f(x₀+H,y_n)
 

Data Registers:           •  R00 = f name                                         ( Registers R00 thru R04 are to be initialized before executing "STM" )

                                      •  R01 = x₀      R05 = H                              R09 to R14: temp
                                      •  R02 = y₀      R06 = x₁                              R15 = next to last H-value
                                      •  R03 = y'₀     R07 = n = 2 , 4 , .... , 16      R16, R17, ........ , R23 = y-approximations
                                      •  R04 = tol     R08 = h = H/n                      R24, R25, ........ , R31 = y'-approximations

                          >>>>   where tol is a small positive number that defines the desired accuracy

Flags:  F09 - F10
Subroutine:  A program that computes  f(x,y)   assuming x is in X-register and y is in Y-register upon entry.
 

-Lines 151-159-167-171 are three-byte GTOs
 
 

 
    ( 236 bytes / SIZE 032 )
 
 

 
Example:        d²y/dx² = -y.( x² + y² )^1/2         y(0) = 1  ,  y'(0) = 0       Evaluate  y(1) and y(PI)
 
 

"T"  ASTO 00

  0  STO 01
  1  STO 02
  0  STO 03    and with  tol = 10 ^-7   STO 04

  1  XEQ "STM"  >>>>        x   =  1                                  ( in 53 seconds )
                            RDN      y(1) =  0.536630616         the 9 decimals are correct
                            RDN     y'(1) = -0.860171925         the last decimal should be a 7

 PI        R/S         >>>>        x   =  3.141592654                  ( in 2mn24s )
                            RDN    y(PI) = -0.411893053         the 9 decimals are exact
                            RDN    y'(PI) =  1.018399901         the last decimal should be a 3
 

Reference:

[1]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4