全域多項(xiàng)式插值指的是在整個插值區(qū)域內(nèi)形成一個多項(xiàng)式函數(shù)作為插值函數(shù)。關(guān)于多項(xiàng)式插值的基本知識，見“計(jì)算基本理論”。

在單項(xiàng)式基插值和牛頓插值形成的表達(dá)式中，求該表達(dá)式在某一點(diǎn)處的值使用的Horner嵌套算法啊，見"Horner嵌套算法"。

1. 單項(xiàng)式(Monomial)基插值

1)插值函數(shù)基　　單項(xiàng)式基插值采用的函數(shù)基是最簡單的單項(xiàng)式：$$\phi_j(t)=t^{j-1}, j=1,2,...n;\quad f(t)=p_{n-1}(t)=x_1+x_2t+x_3t^2+...x_nt^{n-1}=\sum\limits_{j=1}^nx_jt^{j-1}$$　　所要求解的系數(shù)即為單項(xiàng)式系數(shù) $x_1,x_2,...x_n$ ，在這里仍然采用1,2,...n的下標(biāo)記號而不采用和單項(xiàng)式指數(shù)對應(yīng)的0,1,2,...,n-1的下標(biāo)僅僅是出于和前后討論一致的需要。

2)疊加系數(shù)

單項(xiàng)式基插值采用單項(xiàng)式函數(shù)基，若有m個離散數(shù)據(jù)點(diǎn)需要插值，設(shè)使用n項(xiàng)單項(xiàng)式基底：

$$x_1+t_1x_2+t_1^2x_3+...+t_1^{n-1}x_n=y_1\\ x_1+t_2x_2+t_2^2x_3+...+t_2^{n-1}x_n=y_2\\ ......? ?......? ?......? ?......? ?......? ?......\\ x_1+t_mx_2+t_m^2x_3+...+t_m^{n-1}x_n=y_m$$　　系數(shù)矩陣為一 $m\times n$ 的矩陣($m\leq n$)，范德蒙(Vandermonde)矩陣：

$$\begin{bmatrix}1&t_1&t_1^2&...&t_1^{n-1}\\1&t_2&t_2^2&...&t_2^{n-1}\\...&...&...&...&...\\1&t_n&t_n^2&...&t_n^{n-1}\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$　　根據(jù)計(jì)算基本理論中的討論，多項(xiàng)式插值的函數(shù)基一定線性無關(guān)，且只要離散數(shù)據(jù)點(diǎn)兩兩不同，所構(gòu)成的矩陣行也一定線性無關(guān)，這保證了矩陣一定行滿秩。此時當(dāng)且僅當(dāng)m=n時，疊加系數(shù)有且僅有一組解。因此，n=插值基底的個數(shù)=離散數(shù)據(jù)點(diǎn)的個數(shù)=多項(xiàng)式次數(shù)+1。

3)問題條件和算法復(fù)雜度

生成范德蒙矩陣的復(fù)雜度大致在 $O(n^2)$ 量級；由于范德蒙矩陣并沒有什么好的性質(zhì)，既沒有稀疏性，也沒有對稱性，因此只能使用標(biāo)準(zhǔn)的稠密矩陣分解(如LU)來解決，復(fù)雜度在 $O(n^3)$ 量級。因此，問題的復(fù)雜度在 $O(n^3)$ 量級。

范德蒙矩陣存在的問題是，當(dāng)矩陣的維數(shù)越來越高的時候，解線性方程組時問題將越來越病態(tài)，條件數(shù)越來越大；從另一個角度來說，單項(xiàng)式基底本身趨勢相近，次數(shù)增大時將越來越趨于平行(見下圖)。這將造成隨著數(shù)據(jù)點(diǎn)的增加，確定的疊加系數(shù)的不確定度越來越大，因此雖然單項(xiàng)式基很簡單，進(jìn)行插值時卻往往不用這一方法。如果仍然采用單項(xiàng)式基底，有時也會進(jìn)行兩種可以改善以上問題的操作：平移(shifting)和縮放(scaling)，即將 $((t-c)/d)^{j-1}$ 作為基底。常見的平移和縮放將所有數(shù)據(jù)點(diǎn)通過線性變換轉(zhuǎn)移到區(qū)間[-1,1]中，即：$c=(t_1+t_n)/2,d=(t_n-t_1)/2$ 。

4)算法實(shí)現(xiàn)

使用MATLAB實(shí)現(xiàn)單項(xiàng)式插值代碼如下：

function [ polyfunc, vecx, condition ] = MonoPI( vect, vecy, shift, scale )

% 計(jì)算單項(xiàng)式型插值多項(xiàng)式系數(shù)

% 輸入四個參數(shù)：插值點(diǎn)行向量vect，插值點(diǎn)函數(shù)值行向量vecy，平移shift，壓限scale；

% 輸出兩個參數(shù)：插值多項(xiàng)式各項(xiàng)系數(shù)行向量vecx，矩陣條件數(shù)condition；

% 設(shè)置缺省值：若只輸入兩個參數(shù)，則不平移不縮放

if nargin==2

shift = 0; scale = 1;

end

% 求解系數(shù)

vecsize = size(vect, 2);

basis = (vect - shift * ones(1, vecsize))/scale; % 確定基底在各個數(shù)據(jù)點(diǎn)的取值向量basis

Mat = vander(basis); condition = cond(Mat); % 用vander命令生成basis的范德蒙矩陣并求條件數(shù)

[L, U] = lu(Mat); vecx = (U\(L\vecy.‘)).‘; vecx = fliplr(vecx); % 標(biāo)準(zhǔn)lu分解解矩陣方程

% 生成句柄函數(shù)polyfunc

syms t;

monomial = (t - shift)/scale; vecsize = size(vecx, 2); funeval = vecx(vecsize);

for i = vecsize:-1:2 % 生成函數(shù)的算法采用Horner算法提高效率

funeval = vecx(i - 1) + monomial*funeval;

end

polyfunc = matlabFunction(funeval, ‘Vars‘, t);

end

比如對于點(diǎn)：$(-2,-27),(0,-1),(1,0)$ 它具有唯一的二次插值多項(xiàng)式：$p_2(t)=-1+5t-4t^2$ 。調(diào)用以上代碼：

% 命令行輸入

[polyfunc, vecx, condition] = MonoPI(vect, vecy)

% 命令行輸出

polyfunc =

包含以下值的 function_handle:

@(t)-t.*(t.*4.0-5.0)-1.0

vecx =

-1 5 -4

condition =

6.0809

和預(yù)期完全一致。

2. 拉格朗日(Lagrange)插值

1)插值函數(shù)基

拉格朗日插值采用的是一種設(shè)計(jì)巧妙的多項(xiàng)式基，每個基底都是n-1次多項(xiàng)式，而每個基底函數(shù)當(dāng)且僅當(dāng)在第i個數(shù)據(jù)點(diǎn)處取1，在其余數(shù)據(jù)點(diǎn)均為零。這個多項(xiàng)式基是這樣設(shè)計(jì)的：

$$l_j(t)=\frac{(t-t_1)(t-t_2)...(t-t_{j-1})(t-t_{j+1})...(t-t_n)}{(t_j-t_1)(t_j-t_2)...(t_j-t_{j-1})(t_j-t_{j+1})...(t_j-t_n)}=\frac{\prod\limits_{k=1,k\neq j}^n(t-t_k)}{\prod\limits_{k=1, k\neq j}^n(t_j-t_k)}$$　　因此就有：

$$l_j(t_i)=\delta_{ij}, i,j=1,2,...n $$　　其中，$\delta$ 為克羅內(nèi)克(Kronecker)記號，當(dāng)兩個下標(biāo)相等時為1，否則為零；也可以將 $\delta_{ij}$ 理解為一個二階張量，即單位矩陣。只要將各個$t_i$ 帶入定義式，上式是很容易驗(yàn)證的。這意味著拉格朗日插值的疊加系數(shù)的求解將會產(chǎn)生很好的性質(zhì)，即：

2)疊加系數(shù)

需要求解的插值函數(shù)即：$f(t)=\sum\limits_{k=1}^nx_kl_k(t)$ ，而又已知：

$$l_1(t_1)x_1+l_2(t_1)x_2+...+l_n(t_1)x_n=y_1$$

$$l_1(t_2)x_1+l_2(t_2)x_2+...+l_n(t_2)x_n=y_2$$

$$... ... ... ... ...$$

$$l_1(t_n)x_1+l_2(t_n)x_2+...+l_n(t_n)x_n=y_n$$　　寫成矩陣形式就是：

$$\begin{bmatrix}l_1(t_1)&l_2(t_1)&...&l_n(t_1)\\l_1(t_2)&l_2(t_2)&...&l_n(t_2)\\...&...&...&...\\l_1(t_n)&l_2(t_n)&...&l_n(t_n)\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}1&0&..&0\\0&1&..&0\\..&..&..&..\\0&0&..&1\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=I\boldsymbol{x}=\begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$

其矩陣即單位矩陣，因此直接得出 $x_i=y_i$ ，$f(t)=p_{n-1}(t)=y_1l_1(t)+y_2l_2(t)+...+y_nl_n(t)=\sum\limits_{k=1}^ny_kl_k(t)$

3)問題條件和算法復(fù)雜度

拉格朗日插值生成的系數(shù)矩陣為單位矩陣，完全不存在條件病態(tài)的問題，只需要將各個數(shù)據(jù)點(diǎn)的取值作為系數(shù)即可。同樣地，求解系數(shù)也將不存在任何復(fù)雜度。

但是，作為代價的是，函數(shù)求值開銷較大。Horner嵌套算法可以用于單項(xiàng)式和牛頓插值表達(dá)式的求值，將總運(yùn)算量大致控制在n次浮點(diǎn)數(shù)加法和n次浮點(diǎn)數(shù)乘法，但該算法不適用于拉格朗日插值的表達(dá)式。拉格朗日插值的求值的復(fù)雜度至少也要3n次浮點(diǎn)乘(除)法和2n次浮點(diǎn)加法以上，這還是在所有的系數(shù)(將插值系數(shù)和基底的分母合并以后的系數(shù))都已經(jīng)處理完成之后，而處理系數(shù)本身可能就需要 $n^2$ 級別的復(fù)雜度。此外，拉格朗日插值表達(dá)式也不利于求微分和積分；和牛頓插值相比，當(dāng)新增數(shù)據(jù)點(diǎn)時不得不將所有的基底都改寫，很不方便。總而言之，拉格朗日插值屬于非常容易構(gòu)造的一種插值，很適用于在理論上討論某些問題，但在數(shù)值計(jì)算上仍具有很多劣勢。

4)算法實(shí)現(xiàn)

實(shí)現(xiàn)拉格朗日多項(xiàng)式插值一種途徑的MATLAB代碼如下。此處的輸出為多項(xiàng)式函數(shù)句柄。這些函數(shù)(句柄)只需要在函數(shù)名后面加括號變量即可調(diào)用，即polyfun(3)這樣的形式。

function [ polyfun ] = LagrangePI( vect, vecy )

% 生成Lagrange插值多項(xiàng)式

% 輸入兩個參數(shù)：插值點(diǎn)行向量vect，函數(shù)行向量vecy；輸出一個參數(shù)：插值多項(xiàng)式函數(shù)句柄polyfun

vecsize = size(vect, 2);

syms t f term;

f = 0;

for i = 1:vecsize

term = vecy(i);

for j = 1:vecsize

if (j ~= i)

term = term*(t - vect(j))/(vect(i) - vect(j));

end

end

f = f + term;

end

polyfun = matlabFunction(f);

end

但是，由于多項(xiàng)式形式的函數(shù)表達(dá)式帶入后為符號型變量，這意味著每一項(xiàng)的系數(shù)都經(jīng)歷了單獨(dú)計(jì)算，每一項(xiàng)的分子也需要單獨(dú)計(jì)算，這將使得拉格朗日插值表達(dá)式的函數(shù)求值(function evaluation)的復(fù)雜度達(dá)到 $O(n^2)$ 量級；如果想要使得每次求值能夠控制在 $O(n)$ 量級，就必須實(shí)現(xiàn)計(jì)算出除了含有未知量的函數(shù)基分子以外的全部系數(shù)，同時在求值時也需要一些技巧。按照如下的書寫方法可以實(shí)現(xiàn)這一目的：

function [ coefficients ] = newLagrangePI( vect, vecy )

% 生成Lagrange插值多項(xiàng)式的系數(shù)(計(jì)算分母)

% 輸入兩個參數(shù)：插值點(diǎn)行向量vect，函數(shù)行向量vecy；

% 輸出一個參數(shù)：插值多項(xiàng)式的系數(shù)行向量coefficients；

vecsize = size(vect, 2);

coefficients = zeros(1, vecsize);

for i = 1:vecsize

tmp = vecy(i); % 取得基底函數(shù)對應(yīng)的系數(shù)y_i

for j = 1:vecsize % 將其除以函數(shù)基底的分母

if (j~=i)

tmp = tmp/(vect(i) - vect(j));

end

end

coefficients(i) = tmp;

end

end

除了求系數(shù)的函數(shù)還需要一個特別的求值函數(shù)：

function [ funeval ] = evaLagrangePI( coefficients, vect, vecy, t )

% Lagrange插值多項(xiàng)式估值

% 輸入四個參數(shù)：Lagrange插值的完整系數(shù)行向量coefficients，插值點(diǎn)行向量vect，函數(shù)行向量vecy，求值點(diǎn)t；

% 輸出一個參數(shù)：函數(shù)在t處取值funeval

vecsize = size(vect, 2);

[found, index] = ismember(t, vect);

if found % 如果求值點(diǎn)是原數(shù)據(jù)點(diǎn)，則直接返回原始信息中數(shù)據(jù)點(diǎn)的函數(shù)值

funeval = vecy(index);

else % 否則，先計(jì)算全部(t-t_i)的乘積

funeval = 0; product = 1;

for i = 1:vecsize

product = product*(t - vect(i));

end

for i = 1:vecsize % 然后，計(jì)算每一項(xiàng)的值，乘以該項(xiàng)的系數(shù)并且除以該項(xiàng)分子不存在的那項(xiàng)(t-t_i)

funeval = funeval + coefficients(i)*product/(t - vect(i));

end

end

end

同樣是對于三點(diǎn) $(-2,-27),(0,-1),(1,0)$ ，調(diào)用Lagrange插值方法：

vect = [-2, 0, 1]; vecy = [-27, -1, 0];

% 命令行輸入

coefficients = newLagrangePI(vect, vecy)

% 命令行輸出

coefficients =

-4.5000 0.5000 0

% 命令行輸入

val = evaLagrangePI(coefficients, vect, vecy, -2)

% 命令行輸出

val =

-27

% 命令行輸入

val = evaLagrangePI(coefficients, vect, vecy, 0.5)

% 命令行輸出

val =

0.5000

所有的輸出均和實(shí)際的多項(xiàng)式插值 $f(t)=p_2(t)=-1+5t-4t^2$ 吻合。

3. 牛頓(Newton)插值

1)插值函數(shù)基底

單項(xiàng)式基底非常簡潔，缺點(diǎn)是求解方程組所用的是稠密的范德蒙矩陣，可能非常病態(tài)，復(fù)雜度也很高；拉格朗日基底比較精巧復(fù)雜，因?yàn)榍蠼獾南禂?shù)矩陣是單位矩陣，求解很簡單很準(zhǔn)確，缺點(diǎn)是生成表達(dá)式和函數(shù)求值復(fù)雜度很高。牛頓插值方法在二者之間提供了一個折衷選項(xiàng)：基底不如拉格朗日的函數(shù)基那么復(fù)雜，而求解又比單項(xiàng)式基底大大簡化，這來源于牛頓插值選取的基底：$$\pi_j(t)=\prod\limits_{k=1}^{j-1}(t-t_k)=(t-t_1)(t-t_2)...(t-t_{j-1}), j=1,...,n$$　　相對于拉格朗日基底的特殊性($l_j(t_i)=\delta_{ij}$)，牛頓插值基底具有一個弱一點(diǎn)的性質(zhì)：$$\pi_j(t_i)=0,\forall i

2)疊加系數(shù)

$$\pi_1(t_1)x_1+\pi_2(t_1)x_2+...+\pi_n(t_1)x_n=y_1$$

$$\pi_1(t_2)x_1+\pi_2(t_2)x_2+...+\pi_n(t_2)x_n=y_2$$

$$............$$

$$\pi_1(t_n)x_1+\pi_2(t_n)x_2+...+\pi_n(t_n)x_n=y_n$$　　寫成矩陣形式：

$$\begin{bmatrix}\pi_1(t_1)&\pi_2(t_1)&...&\pi_n(t_1)\\ \pi_1(t_2)&\pi_2(t_2)&...&\pi_n(t_2)\\...&...&...&...\\ \pi_1(t_n)&\pi_2(t_n)&...&\pi_n(t_n)\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}1&0&...&0\\1&t_2-t_1&...&0\\...&...&...&...\\1&t_n-t_1&...&(t_n-t_1)..(t_n-t_{n-1})\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$　　也就是說，牛頓插值的系數(shù)求解矩陣為一個下三角矩陣。

3)算法性質(zhì)和算法復(fù)雜度

我們知道對于下三角矩陣，利用向前代入法可以比較便利地解出，其時間復(fù)雜度為 $O(n^2)$ 。再來看生成這個下三角矩陣，復(fù)雜度也是 $O(n^2)$ 的運(yùn)算量。因此求解插值系數(shù)的總復(fù)雜度即 $O(n^2)$ 量級，比稠密矩陣的求解少一個量級。當(dāng)然，求解牛頓插值系數(shù)不一定需要顯示地生成矩陣，然后采用矩陣分解的標(biāo)準(zhǔn)套路；牛頓插值有好幾種生成系數(shù)的方法可供選擇，包括差商方法等，采用遞歸或者迭代都可以獲得良好的效果，還能避免上溢出。

此外，牛頓插值的表達(dá)式在求值時適用Horner嵌套算法(太棒了！)，這將把求值的復(fù)雜度控制在 $O(n)$ 的量級內(nèi)，如果帶上系數(shù)比拉格朗日插值表達(dá)式的求值要更高效。

牛頓插值方法有如下優(yōu)越的性質(zhì)：

3.1 當(dāng)增加數(shù)據(jù)點(diǎn)時，可以僅僅通過添加一項(xiàng)函數(shù)基和增加一個系數(shù)，生成新的牛頓插值多項(xiàng)式。這其實(shí)是可以理解的，因?yàn)楫?dāng)新增點(diǎn) $(t_{n+1},y_{n+1})$ 時，下三角系數(shù)矩陣所有的元素都沒有發(fā)生任何變化，僅僅需要新增一列和一行即可，而在該矩陣右側(cè)新增的一列全為零。這意味著所有的系數(shù) $x_1,x_2,...x_n$ 不僅滿足原線性方程組，也因此必定是新線性方程組解的部分。而基底部分也只是新增了一個基，所以新的插值多項(xiàng)式僅僅是老的插值多項(xiàng)式加一項(xiàng)而已，即 $f(t)^*=p_n(t)=p_{n-1}(t)+x_{n+1}\pi_{n+1}(t)$ 。對于新的這一項(xiàng) $x_{n+1}\pi_{n+1}(t)$ 它的系數(shù)究竟如何取，只需要將新增的數(shù)據(jù)點(diǎn)帶入即可求得：$$f(t_{n+1})^*=p_{n-1}(t_{n+1})+x_{n+1}\pi_{n+1}(t_{n+1})=y_{n+1}\quad \Rightarrow \quad x_{n+1}=\frac{y_{n+1}-p_{n-1}(t_{n+1})}{\pi_{n+1}(t_{n+1})}$$　　生成新系數(shù)的復(fù)雜度大致需要一次函數(shù)求值和一次基底求值，大致為 $O(n)$ 復(fù)雜度。如果迭代地使用這一公式，就可以生成全部牛頓插值多項(xiàng)式系數(shù)，復(fù)雜度 $O(n^2)$ ，和矩陣解法也大致是持平的。

3.2 差商法也可以實(shí)現(xiàn)生成牛頓插值多項(xiàng)式的系數(shù)。其中，差商 $f[]$ 的定義為：

$$f[t_1, t_2,...t_k]=\frac{f[t_2, t_3, ... t_{k}-f[t_1, t_2,...t_{k-1}]}{t_k-t_1}$$　　而牛頓多項(xiàng)式的系數(shù)則取自：$$x_j=f[t_1, t_2... t_j]$$　　對于這個可以有證明：

$$f[t_1]=y_1, x_1=y_1=f[t_1];\quad f[t_1, t_2]=\frac{f[t_2]-f[t_1]}{t_2-t_1}=\frac{y_2-y_1}{t_2-t_1},x_2=\frac{y_2-y_1}{t_2-t_1}=f[t_1, t_2]

$$

若$x_j=f[t_1, t_2, ...t_j]=\frac{f[t_2, t_3,...t_j]-f[t_1, t_2,...t_{j-1}]}{t_j-t_1}$ 對于任意 $j\leq k-1$ 成立，當(dāng) $j=k$ 時：

?考慮對于點(diǎn) $(t_1, y_1), (t_2,y_2),...(t_{k-1},y_{k-1})$ 的 Newton?插值多項(xiàng)式 $p_1(t)$ ；對于點(diǎn) $(t_2, y_2),(t_3, y_3),...$$(t_k, y_k)$ 的 Newton?插值多項(xiàng)式 $p_2(t)$ ，并且已知分差插值系數(shù)對任意 $j\leq k-1$ 均成立，因而有：

$$

p_1(t)=\sum\limits_{j=1}^{k-1}a_j\pi_j(t), \quad p_2(t)=\sum\limits_{j=2}^{k}b_j\pi_j(t),\qquad a_j=f[t_1,...t_{j}],b_j=f[t_2,...t_j]

$$

由 $p_1(t)$ 過點(diǎn) $(t_1, y_1)$ 到 $(t_{k-1},y_{k-1})$ ，$p_2(t)$ 過點(diǎn) $(t_2,y_2)$ 到 $(t_k,y_k)$ ，構(gòu)造插值多項(xiàng)式：

$$

p(t)=\frac{t_k-t}{t_k-t_1}p_1(t)+\frac{t-t_1}{t_k-t_1}p_2(t)

$$

就有該多項(xiàng)式通過點(diǎn) $(t_1, y_1)$ 到 $(t_k,y_k)$ ，因此即為所求的 Newton 插值多項(xiàng)式。帶入 $p_1(t),p_2(t)$ 表達(dá)式，并比較等式兩端最高次項(xiàng)系數(shù)即得：

$$

p(t)=\sum\limits_{j=1}^kx_j\pi_j(t)=\frac{t_k-t}{t_k-t_1}\sum\limits_{j=1}^{k-1}a_j\pi_j(t)+\frac{t-t_1}{t_k-t_1}\sum\limits_{j=2}^{k}b_j\pi_j‘(t)\\

x_k=\frac{-1}{t_k-t_1}a_{k-1}+\frac{1}{t_k-t_1}b_k=\frac{f[t_2,...t_k]-f[t_1,...t_{k-1}]}{t_k-t_1}=f[t_1, ...t_k]\qquad \square

$$

這個證明我摘錄自奧斯陸大學(xué)數(shù)學(xué)系的 Michael S. Floater 在 Newton Interpolation 講義里面寫的證明。

4)算法實(shí)現(xiàn)

根據(jù)3.1，可以通過新增節(jié)點(diǎn)的方法迭代地生成插值系數(shù)。利用這種思路的實(shí)現(xiàn)代碼如下：

function [ vecx_new, vect_new ] = newNPI( vect, vecx, newPoint )

%Newton插值算法新增節(jié)點(diǎn)函數(shù)；

% 輸入三個參數(shù)：原插值點(diǎn)行向量vect，原插值系數(shù)行向量vecx，新增節(jié)點(diǎn)newPoint；

% 輸入兩個參數(shù)：新系數(shù)行向量vecx_new，新插值點(diǎn)行向量vect_new；

vecsize = size(vecx, 2);

vecx_new = zeros(1, vecsize + 1); vecx_new(1:vecsize) = vecx;

vect_new = zeros(1, vecsize + 1); vect_new(1:vecsize) = vect; vect_new(vecsize + 1) = newPoint(1);

p_new = HornerNPI(vect, vecx, newPoint(1)); w_new = 1;

for i = 1:vecsize

w_new = w_new * (newPoint(1) - vect(i));

end

vecx_new(vecsize + 1) = (newPoint(2) - p_new) / w_new;

end

新增節(jié)點(diǎn)函數(shù)newNPI可以單獨(dú)使用；同時也可以反復(fù)調(diào)用生成牛頓插值系數(shù)，如下：

function [ polyfun, vecx ] = newNewtonPI( cvect, cvecy )

% 使用新增節(jié)點(diǎn)函數(shù)逐漸更新產(chǎn)生Newton插值多項(xiàng)式系數(shù)；

% 輸入兩個參數(shù)：插值點(diǎn)行向量cvect，插值系數(shù)行向量cvecx；

% 輸出兩個參數(shù)：多項(xiàng)式函數(shù)句柄polyfun，系數(shù)行向量vecx；

% 迭代生成系數(shù)行向量

vect = cvect(1); vecx = cvecy(1);

vecsize = size(cvect, 2);

for i=2:vecsize

[vecx, vect] = newNPI(vect, vecx, [cvect(i), cvecy(i)]);

end

% 采用Horner嵌套算法生成多項(xiàng)式函數(shù)句柄

syms f t; f = vecx(vecsize);

for i = vecsize-1:-1:1

f = vecx(i) + (t - cvect(i)) * f;

end

polyfun = matlabFunction(f);

end

另一種方法是采用差商。以下是實(shí)現(xiàn)的代碼。和之前的說法不同的是，本代碼使用的并非遞歸，而是正向的類似函數(shù)值緩存的算法。

function [ polyfun, vecx ] = recNewtonPI( vect, vecy )

% 使用差商產(chǎn)生Newton插值多項(xiàng)式系數(shù)；

% 輸入兩個參數(shù)：插值點(diǎn)行向量vect，函數(shù)取值cvecy；

% 輸出兩個參數(shù)：多項(xiàng)式函數(shù)polyfun，系數(shù)行向量vecx；

vecsize = size(vect, 2);

Div = diag(vecy);

% 差商生成系數(shù)行向量vecx

for k = 1:vecsize-1

for i = 1:vecsize-k

Div(i, i+k) = (Div(i+1, i+k) - Div(i, i+k-1))/(vect(i+k) - vect(i));

end

end

vecx = Div(1, :);

% 生成多項(xiàng)式函數(shù)polyfun

syms f t; f = vecx(vecsize);

for i = vecsize-1:-1:1

f = vecx(i) + (t - vect(i)) * f;

end

polyfun = matlabFunction(f);

end

但不論如何，產(chǎn)生的結(jié)果完全一致。用同樣的例子：

vect=[-2, 0, 1]; vecy=[-27, -1, 0];

% 命令行輸入1，調(diào)用新增節(jié)點(diǎn)方法

[polyfun, vecx] = newNewtonPI(vect, vecy)

% 命令行輸入2，調(diào)用差商方法

[polyfun, vecx] = recNewtonPI(vect, vecy)

% 命令行輸出1/2，完全相同

polyfun =

包含以下值的 function_handle:

@(t)-(t.*4.0-1.3e1).*(t+2.0)-2.7e1

vecx =

-27 13 -4

容易檢驗(yàn)，該多項(xiàng)式函數(shù)正是原數(shù)據(jù)點(diǎn)的多項(xiàng)式插值函數(shù)。

總結(jié)

以上是生活随笔為你收集整理的matlab全域基函数,多项式函数插值：全域多项式插值（一）单项式基插值、拉格朗日插值、牛顿插值 [MATLAB]...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。