牛頓插值法C語言實作
請注意：在數學推導過程和解釋中的插值法都是用的 n + 1 n+1 n+1個結點來推導，但是程式中是以 n n n個結點來編寫，數學推導中的 n n n是多項式的次數，程式中的 n n n是結點的個數
一、插值法
設函數 y = f ( x ) y=f(x) y=f(x)在區間 [ a , b ] \left[a,b\right] [a,b]上有定義，且已知在點 a ≤ x 0 < x 1 < . . . < x n ≤ b a\le x_0<x_1<...<x_n\le b a≤x0<x1<...<xn≤b上的值 y 0 , y 1 , . . . , y n y_0,y_1,...,y_n y0,y1,...,yn,如果構造出來一個簡單函數 P ( x ) P(x) P(x)，使得 P ( x i ) = y i , i = 0 , 1 , 2 , . . . n P(x_i)=y_i,i=0,1,2,...n P(xi)=yi,i=0,1,2,...n那麼稱 P ( x ) P(x) P(x)為 f ( x ) f(x) f(x)的插值函數，點 x 0 , x 1 . . . x n x_0,x_1...x_n x0,x1...xn為插值節點， [ a , b ] \left[a,b\right] [a,b]為內插補點區間，求內插補點函數 P ( x ) P(x) P(x)的方法稱為插值法。
如果令 P ( x ) P(x) P(x)為多項式，即
P ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n x n P(x)=a_0+a_1x+a_2x^2+...+a_nx^n P(x)=a0+a1x+a2x2+...+anxn
其中 a i a_i ai為實數，就稱 P ( x ) P(x) P(x)為插值多項式，相應的插值法稱為多項式插值。
二、均差
2.1定義
稱 f [ x 0 , x k ] = f ( x k ) − f ( x 0 ) x k − x 0 f \left [x_0,x_k \right] = \frac{f \left( x_k \right) - f \left( x_0 \right)}{x_k-x_0} f[x0,xk]=xk−x0f(xk)−f(x0)為函數 f ( x ) f(x) f(x)關于點 x 0 , x k x_0,x_k x0,xk的 一階均差 ， f [ x 0 , x 1 , x k ] = f [ x 0 , x k ] − f [ x 1 , x k ] x k − x 1 f\left[x_0,x_1,x_k\right]=\frac {f\left[x_0,x_k\right] - f\left[x_1,x_k\right]} {x_k-x_1} f[x0,x1,xk]=xk−x1f[x0,xk]−f[x1,xk]為 f ( x ) f(x) f(x)的二階均差，一般的，稱
f [ x 0 , x 1 , . . . , x k ] = f [ x 0 , x 1 , . . . , x k ] − f [ x 0 , x 1 , . . . , x k − 1 ] x k − x k − 1 ( ∗ ) f \left[ x_0,x_1,...,x_k \right]=\frac {f \left[ x_0,x_1,...,x_k \right] - f \left[ x_0,x_1,...,x_{k-1} \right]} {x_k-x_{k-1}} \ \ \ \ \ \ \ \ \ \ (*) f[x0,x1,...,xk]=xk−xk−1f[x0,x1,...,xk]−f[x0,x1,...,xk−1] (∗)
為 f ( x ) f(x) f(x)的 k k k階均差
2.2 性質
(1) k k k階均差克表示為函數值 f ( x 0 ) , f ( x 1 ) , . . . , f ( x k ) f(x_0),f(x_1),...,f(x_k) f(x0),f(x1),...,f(xk)的線性組合，即
f [ x 0 , x 1 , . . . , x k ] = ∑ j = 0 k f ( x j ) ( x j − x 0 ) . . . ( x j − x j − 1 ) ( x j − x j + 1 ) . . . ( x j − x k ) f \left[ x_0,x_1,...,x_k \right]=\sum_{j=0} ^k \frac {f(x_j)} {(x_j-x_0)...(x_j-x_{j-1})(x_j-x_{j+1})...(x_j-x_k)} f[x0,x1,...,xk]=j=0∑k(xj−x0)...(xj−xj−1)(xj−xj+1)...(xj−xk)f(xj)
這個性質表明均差與結點的排列次序無關，稱為均差的對稱性，即
f [ x 0 , x 1 , . . . , x k ] = f [ x 1 , x 0 , . . . , x k ] = . . . = f [ x 1 , x 2 , . . . , x k , x 0 ] f \left[ x_0,x_1,...,x_k \right]=f \left[ x_1,x_0,...,x_k \right]=...=f \left[ x_1,x_2,...,x_k,x_0 \right] f[x0,x1,...,xk]=f[x1,x0,...,xk]=...=f[x1,x2,...,xk,x0]
(2)由性質1和 ( ∗ ) (*) (∗)式可得
f [ x 0 , x 1 , . . . , x k ] = f [ x 1 , x 2 , . . . , x k ] − f [ x 0 , x 1 , . . . , x k − 1 ] x k − x 0 f \left[ x_0,x_1,...,x_k \right]=\frac {f\left[ x_1,x_2,...,x_k \right] - f\left[ x_0,x_1,...,x_{k-1} \right]} {x_k-x_0} f[x0,x1,...,xk]=xk−x0f[x1,x2,...,xk]−f[x0,x1,...,xk−1]
(3)若 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上存在 n n n階導數，且結點 x 0 , x 1 , . . . , x n ∈ [ a , b ] x_0,x_1,...,x_n \in [a,b] x0,x1,...,xn∈[a,b],則 n n n階均差與導數的關系為
f [ x 0 , x 1 , . . . , x n ] = f ( n ) ( ξ ) n ! , ξ ∈ [ a , b ] . f\left[ x_0,x_1,...,x_n \right]=\frac {f^{(n)}(\xi)}{n!},\ \ \ \xi\in[a,b]. f[x0,x1,...,xn]=n!f(n)(ξ), ξ∈[a,b].
2.3均差表
均差的計算可以使用均差表
| x k x_k xk | f ( x k ) f(x_k) f(xk) | 一階均差 | 二階均差 | 三階均差 | 四階均差 |
|---|---|---|---|---|---|
| x 0 x_0 x0 | f ( x 0 ) f(x_0) f(x0) | ||||
| x 1 x_1 x1 | f ( x 1 ) f(x_1) f(x1) | f [ x 0 , x 1 ] f\left[x_0,x_1 \right] f[x0,x1] | |||
| x 2 x_2 x2 | f ( x 2 ) f(x_2) f(x2) | f [ x 1 , x 2 ] f\left[x_1,x_2 \right] f[x1,x2] | f [ x 0 , x 1 , x 2 ] f\left[ x_0,x_1,x_2 \right] f[x0,x1,x2] | ||
| x 3 x_3 x3 | f ( x 3 ) f(x_3) f(x3) | f [ x 2 , x 3 ] f\left[x_2,x_3 \right] f[x2,x3] | f [ x 1 , x 2 , x 3 ] f\left[ x_1,x_2,x_3 \right] f[x1,x2,x3] | f [ x 0 , x 1 , x 2 , x 3 ] f\left[ x_0,x_1,x_2,x_3 \right] f[x0,x1,x2,x3] | |
| x 4 x_4 x4 | f ( x 4 ) f(x_4) f(x4) | f [ x 3 , x 4 ] f\left[x_3,x_4 \right] f[x3,x4] | f [ x 2 , x 3 , x 4 ] f\left[ x_2,x_3,x_4 \right] f[x2,x3,x4] | f [ x 1 , x 2 , x 3 , x 4 ] f\left[ x_1,x_2,x_3,x_4 \right] f[x1,x2,x3,x4] | f [ x 0 , x 1 , x 2 , x 3 , x 4 ] f\left[x_0,x_1,x_2,x_3,x_4\right] f[x0,x1,x2,x3,x4] |
| … | … | … | … | … | … |
三、牛頓插值法
3.1 內插補點多項式的逐次生成
當 n = 1 n=1 n=1時，有兩個點，這時的插值多項式 P ( x ) P(x) P(x)滿足 P ( x 0 ) = f ( x 0 ) , P ( x 1 ) = f ( x 1 ) P(x_0)=f(x_0),P(x_1)=f(x_1) P(x0)=f(x0),P(x1)=f(x1)，且 P ( x ) = a 0 + a 1 x P(x)=a_0+a_1x P(x)=a0+a1x，帶入得
P 1 ( x ) = f ( x 0 ) + f ( x 1 ) − f ( x 0 ) x 1 − x 0 ( x − x 0 ) P_1(x)=f(x_0)+\frac{f(x_1)-f(x_0)} {x_1-x_0}\left(x-x_0\right) P1(x)=f(x0)+x1−x0f(x1)−f(x0)(x−x0)
這可以看成 P 0 ( x ) = f ( x 0 ) P_0(x)=f(x_0) P0(x)=f(x0)的修正，即
P 1 ( x ) = P 0 ( x ) + a 1 ( x − x 0 ) P_1(x)=P_0(x)+a_1(x-x_0) P1(x)=P0(x)+a1(x−x0)
其中 a 1 = f ( x 1 ) − f ( x 0 ) x 1 − x 0 a_1=\frac{f(x_1)-f(x_0)}{x_1-x_0} a1=x1−x0f(x1)−f(x0)是函數 f ( x ) f(x) f(x)的一階差商。
當 n = 2 n=2 n=2時，有三個已知點，插值多項式 P ( x ) P(x) P(x)滿足 P ( x 0 ) = f ( x 0 ) , P ( x 1 ) = f ( x 1 ) , P ( x 2 ) = f ( x 2 ) P(x_0)=f(x_0),P(x_1)=f(x_1),P(x_2)=f(x_2) P(x0)=f(x0),P(x1)=f(x1),P(x2)=f(x2)
可表示為
P 2 ( x ) = P 1 ( x ) + a 2 ( x − x 0 ) ( x − x 1 ) P_2(x)=P_1(x)+a_2(x-x_0)(x-x_1) P2(x)=P1(x)+a2(x−x0)(x−x1)
顯然 P 2 ( x ) P_2(x) P2(x)滿足上述條件，将 P ( x 2 ) = f ( x 2 ) P(x_2)=f(x_2) P(x2)=f(x2)帶入得到
a 2 = P 2 ( x 2 ) − P 1 ( x 2 ) ( x 2 − x 0 ) ( x 2 − x 1 ) = f ( x 2 ) − f ( x 0 ) x 2 − x 0 − f ( x 1 ) − f ( x 0 ) x 1 − x 0 x 2 − x 1 a_2=\frac {P_2(x_2)-P_1(x_2)}{(x_2-x_0)(x_2-x_1)}=\frac { \frac {f(x_2)-f(x_0)}{x_2-x_0}-\frac{f(x_1)-f(x_0)}{x_1-x_0}} {x_2-x_1} a2=(x2−x0)(x2−x1)P2(x2)−P1(x2)=x2−x1x2−x0f(x2)−f(x0)−x1−x0f(x1)−f(x0)
是以 a 2 = f [ x 2 , x 0 ] − f [ x 1 , x 0 ] x 2 − x 1 a_2=\frac{f[x_2,x_0]-f[x_1,x_0]}{x_2-x_1} a2=x2−x1f[x2,x0]−f[x1,x0]是函數 f ( x ) f(x) f(x)的二階差商
實際上，根據均差定義，将 x x x看成 [ a , b ] [a,b] [a,b]上的一個點，可得
f ( x ) = f ( x 0 ) + f [ x , x 0 ] ( x − x 0 ) , f(x)=f(x_0)+f[x,x_0](x-x_0), f(x)=f(x0)+f[x,x0](x−x0),
f [ x , x 0 ] = f [ x 0 , x 1 ] + f [ x 0 , x 1 , x 2 ] ( x − x 1 ) , f[x,x_0]=f[x_0,x_1]+f[x_0,x_1,x_2](x-x_1), f[x,x0]=f[x0,x1]+f[x0,x1,x2](x−x1),
. . . ... ...
f [ x , x 0 , . . . , x n − 1 ] = f [ x 0 , x 1 , . . . , x n ] + f [ x , x 0 , . . . , x n ] ( x − x n ) . f[x,x_0,...,x_{n-1}]=f[x_0,x_1,...,x_{n}]+f[x,x_0,...,x_n](x-x_n). f[x,x0,...,xn−1]=f[x0,x1,...,xn]+f[x,x0,...,xn](x−xn).
通過疊代得到
f ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) + . . . + f [ x 0 , x 1 , . . . , x n ] ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) + f [ x , x 0 , . . . , x n ] ω n + 1 ( x ) f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...\\+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})+f[x,x_0,...,x_n]\omega_{n+1}(x) f(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x,x0,...,xn]ωn+1(x)
其中
P n ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) + . . . + f [ x 0 , x 1 , . . . , x n ] ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) R ( x ) = f ( x ) − P n ( x ) = f [ x , x 0 , . . . , x n ] ω n + 1 ( x ) \begin{aligned} P_n(x) &= f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1}) \\ R_{\ }(x) &= f(x)-P_n(x)=f[x,x_0,...,x_n]\omega_{n+1}(x) \end{aligned} Pn(x)R (x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)=f(x)−Pn(x)=f[x,x0,...,xn]ωn+1(x)
四、牛頓插值法的C語言實作
4.1 整體思路
我們已知了 n + 1 n+1 n+1個點 ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) . . . ( x n , f ( x n ) ) (x_0,f(x_0)),(x_1,f(x_1))...(x_n,f(x_n)) (x0,f(x0)),(x1,f(x1))...(xn,f(xn))通過牛頓內插補點多項式的形式 P n ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) + . . . + f [ x 0 , x 1 , . . . , x n ] ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) P_n(x) = f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1}) Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)
我們未知的資料隻有差商值，是以我們輸入資料之後将各階均差算出來就可以得到牛頓內插補點多項式。
4.2 初始資料
4.2.1變量說明
double* x; //自變量值的集合
double* y; //函數值的集合
double* f; //插值系數，
int n; //插值節點個數
4.2.2 輸入各個變量
首先我們要确定插值節點的個數，是以先在main函數中輸入變量n的值
printf("輸入要插值節點的個數:");
scanf("%d", &n);
然後我們要為x指針和y指針開辟記憶體空間，x和y指針開n個記憶體空間
x = (double*)malloc(sizeof(double) * n );
y = (double*)malloc(sizeof(double) * n );
再輸入各個節點的值，在這裡我們建立一個data函數，裡面采用for循環
void data(double* x, double* y, int n){
for(int i=0;i<n+1;i++) {
cout << "x[" << i << "]:";
cin >> x[i];
cout << "y[" << i << "]:";
cin >> y[i];
}
}
現在我們就得到了三個變量（x,y,n）并且都已指派
4.3 計算各階均差和資料的存儲
這裡先介紹差商的存儲，因為差商可以表示成一個表，是以用矩陣或二維數組存儲是最友善的，但是考慮到如果有n個結點，生成的矩陣就有 1 + 2 + . . . + ( n − 1 ) = n ( n − 1 ) 2 1+2+...+(n-1)=\frac{n(n-1)}2 1+2+...+(n−1)=2n(n−1)個空間浪費，是以我們使用一維數組來存儲。我們以上表為例
| x k x_k xk | f ( x k ) f(x_k) f(xk) | 一階均差 | 二階均差 | 三階均差 | 四階均差 |
|---|---|---|---|---|---|
| x 0 x_0 x0 | f ( x 0 ) f(x_0) f(x0) | ||||
| x 1 x_1 x1 | f ( x 1 ) f(x_1) f(x1) | f [ x 0 , x 1 ] f\left[x_0,x_1 \right] f[x0,x1] | |||
| x 2 x_2 x2 | f ( x 2 ) f(x_2) f(x2) | f [ x 1 , x 2 ] f\left[x_1,x_2 \right] f[x1,x2] | f [ x 0 , x 1 , x 2 ] f\left[ x_0,x_1,x_2 \right] f[x0,x1,x2] | ||
| x 3 x_3 x3 | f ( x 3 ) f(x_3) f(x3) | f [ x 2 , x 3 ] f\left[x_2,x_3 \right] f[x2,x3] | f [ x 1 , x 2 , x 3 ] f\left[ x_1,x_2,x_3 \right] f[x1,x2,x3] | f [ x 0 , x 1 , x 2 , x 3 ] f\left[ x_0,x_1,x_2,x_3 \right] f[x0,x1,x2,x3] | |
| x 4 x_4 x4 | f ( x 4 ) f(x_4) f(x4) | f [ x 3 , x 4 ] f\left[x_3,x_4 \right] f[x3,x4] | f [ x 2 , x 3 , x 4 ] f\left[ x_2,x_3,x_4 \right] f[x2,x3,x4] | f [ x 1 , x 2 , x 3 , x 4 ] f\left[ x_1,x_2,x_3,x_4 \right] f[x1,x2,x3,x4] | f [ x 0 , x 1 , x 2 , x 3 , x 4 ] f\left[x_0,x_1,x_2,x_3,x_4\right] f[x0,x1,x2,x3,x4] |
| … | … | … | … | … | … |
從一階均差開始， f [ x 0 , x 1 ] f\left[x_0,x_1 \right] f[x0,x1]作為 f [ 0 ] f[\ 0\ ] f[ 0 ]， f [ x 1 , x 2 ] f\left[x_1,x_2 \right] f[x1,x2]作為 f [ 1 ] f[\ 1\ ] f[ 1 ] …一階均差有n-1個數，二階均差有n-2個數 … n-1階均差有1個數。是以一共需要 1 + 2 + . . . + ( n − 1 ) = n ( n − 1 ) 2 1+2+...+(n-1)=\frac {n(n-1)}2 1+2+...+(n−1)=2n(n−1)個空間，是以我們申請 n ( n − 1 ) 2 \frac{n(n-1)}2 2n(n−1)個空間
我們現在有了 n ( n − 1 ) 2 \frac {n(n-1)}2 2n(n−1)個記憶體空間，前n-1個為一階差商，是以我們可以通過一個for循環來算出來所有一階差商
for(j=0,k=0;i<n;j++){
f[k++]=( y[j + 1] - y[j] ) / ( x[j + 1] - x[j])
}
我們根據差商定義有
f [ x k , x k + 1 , x k + 2 ] = f [ x k + 2 , x k + 1 ] − f [ x k + 1 , x k ] x k + 2 − x k f[x_k,x_{k+1},x_{k+2}]=\frac {f[x_{k+2},x_{k+1}]-f[x_{k+1},x_k]} {x_{k+2}-x_k} f[xk,xk+1,xk+2]=xk+2−xkf[xk+2,xk+1]−f[xk+1,xk]
假設 f [ x k , x k + 1 , x k + 2 ] f[x_k,x_{k+1},x_{k+2}] f[xk,xk+1,xk+2]是第 k k k個元素，則 f [ x k + 1 , x k ] f[x_{k+1},x_{k}] f[xk+1,xk]是第 k − n k-n k−n個元素， f [ x k + 2 , x k + 1 ] f[x_{k+2},x_{k+1}] f[xk+2,xk+1]是第 k − n + 1 k-n+1 k−n+1個元素；
假設 f [ x k , x k + 1 , x k + 2 , x k + 3 ] f[x_k,x_{k+1},x_{k+2},x_{k+3}] f[xk,xk+1,xk+2,xk+3]是第 k k k個元素，則 f [ x k , x k + 1 , x k + 2 ] f[x_k,x_{k+1},x_{k+2}] f[xk,xk+1,xk+2]是第 k − n + 1 k-n+1 k−n+1個元素 f [ x k + 1 , x k + 2 , x k + 3 ] f[x_{k+1},x_{k+2},x_{k+3}] f[xk+1,xk+2,xk+3]是第 k − n + 2 k-n+2 k−n+2個元素
歸納得
假設 i + 1 i+1 i+1階差商第一個元素 f [ x k , x k + 1 , . . . x k + i ] f[x_k,x_{k+1},...x_{k+i}] f[xk,xk+1,...xk+i]是第 k k k個元素，則 f [ x k , x k + 1 , . . . x k + i − 1 ] f[x_k,x_{k+1},...x_{k+i-1}] f[xk,xk+1,...xk+i−1]是第 k − n + ( i + 1 ) − 2 = k − n + i − 1 k-n+(i+1)-2=k-n+i-1 k−n+(i+1)−2=k−n+i−1個元素， f [ x k + 1 , x k + 2 , . . . x k + i ] f[x_{k+1},x_{k+2},...x_{k+i}] f[xk+1,xk+2,...xk+i]是第 k − n + i k-n+i k−n+i個元素
是以我們可以通過一個二重for循環來得出所有差商表的數，循環變量 i i i代表差商階數減1，循環變量 j j j代表本階差商的第幾個數， i + 1 i+1 i+1階差商有 n − ( i + 1 ) = n − i − 1 n-(i+1)=n-i-1 n−(i+1)=n−i−1個數，是以循環可以寫出。
void newton(double* x, double* y, double* f, int n)
{
int i = 0, j, k = 0;
for (i = 0; i < n; i++)
for (j = 0; j < n - i; j++) {
if (i == 0)
f[k++] = (y[j + 1] - y[j]) / (x[j + 1] - x[j]);
else {
f[k] = (f[k + i - n] - f[k + i - n - 1]) / (x[j + i + 1] - x[j]);
k++;
}
}
}
4.4 輸出均差表
f f f的前 n − 1 n-1 n−1個元素是一階差商，接下來的 n − 2 n-2 n−2個元素是二階差商，…
是以我們隻需要将 x x x指針中的元素循環一次，每次循環輸出一行資料。
第一行的資料有 x [ 0 ] , y [ 0 ] x[0],y[0] x[0],y[0]，第二行資料有 x [ 1 ] , y [ 1 ] , f [ 0 ] x[1],y[1],f[0] x[1],y[1],f[0]。
是以第 i + 1 i+1 i+1行資料有 i + 2 i+2 i+2個數，有 i i i個差商的資料，是以我們可以用一個 n n n次的循環，先輸出 x [ i ] x[i] x[i]和 y [ i ] y[i] y[i]，再用一個 i i i次循環來輸出差商資料。
void printnew(double* x, double* y, double* f, int n)
{
/*表頭*/
int i, j, k = 0;
printf("差商表:\n");
printf("x\t ");
for (i = 0; i < n; i++)
printf("f(x%d)\t\t", i);
printf("\n");
for (i = 0; i < n; i++)
printf("----------------");
printf("\n");
/*表頭部分結束*/
for (i = 0; i < n; i++) {
printf("%-10g %-10g", x[i], y[i]);
k = i - 1;
for (j = 0; j < i; j++) {
printf(" %-10g", f[k]);
k += n - 2 - j;
}
if (j == i)
printf("\n");
}
}
4.5 根據牛頓插值多項式預測下一個節點的值
先在程式裡得到預測的橫坐标 a a a，然後我們在建立一個有 n n n個元素的數組來存 P n ( x ) P_n(x) Pn(x)除了差商和 f ( x 0 ) f(x_0) f(x0)的部分
P n ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) + . . . + f [ x 0 , x 1 , . . . , x n ] ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) P_n(x) = f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1}) Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)
先算 ( a − x 0 ) (a-x_0) (a−x0)， ( a − x 0 ) ( a − x 1 ) (a-x_0)(a-x_1) (a−x0)(a−x1)，…， ( a − x 0 ) ( a − x 1 ) . . . ( a − x n − 2 ) (a-x_0)(a-x_1)...(a-x_{n-2}) (a−x0)(a−x1)...(a−xn−2)
b[0]=1;
for (i = 0; i < n - 1; i++)
b[i + 1] = b[i] * (a - x[i]);
然後再算 P n ( a ) P_n(a) Pn(a)
for (i = 0; i < n; i++) {
if (i == 0)
a = y[0];
else {
a += b[i] * f[k];
k += n - i;
}
}
最後輸出
是以我們得到了完整的函數體
void getValue(double* x, double* y, double* f, int n)
{
double a, * b;
int i, k = 0;
b = (double*)malloc(sizeof(double) * n);
printf("輸入要要插入節點的X的值:");
cin >> a;
b[0] = 1.0;
for (i = 0; i < n - 1; i++)
b[i + 1] = b[i] * (a - x[i]);
for (i = 0; i < n; i++) {
if (i == 0)
a = y[0];
else {
a += b[i] * f[k];
k += n - i;
}
}
printf("插值節點對應的Y的值:%g\n", a);
}
4.6編寫main函數，完成程式
我們現在有四個函數
void data(double* x, double* y, int n);
void newton(double* x, double* y, double* f, int n);
void printnew(double* x, double* y, double* f, int n);
void newvalue(double* x, double* y, double* f, int n);
他們的功能分别為：輸入資料，計算系數，輸出差商表，預測結點的值。
是以我們可以通過這幾個函數來編寫main函數
int main()
{
int n;
double* x, * y, * f;
printf("輸入要插值節點的個數:");
scanf("%d", &n);
x = (double*)malloc(sizeof(double) * n);
y = (double*)malloc(sizeof(double) * n);
f = (double*)malloc(sizeof(double) * (n - 1) * n / 2);
data(x, y, n);
newton(x, y, f, n - 1);
printnew(x, y, f, n);
do {
newvalue(x, y, f, n);
} while (1);
return 0;
}
4.7測試程式
我們運作程式，用 y = ln x y=\ln x y=lnx為例子
得到的 ln 0.54 ≈ − 0.616143 \ln 0.54\approx-0.616143 ln0.54≈−0.616143
![241 Different Ways to Add Parentheses（C代碼版）[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)