數(shù)值分析（擬合、插值和逼近）之?dāng)?shù)據(jù)插值方法（線性插值、二次插值、Cubic插值、埃米爾特、拉格朗日多項式插值、牛頓插值、樣條插值）（含opengl程序）

戰(zhàn)神之家 2019-02-13

展開全文

插值、擬合和逼近的區(qū)別

據(jù)維基百科，科學(xué)和工程問題可以通過諸如采樣、實驗等方法獲得若干離散的數(shù)據(jù)，根據(jù)這些數(shù)據(jù)，我們往往希望得到一個連續(xù)的函數(shù)（也就是曲線）或者更加密集的離散方程與已知數(shù)據(jù)相吻合，這過程就叫做擬合。通過擬合得到的函數(shù)獲得未知點的數(shù)據(jù)的方法，叫做插值。其中，擬合函數(shù)經(jīng)過所有已知點的插值方法，叫做內(nèi)插。
擬合是已知點列，從整體上靠近它們；插值是已知點列并且完全經(jīng)過點列；逼近是已知曲線，或者點列，通過逼近使得構(gòu)造的函數(shù)無限靠近它們。

最小二乘意義下的擬合，是要求擬合函數(shù)與原始數(shù)據(jù)的均方誤差達(dá)到極小，是一種整體意義的逼近，對局部性質(zhì)沒有要求。而所謂“插值”，就是要在原有離散數(shù)據(jù)之間“插入”一些值，這就要求插值函數(shù)必須通過所有的離散點，插值函數(shù)在離散點之外的那些點都相當(dāng)于“插入”的值。插值有全局插值，也有局部插值（比如分段線性插值），插值誤差通?？紤]的是逐點誤差或最大模誤差，插值的好壞往往通過某些局部的性質(zhì)來體現(xiàn)，比如龍格現(xiàn)象或吉布斯振蕩。

[知乎擬合與插值的區(qū)別？]

插值方法

多項式插值

對于大部分多項式插值函數(shù)，插值點的高度值可以視為所有（或某些）節(jié)點高度值的線性組合，而線性組合的系數(shù)一般是x坐標(biāo)的多項式函數(shù)，稱作基函數(shù)。對于一個節(jié)點的基函數(shù)，它在x等于該節(jié)點的x時等于1，在x等于其他節(jié)點的x時等于0。這就保證曲線必定經(jīng)過所有節(jié)點，所以屬于內(nèi)插方法。

在本小節(jié)，均以一組隨機數(shù)作為已知的高度值，使它們對應(yīng)于間隔固定的x坐標(biāo)，使用不同的插值函數(shù)獲得各已知點（稱為插值函數(shù)的節(jié)點）之外其它x坐標(biāo)所對應(yīng)的高度值，畫出這些點所對應(yīng)的曲線。再把所有高度值轉(zhuǎn)換成灰度值，以顏色的變化比較各插值函數(shù)。

原點列如圖：（假定橫向為x，縱向為y。各點x坐標(biāo)的間隔是固定的，但y坐標(biāo)是隨機的）

線性插值

線性插值是用一系列首尾相連的線段依次連接相鄰各點，每條線段內(nèi)的點的高度作為插值獲得的高度值。

以(xi,yi)表示某條線段的前一個端點，(x(i+1),y(i+1))表示該線段的后一個端點，則對于在[xi,x(i+1)]范圍內(nèi)的橫坐標(biāo)為x的點，其高度y為：

$p(x) = f(x_0) + \frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0). \,\!$

為便于與后面各函數(shù)比較，寫成比較對稱的形式：

其中，yi和y(i+1)的兩個參數(shù)稱為基函數(shù)，二者之和為1，分別代表yi和y(i+1)對插值點高度的權(quán)值。

[wikipedia線性插值]

插值圖像如下：

將高度轉(zhuǎn)化為灰度，得到如下條帶：

線性插值的特點是計算簡便，但光滑性很差。如果用線性插值擬合一條光滑曲線，對每一段線段，原曲線在該段內(nèi)二階導(dǎo)數(shù)絕對值的最大值越大，擬合的誤差越大。

二次插值

如果按照線性插值的形式，以每3個相鄰點做插值，就得到了二次插值：

OpenGL實現(xiàn)代碼如下：

void quadratic(float p[20][2])
{
float x,y;
int i;
float x01,x02,x12;
glColor3f(0.0,0.0,1.0);
glBegin(GL_LINE_STRIP);
for(i=0;i<20;i+=2)
{
x01=p[i][0]-p[i+1][0];
x02=p[i][0]-p[i+2][0];
x12=p[i+1][0]-p[i+2][0];
for(x=p[i][0];x<=p[i+2][0];x+=1.0)
{
y=(x-p[i+1][0])*(x-p[i+2][0])/x01/x02*p[i][1]-(x-p[i][0])*(x-p[i+2][0])/x01/x12*p[i+1][1]+(x-p[i][0])*(x-p[i+1][0])/x02/x12*p[i+2][1];
glVertex2f(x,y);
}
}
glEnd();
}

二次（分段）插值圖像如下：

轉(zhuǎn)換成灰度值如圖：

二次插值在每段二次曲線內(nèi)是光滑的，但在每條曲線的連接處其光滑性可能甚至比線性插值還差。二次插值只適合3個節(jié)點的情形，當(dāng)節(jié)點數(shù)超過3個時，就需要分段插值了。

Cubic插值

如果想要保證各段曲線連接處光滑（一階導(dǎo)數(shù)相同），并且不想使用除法運算，可以考慮Cubic插值函數(shù)：

其中，v代表插值點，v0、v1、v2、v3代表4個連續(xù)的節(jié)點。t取值為[0,1]，將會產(chǎn)生一段連接v1和v2的曲線。也就是說，如果有n個節(jié)點，Cubic插值函數(shù)將會產(chǎn)生(n-2)段曲線，位于首尾兩端的節(jié)點不會納入曲線。

實現(xiàn)代碼如下：

float cubic(float v0,float v1,float v2,float v3,float x)
{
float v32=v3-v2;
float v01=v0-v1;
float v20=v2-v0;
float temp=(v32-v01)*x;
temp=(temp+v01+v01-v32)*x;
temp=(temp+v20)*x+v1;
return temp;
}
void drawCubic(float p[20])
{
float x,y;
int step;
float delta;
glColor3f(0.0,0.0,1.0);
glBegin(GL_LINE_STRIP);
for(step=0;step<17;step++)
{
for(delta=0.0;delta<=1.0;delta+=0.05)
{
x=delta*(p[step+1][0]-p[step][0])+p[step][0];
y=cubic(p[step],p[step+1],p[step+2],p[step+3],delta);
glVertex2f(x,y);
}
}
glEnd();
}

Cuibc插值圖像如下：

轉(zhuǎn)化成灰度如圖：

Cubic插值可以產(chǎn)生整體上光滑的曲線，但容易產(chǎn)生較劇烈的波動，使得曲線的最高點比最高的節(jié)點還高、曲線的最低點比最低的節(jié)點還低。所以對顏色等取值有嚴(yán)格的上下界的數(shù)據(jù)進(jìn)行插值時，會造成曲線的截取，破壞其光滑性。比如顏色（RGB三個分量取值范圍都是[0,255]），如果最高的節(jié)點是255，最低的節(jié)點是0，那么插值后負(fù)數(shù)會被截取成0，大于255的數(shù)會被截取成255。

拉格朗日多項式插值

解決插值函數(shù)的直接構(gòu)造問題以及插值誤差的估計問題，但是同時也使得插值函數(shù)值的計算變得復(fù)雜。

依照線性插值和二次插值的思路，可以增加基函數(shù)分子和分母的階數(shù)，構(gòu)造拉格朗日插值多項式：

一個n次的拉格朗日插值函數(shù)可以繪制經(jīng)過(n+1)個節(jié)點的曲線，但運算量非常大。而且在次數(shù)比較高時，容易產(chǎn)生劇烈的震蕩（龍格現(xiàn)象）。所以要選擇位置特殊的節(jié)點（比如切比雪夫多項式的零點）進(jìn)行插值，或使用多個次數(shù)較低的拉格朗日函數(shù)分段插值。（關(guān)于拉格朗日多項式和龍格現(xiàn)象，詳見維基百科鏈接）

分段插值實現(xiàn)代碼如下：

//n為次數(shù)，nmax為節(jié)點的總數(shù)。n不大于nmax
void lagrange(float p[][2],int n,int nmax)
{
float x,y;
int i,j,t;
float temp;
glColor3f(0.0,0.0,1.0);
for(i=0;i<=(nmax-1);i+=(n-1))
{
glBegin(GL_LINE_STRIP);
for(x=p[i][0];x<=p[i+n-1][0];x+=1.0)
{
y=0.0;
for(j=0;j<n;j++)
{
temp=1.0;
for(t=0;t<n;t++)
{
if(t==j)
continue;
temp*=(x-p[i+t][0])/(p[i+j][0]-p[i+t][0]);
}
y+=temp*p[i+j][1];
}
glVertex2f(x,y);
}
glEnd();
}
}

使用4次拉格朗日多項式分段插值：

轉(zhuǎn)化為灰度：

拉格朗日多項式插值也存在連接處不光滑的問題。

如果直接使用20次的拉格朗日插值，得到的圖像如下：

這樣的插值曲線顯然是不能容忍的~

牛頓插值

從算法角度考慮，提出便于計算的插值方法，也就是牛頓插值公式。拉格朗日插值每增加一個節(jié)點，整個函數(shù)要重新計算，計算量巨大。而牛頓插值每增加一個點只需要在多項式的最后增加一項，而且各基函數(shù)的系數(shù)可以遞歸計算，減少了很多計算量。

均差（Divided differences）

是遞歸除法過程。在數(shù)值分析中，也稱差商（Difference quotient），可用于計算牛頓多項式形式的多項式插值的系數(shù)。

定義

[數(shù)值分析鐘爾杰電子科大]

差商表（高階差商是兩個低一階差商的差商）
	$0$ 階差商	階差商	階差商	階差商	$\ldots$	階差商
$x_{0}$	$f[x_{{0}}]$
$x_{{1}}$	$f[x_{{1}}]$	$f[x_{{0}},x_{{1}}]$
$x_{{2}}$	$f[x_{{2}}]$	$f[x_{{1}},x_{{2}}]$	$f[x_{{0}},x_{{1}},x_{{2}}]$
$x_{{3}}$	$f[x_{{3}}]$	$f[x_{{2}},x_{{3}}]$	$f[x_{{1}},x_{{2}},x_{{3}}]$	$f[x_{{0}},x_{{1}},x_{{2}},x_{{3}}]$
$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
$x_{{k}}$	$f[x_{{k}}]$	$f[x_{{k-1}},x_{{k}}]$	$f[x_{{k-2}},x_{{k-1}},x_{{k}}]$	$f[x_{{k-3}},x_{{k-2}},x_{{k-1}},x_{{k}}]$	$\ldots$	$f[x_{{0}},\ldots ,x_{{k}}]$

由函數(shù)表出發(fā)，從上往下從左往右依次計算出一階，二階，。。。，各階均差。

[wikipedia均差]

牛頓多項式

將n次插值多項式寫作如下形式：

$N(x)=[y_{0}]+[y_{0},y_{1}](x-x_{0})+\cdots +[y_{0},\ldots ,y_{k}](x-x_{0})(x-x_{1})\cdots (x-x_{{k-1}})$

也就是牛頓插值公式中的各項系數(shù)就是均差表中計算出的各階均差的第一個數(shù)據(jù)！

即牛頓插值公式為：

[wikipedia牛頓多項式]

[wikipedia均差小節(jié)：牛頓插值法]

實現(xiàn)代碼如下：

void newton(float p[][2],int n,int nmax)
{
float x,y;
int i,j,t;
float temp;
float f[20][20];
glColor3f(0.0,0.0,1.0);
for(i=0;i<=(nmax-1);i+=(n-1))
{
for(t=0;t<n;t++)
{
f[t][0]=(p[i+t][1]-p[i+t+1][1])/(p[i+t][0]-p[i+t+1][0]);
for(j=1;j<=t;j++)
f[t][j]=(f[t-1][j-1]-f[t][j-1])/(p[i+t-j][0]-p[i+t+1][0]);
}
glBegin(GL_LINE_STRIP);
for(x=p[i][0];x<=p[i+n-1][0];x+=1.0)
{
y=p[i][1];
temp=1.0;
for(j=0;j<n;j++)
{
temp*=x-p[i+j][0];
y+=temp*f[j][j];
}
glVertex2f(x,y);
}
glEnd();
}
}

可能由于計算精度的原因，牛頓插值繪制的曲線與拉格朗日插值的曲線略有不同。但次數(shù)較高時，牛頓插值也會產(chǎn)生劇烈的震蕩。分段4次牛頓插值圖像如下：

轉(zhuǎn)化成灰度如下：

牛頓插值也存在連接處不光滑的缺陷。

埃爾米特插值

以上各多項式插值方法的插值條件都是各節(jié)點的坐標(biāo)，在以低階函數(shù)分段插值時雖然可以保持曲線的穩(wěn)定（比較平緩），但在各分段曲線的連接處無法保持光滑（一階導(dǎo)數(shù)相等）。埃爾米特插值方法不但規(guī)定了各節(jié)點的坐標(biāo)值，還規(guī)定了曲線在每個節(jié)點的各階導(dǎo)數(shù)，這樣就可以既保持曲線的穩(wěn)定，又保證在連接處足夠光滑。

以3次二重（"m重"就是規(guī)定坐標(biāo)和曲線在所有節(jié)點處1到m-1階導(dǎo)數(shù)的值）埃爾米特插值為例：

4個基函數(shù)滿足分別只在y0，y1，y0的導(dǎo)數(shù)，y1的導(dǎo)數(shù)處等于1，而在其他3個條件下等于0?？梢园寻柮滋夭逯悼醋鲗ψ鴺?biāo)和導(dǎo)數(shù)的插值的組合。

曲線在每個節(jié)點的導(dǎo)數(shù)可以根據(jù)相鄰節(jié)點和它的相對位置確定，也可以完全隨機生成。只要位置比較高和比較低的節(jié)點的導(dǎo)數(shù)絕對值不是很大，就可以使整條曲線落在最高與最低節(jié)點定義的帶狀區(qū)域內(nèi)，這樣就可以避免對插值的截取。

對于本節(jié)的示例節(jié)點，一種可能的導(dǎo)數(shù)值如下：grad[20]={-4.0,4.0,0.5,-2.0,1.0,-2.0,-2.0,2.0,1.0,1.0,0.0,-1.0,-4.0,3.0,0.0,-3.0,3.0,0.0,-4.0,-4.0};

實現(xiàn)代碼如下：

//p[i][0],p[i][1]為點i的坐標(biāo)，p[i][2]為曲線在點i的導(dǎo)數(shù)
//nmax為節(jié)點的總數(shù)
void hermite(float p[][3],int nmax)
{
float x,y;
float a1,a0,b1,b0;
float x00,x11;
int i;
float temp;
glColor3f(0.0,0.0,1.0);
for(i=0;i<nmax;i++)
{
glBegin(GL_LINE_STRIP);
x00=p[i][0];
x11=p[i+1][0];
for(x=p[i][0];x<=p[i+1][0];x+=1.0)
{
temp=(x11-x00)*(x11-x00);
a0=(x11-3*x00+2*x)*(x11-x)*(x11-x)/temp/(x11-x00);
a1=(-x00+3*x11-2*x)*(x-x00)*(x-x00)/temp/(x11-x00);
b0=(x-x00)*(x-x11)*(x-x11)/temp;
b1=(x-x00)*(x-x00)*(x-x11)/temp;
y=a0*p[i][1]+a1*p[i+1][1]+b0*p[i][2]+b1*p[i+1][2];
glVertex2f(x,y);
}
glEnd();
}
}

分段3次埃爾米特插值圖像如下：

轉(zhuǎn)化為灰度如下：

可見雖然n節(jié)點的埃爾米特插值是由(n-1)段曲線構(gòu)成，但在每一個連接處都是光滑的。

樣條插值

B-樣條

B-樣條是Bezier樣條的一般化，也可推廣為NURBS樣條。先來介紹一下B-樣條：

該式定義了一個n次的B-樣條，它有(m+1)個控制點（樣條曲線的“節(jié)點”稱作控制點，因為這些點控制曲線的彎曲方向和程度，但曲線不一定經(jīng)過這些點），也就有(m+1)個基函數(shù)。一般繪制的是完整的曲線，u(min)取0，u(max)取1。當(dāng)u取值均勻時，該樣條稱作均勻B-樣條。當(dāng)m=n時，B-樣條退化為Bezier樣條。

函數(shù)繪制的曲線始終落在其控制點的凸包（包含一個點集所有點的凸多邊形，而且該凸多邊形所有頂點都來自這個點集）中。對于一個n次的B-樣條，改變一個控制點的位置，只改變它所在的n段曲線（由n+1個控制點定義，且以該點起始）的形狀，而不對其余的(m-n)段曲線產(chǎn)生影響。

Bezier樣條

Bezier（貝塞爾）樣條是法國工程師皮埃爾·貝塞爾（Pierre Bézier）在1962年為了設(shè)計汽車主體外形曲線而提出的。其一般式表達(dá)式為：

其中，u取值為[0,1]，pk(k=0,...,n)為(n+1)個節(jié)點，n稱為階數(shù)。

Bezier樣條還可以遞歸定義為：

含義是n階Bezier樣條是兩條(n-1)階Bezier樣條的插值。

當(dāng)階數(shù)降為1時，Bezier插值退化成線性插值。改變?nèi)我庖粋€控制點的位置，整條曲線的形狀都會發(fā)生變化。

比較常用的Bezier樣條是3次Bezier：

Beizer樣條在首尾端的切線是前兩個點和最后兩個點的連線。除了第一個點和最后一個點，其他控制點一般都不在曲線上。

3階（4控制點）的Bezier函數(shù)圖像如下：（黑色曲線為Bezier曲線，藍(lán)色折線為控制點的連線）

Bezier樣條可以實現(xiàn)非?？焖俚倪\算。為了方便說明，將3次Bezier樣條表示成如下形式：

對于任意給定的u，可以通過合并的方式將原來的7次乘法、4次加法減少為3次乘法、3次加法：

一般情況下，應(yīng)用Bezier樣條繪制曲線時，都是先給定一個很小的步長t（步長足夠小才能保證Bezier曲線的精確），讓t從0取到1，從頭到尾繪制整條曲線。在t不變的條件下，可以使用更快速的差分方法，利用前一個點計算出下一個點的值，將每步的計算量減小到只有3次加法：

只需要在繪制曲線之前計算4個常數(shù)，就可以很快地計算出曲線上的所有點：

采用這種方式，Bezier樣條的運算量只隨階數(shù)線性增長。

實現(xiàn)代碼如下：

void bezier3(float xx0,float yy0,float xx1,float yy1,float xx2,float yy2,float xx3,float yy3)
{
GLint i;
GLfloat x,y;
GLdouble ax,bx,cx,ay,by,cy;
GLdouble t;
GLdouble dx,d2x,dy,d2y,delta;
GLdouble d3x,d3y;
ax=-xx0+3.0*xx1-3.0*xx2+xx3;
bx=3.0*xx0-6.0*xx1+3.0*xx2;
cx=-3.0*xx0+3.0*xx1;
ay=-yy0+3.0*yy1-3.0*yy2+yy3;
by=3.0*yy0-6.0*yy1+3.0*yy2;
cy=-3.0*yy0+3.0*yy1;
delta=0.0005;
d3x=6.0*ax*delta*delta*delta;
d3y=6.0*ay*delta*delta*delta;
x=xx0;
y=yy0;
glBegin(GL_LINE_STRIP);
glVertex2f(x,y+200);
d2x=6.0*ax*delta*delta*delta+2.0*bx*delta*delta;
dx=ax*delta*delta*delta+bx*delta*delta+cx*delta;
d2y=6.0*ay*delta*delta*delta+2.0*by*delta*delta;
dy=ay*delta*delta*delta+by*delta*delta+cy*delta;
for(t=0.0;t<=1.0;t+=delta)
{
d2x+=d3x;
dx+=d2x;
x+=(float)dx;
d2y+=d3y;
dy+=d2y;
y+=(float)dy;
glVertex2f(x,y);
}
glEnd();
}