Lagrangian Interpolation

[Namir]

This article presents an HP-41C program for Langrangian Interpolation of N points, where N > 1

Program XLGINT


Usage
XEQ XLGINT (same as pressing A)

A Program prompts you to enter N, N points, and the interpolated value of x
B Program prompts you to enter the interpolated value of x

Program calculates and displays the value of the interpolated Y.

Example
Consider the data in the following table:

Code:

X    Y
1    1
5    25
10    100

Using the above data, calculate Y for X = 4. The Steps involved are:

Code:

Step     Task    Command/Input                            Output
1    Start the program.               XEQ "XLGINT"           N?
2    Enter the number of points.       3 [R/S]                       Y1/^X1?
3    Enter the first data point.       1 [ENTER] 1 [R/S]      Y2/^X2?
4    Enter the second data point.  25 [ENTER] 5 [R/S]  Y3X3?
5    Enter the third data point.       100 [ENTER] 10[R/S]  READY
6    Start the interpolation.       [B]                             X?
7    Enter the interpolated value of X      4 [R/S]                Y=16.00000

Algorithm

Code:

INPUT N, array X(1..N), Y(1..N), and Xint
Yint = 0
FOR I = 1 TO N
  Product = Y(I)
  FOR J = 1 to N
    IF I <> J THEN
      Product = Product * (Xint - X(J)) / (X(I) - X(J))
    ENDIF
  NEXT J
  Yint = Yint + Product
NEXT I
Show Yint

Memory Map
R00 = N
R01 = Xint
R02 = Yint
R03 = Product
R05 = J for X(J)
R06 = I for X(I)
R07 = I for Y(I)
R10 = X(1)
R11 = X(2)
...
R10+N-1 = X(N)
R10+N = Y(1)
R11+N = Y(2)
...
R10+2N-1 = Y(N)

Source Code
The source code for the HP-41C program appears below. Please note the following:
1) Text appearing in a pair of double quotes represents characters in the Alpha register.
2_ The |- characters represent the single append character for the Alpha register.
3)
The blank lines are intentionally inserted to separate logical blocks of commands:

Code:

LBL XLGINT    
1 LBL XLGINT    
2 LBL A    
3 "N?"    
4 PROMPT
5 INT    
6 STO 00    
7 2    
8 X>Y?    # is N < 2? then re-prompt N
9 GTO A    
10 X<>Y    
11 1E3    
12 /    
13 1    
14 +    
15 STO 05    # Initializes indices for storage
16 10    
17 STO 06    # I for X(I)
18 RCL 00    
19 +    
20 STO 07    # I for Y(I)  
21 LBL 00    # start input loop
22 FIX 0
23 CF 29    
24 "Y"    
25 ARCL 05    
26 "|-^X"    
27 ARCL 05    
28 "|-?"    
29 FIX 5
30 SF 29    
31 PROMPT    
32 STO IND 06    
33 X<>Y    
34 STO IND 07    
35 1    
36 ST+ 06    
37 ST+ 07    
38 ISG 05    
39 GTO 00    # end input loop
40 LBL B    
41 "X?"    
42 PROMPT    
43 STO 01     
44 XEQ 09    
45 STO 06    # set I for X(I) 
46 XEQ 10    
47 STO 07    # set I for Y(I)  
48 0    
49 STO 02    # y = 0  
50 LBL 01    # outer loop starts 
51 RCL IND 07    
52 STO 03    
53 XEQ 09    
54 STO 05    # set J for X(J)
55 LBL 02    # inner loop starts
56 RCL 06    
57 RCL 05    
58 X=Y?    # if I = J the skip iteration
59 GTO 03      
60 RCL 01    
61 RCL IND 05    
62 -    
63 RCL IND 06    
64 RCL IND 05    
65 -    
66 /    
67 ST* 03    # update product   
68 LBL 03    
69 ISG 05    
70 GTO 02    # inner loop ends 
71 RCL 03    
72 ST+ 02     
73 1    
74 ST+ 07    # update index of Y() by simply adding 1
75 ISG 06    
76 GTO 01    # outer loop ends
77 "Y = "    
78 ARCL 02    
79 PROMPT    
80 GTO B        
81 LBL 09    # subroutine to calculate index range for
82 10    # accessing array X()
83 STO Y    
84 RCL 00    
85 +    
86 1    
87 -    
88 1E3    
89 /    
90 +    
91 RTN    
92 LBL 10    # subroutine to calculate starting index # for accessing array Y()
93 10    
94 RCL 00    
95 +    
96 RTN