(42S) Newton Polynomial Interpolation
03-09-2019, 04:49 PM
Post: #1
 Thomas Klemm Senior Member Posts: 2,068 Joined: Dec 2013
(42S) Newton Polynomial Interpolation
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

$$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

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

30 ENTER 0.5
XEQ "P+"

23
XEQ "f(x)"
0.3913

23
SIN
%CH
-0.1397

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.
05-11-2022, 09:32 PM (This post was last modified: 05-11-2022 09:42 PM by Thomas Klemm.)
Post: #2
 Thomas Klemm Senior Member Posts: 2,068 Joined: Dec 2013
(67) Newton Polynomial Interpolation
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
 « Next Oldest | Next Newest »

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