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

時間 2019-11-10

標籤數值分析擬合插值數據方法線性二次 cubic 埃米爾拉格朗日多項式牛頓 opengl 程序欄目應用數學简体版

原文原文鏈接

插值、擬合和逼近的區別 算法

據維基百科，科學和工程問題能夠經過諸如採樣、實驗等方法得到若干離散的數據，根據這些數據，咱們每每但願獲得一個連續的函數（也就是曲線）或者更加密集的離散方程與已知數據相吻合，這過程就叫作擬合。經過擬合獲得的函數得到未知點的數據的方法，叫作插值。其中，擬合函數通過全部已知點的插值方法，叫作內插。 ide

擬合是已知點列，從總體上靠近它們；插值是已知點列而且徹底通過點列；逼近是已知曲線，或者點列，經過逼近使得構造的函數無限靠近它們。函數

最小二乘意義下的擬合，是要求擬合函數與原始數據的均方偏差達到極小，是一種總體意義的逼近，對局部性質沒有要求。而所謂"插值"，就是要在原有離散數據之間"插入"一些值，這就要求插值函數必須經過全部的離散點，插值函數在離散點以外的那些點都至關於"插入"的值。插值有全局插值，也有局部插值（好比分段線性插值），插值偏差一般考慮的是逐點偏差或最大模偏差，插值的好壞每每經過某些局部的性質來體現，好比龍格現象或吉布斯振盪。 ui

[知乎擬合與插值的區別？] 加密

插值方法 spa

多項式插值 .net

對於大部分多項式插值函數，插值點的高度值能夠視爲全部（或某些）節點高度值的線性組合，而線性組合的係數通常是x座標的多項式函數，稱做基函數。對於一個節點的基函數，它在x等於該節點的x時等於1，在x等於其餘節點的x時等於0。這就保證曲線一定通過全部節點，因此屬於內插方法。設計

在本小節，均以一組隨機數做爲已知的高度值，使它們對應於間隔固定的x座標，使用不一樣的插值函數得到各已知點（稱爲插值函數的節點）以外其它x座標所對應的高度值，畫出這些點所對應的曲線。再把全部高度值轉換成灰度值，以顏色的變化比較各插值函數。 3d

原點列如圖：（假定橫向爲x，縱向爲y。各點x座標的間隔是固定的，但y座標是隨機的） orm

線性插值

線性插值是用一系列首尾相連的線段依次鏈接相鄰各點，每條線段內的點的高度做爲插值得到的高度值。

以(xi,yi)表示某條線段的前一個端點，(x(i+1),y(i+1))表示該線段的後一個端點，則對於在[xi,x(i+1)]範圍內的橫座標爲x的點，其高度y爲：

爲便於與後面各函數比較，寫成比較對稱的形式：

其中，yi和y(i+1)的兩個參數稱爲基函數，兩者之和爲1，分別表明yi和y(i+1)對插值點高度的權值。

[wikipedia線性插值]

插值圖像以下：

將高度轉化爲灰度，獲得以下條帶：

線性插值的特色是計算簡便，但光滑性不好。若是用線性插值擬合一條光滑曲線，對每一段線段，原曲線在該段內二階導數絕對值的最大值越大，擬合的偏差越大。

二次插值

若是按照線性插值的形式，以每3個相鄰點作插值，就獲得了二次插值：

OpenGL實現代碼以下：

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();
}
二次（分段）插值圖像以下：

轉換成灰度值如圖：

二次插值在每段二次曲線內是光滑的，但在每條曲線的鏈接處其光滑性可能甚至比線性插值還差。二次插值只適合3個節點的情形，當節點數超過3個時，就須要分段插值了。
Cubic插值
若是想要保證各段曲線鏈接處光滑（一階導數相同），而且不想使用除法運算，能夠考慮Cubic插值函數：

其中，v表明插值點，v0、v1、v2、v3表明4個連續的節點。t取值爲[0,1]，將會產生一段鏈接v1和v2的曲線。也就是說，若是有n個節點，Cubic插值函數將會產生(n-2)段曲線，位於首尾兩端的節點不會歸入曲線。
實現代碼以下：
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插值圖像以下：

轉化成灰度如圖：

Cubic插值能夠產生總體上光滑的曲線，但容易產生較劇烈的波動，使得曲線的最高點比最高的節點還高、曲線的最低點比最低的節點還低。因此對顏色等取值有嚴格的上下界的數據進行插值時，會形成曲線的截取，破壞其光滑性。好比顏色（RGB三個份量取值範圍都是[0,255]），若是最高的節點是255，最低的節點是0，那麼插值後負數會被截取成0，大於255的數會被截取成255。
拉格朗日多項式插值
解決插值函數的直接構造問題以及插值偏差的估計問題，可是同時也使得插值函數值的計算變得複雜。
依照線性插值和二次插值的思路，能夠增長基函數分子和分母的階數，構造拉格朗日插值多項式：

一個n次的拉格朗日插值函數能夠繪製通過(n+1)個節點的曲線，但運算量很是大。並且在次數比較高時，容易產生劇烈的震盪（龍格現象）。因此要選擇位置特殊的節點（好比切比雪夫多項式的零點）進行插值，或使用多個次數較低的拉格朗日函數分段插值。（關於拉格朗日多項式和龍格現象，詳見維基百科連接）
分段插值實現代碼以下：
//n爲次數，nmax爲節點的總數。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次拉格朗日多項式分段插值：

轉化爲灰度：

拉格朗日多項式插值也存在鏈接處不光滑的問題。
若是直接使用20次的拉格朗日插值，獲得的圖像以下：

這樣的插值曲線顯然是不能容忍的~
牛頓插值
從算法角度考慮，提出便於計算的插值方法，也就是牛頓插值公式。拉格朗日插值每增長一個節點，整個函數要從新計算，計算量巨大。而牛頓插值每增長一個點只須要在多項式的最後增長一項，並且各基函數的係數能夠遞歸計算，減小了不少計算量。
均差（Divided differences）
是遞歸除法過程。在數值分析中，也稱差商（Difference quotient），可用於計算牛頓多項式形式的多項式插值的係數。
定義

[數值分析鍾爾傑電子科大]

階差商	階差商	階差商	階差商	階差商

差商表（高階差商是兩個低一階差商的差商）

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

[wikipedia均差]

牛頓多項式

將n次插值多項式寫做以下形式：

也就是牛頓插值公式中的各項係數就是均差表中計算出的各階均差的第一個數據！

即牛頓插值公式爲：

[wikipedia牛頓多項式]

[wikipedia均差小節：牛頓插值法]

實現代碼以下：

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();
}
}
可能因爲計算精度的緣由，牛頓插值繪製的曲線與拉格朗日插值的曲線略有不一樣。但次數較高時，牛頓插值也會產生劇烈的震盪。分段4次牛頓插值圖像以下：

轉化成灰度以下：

  牛頓插值也存在鏈接處不光滑的缺陷。
埃爾米特插值
以上各多項式插值方法的插值條件都是各節點的座標，在以低階函數分段插值時雖然能夠保持曲線的穩定（比較平緩），但在各分段曲線的鏈接處沒法保持光滑（一階導數相等）。埃爾米特插值方法不但規定了各節點的座標值，還規定了曲線在每一個節點的各階導數，這樣就能夠既保持曲線的穩定，又保證在鏈接處足夠光滑。
以3次二重（"m重"就是規定座標和曲線在全部節點處1到m-1階導數的值）埃爾米特插值爲例：

4個基函數知足分別只在y0，y1，y0的導數，y1的導數處等於1，而在其餘3個條件下等於0。能夠把埃爾米特插值看做對座標和導數的插值的組合。
曲線在每一個節點的導數能夠根據相鄰節點和它的相對位置肯定，也能夠徹底隨機生成。只要位置比較高和比較低的節點的導數絕對值不是很大，就可使整條曲線落在最高與最低節點定義的帶狀區域內，這樣就能夠避免對插值的截取。
對於本節的示例節點，一種可能的導數值以下：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};
實現代碼以下：
//p[i][0],p[i][1]爲點i的座標，p[i][2]爲曲線在點i的導數
//nmax爲節點的總數
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次埃爾米特插值圖像以下：

轉化爲灰度以下：

可見雖然n節點的埃爾米特插值是由(n-1)段曲線構成，但在每個鏈接處都是光滑的。
樣條插值
B-樣條
B-樣條是Bezier樣條的通常化，也可推廣爲NURBS樣條。先來介紹一下B-樣條：

該式定義了一個n次的B-樣條，它有(m+1)個控制點（樣條曲線的"節點"稱做控制點，由於這些點控制曲線的彎曲方向和程度，但曲線不必定通過這些點），也就有(m+1)個基函數。通常繪製的是完整的曲線，u(min)取0，u(max)取1。當u取值均勻時，該樣條稱做均勻B-樣條。當m=n時，B-樣條退化爲Bezier樣條。
函數繪製的曲線始終落在其控制點的凸包（包含一個點集全部點的凸多邊形，並且該凸多邊形全部頂點都來自這個點集）中。對於一個n次的B-樣條，改變一個控制點的位置，只改變它所在的n段曲線（由n+1個控制點定義，且以該點起始）的形狀，而不對其他的(m-n)段曲線產生影響。
Bezier樣條
Bezier（貝塞爾）樣條是法國工程師皮埃爾·貝塞爾（Pierre Bézier）在1962年爲了設計汽車主體外形曲線而提出的。其通常式表達式爲：

其中，u取值爲[0,1]，pk(k=0,...,n)爲(n+1)個節點，n稱爲階數。
Bezier樣條還能夠遞歸定義爲：

含義是n階Bezier樣條是兩條(n-1)階Bezier樣條的插值。
當階數降爲1時，Bezier插值退化成線性插值。改變任意一個控制點的位置，整條曲線的形狀都會發生變化。
比較經常使用的Bezier樣條是3次Bezier：

Beizer樣條在首尾端的切線是前兩個點和最後兩個點的連線。除了第一個點和最後一個點，其餘控制點通常都不在曲線上。
3階（4控制點）的Bezier函數圖像以下：（黑色曲線爲Bezier曲線，藍色折線爲控制點的連線）

Bezier樣條能夠實現很是快速的運算。爲了方便說明，將3次Bezier樣條表示成以下形式：

對於任意給定的u，能夠經過合併的方式將原來的7次乘法、4次加法減小爲3次乘法、3次加法：

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

只須要在繪製曲線以前計算4個常數，就能夠很快地計算出曲線上的全部點：

採用這種方式，Bezier樣條的運算量只隨階數線性增加。
實現代碼以下：
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();
}
NURBS樣條
NURBS（Non-Uniform Rational B-Splines 非均勻有理B-樣條），是貝塞爾樣條的推廣。"非均勻"的意思是控制點的間隔能夠是不均勻的，"有理"的意思是各控制點帶有不一樣的權值。和Bezier樣條相比，它對曲線形狀的控制更自由：

其基函數B(k,d)與B-樣條的基函數相同，w(k)爲各點的權因子。和B-樣條同樣，改變一個控制點的位置，只改變它所在的n段曲線的形狀，而不對其他的(m-n)段曲線產生影響。
插值的一些應用
噪聲
噪聲圖像通常須要二維或三維插值函數，不過了解了各一維插值函數，就很容易擴展到二維和三維了。關於二維噪聲圖像的產生，請參見《二維噪聲圖像》。

水波
下面這個圖像就是用插值函數繪製的模擬水波的三維網格（靜態圖）：

雖然這張網格看起來彷佛須要二維插值，但這張網格是由6個不一樣波長和頻率的Gestner合成的，每一個獨立的波形只須要一維插值生成。
[經常使用的插值函數 ]
from: http://blog.csdn.net/pipisorry/article/details/62227459
轉：http://blog.csdn.net/pipisorry/article/details/62227459

來自 <http://www.javashuo.com/article/p-alxwbwyu-ks.html>

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。