hp41programs

Hermite Formula

Osculating Polynomial for the HP-41


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