Post Reply 
Lagrangian Interpolation
12-18-2013, 06:04 AM (This post was last modified: 01-17-2014 05:18 AM by Namir.)
Post: #1
Lagrangian Interpolation
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 thir dprompt



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
Find all posts by this user
Quote this message in a reply
03-05-2015, 05:17 AM
Post: #2
RE: Lagrangian Interpolation
Very nice program. Just loaded it and it works very well.

Thank you for posting !
Bob
Find all posts by this user
Quote this message in a reply
03-05-2015, 09:33 PM
Post: #3
RE: Lagrangian Interpolation
(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
Find all posts by this user
Quote this message in a reply
03-07-2015, 11:49 PM
Post: #4
RE: Lagrangian Interpolation
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
Find all posts by this user
Quote this message in a reply
03-08-2015, 06:02 AM (This post was last modified: 03-08-2015 07:59 PM by Thomas Klemm.)
Post: #5
RE: Lagrangian Interpolation
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
Find all posts by this user
Quote this message in a reply
03-08-2015, 07:58 PM (This post was last modified: 03-08-2015 08:22 PM by Thomas Klemm.)
Post: #6
RE: Lagrangian Interpolation
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.


Attached File(s)
.zip  lagrange-interpolation.zip (Size: 370 bytes / Downloads: 1)
Find all posts by this user
Quote this message in a reply
03-09-2015, 03:30 AM
Post: #7
RE: Lagrangian Interpolation
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
Find all posts by this user
Quote this message in a reply
03-09-2015, 03:37 AM
Post: #8
RE: Lagrangian Interpolation
(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
Find all posts by this user
Quote this message in a reply
03-09-2015, 09:36 AM
Post: #9
RE: Lagrangian Interpolation
(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
Find all posts by this user
Quote this message in a reply
03-09-2015, 09:50 PM
Post: #10
RE: Lagrangian Interpolation
Thomas,

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

Regards,
Bob
Find all posts by this user
Quote this message in a reply
03-11-2015, 08:04 AM
Post: #11
RE: Lagrangian Interpolation
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
Find all posts by this user
Quote this message in a reply
03-13-2015, 05:33 AM
Post: #12
RE: Lagrangian Interpolation
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
Find all posts by this user
Quote this message in a reply
04-04-2015, 07:40 PM
Post: #13
RE: Lagrangian Interpolation - Newton Polynomial
Newton Polynomial Version works just great. Thanks Thomas for posting !
Bob
Find all posts by this user
Quote this message in a reply
Post Reply 




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