Bulirsch-Stoer

# 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                  d2y/dx2 =  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(x0) = y0
-It uses an extrapolation to the limit and the modified midpoint formula to compute y(x0+H)
-Let

z0 = y0
z1 = z0 + h.f(x0,z0)
z2 = z0 + 2h.f(x0+h,z1)
z3 = z1 + 2h.f(x0+2h,z2)                                     where  h = H/n
..................................
zn = zn-2 + 2h.f(x0+(n-1).h,zn-1)

y(x0+H) ~  yn = (1/2).[ zn + zn-1 + h.f(x0+H,zn)

-The numbers yn 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 = x0      R04 = x1     R07 = h = H/n             R13 = next to last H
•  R02 = y0      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

 01  LBL "BST"   02  STO 04   03  9   04  STO 06   05  SIGN   06  STO 05   07  GTO 01   08  LBL 00   09  2   10  ST/ 05   11  GTO 02   12  LBL 01   13  RCL 02   14  RCL 01   15  XEQ IND 00   16  STO 11   17  6   18  RCL 06   19  X>Y?   20  GTO 02   21  RCL 05   22  ST+ 05   23  LBL 02   24  SF 09   25  RCL 05   26  STO 13   27  ABS   28  RCL 04   29  RCL 01   30  - 31  ABS   32  X<=Y?   33  CF 09   34  X>Y?   35  X<>Y   36  LASTX   37  SIGN   38  *   39  STO 05   40  SF 10   41  CLX   42  STO 06          43  13   44  STO 12   45  LBL 03   46  2   47  ST+ 06   48  18   49  RCL 06   50  X=Y?   51  GTO 00   52  DSE X   53   E3   54  /   55  ISG X   56  STO 08   57  RCL 05   58  RCL 06   59  /   60  STO 07 61  RCL 11   62  *   63  RCL 02   64  STO 09   65  +   66  ISG 12   67  LBL 04   68  STO 10   69  RCL 08   70  INT   71  RCL 07   72  *   73  RCL 01   74  +   75  XEQ IND 00   76  RCL 07   77  ST+ X   78  *   79  RCL 10   80  X<> 09   81  +   82  ISG 08   83  GTO 04   84  STO 10   85  RCL 01   86  RCL 05   87  +   88  XEQ IND 00   89  RCL 07   90  * 91  RCL 09   92  +   93  RCL 10   94  +   95  2   96  /   97  STO IND 12   98  FS?C 10   99  GTO 03 100  RCL 06 101  LASTX 102  ST- Y 103  / 104  STO 10 105  LBL 05 106  RCL 12 107  RCL 10 108  - 109  RCL IND 12 110  ENTER^ 111  X<> IND Z 112  STO Z 113  - 114  RCL 06 115  X^2 116  ST* Y 117  RCL 10 118  ST+ X 119  X^2 120  - 121  / 122  + 123  STO IND 12 124  DSE 10 125  GTO 05 126  LASTX 127  ABS 128  RCL 03 129  X Z 132  STO 02 133  RCL 05 134  ST+ 01 135  FS? 09 136  GTO 01 137  CLX 138  RCL 01 139  RTN 140  STO 04 141  RCL 05 142  RCL 13 143  X=Y? 144  GTO 01 145  STO 05 146  9 147  STO 06 148  GTO 01 149  END

( 203 bytes / SIZE 022 )

 STACK INPUTS OUTPUTS Y / y(x1) X x1 x1

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

 LBL "T"   X<>Y    2   /   X^2   *   RTN

"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-x2/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  x1 <  x0 ) - 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  x1-x0 is smaller than H , H is replaced by  x1-x0  ( 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  Hinitial = 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(x0) = y0  ,    z(x0) = z0

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

•  R01 = x0      R05 = H                              R09 to R17: temp
•  R02 = y0      R06 = x1                              R18 = next to last H-value
•  R03 = z0      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

 01  LBL "BST2"   02  STO 06   03  9   04  STO 07   05  SIGN   06  STO 05   07  GTO 01   08  LBL 00   09  2   10  ST/ 05   11  GTO 02   12  LBL 01   13  RCL 03   14  RCL 02   15  RCL 01   16  XEQ IND 00   17  STO 15   18  X<>Y   19  STO 14   20  6   21  RCL 07   22  X>Y?   23  GTO 02   24  RCL 05   25  ST+ 05   26  LBL 02   27  SF 09   28  RCL 05   29  STO 18   30  ABS   31  RCL 06   32  RCL 01   33  -   34  ABS   35  X<=Y?   36  CF 09   37  X>Y?   38  X<>Y   39  LASTX 40  SIGN   41  *   42  STO 05   43  SF 10   44  CLX   45  STO 07   46  18   47  STO 16   48  26   49  STO 17   50  LBL 03   51  2   52  ST+ 07   53  18   54  RCL 07           55  X=Y?   56  GTO 00   57  1   58  ST+ 16   59  ST+ 17   60  -   61   E3   62  /   63  ISG X   64  STO 09   65  RCL 05   66  RCL 07   67  /   68  STO 08   69  RCL 14   70  *   71  RCL 02   72  STO 10   73  +   74  STO 11   75  RCL 08   76  RCL 15   77  *   78  RCL 03 79  STO 12   80  +   81  STO 13   82  LBL 04   83  RCL 13   84  RCL 11   85  RCL 09   86  INT   87  RCL 08   88  *   89  RCL 01   90  +   91  XEQ IND 00   92  RCL 08   93  ST+ X   94  ST* Z   95  *   96  RCL 12   97  +   98  X<> 13   99  STO 12 100  CLX 101  RCL 10 102  + 103  X<> 11 104  STO 10 105  ISG 09 106  GTO 04 107  RCL 13 108  RCL 11 109  RCL 01 110  RCL 05 111  + 112  XEQ IND 00 113  RCL 08 114  ST* Z 115  * 116  RCL 12 117  + 118  RCL 13 119  + 120  STO IND 17 121  CLX 122  RCL 10 123  + 124  RCL 11 125  + 126  2 127  ST/ IND 17 128  / 129  STO IND 16 130  FS?C 10 131  GTO 03 132  16.017 133  STO 09 134  LBL 05 135  RCL 07 136  2 137  ST- Y 138  / 139  STO 10 140  LBL 06 141  RCL IND 09 142  STO 11 143  RCL 10 144  - 145  RCL IND 11 146  ENTER^ 147  X<> IND Z 148  STO Z 149  - 150  RCL 07 151  X^2 152  ST* Y 153  RCL 10 154  ST+ X 155  X^2 156  - 157  / 158  + 159  STO IND 11 160  DSE 10 161  GTO 06 162  LASTX 163  ABS 164  X<> 12 165  ISG 09 166  GTO 05 167  RCL 12 168  XY 170  RCL 04 171  X

( 267 bytes / SIZE 035 )

 STACK INPUTS OUTPUTS Z / z(x1) Y / y(x1) X x1 x1

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

This second-order equation is equivalent to the system

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

 LBL "T"    RCL Z    *    +   ST+ X   CHS   RTN

"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     d2y/dx2 = f(x,y)     with initial values    y(x0) = y0  ,   y'(x0) = y'0
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

d0 = h.[ y'0 + (h/2).f(x0,y0) ]
y1 = y0 + d0
dk = dk-1 + h2f(x0+k.h,yk)                        dk = yk+1 - yk          k = 1 , 2 , ..... , n-1         h = H/n
y'n = dn-1/h + (h/2).f(x0+H,yn)

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

•  R01 = x0      R05 = H                              R09 to R14: temp
•  R02 = y0      R06 = x1                              R15 = next to last H-value
•  R03 = y'0     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

 01  LBL "STM"   02  STO 06   03  9   04  STO 07   05  SIGN   06  STO 05   07  GTO 01   08  LBL 00   09  2   10  ST/ 05   11  GTO 02   12  LBL 01   13  RCL 02   14  RCL 01   15  XEQ IND 00   16  2   17  /   18  STO 12   19  6   20  RCL 07   21  X>Y?   22  GTO 02   23  RCL 05   24  ST+ 05   25  LBL 02   26  SF 09   27  RCL 05   28  STO 15   29  ABS   30  RCL 06   31  RCL 01   32  -   33  ABS   34  X<=Y?   35  CF 09 36  X>Y?   37  X<>Y   38  LASTX   39  SIGN   40  *   41  STO 05   42  SF 10   43  CLX   44  STO 07   45  15   46  STO 13           47  23   48  STO 14   49  LBL 03   50  2   51  ST+ 07   52  18   53  RCL 07   54  X=Y?   55  GTO 00   56  1   57  ST+ 13   58  ST+ 14   59  -   60   E3   61  /   62  ISG X   63  STO 09   64  RCL 05   65  RCL 07   66  /   67  STO 08   68  RCL 12   69  *   70  RCL 03 71  +   72  RCL 08   73  *   74  STO 11   75  RCL 02   76  +   77  LBL 04   78  STO 10   79  RCL 09   80  INT   81  RCL 08   82  *   83  RCL 01   84  +   85  XEQ IND 00   86  RCL 08   87  X^2   88  *   89  RCL 11   90  +   91  STO 11   92  RCL 10   93  +   94  ISG 09   95  GTO 04   96  STO IND 13   97  RCL 01   98  RCL 05   99  + 100  XEQ IND 00 101  2 102  / 103  RCL 11 104  RCL 08 105  ST* Z 106  / 107  + 108  STO IND 14 109  FS?C 10 110  GTO 03 111  13.014 112  STO 09 113  LBL 05 114  RCL 07 115  2 116  ST- Y 117  / 118  STO 10 119  LBL 06 120  RCL IND 09 121  STO 11 122  RCL 10 123  - 124  RCL IND 11 125  ENTER^ 126  X<> IND Z 127  STO Z 128  - 129  RCL 07 130  X^2 131  ST* Y 132  RCL 10 133  ST+ X 134  X^2 135  - 136  / 137  + 138  STO IND 11 139  DSE 10 140  GTO 06 141  LASTX 142  ABS 143  X<> 08 144  ISG 09 145  GTO 05 146  RCL 08 147  XY 149  RCL 04 150  X

( 236 bytes / SIZE 032 )

 STACK INPUTS OUTPUTS Z / y'(x1) Y / y(x1) X x1 x1

Example:        d2y/dx2 = -y.( x2 + y2 )1/2         y(0) = 1  ,  y'(0) = 0       Evaluate  y(1) and y(PI)

 LBL "T"   X^2    RCL Y   X^2   +   SQRT   *   CHS   RTN

"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:

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