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: 824 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:
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
 « Next Oldest | Next Newest »

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