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
 
 

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