(42S) Newton Polynomial Interpolation

[Thomas Klemm]

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 END

Initialisation

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 END

Example:

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.

[Thomas Klemm]

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        ; RTN

Instead of P+ and f(x) use the labels A and B.

Registers

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