全域多項式插值指的是在整個插值區域內形成一個多項式函數作為插值函數。關于多項式插值的基本知識，見“計算基本理論”。

在單項式基插值和牛頓插值形成的表達式中，求該表達式在某一點處的值使用的Horner嵌套算法啊，見"Horner嵌套算法"。

1. 單項式(Monomial)基插值

1)插值函數基　　單項式基插值采用的函數基是最簡單的單項式：$$\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}$$　　所要求解的系數即為單項式系數 $x_1,x_2,...x_n$ ，在這里仍然采用1,2,...n的下標記號而不采用和單項式指數對應的0,1,2,...,n-1的下標僅僅是出于和前后討論一致的需要。

2)疊加系數

單項式基插值采用單項式函數基，若有m個離散數據點需要插值，設使用n項單項式基底：

$$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$$　　系數矩陣為一 $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}$$　　根據計算基本理論中的討論，多項式插值的函數基一定線性無關，且只要離散數據點兩兩不同，所構成的矩陣行也一定線性無關，這保證了矩陣一定行滿秩。此時當且僅當m=n時，疊加系數有且僅有一組解。因此，n=插值基底的個數=離散數據點的個數=多項式次數+1。

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

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

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

4)算法實現

使用MATLAB實現單項式插值代碼如下：

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

% 計算單項式型插值多項式系數

% 輸入四個參數：插值點行向量vect，插值點函數值行向量vecy，平移shift，壓限scale；

% 輸出兩個參數：插值多項式各項系數行向量vecx，矩陣條件數condition；

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

if nargin==2

shift = 0; scale = 1;

end

% 求解系數

vecsize = size(vect, 2);

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

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

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

% 生成句柄函數polyfunc

syms t;

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

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

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

end

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

end

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

% 命令行輸入

[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

和預期完全一致。

2. 拉格朗日(Lagrange)插值

1)插值函數基

拉格朗日插值采用的是一種設計巧妙的多項式基，每個基底都是n-1次多項式，而每個基底函數當且僅當在第i個數據點處取1，在其余數據點均為零。這個多項式基是這樣設計的：

$$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$ 為克羅內克(Kronecker)記號，當兩個下標相等時為1，否則為零；也可以將 $\delta_{ij}$ 理解為一個二階張量，即單位矩陣。只要將各個$t_i$ 帶入定義式，上式是很容易驗證的。這意味著拉格朗日插值的疊加系數的求解將會產生很好的性質，即：

2)疊加系數

需要求解的插值函數即：$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)問題條件和算法復雜度

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

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

4)算法實現

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

function [ polyfun ] = LagrangePI( vect, vecy )

% 生成Lagrange插值多項式

% 輸入兩個參數：插值點行向量vect，函數行向量vecy；輸出一個參數：插值多項式函數句柄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

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

function [ coefficients ] = newLagrangePI( vect, vecy )

% 生成Lagrange插值多項式的系數(計算分母)

% 輸入兩個參數：插值點行向量vect，函數行向量vecy；

% 輸出一個參數：插值多項式的系數行向量coefficients；

vecsize = size(vect, 2);

coefficients = zeros(1, vecsize);

for i = 1:vecsize

tmp = vecy(i); % 取得基底函數對應的系數y_i

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

if (j~=i)

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

end

end

coefficients(i) = tmp;

end

end

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

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

% Lagrange插值多項式估值

% 輸入四個參數：Lagrange插值的完整系數行向量coefficients，插值點行向量vect，函數行向量vecy，求值點t；

% 輸出一個參數：函數在t處取值funeval

vecsize = size(vect, 2);

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

if found % 如果求值點是原數據點，則直接返回原始信息中數據點的函數值

funeval = vecy(index);

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

funeval = 0; product = 1;

for i = 1:vecsize

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

end

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

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

end

end

end

同樣是對于三點 $(-2,-27),(0,-1),(1,0)$ ，調用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

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

3. 牛頓(Newton)插值

1)插值函數基底

單項式基底非常簡潔，缺點是求解方程組所用的是稠密的范德蒙矩陣，可能非常病態，復雜度也很高；拉格朗日基底比較精巧復雜，因為求解的系數矩陣是單位矩陣，求解很簡單很準確，缺點是生成表達式和函數求值復雜度很高。牛頓插值方法在二者之間提供了一個折衷選項：基底不如拉格朗日的函數基那么復雜，而求解又比單項式基底大大簡化，這來源于牛頓插值選取的基底：$$\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}$)，牛頓插值基底具有一個弱一點的性質：$$\pi_j(t_i)=0,\forall i

2)疊加系數

$$\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}$$　　也就是說，牛頓插值的系數求解矩陣為一個下三角矩陣。

3)算法性質和算法復雜度

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

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

牛頓插值方法有如下優越的性質：

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

3.2 差商法也可以實現生成牛頓插值多項式的系數。其中，差商 $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}$$　　而牛頓多項式的系數則取自：$$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$ 成立，當 $j=k$ 時：

?考慮對于點 $(t_1, y_1), (t_2,y_2),...(t_{k-1},y_{k-1})$ 的 Newton?插值多項式 $p_1(t)$ ；對于點 $(t_2, y_2),(t_3, y_3),...$$(t_k, y_k)$ 的 Newton?插值多項式 $p_2(t)$ ，并且已知分差插值系數對任意 $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)$ 過點 $(t_1, y_1)$ 到 $(t_{k-1},y_{k-1})$ ，$p_2(t)$ 過點 $(t_2,y_2)$ 到 $(t_k,y_k)$ ，構造插值多項式：

$$

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

$$

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

$$

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

$$

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

4)算法實現

根據3.1，可以通過新增節點的方法迭代地生成插值系數。利用這種思路的實現代碼如下：

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

%Newton插值算法新增節點函數；

% 輸入三個參數：原插值點行向量vect，原插值系數行向量vecx，新增節點newPoint；

% 輸入兩個參數：新系數行向量vecx_new，新插值點行向量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

新增節點函數newNPI可以單獨使用；同時也可以反復調用生成牛頓插值系數，如下：

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

% 使用新增節點函數逐漸更新產生Newton插值多項式系數；

% 輸入兩個參數：插值點行向量cvect，插值系數行向量cvecx；

% 輸出兩個參數：多項式函數句柄polyfun，系數行向量vecx；

% 迭代生成系數行向量

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嵌套算法生成多項式函數句柄

syms f t; f = vecx(vecsize);

for i = vecsize-1:-1:1

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

end

polyfun = matlabFunction(f);

end

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

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

% 使用差商產生Newton插值多項式系數；

% 輸入兩個參數：插值點行向量vect，函數取值cvecy；

% 輸出兩個參數：多項式函數polyfun，系數行向量vecx；

vecsize = size(vect, 2);

Div = diag(vecy);

% 差商生成系數行向量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, :);

% 生成多項式函數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

但不論如何，產生的結果完全一致。用同樣的例子：

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

% 命令行輸入1，調用新增節點方法

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

% 命令行輸入2，調用差商方法

[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

容易檢驗，該多項式函數正是原數據點的多項式插值函數。

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。