Lagrangian Interpolation
12-18-2013, 09:10 PM (This post was last modified: 01-14-2014 01:30 PM by Namir.)
Post: #1
 Namir Senior Member Posts: 950 Joined: Dec 2013
Lagrangian Interpolation
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:
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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)