插值和擬合

本文轉載自查看原文 2020-04-05 22:13 1760 數學

2020.4.8

matlab插值

預先選用什么插值方法，可以關注曲線的實際形狀。比如道路標線的話，可以選用分段線性插值，因為道路標線一般是直的，衛星軌道的話可以選擇三次多項式插值。因為衛星軌道是橢圓。

clc;clear;close all;
x=0:2*pi;
y=sin(x);
xx=0:0.5:2*pi;

%interp1對sin函數進行分段線性插值，調用interp1的時候，默認的是分段線性插值
y1=interp1(x,y,xx);
subplot(2,2,1);plot(x,y,'o',xx,y1,'r')
title('分段線性插值')

%臨近插值
y2=interp1(x,y,xx,'nearest');
subplot(2,2,2);
plot(x,y,'o',xx,y2,'r');
title('臨近插值')

%球面線性插值
y3=interp1(x,y,xx,'spline');
subplot(2,2,3);
plot(x,y,'o',xx,y3,'r')
title('球面插值')

%三次多項式插值法
y4=interp1(x,y,xx,'cubic');
subplot(2,2,4);
plot(x,y,'o',xx,y4,'r');
title('三次多項式插值')

(1) Nearest方法速度最快，占用內存最小，但一般來說誤差最大，插值結果最不光滑。

(2) Spline三次樣條插值是所有插值方法中運行耗時最長的，插值函數及其一二階導函數都連續，是最光滑的插值方法。占用內存比cubic方法小，但是已知數據分布不均勻的時候可能出現異常結果。

(3) Cubic三次多項式插值法中，插值函數及其一階導數都是連續的，所以插值結果比較光滑，速度比Spline快，但是占用內存最多。
————————————————
版權聲明：本文為CSDN博主「風翼冰舟」的原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/zb1165048017/article/details/48343507

2020.4.6

一、插值和擬合的區別

什么叫插值？插值是數學領域數值分析中的通過已知的離散數據求未知數據的過程或方法。

首先插值和擬合都是根據某個未知函數（或已知但難於求解的函數）的幾個已知數據點求出變化規律和特征相似的近似曲線的過程。但是插值法要求的是近似的曲線需要完全經過數據點，而擬合則是得到最接近的結果，強調最小方差的概念。插值和擬合的區別如下圖所示^[1]（其中左邊為插值，右邊為擬合）：

二、常見插值法

1.基本概念

設函數 $y=f(x)$ 有n個已知數據點 $[x_i,y_i],i=0,1,…,n$ ，若存在一簡單函數 $P(x)$ ，使 $P(x_i)=y_i,i=1,2,…,n$ 成立。則 $P(x)$ 為 $f(x)$ 的插值函數，點 $x_i$ 為插值節點，求 $P(x)$ 的方法稱為插值法。

2.拉格朗日插值（可以參考：https://www.zhihu.com/question/58333118/answer/262507694）

在數值分析中，拉格朗日插值法是以法國十八世紀數學家約瑟夫·路易斯·拉格朗日命名的一種多項式插值方法。許多實際問題中都用函數來表示某種內在聯系或規律，而不少函數都只能通過實驗和觀測來了解。如對實踐中的某個物理量進行觀測，在若干個不同的地方得到相應的觀測值，拉格朗日插值法可以找到一個多項式，其恰好在各個觀測的點取到觀測到的值。這樣的多項式稱為拉格朗日（插值）多項式。數學上來說，拉格朗日插值法可以給出一個恰好穿過二維平面上若干個已知點的多項式函數。拉格朗日插值法最早被英國數學家愛德華·華林於1779年發現^[1]，不久后（1783年）由萊昂哈德·歐拉再次發現。1795年，拉格朗日在其著作《師范學校數學基礎教程》中發表了這個插值方法，從此他的名字就和這個方法聯系在一起^[2]。

對於給定的若n+1個點 $(x_0, y_0),(x_1, y_1),\ldots,(x_n, y_n)$ ，對應於它們的次數不超過n的拉格朗日多項式 $\scriptstyle L$ 只有一個。如果計入次數更高的多項式，則有無窮個，因為所有與 $\scriptstyle L$ 相差 $\lambda (x-x_0)(x-x_1)\ldots(x-x_n)$ 的多項式都滿足條件。

Lagrange插值多項式的公式如下：

基函數：

$\begin{aligned} l_{i}(x) &=\frac{\left(x-x_{0}\right) \cdots\left(x-x_{i-1}\right)\left(x-x_{i+1}\right) \cdots\left(x-x_{n}\right)}{\left(x_{i}-x_{0}\right) \cdots\left(x_{i}-x_{i-1}\right)\left(x_{i}-x_{i+1}\right) \cdots\left(x_{i}-x_{n}\right)} \\ &=\prod_{j=0 \atop j \neq i}^{n} \frac{x-x_{j}}{x_{i}-x_{j}}, \quad(i=0,1, \cdots, n) \end{aligned}$

插值多項式：

$L_{n}(x)=\sum_{i=0}^{n} y_{i} l_{i}(x)=\sum_{i=0}^{n} y_{i}\left(\prod_{j=0 \atop j \neq i}^{n} \frac{x-x_{j}}{x_{i}-x_{j}}\right)$

插值余項：

$R_{n}(x)=f(x)-L_{n}(x)=\frac{f^{(n+1)}(\xi)}{(n+1) !} \omega_{n+1}(x), \xi \in(a, b)$

當 $\max\limits_{a\leq{x}\leq{b}}|f^{n+1}(x)|=M_{n+1}$ 時，截斷誤差界是：

$R_{n}(x)=\frac{M^{(n+1)}(x)}{(n+1) !} \omega_{n+1}(x)$

其中： $\omega_{n+1}(x)=\prod_{i=0}^{n}\left(x-x_{j}\right)$

范例

假設有某個多項式函數 $f ，已知它在三個點上的取值為：$

$f (4) = 10$
$f (5) = 5.25$
$f (6) = 1$

要求 $f (18)的值。$

首先寫出每個拉格朗日基本多項式：

matlab代碼：

function y=Lagrange(x0,y0,x) %輸入：x0:節點變量數據 % y0:節點函數值 % x:插值數據 %輸出：y:插值函數值 n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end end

例1：設 $f(x)=ln(x)$ ,並給出如下節點數據

x0=[0.4,0.5,0.7,0.8]

y0=[-0.916291,-0.693147,-0.356675,-0.223144]

估計x=0.6的值

帶入Lagrange()函數，得到 $y=-0.5100$ ，實際值： $ln(0.6)= -0.5108$

優點與缺點

拉格朗日插值法的公式結構整齊緊湊，在理論分析中十分方便，然而在計算中，當插值點增加或減少一個時，所對應的基本多項式就需要全部重新計算，於是整個公式都會變化，非常繁瑣^[5]。這時可以用重心拉格朗日插值法或牛頓插值法來代替。此外，當插值點比較多的時候，拉格朗日插值多項式的次數可能會很高，因此具有數值不穩定的特點，也就是說盡管在已知的幾個點取到給定的數值，但在附近卻會和“實際上”的值之間有很大的偏差（如右下圖）^[6]。這類現象也被稱為龍格現象，解決的辦法是分段用較低次數的插值多項式。

2.分段線性插值：

2.1基本原理：

將每兩個相鄰的節點用直線連起來，如此形成的一條折線就是分段線性插值函數。計算點的插值時，只用到左右的兩個節點，計算量與節點個數n（初始值x0，y0的長度,n=length(x0)）無關，而拉格朗日插值與n值有關。分段線性插值中n越大，分段越多，插值誤差越小。

2.2 Matlab實現分段線性插值：

用matlab實現分段線性插值不需要自己手動編寫函數，matlab有現成的一維插值函數interp1

y=interp1(x0,y0,x,'method')

該式可以根據X，Y的值來計算函數在X1處的值。其中X，Y是兩個等長的已知向量，分別表示采樣點和采樣值。X1是一個向量或標量，表示要插值的點。

method指定插值方法，其值可為：

linear：線性插值（默認）線形插值，默認方法。將與插值點靠近的兩個數據點用直線連接，然后在直線上選取對應插值點的數據。

nearest：最近項插值，最近點插值。選擇最近樣本點的值作為插值數據。

spline：逐次3次樣條插值，每一個分段內構造一個三次多項式，使其插值函數除滿足插值條件外，還要求在各節點處具有連續的一階和二階導數。

cubic：保凹凸性 3 次插值

pchip: 分段3次埃爾米特插值。采用分段三次多項式，除滿足插值條件，還需滿足在若干節點處相鄰段插值函數的一階導數相等，使得曲線光滑的同時，還具有保形性。

所有插值方法都要求x0單調。

MATLAB的二維插值函數為interp2(),

調用格式：y1 = interp2(X,Y,Z,X1,Y1,method)

其中X，Y兩個向量，表示兩個參數的采樣點，Z是采樣點對應的函數值。X1，Y1是兩個標量或向量，表示要插值的點。指定的算法method計算二維插值。

’linear’：雙線性插值算法(缺省算法);
’nearest’：最臨近插值;
’spline’：三次樣條插值;
’cubic’：雙三次插值。

3.三次樣條插值：

使用三次樣條插值有兩種方法：其中一種就是第二種插值方式（分段線性插值），只需要將method修改為spline即可實現。

還有一種實現方式如下：

pp=csape(x0,y0);

y=ppval(pp,x);

其中x0，y0，x與前面含義相同，返回值y即插值結果。

4.彩蛋（例題）：

表1給出的 x, y 數據位於機翼斷面的下輪廓線上，假設需要得到 x 坐標每改變0.1 時的 y 坐標。試完成加工所需數據，畫出曲線。要求用 Lagrange、分段線性和三次樣條三種插值方法計算。

表1

x 0 3 5 7 9 11 12 13 14 15
y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6
解：編寫代碼如下：

clear,clc;close all;
x0=[0,3,5,7,9,11,12,13,14,15];
y0=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
x1=[0:0.1:15];

%拉格朗日插值
y1=lagrange(x0,y0,x1);
subplot(1,3,1);hold on;plot(x0,y0,'ro');
subplot(1,3,1);hold on;plot(x1,y1,'b.');%axis equal;axis off;
title('拉格朗日插值')
%分段線性插值
y2=interp1(x0,y0,x1);
subplot(1,3,2);hold on;plot(x0,y0,'ro');
subplot(1,3,2);hold on;plot(x1,y2,'b.');%axis equal;axis off;
title('分段線性插值')
%三次樣條插值
y3=interp1(x0,y0,x1,'spline');
subplot(1,3,3);hold on;plot(x0,y0,'ro');
subplot(1,3,3);hold on;plot(x1,y3,'b.');%axis equal;axis off;
title('三次樣條插值');
function y=lagrange(x0,y0,x)
%輸入：x0:節點變量數據
%      y0:節點函數值
%      x:插值數據
%輸出：y:插值函數值
    n=length(x0);m=length(x);
    for i=1:m
        z=x(i);
        s=0.0;
        for k=1:n
            p=1.0;
            for j=1:n
                if j~=k
                    p=p*(z-x0(j))/(x0(k)-x0(j));
                end
            end
            s=p*y0(k)+s;
        end
        y(i)=s;
    end
end

結果圖為:

紅色圓圈是原始數據點，藍色點點是插值點（包含插值位置和插值值）。

由圖可知：拉格朗日插值誤差較大，插值大量點的時候不建議使用。分段線性插值結果最好，使用的是method是linear（默認）。三次樣條插值有點誤差，但不是很大，得到的曲線比較光滑。

————————————————
版權聲明：本文為CSDN博主「Da_wan」的原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/Da_wan/article/details/82223572

3.Newton插值

上面介紹的拉格朗日插值多項式，當插值節點增減時，計算要全部進行，很不方便，所以提出一種Newton插值。

（1）差商和差分的性質

一階差商（均差）：

$f\left[x_{i}, x_{j}\right]=\frac{f\left(x_{i}\right)-f\left(x_{j}\right)}{x_{i}-x_{j}}$

二階差商（均差）：一階差商的差商

$f\left[x_{i}, x_{j},x_{k}\right]=\frac{f\left[x_{i}, x_{j}\right]-f\left[x_{j}, x_{k}\right]}{x_{i}-x_{k}}$

$k$ 階差商（均差）：

$f\left[x_{0}, x_{1}, \cdots, x_{k}\right]=\frac{f\left[x_{0}, x_{1}, \cdots, x_{k-1}\right]-f\left[x_{1}, x_{2}, \cdots, x_{k}\right]}{x_{0}-x_{k}}$

（2）Newton插值多項式

$\begin{aligned} N_{n}(x) &=f\left(x_{0}\right)+\left(x-x_{0}\right) f\left[x_{0}, x_{1}\right]+\cdots \\ &+\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-1}\right) f\left[x_{0}, x_{1}, \cdots, x_{n}\right] \end{aligned}$

（3）誤差

$\begin{aligned} R_{n}(x) &=\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n}\right) f\left[x_{0}, x_{1}, \cdots, x_{n}\right] \end{aligned}\\=w_{n+1}(x)f[x,x_{0},x_{1},…,x_{n}]\\=\frac{M^{(n+1)}(x)}{(n+1) !} \omega_{n+1}(x)$

(4)差商與導數的關系

$f\left[x_{0}, x_{1}, \cdots, x_{k}\right]=\frac{f^{n}(\xi)}{n!}$

matlab代碼如下：

function [YY,y]=newton_chazhi(X,Y,x,M) %輸入為：X-插值點的x軸向量 %Y-插值點的y軸向量 %需要求解的x變量 %M為多項式次數 %輸出YY為差分表 %y是x對應的因變量 m=length(X); YY=zeros(m); YY(:,1)=Y; %求查分表 for i=2:m for j=i:m YY(j,i)=(YY(j,i-1)-YY(j-1,i-1))/(X(j)-X(j-i+1)); end end y=Y(1); %計算newton插值公式 for i=1:M xl=1; for j=1:i xl=xl*(x-X(j)); end y=y+xl*YY(i+1,i+1); end end function [YY,y]=main() X=[0.40,0.55,0.65,0.80,0.90,1.05]; Y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382]; x=0.596; M=4; [YY,y]=newton_chazhi(X,Y,x,M); end

4.三次Hermite插值

$P(x)=f(x_0)+f[x_0,x_1](x_1-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+A(x-x_0)(x-x_1)(x-x_2)$

其中： $A=\frac{f'(x_1)-f[x_0,x_1]-(x_1-x_0)f[x_0,x_1,x_2]}{(x_1-x_0)(x_1-x_2)}$ ，余項表達式：

$R(x)=\frac{1}{4!}f^{(4)}(\xi)(x-x_0)(x-x_1)^2(x-x_2)$

用Hermite進行插值時需要有確定的三個數據點以及中間點的一階導數。

例：給定 $f(x)=x^{3/2},x_0=\frac{1}{4},x_1=1,x_2=\frac{9}{4},$ 求三次樣條插值多項式 $P(x)$ ，及余項表達式。

matlab代碼：

function y=hermiter_chazhi(X,Y,x1,x) %求差分表 m=length(X); YY=zeros(m); YY(:,1)=Y; for i=2:m for j=i:m YY(j,i)=(YY(j,i-1)-YY(j-1,i-1))/(X(j)-X(j-i+1)); end end %求A A=x1-YY(2,2)-(X(2)-X(1))*YY(3,3); %求插值 y=Y(1)+YY(2,2).*(x-X(1))+YY(3,3).*(x-X(1)).*(x-X(2))+A.*(x-X(1)).*(x-X(2)).*(x-X(3)); end function y=main_her() x=[1/4:0.01:9/4]; f=x.^(3/2); X=[1/4,1,9/4]; Y=[1/8,1,27/8]; x1=3/2; y=hermiter_chazhi(X,Y,x1,x); plot(x,y,"r") grid on hold on plot(x,f,"b") scatter(X,Y) legend("插值曲線","實際曲線") end

5.三次樣條插值

三次樣條插值參考下面的博客：

HappyWang：軌跡生成--三次樣條插值zhuanlan.zhihu.com

三、常見的擬合方法

HappyWang：回歸預測（1）--線性回歸和多項式擬合zhuanlan.zhihu.com

參考

^https://www.zhihu.com/question/24276013

作者：HappyWang
鏈接：https://zhuanlan.zhihu.com/p/98431641
來源：知乎
著作權歸作者所有。商業轉載請聯系作者獲得授權，非商業轉載請注明出處。

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 Matlab隨筆之插值與擬合（上） MATLAB實現插值和擬合 matlab——插值與擬合 MATLAB擬合和插值 matlab學習——05插值和擬合(一維二維插值,擬合) Matlab隨筆之插值與擬合（下）插值回歸擬合逼近的區別 Matlab曲面擬合和插值 Python SciPy庫——插值與擬合 Matlab數學建模學習筆記——插值與擬合