Lagrangian Interpolation

[Namir]

Langrangian Interpolation of 3 points

Usage

A program prompts you to enter 3 points and the
interpolated x

Enter y1 and the x1 at first prompt
Enter y2 and the x2 at second prompt
Enter y3 and the x3 at third prompt
Enter xint at fourth prompt
Note calculator displays 1, 2, 3 at the first, second,
and third prompt



B program prompts you to enter the interpolated x

Program calculates and display the value of the interpolated Y.

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 = x1
R01 = x2
R02 = X3
R03 = y1
R04 = y2
R05 = y3
R06 = x
R07 = y
R08 = Product
R09 = used
A = I for X(I)
B = J for X(J)
C = I for Y(I)


Source Code

Code:


1 LBL A
2 1
3 R/S        # enter Y1 and X1
4 STO 0
5 x<>y
6 STO 3
7 2
8 R/S        # enter Y2 and X2
9 STO 1
10 x<>y
11 STO 4
12 3
13 R/S        # enter Y2 and X3
14 STO 2
15 x<>y
16 STO 5
17 R/S        # enter Xint
18 LBL B
19 STO 6
20 2
21 STO B        # set J for X(J)
22 5
23 STO C        # set I for Y(I)
24 0
25 STO 7        # y = 0
26 LBL 0        # outer loop starts
27 RCL C
28 GSB 9
29 STO 8
30 2
31 STO A        # set I fo X(I)
32 LBL 1        # inner loop starts
33 RCL A
34 RCL B
35 x=y?        if I = J the skip iteration
36 GTO 2
37 RCL 6
38 RCL A
39 GSB 9
40 STO 9
41 -
42 RCL B 
43 GSB 9
44 RCL 9
45 -
46 /
47 STO* 8        # update product
48 LBL 2
49 RCL A
50 1
51 -
52 STO A
53 x>0?
54 GTO 1        # inner loop ends
55 x=0?
56 GTO 1
57 RCL 8
58 STO+ 7
59 RCL B
60 1
61 -
62 STO B
63 RCL C
64 1
65 -
66 STO C
67 RCL B
68 x>0?
69 GTO 0        # outer loop ends
70 x=0?
71 GTO 0
72 RCL 7
73 R/S
74 GTO B
75 LBL 9        # INDIRECT ADDRESS
76 STO I
77 RDN
78 RCL (I)
79 RTN

[bshoring]

Very nice program. Just loaded it and it works very well.

Thank you for posting !
Bob

[PedroLeiva]

(03-05-2015 05:17 AM)bshoring Wrote:  Very nice program. Just loaded it and it works very well.

Thank you for posting !
Bob

Namir, thank you very much for your contribution. I guess it is coded for HP- 67. Could you please include a sample calculation for testing the correct entry of the programming steps.
From already thank you for your time. Pedro

[bshoring]

A program for the HP-25 gave this example:

X1=-11 Y1=121 X2=3 Y2=9 X3=2 Y3=4
X=2.5 Y=6.25 etc.

You should get the same results on this program.

Bob

[Thomas Klemm]

You could use the barycentric interpolation formula:
\[
L(x) = \frac{\sum_{j=0}^k \frac{w_j}{x-x_j}y_j}{\sum_{j=0}^k \frac{w_j}{x-x_j}}
\]

This avoids the nested loop at the cost of 3 additional registers for the weights \(w_j\). These weights have to be computed only once for the given data-set.
With only 3 points you could even unroll the loop which would probably speed up the calculation.

Cheers
Thomas

[Thomas Klemm]

This program uses the method mentioned in my previous post.

Memory Map
\(
\begin{matrix}
I & : & x \\
R_0 & : & x_1 \\
R_1 & : & x_2 \\
R_2 & : & x_3 \\
R_3 & : & y_1 \\
R_4 & : & y_2 \\
R_5 & : & y_3 \\
R_6 & : & w_1 \\
R_7 & : & w_2 \\
R_8 & : & w_3 \\
R_{S4} & : & \Sigma x \\
R_{S8} & : & \Sigma xy \\
\end{matrix}
\)

Code:

001  31 25 11  : LBL A
002        01  : 1
003        84  : R/S
004     33 00  : STO 0
005     35 52  : x<>y
006     33 03  : STO 3
007        02  : 2
008        84  : R/S
009     33 01  : STO 1
010     35 52  : x<>y
011     33 04  : STO 4
012        03  : 3
013        84  : R/S
014     33 02  : STO 2
015     35 52  : x<>Y
016     33 05  : STO 5
017        01  : 1
018        42  : CHS
019     33 06  : STO 6
020     33 07  : STO 7
021     33 08  : STO 8
022     34 00  : RCL 0
023     34 01  : RCL 1
024        51  : -
025  33 81 06  : STO/ 6
026  33 81 07  : STO/ 7
027     34 01  : RCL 1
028     34 02  : RCL 2
029        51  : -
030  33 81 07  : STO/ 7
031  33 81 08  : STO/ 8
032     34 02  : RCL 2
033     35 00  : RCL 0
034        51  : -
035  33 81 08  : STO/ 8
036  33 81 06  : STO/ 6
037     31 42  : P<>S
038        84  : R/S
039  31 25 12  : LBL B
040     31 43  : CL REG
041     35 33  : ST I
042     31 42  : P<>S
043     34 03  : RCL 3
044     34 06  : RCL 6
045     35 34  : RC I
046     34 00  : RCL 0
047        51  : -
048        81  : /
049        21  : Σ+
050     34 04  : RCL 4
051     34 07  : RCL 7
052     35 34  : RC I
053     34 01  : RCL 1
054        51  : -
055        81  : /
056        21  : Σ+
057     34 05  : RCL 5
058     34 08  : RCL 8
059     35 34  : RC I
060     34 02  : RCL 2
061        51  : -
062        81  : /
063        21  : Σ+
064     31 42  : P<>S
065     34 08  : RCL 8
066     34 04  : RCL 4
067        81  : /
068     35 22  : RTN

Explanation

Lines 001-016 are the same as in Namir's orginal program.

Lines 017-036 calculate the weights

\(w_j = \frac{1}{\prod_{i=0,i \neq j}^k(x_j-x_i)}\).

However this formula is slightly changed so that we only need the differences \(x_1-x_2\), \(x_2-x_3\) and \(x_3-x_1\):

\(w_1=\frac{-1}{(x_1-x_2)(x_3-x_1)}\)

\(w_2=\frac{-1}{(x_2-x_3)(x_1-x_2)}\)

\(w_3=\frac{-1}{(x_3-x_1)(x_2-x_3)}\)

As you can see the same factor can be used in two of the weights. Nothing fancy it's just to get the sign correct. That's why we start with \(-1\) in lines 017-021.

Lines 043-063 recall \(y_j\) and calculate \(\frac{w_j}{x-x_j}\) and then use \(\Sigma +\) to calculate both \(\sum_{j=0}^k \frac{w_j}{x-x_j}\) and \(\sum_{j=0}^k \frac{w_j}{x-x_j}y_j\) in one single step.

Lines 065-067 finally recall these values and calculate \(\frac{\Sigma xy}{\Sigma x}\).

Cheers
Thomas


For those with a HP-15C here's the corresponding program:

Code:

001 - 42,21,11  LBL A
002 -        1  1
003 -       31  R/S
004 -    44  8  STO 8
005 -       34  x<>y
006 -    44 .1  STO .1
007 -        2  2
008 -       31  R/S
009 -    44  9  STO 9
010 -       34  x<>y
011 -    44 .2  STO .2
012 -        3  3
013 -       31  R/S
014 -    44 .0  STO .0
015 -       34  x<>y
016 -    44 .3  STO .3
017 -        1  1
018 -       16  CHS
019 -    44 .4  STO .4
020 -    44 .5  STO .5
021 -    44 .6  STO .6
022 -    45  8  RCL 8
023 - 45,30, 9  RCL- 9
024 - 44,10,.4  STO/ .4
025 - 44,10,.5  STO/ .5
026 -    45  9  RCL 9
027 - 45,30,.0  RCL- .0
028 - 44,10,.5  STO/ .5
029 - 44,10,.6  STO/ .6
030 -    45 .0  RCL .0
031 - 45,30, 8  RCL- 8
032 - 44,10,.6  STO/ .6
033 - 44,10,.4  STO/ .4
034 -       31  R/S
035 - 42,21,12  LBL B
036 -    44  0  STO 0
037 -    42 32  CLEAR Σ
038 -    45 .1  RCL .1
039 -    45 .4  RCL .4
040 -    45  0  RCL 0
041 - 45,30, 8  RCL- 8
042 -       10  /
043 -       49  Σ+
044 -    45 .2  RCL .2
045 -    45 .5  RCL .5
046 -    45  0  RCL 0
047 - 45,30, 9  RCL- 9
048 -       10  /
049 -       49  Σ+
050 -    45 .3  RCL .3
051 -    45 .6  RCL .6
052 -    45  0  RCL 0
053 - 45,30,.0  RCL- .0
054 -       10  /
055 -       49  Σ+
056 -    45  7  RCL 7
057 - 45,10, 3  RCL/ 3
058 -    43 32  RTN

Added a card of the program that can be used with Jacques Laporte's HP-67 Microcode Simulator.

[bshoring]

Thomas,

I assume pressing Keys A & B are the same as in Namir's, am I right?
However after entering three data pairs, when I enter an X value and press B I am getting an Error message.

For example, if I enter these X & Y values:
100 Enter 1
105 Enter 2
110.25 Enter 3

and then press 1 B, I either get an error message or get a number like 0.49 that doesn't correspond. With Namir's program, pressing 1 B would give me 100.

I have checked and rechecked the program listing and can't find any errors. Am I entering the numbers wrong ?

Thanks,
Bob

[PedroLeiva]

(03-07-2015 11:49 PM)bshoring Wrote:  A program for the HP-25 gave this example:

X1=-11 Y1=121 X2=3 Y2=9 X3=2 Y3=4
X=2.5 Y=6.25 etc.

You should get the same results on this program.

Bob

Thank you very much Bob, I'll try the values in my HP 35s. Pedro

[Thomas Klemm]

(03-09-2015 03:30 AM)bshoring Wrote:  I have checked and rechecked the program listing and can't find any errors. Am I entering the numbers wrong ?

No. That's the problem with this formula:

\(\ell_j(x) = \ell(x)\frac{w_j}{x-x_j}\)

This leads to a division by 0 if \(x=x_j\). The reasoning is that you usually don't want to evaluate the interpolation at the known points.
You could try a value very close to 1 instead.

But I think your point is valid and it made me wonder if we could just use the Lagrange polynomials as Namir did but still unroll the loops. Thus the weights are now calculated in lines 017-031:

\(w_1=\frac{y_1}{(x_1-x_2)(x_3-x_1)}\)

\(w_2=\frac{y_2}{(x_2-x_3)(x_1-x_2)}\)

\(w_3=\frac{y_3}{(x_3-x_1)(x_2-x_3)}\)

Memory Map
\(
\begin{matrix}
I & : & x \\
R_0 & : & x_1 \\
R_1 & : & x_2 \\
R_2 & : & x_3 \\
R_3 & : & w_1 \\
R_4 & : & w_2 \\
R_5 & : & w_3 \\
\end{matrix}
\)

By factoring the sum I was able to calculate it using only the stack in lines 035-057:

\((x-x_2)(x-x_3)w_1+(x-x_1)(x-x_3)w_2+(x-x_1)(x-x_2)w_3=((x-x_2)w_1+(x-x_1)w_2)(x-x_3)+(x-x_1)(x-x_2)w_3\)

Code:

001  31 25 11  : LBL A
002        01  : 1
003        84  : R/S
004     33 00  : STO 0
005     35 52  : x<>y
006     33 03  : STO 3
007        02  : 2
008        84  : R/S
009     33 01  : STO 1
010     35 52  : x<>y
011     33 04  : STO 4
012        03  : 3
013        84  : R/S
014     33 02  : STO 2
015     35 52  : x<>Y
016     33 05  : STO 5
017     34 00  : RCL 0
018     34 01  : RCL 1
019        51  : -
020  33 81 03  : STO/ 3
021  33 81 04  : STO/ 4
022     34 01  : RCL 1
023     34 02  : RCL 2
024        51  : -
025  33 81 04  : STO/ 4
026  33 81 05  : STO/ 5
027     34 02  : RCL 2
028     34 00  : RCL 0
029        51  : -
030  33 81 05  : STO/ 5
031  33 81 03  : STO/ 3
032        84  : R/S
033  31 25 12  : LBL B
034     35 33  : ST I
035     34 04  : RCL 4
036     35 34  : RC I
037     34 00  : RCL 0
038        51  : -
039        71  : *
040     34 05  : RCL 5
041     35 82  : LST x
042        71  : *
043     35 34  : RC I
044     34 01  : RCL 1
045        51  : -
046        71  : *
047     35 52  : x<>y
048     34 03  : RCL 3
049     35 82  : LST x
050        71  : *
051        61  : +
052     35 34  : RC I
053     34 02  : RCL 2
054        51  : -
055        71  : *
056        61  : +
057        42  : CHS
058     35 22  : RTN

It made the program even a little shorter.

Cheers
Thomas

[bshoring]

Thomas,

Thanks for the updated code. This works just great ! And it's fast, too.

Regards,
Bob

[Thomas Klemm]

For the sake of completeness I wrote a program using the Newton Polynomial

\(N(x) = [y_0] + [y_0,y_1](x-x_0) + \cdots + [y_0,\ldots,y_k](x-x_0)(x-x_1)\cdots(x-x_{k-1})\)

where

\([y_0,\ldots,y_j]\)

is the notation for divided differences.

Code:

001  31 25 11  : LBL A
002        01  : 1
003        84  : R/S
004     33 00  : STO 0
005     35 52  : x<>y
006     33 03  : STO 3
007        02  : 2
008        84  : R/S
009     33 01  : STO 1
010     35 52  : x<>y
011     33 04  : STO 4
012        03  : 3
013        84  : R/S
014     33 02  : STO 2
015     35 52  : x<>y
016     34 04  : RCL 4
017        51  : -
018     34 02  : RCL 2
019     34 01  : RCL 1
020        51  : -
021        81  : /
022     34 04  : RCL 4
023     34 03  : RCL 3
024        51  : -
025     34 01  : RCL 1
026     34 00  : RCL 0
027        51  : -
028        81  : /
029     33 04  : STO 4
030        51  : -
031     34 02  : RCL 2
032     34 00  : RCL 0
033        51  : -
034        81  : /
035     34 05  : STO 5
036        84  : R/S
037  31 25 12  : LBL B
038     35 33  : ST I
039     34 01  : RCL 1
040        51  : -
041     34 05  : RCL 5
042        71  : *
043     34 04  : RCL 4
044        61  : +
045     35 34  : RC I
046     34 00  : RCL 0
047        51  : -
048        71  : *
049     34 03  : RCL 3
050        61  : +
051     35 22  : RTN

It is a little shorter and maybe a little faster as well.

Cheers
Thomas


Thanks to its extended register arithmetic the program for the HP-15C is even 10 lines shorter:

Code:

001 - 42,21,11  LBL A       
002 -        1  1           
003 -       31  R/S         
004 -    44  0  STO 0       
005 -       34  x<>y        
006 -    44  3  STO 3       
007 -        2  2           
008 -       31  R/S         
009 -    44  1  STO 1       
010 -       34  x<>y        
011 -    44  4  STO 4       
012 -        3  3           
013 -       31  R/S         
014 -    44  2  STO 2       
015 -       34  x<>y        
016 - 45,30, 4  RCL- 4      
017 -    45  2  RCL 2       
018 - 45,30, 1  RCL- 1      
019 -       10  /           
020 -    45  4  RCL 4       
021 - 45,30, 3  RCL- 3      
022 -    45  1  RCL 1       
023 - 45,30, 0  RCL- 0      
024 -       10  /           
025 -    44  4  STO 4       
026 -       30  -           
027 -    45  2  RCL 2       
028 - 45,30, 0  RCL- 0      
029 -       10  /           
030 -    44  5  STO 5       
031 -       31  R/S         
032 - 42,21,12  LBL B       
033 -    44 25  STO I       
034 - 45,30, 1  RCL- 1      
035 - 45,20, 5  RCLx 5      
036 - 45,40, 4  RCL+ 4      
037 -    45 25  RCL I       
038 - 45,30, 0  RCL- 0      
039 -       20  x           
040 - 45,40, 3  RCL+ 3      
041 -    43 32  RTN

[bshoring]

Thanks, Thomas for posting the listing for the Newton polynomial.
Looks nice. Haven't had time to try it yet. It also looks like it will work on an HP-25 as well as the HP-67.

Regards,
Bob

[bshoring]

Newton Polynomial Version works just great. Thanks Thomas for posting !
Bob

[Thomas Klemm]

(03-13-2015 05:33 AM)bshoring Wrote:  It also looks like it will work on an HP-25 as well as the HP-67.

This program just fits in the 49 steps available for the HP-25:

Code:

01: 01    :    1
02: 74    :    R/S
03: 23 00 :    STO 0
04: 21    :    x<>y
05: 23 03 :    STO 3
06: 02    :    2
07: 74    :    R/S
08: 23 01 :    STO 1
09: 21    :    x<>y
10: 23 04 :    STO 4
11: 03    :    3
12: 74    :    R/S
13: 23 02 :    STO 2
14: 21    :    x<>y
15: 24 04 :    RCL 4
16: 41    :    -
17: 24 02 :    RCL 2
18: 24 01 :    RCL 1
19: 41    :    -
20: 71    :    ÷
21: 24 04 :    RCL 4
22: 24 03 :    RCL 3
23: 41    :    -
24: 24 01 :    RCL 1
25: 24 00 :    RCL 0
26: 41    :    -
27: 71    :    ÷
28: 23 04 :    STO 4
29: 41    :    -
30: 24 02 :    RCL 2
31: 24 00 :    RCL 0
32: 41    :    -
33: 71    :    ÷
34: 23 05 :    STO 5
35: 74    :    R/S
36: 23 06 :    STO 6
37: 24 01 :    RCL 1
38: 41    :    -
39: 24 05 :    RCL 5
40: 61    :    ×
41: 24 04 :    RCL 4
42: 51    :    +
43: 24 06 :    RCL 6
44: 24 00 :    RCL 0
45: 41    :    -
46: 61    :    ×
47: 24 03 :    RCL 3
48: 51    :    +
49: 13 35 :    GTO 35

Definition of the Polynomial

Example:

Find a quadratic polynomial given these 3 points: \(P_1(-5, 12)\), \(P_2(1, 13)\) and \(P_3(2, 11)\)

CLEAR PRGM
R/S
1.0000

12 ENTER -5
R/S
2.0000

13 ENTER 1
R/S
3.0000

11 ENTER 2
R/S


These are the coefficients of the Newton polynomial:

R3: 12.000000
R4:  0.166667
R5: -0.309524

This leads to the formula:

\(f(x) = 12 + (x+5)(\frac{1}{6} - (x-1)\frac{13}{42})\)

Interpolation of the Polynomial

Example:

Evaluate the polynomial at \(x=0.5\).

0.5
R/S

13.7679

Hint: You can just enter new values and hit R/S to evaluate the polynomial.

Cheers
Thomas

[PedroLeiva]

Here is my compilation for HP 25, I have not modified anything, just for your file

To make things clear, in the HP 25 Application Programs manual, Chap 5 Numerical Methods: Linear interpolation, page 85. cites a program that works with two points instead of three. This means that when the function is a straight line I have to use this and when it is a curve the Lagrangian method?

I apologize for my little knowledge in mathematics
Pedro

[Thomas Klemm]

(03-14-2019 03:55 PM)PedroLeiva Wrote:  This means that when the function is a straight line I have to use this and when it is a curve the Lagrangian method?

You could enter a dummy point (e.g. 0 ENTER) as third point.
Just make sure that the x-value is different from the other 2 points.
And then clear register 5:

CLx
STO 5

If you evaluate now the polynomial of the previous example at \(x=2\) you get:

2
R/S

13.1667

Or then you can just shorten the existing program to:

Code:

01: 01    :    1
02: 74    :    R/S
03: 23 00 :    STO 0
04: 21    :    x<>y
05: 23 01 :    STO 1
06: 02    :    2
07: 74    :    R/S
08: 21    :    x<>y
09: 24 01 :    RCL 1
10: 41    :    -
11: 21    :    x<>y
12: 24 00 :    RCL 0
13: 41    :    -
14: 71    :    ÷
15: 23 02 :    STO 2
16: 74    :    R/S
17: 24 00 :    RCL 0
18: 41    :    -
19: 24 02 :    RCL 2
20: 61    :    ×
21: 24 01 :    RCL 1
22: 51    :    +
23: 13 16 :    GTO 16

HTH
Thomas

PS: You could also use the mid-point of the first two points as a third point.
This would lead to a "natural" way to make the entry in register 5 equal to 0.

[PedroLeiva]

Both altrnatives (0 ENTER Xn) plus Clx STO 5, or (12,5 ENTER -2) gives me the same results as yours: 13.1667. So it works as you suggest. So this is a way to use Lagrangian for linear interpolation (1st. order polynomials), only two points.

For 2nd order polynomials we can use your program, just loading three points; so we cannot use it for 3rd. or 4th. order, considering we will need 4 or 5 points. Am I right?
Pedro

[Thomas Klemm]

(03-14-2019 07:22 PM)PedroLeiva Wrote:  so we cannot use it for 3rd. or 4th. order, considering we will need 4 or 5 points. Am I right?

Correct.
In this case you could use (42S) Newton Polynomial Interpolation.
It allows to enter as many points as you like and memory allows.
Given the restrictions of the HP-25 I'm afraid we can't go further than 3 points.

Cheers
Thomas

[PedroLeiva]

I think I was asking HP25 too much. Anyway, I found that the example of the manual for linear interpolation works very well in both cases………….
Example from HP 25 Application Programs, Chapter 5 Numerical Methods/Linear Interpolation: p86
Given
f(7.3)= 1.9879
f(7.4)= 2.0015
find by linear interpolation f(7.37)

A. CLEAR PRGM
R/S
1.0000
1.9879 ENTER 7.3
R/S
2.0000
2.0015 ENTER 7.4
R/S
3.0000
0 ENTER 7
[R/S]
CLx STO 5
7.37 [R/S]  1.9974

B. CLEAR PRGM
R/S
1.0000
1.9879 ENTER 7.3
R/S
2.0000
2.0015 ENTER 7.4
R/S
3.0000
1.9947 [ENTER] 7.35
[R/S]
7.37 [R/S]  1.9974

I appreciate your attention, Pedro

[Thomas Klemm]

(03-14-2019 08:12 PM)PedroLeiva Wrote:  Example from HP 25 Application Programs, Chapter 5 Numerical Methods/Linear Interpolation: p86
Given
f(7.3)= 1.9879
f(7.4)= 2.0015
find by linear interpolation f(7.37)

For those who want to give it a try:

Code:

01: 23 04    : STO 4
02: 24 00    : RCL 0
03: 41       : -
04: 24 03    : RCL 3
05: 61       : *
06: 24 02    : RCL 2
07: 24 04    : RCL 4
08: 41       : -
09: 24 01    : RCL 1
10: 61       : *
11: 51       : +
12: 24 02    : RCL 2
13: 24 00    : RCL 0
14: 41       : -
15: 71       : /
16: 13 00    : GTO 00

References

• HP-25C emulator
• HP-25 Applications Programs
• (25) 55 Applications Programs for Teenix MultiCalculator Emulator and Hardware