多項式函數插值：全域多項式插值（一）單項式基插值、拉格朗日插值、牛頓插值 [MATLAB]

時間  2019-11-07
標籤 多項式 函數 插值 單項式 拉格朗日 牛頓 matlab 欄目 MATLAB 简体版
原文   原文鏈接

　　全域多項式插值指的是在整個插值區域內造成一個多項式函數做爲插值函數。關於多項式插值的基本知識，見「計算基本理論」。html

　　在單項式基插值和牛頓插值造成的表達式中，求該表達式在某一點處的值使用的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）疊加係數lua

　　單項式基插值採用單項式函數基，如有m個離散數據點須要插值，設使用n項單項式基底：spa

$$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）問題條件和算法複雜度htm

　　生成範德蒙矩陣的複雜度大體在 $O(n^2)$ 量級；因爲範德蒙矩陣並無什麼好的性質，既沒有稀疏性，也沒有對稱性，所以只能使用標準的稠密矩陣分解（如LU）來解決，複雜度在 $O(n^3)$ 量級。所以，問題的複雜度在 $O(n^3)$ 量級。 blog

　　範德蒙矩陣存在的問題是，當矩陣的維數愈來愈高的時候，解線性方程組時問題將愈來愈病態，條件數愈來愈大；從另外一個角度來講，單項式基底自己趨勢相近，次數增大時將愈來愈趨於平行（見下圖）。這將形成隨着數據點的增長，肯定的疊加係數的不肯定度愈來愈大，所以雖然單項式基很簡單，進行插值時卻每每不用這一方法。若是仍然採用單項式基底，有時也會進行兩種能夠改善以上問題的操做：平移(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<j$$　　求出的多項式形如：$f(t)=p_{n-1}(t)=\sum\limits_{j=1}^nx_j\pi_j(t)=x_1+x_2(t-t_1)+...+x_n(t-t_1)(t-t_2)...(t-t_{n-1})$

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

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

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程