The Museum of HP Calculators

Spectral Analysis for the HP-41

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

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

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

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                 a three-byte GTO
131 X<> Z
132 STO 02
133 RCL 05
134 ST+ 01
135 FS? 09
136 GTO 01                 a three-byte GTO
137 CLX
138 RCL 01
139 RTN
140 STO 04
141 RCL 05
142 RCL 13
143 X=Y?
144 GTO 01                 a three-byte GTO
145 STO 05
146 9
147 STO 06
148 GTO 01                 a three-byte GTO
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"        /          RTN              "T" ASTO 00
X<>Y          X^2
2                  *

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.

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                 a three-byte GTO
173 R^
174 STO 03
175 RCL IND 16
176 STO 02
177 RCL 05
178 ST+ 01
179 FS? 09
180 GTO 01                 a three-byte GTO
181 CLX
182 RCL 01
183 RTN
184 STO 06
185 RCL 05
186 RCL 18
187 X=Y?
188 GTO 01                 a three-byte GTO
189 STO 05
190 9
191 STO 07
192 GTO 01                 a three-byte GTO
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"        +                  RTN             "T" ASTO 00
RCL Z          ST+ X
*                   CHS

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:

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

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                 a three-byte GTO
152 R^
153 STO 03
154 RCL IND 13
155 STO 02
156 RCL 05
157 ST+ 01
158 FS? 09
159 GTO 01                 a three-byte GTO
160 CLX
161 RCL 01
162 RTN
163 STO 06
164 RCL 05
165 RCL 15
166 X=Y?
167 GTO 01                 a three-byte GTO
168 STO 05
169 9
170 STO 07
171 GTO 01                 a three-byte GTO
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           *              "T" ASTO 00
X^2            +               CHS
RCL Y        SQRT       RTN

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:

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

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall