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

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<Y?
130 GTO 03
131 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(x₁)
X	x₁	x₁

Example: dy/dx = x.(y/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 STO02 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 producean 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 R13to 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

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 X<Y?
169 X<>Y
170 RCL 04
171 X<Y?
172 GTO 03
173 R^
174 STO 03
175 RCL IND 16
176 STO 02
177 RCL 05
178 ST+ 01
179 FS? 09
180 GTO 01
181 CLX
182 RCL 01
183 RTN
184 STO 06
185 RCL 05
186 RCL 18
187 X=Y?
188 GTO 01
189 STO 05
190 9
191 STO 07
192 GTO 01
193 END

( 267 bytes / SIZE 035 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₁)
Y	/	y(x₁)
X	x₁	x₁

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

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

"T" ASTO 00

0 STO 01
1 STO02
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

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 X<Y?
148 X<>Y
149 RCL 04
150 X<Y?
151 GTO 03
152 R^
153 STO 03
154 RCL IND 13
155 STO 02
156 RCL 05
157 ST+ 01
158 FS? 09
159 GTO 01
160 CLX
161 RCL 01
162 RTN
163 STO 06
164 RCL 05
165 RCL 15
166 X=Y?
167 GTO 01
168 STO 05
169 9
170 STO 07
171 GTO 01
172 END

( 236 bytes / SIZE 032 )

STACK	INPUTS	OUTPUTS
Z	/	y'(x₁)
Y	/	y(x₁)
X	x₁	x₁

Example: d²y/dx² = -y.( x² + y² )^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 STO02
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

hp41programs

The Bulirsch-Stoer Method for the HP-41