 Lagrangian Interpolation - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP-41C Software Library (/forum-11.html) +--- Thread: Lagrangian Interpolation (/thread-158.html) Lagrangian Interpolation - Namir - 12-18-2013 09:10 PM 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```