hp41programs

MinMaxParabola Min-Max Parabola for the HP-41
 

Overview
 

-Given a set of n data points ( x1 , y1 ) ..... ( xn , yn )     ( n > 3 )
  the min-max parabola:  y = a.x2 + b.x + c  minimizes the maximum "error"     h = supi  | yi - a xi2 -  b xi - c |    over the entire set of points.

-The following program uses the "exchange method":

  1°) The min-max parabola is found for the first 4 points.
  2°) The greatest error is computed at all data points for this parabola.
  3°) One of the first 4 points is exchanged with a data point at which the greatest error occurs ( so that the errors are still of alternating sign for this new quadruple )

-The process is repeated until the greatest error occurs at one of the first 4 points ( step 3 is skipped if it happens the first time )
 

Program Listing
 

Data Registers:           •  R00 = n = number of points                          ( Registers R00 thru R2n are to be initialized before executing "MMP" )

                                      •  R01 = x1   •  R03 = x2    ..........  •  R2n-1 = xn                    ( it's necessary to have x1 < x2 < x3 < x4 )
                                      •  R02 = y1   •  R04 = y2    ..........  •  R2n = yn

Flag:  F10 is cleared at the end
Subroutines: /

-Lines 122-160-175-186-196-206-214  are three-byte GTOs
 
 

  01  LBL "MMP"
  02  LBL 01
  03  RCL 06
  04  RCL 02
  05  -
  06  RCL 07
  07  RCL 03
  08  -
  09  *
  10  RCL 08
  11  RCL 04
  12  -
  13  RCL 05
  14  RCL 01
  15  -
  16  *
  17  -
  18  RCL 01
  19  RCL 03
  20  -
  21  RCL 07
  22  -
  23  RCL 05
  24  ST+ Y
  25  RCL 01
  26  -
  27  *
  28  RCL 07
  29  RCL 03
  30  -
  31  *
  32  /
  33  STO M       
  34  RCL 01
  35  RCL 05
  36  +
  37  *
  38  RCL 06
  39  RCL 02
  40  -
  41  RCL 05
  42  RCL 01
  43  -
  44  /
  45  X<>Y
  46  -
  47  STO N
  48  RCL 03
  49  RCL 05
  50  +
  51  *
  52  RCL 03
  53  X^2
  54  RCL 05
  55  X^2
  56  +
  57  RCL M       
  58  *
  59  +
  60  RCL 04
  61  RCL 06
  62  +
  63  X<>Y
  64  -
  65  2
  66  /
  67  STO O
  68  RCL 01
  69  RCL M
  70  *
  71  RCL N
  72  +
  73  RCL 01
  74  *
  75  +
  76  RCL 02
  77  -
  78  CF 10
  79  X<0?
  80  SF 10
  81  ENTER^
  82  ABS
  83  RCL 00
  84  ST+ X
  85   E3
  86  /
  87  8
  88  +
  89  STO P 
  90  STO Q
  91  LBL 02
  92  ISG P
  93  X<0?
  94  GTO 03
  95  CLX
  96  RCL IND P
  97  RCL M
  98  *
  99  RCL N
100  +
101  RCL IND P
102  *
103  RCL O
104  +
105  ISG P
106  RCL IND P
107  -
108  ABS
109  X<=Y?
110  GTO 02
111  LASTX
112  X<>Y
113  RCL P
114  STO Q
115  GTO 02
116  LBL 03
117  X<> Q
118  INT
119  STO L
120  8
121  X=Y?
122  GTO 09 
123  R^
124  X<0?
125  FS?C 10
126  FS? 30
127  SF 10
128  LASTX
129  DSE X
130  STO Y
131  RCL IND X
132  RCL 07
133  X<Y?
134  GTO 04
135  CLX
136  RCL 05
137  X<Y?
138  GTO 05
139  CLX
140  RCL 03
141  X<Y?
142  GTO 06
143  CLX
144  RCL 01
145  X<Y?
146  GTO 07
147  FC? 10
148  GTO 08
149  X<> 03
150  X<> 05
151  X<> 07
152  X<> IND T
153  STO 01
154  RCL 02
155  X<> 04
156  X<> 06
157  X<> 08
158  X<> IND L
159  STO 02
160  GTO 01  
161  LBL 04
162  FS? 10
163  GTO 00
164  X<> 05
165  X<> 03
166  X<> 01
167  X<> IND T
168  STO 07
169  RCL 02
170  X<> IND L
171  X<> 08
172  X<> 06
173  X<> 04
174  STO 02
175  GTO 01
176  LBL 05
177  FC? 10
178  GTO 06
179  LBL 00
180  RCL 07
181  X<> IND T
182  STO 07
183  RCL 08
184  X<> IND L
185  STO 08
186  GTO 01  
187  LBL 06
188  FS? 10
189  GTO 00
190  RCL 05
191  X<> IND T
192  STO 05
193  RCL 06
194  X<> IND L
195  STO 06
196  GTO 01   
197  LBL 07
198  FS? 10
199  GTO 00
200  LBL 08
201  X<> IND T
202  STO 01
203  RCL 02
204  X<> IND L
205  STO 02
206  GTO 01   
207  LBL 00
208  RCL 03
209  X<> IND T
210  STO 03
211  RCL 04
212  X<> IND L
213  STO 04
214  GTO 01   
215  LBL 09
216  R^
217  RCL O
218  RCL N
219  RCL M
220  CF 10
221  CLA
222  END

 
   ( 320 bytes / SIZE 2n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /          +/-h
           Z             /             c
           Y             /             b
           X             /             a

Example:     With the following set of n = 11 points    ( 11  STO 00 )
 
 

     x     0     1     2     3     4     5     6     7     8     9    10
     y   1.2     6    9.7   13.7   16.8     20   21.7    24   24.9   25.7    26

-Store these 22 numbers into

               R01      R03     ..............................................................................................      R21
               R02      R04     ..............................................................................................      R22      respectively

             XEQ "MMP"  yields          a = -0.25          ( in 30 seconds )
                                    RDN          b =  5
                                    RDN          c = 0.975
                                    RDN    +/- h = -0.275

   >>>  whence the equation of the min-max parabola is  y = -0.25 x2 + 5 x + 0.975  and the maximum "error" is  0.275

Notes:

-When the program stops, the order of the points will have probably changed.
-The abscissas can be unequally spaced.
-Due to roundoff errors, it's perhaps not impossible that the algorithm oscillates for ever between 2 very close | h |-values
-Add  VIEW Z  after line 120 to check
-If it ever happens, set flag F21 ( the program will stop each time VIEW Z is executed ) and press XEQ 09 to get the results.
 

Reference:

[1]  F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9