The Museum of HP Calculators

Spectral Analysis for the HP-41

This program is Copyright © 2007 by Jean-Marc Baillard and is used here by permission.

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

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

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(x1) X x1 x1

Example:     dy/dx =   x.(y/2)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  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.

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(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"        +                  RTN             "T"  ASTO 00
RCL Z          ST+ X
*                   CHS

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:

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

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'(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           *              "T"  ASTO 00
X^2            +               CHS
RCL Y        SQRT       RTN

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