+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (42S) Newton Polynomial Interpolation (/thread-12595.html)
(42S) Newton Polynomial Interpolation - Thomas Klemm - 03-09-2019 04:49 PM
This program uses divided differences to determine 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}})\)
Definition of the Polynomial
The polynomial of degree \(n\) is defined by \(n+1\) points \(P_i=(x_i, y_i)\) for \(i \in \{0, \cdots, n\}\), where all \(x_i\) are different.
This program can be used to enter the points \(P_i\):
Code:
00 { 58-Byte Prgm }
01▸LBL "P+"
02 RCL 00
03 1
04 +
05 STO 00
06 3.002
07 ×
08 STO 01
09 R↓
10▸LBL 00
11 STO IND 01
12 DSE 01
13 GTO 01
14 R↓
15 STO IND 01
16 RTN
17▸LBL 01
18 RCL 01
19 2
20 -
21 R↓
22 RCL- IND ST T
23 RCL 00
24 2
25 ×
26 RCL- 01
27 2
28 ×
29 X<> ST T
30 RCL- IND ST T
31 ÷
32 GTO 00
33 ENDInitialisation
Before you enter points, make sure the counter is reset:
0
STO 00
Example:
Find a quadratic polynomial given these 3 points: \(P_0(-5, 12)\), \(P_1(1, 13)\) and \(P_2(2, 11)\)
0 STO 00
-5 ENTER 12
XEQ "P+"
1 ENTER 13
XEQ "P+"
2 ENTER 11
XEQ "P+"
These are the coefficients of the Newton polynomial:
R03: 12.000000
R05: 0.166667
R07: -0.309524
This leads to the formula:
\(f(x) = 12 + (x+5)(\frac{1}{6} - (x-1)\frac{13}{42})\)
Interpolation of the Polynomial
This program can be used to evaluate the polynomial:
Code:
00 { 37-Byte Prgm }
01▸LBL "f(x)"
02 1.001
03 RCL 00
04 2
05 ×
06 +
07 STO 01
08 CLX
09▸LBL 00
10 RCL IND 01
11 X<>Y
12 DSE 01
13 R↑
14 RCL- IND 01
15 ×
16 +
17 DSE 01
18 GTO 00
19 ENDExample:
Evaluate the polynomial at \(x=0.5\).
0.5
XEQ "f(x)"
13.7679
Adding more data points
An advantage of this method is that the polynomial can be evaluated before all points have been entered.
So you can add more points if the interpolation is not considered sufficient.
Example:
Interpolation of \(\sin(x)\) using values at 0°, 45° and 90°.
0 STO 00
0 ENTER
XEQ "P+"
45 ENTER 0.5 √x
XEQ "P+"
90 ENTER 1
XEQ "P+"
Compare the interpolation at 23° with the correct value:
23
XEQ "f(x)"
0.4132
23
SIN
0.3907
%CH
-5.4289
Add an additional point at 30°:
30 ENTER 0.5
XEQ "P+"
23
XEQ "f(x)"
0.3913
23
SIN
%CH
-0.1397
Add an additional point at 60°:
60 ENTER 0.75 √x
XEQ "P+"
23
XEQ "f(x)"
0.3906
23
SIN
%CH
0.0224
Registers
This is the situation with 3 points:
R00: 3 … number of points
R01: k … index
R02: \(x_0\)
R03: \([y_0]\)
R04: \(x_1\)
R05: \([y_0, y_1]\)
R06: \(x_2\)
R07: \([y_0, y_1, y_2]\)
R08: \([y_1, y_2]\)
R09: \([y_2]\)
When the next point is entered, some of the values at the end are overwritten:
R00: 4 … number of points
R01: k … index
R02: \(x_0\)
R03: \([y_0]\)
R04: \(x_1\)
R05: \([y_0, y_1]\)
R06: \(x_2\)
R07: \([y_0, y_1, y_2]\)
R08: \(x_3\)
R09: \([y_0, y_1, y_2, y_3]\)
R10: \([y_1, y_2, y_3]\)
R11: \([y_2, y_3]\)
R12: \([y_3]\)
Cheers
Thomas
PS: For 3 points I've already posted programs for the HP-67 and HP-15C.
(67) Newton Polynomial Interpolation - Thomas Klemm - 05-11-2022 09:32 PM
Program
This program for the HP-67 works similar to the one for the HP-42S:
Code:
001: 31 25 11 ; LBL A
002: 34 15 ; RCL E
003: 02 ; 2
004: 71 ; *
005: 33 14 ; STO D
006: 35 53 ; Rv
007: 34 15 ; RCL E
008: 01 ; 1
009: 61 ; +
010: 33 15 ; STO E
011: 34 14 ; RCL D
012: 61 ; +
013: 35 33 ; ST I
014: 35 53 ; Rv
015: 31 25 00 ; LBL 0
016: 33 24 ; STO (i)
017: 34 14 ; RCL D
018: 31 81 ; x>0
019: 22 01 ; GTO 1
020: 31 33 ; Rv
021: 35 84 ; DSZ
022: 35 54 ; SPACE
023: 33 24 ; STO (i)
024: 35 22 ; RTN
025: 31 25 01 ; LBL 1
026: 35 53 ; Rv
027: 31 33 ; DSZ
028: 31 33 ; DSZ
029: 31 33 ; DSZ
030: 34 24 ; RCL (i)
031: 51 ; -
032: 34 14 ; RCL D
033: 02 ; 2
034: 51 ; -
035: 33 14 ; STO D
036: 35 24 ; x<>I
037: 34 24 ; RCL (i)
038: 35 52 ; x<>y
039: 35 33 ; ST I
040: 44 ; CLx
041: 61 ; +
042: 35 54 ; R^
043: 35 52 ; x<>y
044: 51 ; -
045: 81 ; /
046: 31 34 ; ISZ
047: 31 34 ; ISZ
048: 22 00 ; GTO 0
049: 31 25 12 ; LBL B
050: 00 ; 0
051: 34 15 ; RCL E
052: 02 ; 2
053: 71 ; *
054: 01 ; 1
055: 51 ; -
056: 35 33 ; ST I
057: 35 53 ; Rv
058: 22 03 ; GTO 3
059: 31 25 02 ; LBL 2
060: 35 52 ; x<>y
061: 34 24 ; RCL (i)
062: 51 ; -
063: 71 ; *
064: 31 34 ; ISZ
065: 31 25 03 ; LBL 3
066: 34 24 ; RCL (i)
067: 61 ; +
068: 31 33 ; DSZ
069: 31 33 ; DSZ
070: 31 33 ; DSZ
071: 22 02 ; GTO 2
072: 35 52 ; x<>y
073: 34 00 ; RCL 0
074: 51 ; -
075: 71 ; *
076: 34 01 ; RCL 1
077: 61 ; +
078: 35 22 ; RTNRegisters
This is the situation with 3 points:
R0: \(x_0\)
R1: \([y_0]\)
R2: \(x_1\)
R3: \([y_0, y_1]\)
R4: \(x_2\)
R5: \([y_0, y_1, y_2]\)
R6: \([y_1, y_2]\)
R7: \([y_2]\)
…
D: ? … used
E: 3 … number of points
I: k … index
When the next point is entered, some of the values at the end are overwritten:
R0: \(x_0\)
R1: \([y_0]\)
R2: \(x_1\)
R3: \([y_0, y_1]\)
R4: \(x_2\)
R5: \([y_0, y_1, y_2]\)
R6: \(x_3\)
R7: \([y_0, y_1, y_2, y_3]\)
R8: \([y_1, y_2, y_3]\)
R9: \([y_2, y_3]\)
R10: \([y_3]\)
…
D: ? … used
E: 4 … number of points
I: k … index
For \( n \) points \( 3 n - 1 \) registers are used.
Thus we can enter at most \( 8 \) points.
Initialisation
Before you enter points, make sure the counter is reset:
0
STO E
Examples
The examples from the previous post work the same way.
Only the index of the registers is shifted by 2:
Code:
00: -5
01: 12
02: 1
03: 0.1666666667
04: 2
05: -0.3095238096
06: -2
07: 11