The Museum of HP Calculators

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

-An osculating polynomial p(x) takes the same values yi  as a given function f at n specified arguments xi     p(xi) = yi = f(xi)   for   i = 1,2,....,n
  but its first derivative also verifies:   p'(xi) = y'i = f '(xi)   for   i = 1,2,....,n
-The following program employs Hermite's formula:

   p(x) = SUMi=1,2,....,n [ 1 - 2.L'i(xi).(x-xi) ].Li2(x).yi  +  SUMi=1,2,....,n (x-xi).Li2(x).y'i

   where    Li(x) = ( x - x1 ) ...... ( x - xi-1 ) ( x - xi+1 ) ...... ( x - xn ) / [ ( xi - x1 ) ...... ( xi - xi-1 ) ( xi - xi+1 ) ( xi - xn ) ]   are the Lagrange's multiplier functions.

-This is the unique polynomial of degree < 2n  that verifies  p(xi) = yi  and  p'(xi) = y'i   for   i = 1,2,....,n
-So, we can hope a better interpolation than with the simple collocation ( Lagrange ) polynomial.
 

Program Listing
 

Data Registers:         R00 = bbb.eee                                                            ( Registers Rbbb thru Reee are to be initialized before executing "HPI"  )

                                 •  Rbbb   = x1   •  Rbb+3 = x2  ..........   •  Ree-2 = xn
                                 •  Rbb+1 = y1   •  Rbb+4 = y2  ..........   •  Ree-1 = yn
                                 •  Rbb+2 = y'1  •  Rbb+5 = y'2  ..........   •  Reee  = y'n
Flags: /
Subroutines: /
 

01  LBL "HPI"
02  CLA
03  STO M
04  X<>Y
05  STO 00
06  STO O
07  LBL 01
08  CLX
09  STO Q
10  RCL 00
11  STO N
12  SIGN
13  LBL 02
14  RCL M
15  RCL IND O
16  RCL IND N
17  ST- Z
18  -
19  X#0?
20  GTO 03
21  X<> Z
22  GTO 04
23  LBL 03
24  1/X
25  ST- Q
26  *
27  *
28  LBL 04
29  ISG N
30  ISG N
31  ISG N
32  GTO 02
33  X^2
34  STO Y
35  RCL M
36  RCL IND O
37  -
38  ST* Y
39  RCL Q
40  ST+ X
41  *
42  R^
43  ST* Y
44  +
45  ISG O
46  RCL IND O
47  *
48  X<>Y
49  ISG O
50  RCL IND O
51  *
52  +
53  ST+ P
54  ISG O
55  GTO 01
56  RCL 00
57  RCL M
58  SIGN
59  X<> P
60  CLA
61  END

( 101 bytes )
 
 
    STACK      INPUTS    OUTPUTS
         Y       bbb.eee      bbb.eee
         X            x        p(x)
         L            /          x

  where  bbb.eee is the control number of the array  ( bbb > 000 )

Example:

  Evaluate  p(6) & p(8)  for the osculating polynomial of degree < 10  that takes the values prescribed below:
 
    xi     1     2     4     7    10
    yi     1     4     6     7     5
    y'i     3     2     1    -1    -2

 For instance, you store these 15 numbers  into

               R01      R04     R07      R10     R13
               R02      R05     R08      R11     R14
               R03      R06     R09      R12     R15          respectively, then

 1.015  ENTER^
     6     XEQ "HPI"  >>>>    p(6) = 7.505337940    ( in 20 seconds )

  RDN   8   R/S  >>>>  p(8) = 5.746750036
 

Notes:

-If you don't want to use synthetic registers, replace M N O P Q by R01  R02  R03  R04  R05  and choose  bbb > 005
-On the other hand, you can replace R00 by synthetic register a  but in this case, the program must not be called as more than a first level subroutine.

Reference:

 [1]  Francis Scheid   "Numerical Analysis"       McGraw-Hill   ISBN 07-055197-9
 
 

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