Osculating Polynomial for the HP-41

Overview
 

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

   p(x) = SUM_i=1,2,....,n [ 1 - 2.L'_i(x_i).(x-x_i) ].L_i²(x).y_i  +  SUM_i=1,2,....,n (x-x_i).L_i²(x).y'_i

   where    L_i(x) = ( x - x₁ ) ...... ( x - x_i-1 ) ( x - x_i+1 ) ...... ( x - x_n ) / [ ( x_i - x₁ ) ...... ( x_i - x_i-1 ) ( x_i - x_i+1 ) ( x_i - x_n ) ]   are the Lagrange's multiplier functions.

-This is the unique polynomial of degree < 2n  that verifies  p(x_i) = y_i  and  p'(x_i) = 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   = x₁   •  Rbb+3 = x₂  ..........   •  Ree-2 = x_n
                                 •  Rbb+1 = y₁   •  Rbb+4 = y₂  ..........   •  Ree-1 = y_n
                                 •  Rbb+2 = y'₁  •  Rbb+5 = y'₂  ..........   •  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