Matlab數(shù)據(jù)插值與擬合

數(shù)據(jù)插值與擬合在工程實踐與科學(xué)實驗中,常常需要從一組試驗數(shù)據(jù)之中找到自變量與因變量之間的關(guān)系，一般可用一個近似函數(shù)表示。函數(shù)產(chǎn)生的方法因觀測數(shù)據(jù)的要求不同而異，數(shù)據(jù)插值與擬合是兩種常用的方法。4.1MATLAB中的插值函數(shù)4.2拉格朗日插值法4.3利用均差的牛頓插值法4.4利用差分的牛頓插值法4.5Hermite插值4.6spline三次樣條插值4.7多項式曲線擬合4.8最小二乘擬合4.1MATLAB中的插值函數(shù)函數(shù)插值來源于函數(shù)的以下問題：只知道函數(shù)在某區(qū)間有定義且已得到區(qū)間內(nèi)一些離散點的值，希望用簡單的表達式近似給出函數(shù)在此區(qū)間上的整體描述，并能與離散點上的值相等。插值法按插值函數(shù)的形式主要分為以下幾種形式：〔1〕代數(shù)多項式插值；〔2〕三角多項式插值；〔3〕有理分式插值。代數(shù)多項式插值是最常用的插值方式，其內(nèi)容也是最豐富的，它又可分為以下幾種插值方式：〔1〕非等距節(jié)點插值，包括拉格朗日插值、利用均差的牛頓插值和埃特金插值；〔2〕非等距節(jié)點插值，包括利用差分的牛頓插值和高斯插值等；〔3〕在插值中增加了導(dǎo)數(shù)的Hermite〔埃爾米特〕插值；〔4〕分段插值，包括分段線性插值、分段Hermite〔埃爾米特〕插值和樣條函數(shù)插值；〔5〕反插值。按被插值函數(shù)的變量個數(shù)還可把插值法分為一元插值和多元插值。4.1.1一元插值函數(shù)MATLAB中的一元插值函數(shù)為interp1()，它的功能是一維數(shù)據(jù)插值〔表格查找〕。該命令對數(shù)據(jù)點之間進行計算內(nèi)插值，它出一元函數(shù)f(x)在中間點的數(shù)值，其中函數(shù)f(x)由所給數(shù)據(jù)決定。一元插值函數(shù)interp1()的幾種調(diào)用格式如表4-1所示。表4-1一維插值插值函數(shù)interp1的語法格式語法形式說明y=interp1(x,Y,xi)由已知點集(x,Y）插值計算xi上的函數(shù)值y=interp1(x,Y,xi)相當(dāng)于x=1:length(Y)的interp(x,Y,xi)y=interp1(x,Y,xi,method)用指定插值方法計算插值點xi上的函數(shù)值y=interp1(x,Y,xi,method,’extrap’)對xi中超出已知點集的插值點用指定插值方法計算函數(shù)值y=interp1(x,Y,xi,method,’extrap’,extrapval)用指定方法插值xi上的函數(shù)值，超出已知點集處函數(shù)值取extrapvaly=interp1(x,Y,xi,method,’pp’)用指定方法插值，但返回結(jié)果為分段多項式method方法描述‘nearest’最鄰近插值：插值點處函數(shù)值取與插值點最鄰近的已知點的函數(shù)值‘liner’分段線性插值：插值點處函數(shù)值由連接其最鄰近的兩側(cè)點的線性函數(shù)預(yù)測，MATLAB中interp1的默認方法‘spline’樣條插值：默認為三次樣條插值?？捎胹pline函數(shù)代替‘pchip’三次Hermite多項式插值。可用pchip函數(shù)代替‘cubic’同‘pchip’，三次Hermite多項式插值MATLAB中一維插值有多種算法，由interp1函數(shù)中的method指定。MATLAB中一維插值的各種算法如表4-2所示。表4-2一維插值算法(method)1.Linear〔分段線性插值〕它的算法是在每個小區(qū)間[xi,xi+1]上采用簡單的線性插值。在區(qū)間[xi,xi+1]上的子插值多項式為：由此整個區(qū)間[xi,xi+1]上的插值函數(shù)為：

其中定義如下：分段線性插值方法在速度和誤差之間取得了比較好的均衡，其插值函數(shù)具有連續(xù)性，但在數(shù)據(jù)點處的斜率一般不會改變，因此不是光滑的。分段線性插值方法是MATLAB一維插值默認的方法。2.Spline〔樣條插值〕樣條插值是用分段低次多項式去逼近函數(shù)。樣條函數(shù)可以給出光滑的插值曲線，只要在插值區(qū)間端點提供某些導(dǎo)數(shù)信息，樣條插值可以適應(yīng)不同光滑需求。三次樣條是使用最為廣泛的樣條插值，它在每個子區(qū)間[xi,xi+1]上都是有二階連續(xù)導(dǎo)數(shù)的三次多項式，即

其中都是三次多項式。對于給定的離散的測量數(shù)據(jù)經(jīng)x，y〔稱為斷點〕，要尋找一個三次多項式y(tǒng)=p(x)，以逼近每對數(shù)據(jù)(xi,yi)點間曲線。過兩點(xi,yi)和(xi+1,yi+1)只能確定一條直線，而通過一點的三次多項式曲線有無窮多條。為使通過中間斷點的三次多項式曲線具有唯一性，要增加以下的連續(xù)條件和邊界條件〔因為三次多項式有4個系數(shù)〕：〔1〕三次多項式在點(xi,yi)處有：；〔2〕三次多項式在點(xi,yi)處有：；〔3〕三次多項式在點(xi,yi)處有：；〔4〕邊界條件：。表4-2中各種方法中：〔1〕nearest方法速度最快，占用內(nèi)存最小，但一般來說誤差最大，插值結(jié)果最不光滑；〔2〕spline三次樣條插值是所有插值方法中運行耗時最長的，其插值函數(shù)以及插值函數(shù)的一階、二階導(dǎo)函數(shù)都連續(xù)，因此是最光滑的插值方法，占用內(nèi)存上比cubic方法小，但當(dāng)數(shù)據(jù)點不均勻分布時可能出現(xiàn)異常結(jié)果?！?〕cubic三次多項式插值法中插值函數(shù)及其一階導(dǎo)數(shù)都是連續(xù)的，因此其插值結(jié)果也比較光滑，運算速度比spline方法略快，但占用內(nèi)存最多。在實際的使用中，應(yīng)根據(jù)實際需求和運算條件選擇適宜的算法。例4-1用interp1對sin函數(shù)進行分段線性插值。解：在MATLAB命令窗口中輸入以下命令：>>x=0:2*pi;>>y=sin(x);>>xx=0:0.5:2*pi>>yy=interp1(x,y,xx);>>plot(x,y,'s',xx,yy)注：例4-1中用默認的(分段線性插值的linear)對的7個sin函數(shù)的數(shù)據(jù)點進行插值，用plot畫出插值結(jié)果。從圖中可以看出分段線性就是聯(lián)結(jié)兩個鄰近的點的線性函數(shù)插值計算該區(qū)間內(nèi)插值點上的函數(shù)值。例4-2用其他一維插值方法對以下7個離散數(shù)據(jù)點

(1,3.5)、(2,2.1)、(3,1.3)、(4.0.8)、(5,2.9)、(6,4.2)、(7,5.7)進行一維插值方法。解：在MATLAB命令窗口中輸入以下命令：

>>x=[1234567];>>y=[3.52.11.30.82.94.25.7];>>xx=1:0.5:7;>>y1=interp1(x,y,xx,'nearest');>>y2=interp1(x,y,xx,'spline');>>y3=interp1(x,y,xx,'cubic');>>plot(x,y,'o',xx,y1,'-',xx,y2,'-.',xx,y3,':')4.2拉格朗日插值法拉格朗日插值法是基于基函數(shù)的插值方法，插值多項式可表示為其中稱為i次基函數(shù)：在MATLAB中編程實現(xiàn)拉格朗日插值法函數(shù)為：Language。功能：求數(shù)據(jù)點的拉格朗日多項式；調(diào)用格式：f=Language(x,y)或f=Language(x,y,x0)。其中，x為數(shù)據(jù)點的x坐標向量；y為數(shù)據(jù)點的y坐標向量；x0為插值點的x坐標；f為求得的拉格朗日多項式或x0處的插值。functionf=Language(x,y,x0)%求數(shù)據(jù)點的拉格朗日多項式%數(shù)據(jù)點的x坐標向量:x%數(shù)據(jù)點的y坐標向量:y%為插值點的x坐標:x0%求得的拉格朗日多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等！');return;end%檢錯f=0.0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));%計算拉格朗日基函數(shù)end;f=f+l;%計算拉格朗日插值函數(shù)simplify(f);%化簡if(i==n)if(nargin==3)f=subs(f,'t',x0);%計算插值點的函數(shù)值elsef=collect(f);%將插值多項式展開f=vpa(f,6);%將插值多項式的系數(shù)化成6位精度的小數(shù)endendend例4-3根據(jù)下表的數(shù)據(jù)點求出其拉格朗日插值多項式，并計算當(dāng)x=1.6時y的值。解：>>

x=[11.21.82.54];>>y=[0.84150.93200.97380.5985-0.7568];>>f=language(x,y)f=1.05427*t-.145485e-1*t^2-.204917*t^3+.328112e-1*t^4-.261189e-1>>f=language(x,y,1.6)f=0.9992x11.21.82.54y0.84150.93200.97380.5985-0.75684.3利用均差的牛頓插值法函數(shù)f的零階均差定義為，一階定義均差為一般地，函數(shù)f的k階均差定義為：利用均差的牛頓插值法多項式為

系數(shù)的計算過程如表4-3所示表4-3均差計算表格一階均差二階均差三階均差……n階均差……………在MATLAB中編程實現(xiàn)利用均差牛頓插值法函數(shù)為：Newton。功能：求數(shù)據(jù)點的均差形式的牛頓插值多項式；調(diào)用格式：f=Newton(x,y)或f=Newton(x,y,x0)。其中，x為數(shù)據(jù)點的x坐標向量；y為數(shù)據(jù)點的y坐標向量；x0為插值點的x坐標；f為求得的牛頓插值法多項式或x0處的插值。functionf=Newton(x,y,x0)%求數(shù)據(jù)點的均差形式牛頓插值多項式%數(shù)據(jù)點的x坐標向量:x%數(shù)據(jù)點的y坐標向量:y%為插值點的x坐標:x0%求得的均差形式牛頓插值多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;endf=y(1);y1=0;l=1;

for(i=1:n-1)for(j=i+1:n)y1(j)=(y(j)-y(i))/(x(j)-x(i));endc(i)=y1(i+1);l=l*(t-x(i));f=f+c(i)*l;simplify(f);y=y1;

if(i==n-1)if(nargin==3)f=subs(f,'t',x0);elsef=collect(f);%將插值多項式展開

f=vpa(f,6);endendend例4-4根據(jù)下表的數(shù)據(jù)點求出其均差形式牛頓插值多項式，并計算當(dāng)x=2.0時y的值。解：>>x=[11.21.82.54];>>y=[11.443.246.2516];>>f=Newton(x,y)f=.182711e-14-.482154e-14*t+1.00000*t^2-.169177e-14*t^3+.211471e-15*t^4>>f=Newton(x,y,2.0)f=4x11.21.82.54y11.443.246.25164.4利用差分的牛頓插值法差分分為向前差分、后向差分和中心差分三種，它們的記法及定義如下為：

其中：代表向前差分；代表向后差分；代表向后差分。假設(shè)。為了方便，可構(gòu)造如表4-4所示的差分表〔〕。表4-4差分計算表格………………向前牛頓插值向前牛頓插值多項式可表示如下：其中叫做步長，，且的取值范圍為。在MATLAB中編程實現(xiàn)向前牛頓插值法函數(shù)為：Newtonforward。功能：求數(shù)據(jù)點的向前牛頓插值法多項式；調(diào)用格式：f=Newtonforward(x,y)或f=Newtonforward(x,y,x0)。其中，x為數(shù)據(jù)點的x坐標向量；y為數(shù)據(jù)點的y坐標向量；x0為插值點的x坐標；f為求得的向前牛頓插值法多項式或x0處的插值。functionf=Newtonforward(x,y,x0)%求數(shù)據(jù)點的向前差分牛頓插值多項式%數(shù)據(jù)點的x坐標向量:x%數(shù)據(jù)點的y坐標向量:y%為插值點的x坐標:x0%求得的向前差分牛頓插值多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;endf=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');return;endfor(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);endc(i)=y1(1);l=t;for(k=1:i-1)l=l*(t-k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(1))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend向前牛頓插值向后牛頓插值多項式可表示如下：其中叫做步長，，且的取值范圍為。在MATLAB中編程實現(xiàn)向后牛頓插值法函數(shù)為：Newtonback。功能：求數(shù)據(jù)點的向前牛頓插值法多項式；調(diào)用格式：f=Newtonback(x,y)或f=Newtonback(x,y,x0)。其中，x為數(shù)據(jù)點的x坐標向量；y為數(shù)據(jù)點的y坐標向量；x0為插值點的x坐標；f為求得的向前牛頓插值法多項式或x0處的插值。functionf=Newtonback(x,y,x0)%求數(shù)據(jù)點的向后差分牛頓插值多項式%數(shù)據(jù)點的x坐標向量:x%數(shù)據(jù)點的y坐標向量:y%為插值點的x坐標:x0%求得的向前差分牛頓插值多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;endf=y(n);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');return;endfor(i=1:n-1)for(j=i+1:n)y1(j)=y(j)-y(j-1);endc(i)=y1(n);l=t;for(k=1:i-1)l=l*(t+k);end;

f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(n))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend例4-5根據(jù)下表的數(shù)據(jù)點求出其差分形式的牛頓插值多項式，并計算當(dāng)x=1.55時y的值。解：>>x=[11.21.41.61.8];>>y=[0.84150.9320 0.98540.99960.9738];>>f=Newtonforward(x,y)f=.841500+.108025*t-.169042e-1*t^2-.675000e-3*t^3+.541667e-4*t^4>>f=Newtonforward(x,y,1.55)f=0.9998f=Newtonback(x,y)f=.973800-.457417e-1*t-.198042e-1*t^2+.191667e-3*t^3+.541667e-4*t^4>>f=Newtonback(x,y,1.55)f=0.9998x11.21.41.61.8y0.84150.93200.98540.99960.97384.5Hermite插值Hermite插值滿足在節(jié)點上等于給定函數(shù)值，而且在節(jié)點上的導(dǎo)數(shù)值也等于給定的導(dǎo)數(shù)值。對于高階導(dǎo)數(shù)的情況，Hermite插值多項式比較復(fù)雜，在實際中，常常遇到的是函數(shù)值與一階導(dǎo)數(shù)給定的情況。在此情況下，n個節(jié)點x1，x2，…，xn的Hermite插值多項式的表達式如下：其中，，，在MATLAB中編程實現(xiàn)Hermite插值法函數(shù)為：Hermite。功能：求數(shù)據(jù)點的Hermite插值法多項式；調(diào)用格式：f=Hermite(x,y,y_1)或f=Hermite(x,y,y_1,x0)。其中，x為數(shù)據(jù)點的x坐標向量；y為數(shù)據(jù)點的y坐標向量；y_1為數(shù)據(jù)點導(dǎo)數(shù)向量；x0為插值點的x坐標；f為求得的Hermite插值法多項式或x0處的插值。functionf=Hermite(x,y,y_1,x0)%求數(shù)據(jù)點的向后差分牛頓插值多項式%數(shù)據(jù)點的x坐標向量:x%數(shù)據(jù)點的y坐標向量:y%數(shù)據(jù)點的導(dǎo)數(shù)向量:y_1%求得的Hermite插值多項式或x0處的插值:fsymst;f=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的導(dǎo)數(shù)的維數(shù)不相等！');return;endelsedisp('x和y的維數(shù)不相等！');return;endfori=1:nh=1.0;a=0.0;forj=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endend

f=f+h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));

if(i==n)if(nargin==4)f=subs(f,'t',x0);elsef=vpa(f,6);endendend例4-6根據(jù)下表的數(shù)據(jù)點求出其Hermite插值多項式，并計算當(dāng)x=1.144時y的值。解：

>>x=1:0.2:1.8;>>y=[11.09541.18321.26491.3416];>>y_1=[0.50.45640.42260.39530.3727];>>f=Hermite(x,y,y_1)f=678.168*(t-1.20000)^2*(t-1.40000)^2*(t-1.60000)^2*(t-1.80000)^2*(-20.3333+21.3333*t)+10850.7*(t-1.)^2*(t-1.40000)^2*(t-1.60000)^2*(t-1.80000)^2*(-10.4063+9.58473*t)+24414.1*(t-1.)^2*(t-1.20000)^2*(t-1.60000)^2*(t-1.80000)^2*(.591560+.422600*t)+10850.7*(t-1.)^2*(t-1.20000)^2*(t-1.40000)^2*(t-1.80000)^2*(17.4978-10.1455*t)+678.168*(t-1.)^2*(t-1.20000)^2*(t-1.40000)^2*(t-1.60000)^2*(50.9807-27.5773*t)>>f=Hermite(x,y,y_1,1.44)f=1.2000x11.21.41.61.8y11.09541.18321.26491.3416y’0.50000.45640.42260.39530.37274.6spline三次樣條插值Spline插值為分段三次插值，即：在每一個小區(qū)間上是不超過三次的多項式且具有二階連續(xù)導(dǎo)數(shù)，利用具有一階導(dǎo)數(shù)邊界條件的源程序為：naspline.mfunctionm=naspline(x,y,dy0,dyn,xx)%用途：三階樣條插值〔一階導(dǎo)數(shù)邊界條件〕%格式：x為節(jié)點向量；y為數(shù)據(jù)；dyo，dyn為左右兩端點的一階導(dǎo)數(shù)%如果xx缺省，那么輸出各節(jié)點的一階導(dǎo)數(shù)值，否那么m為xx的三階樣條插值n=length(x)-1;%計算小區(qū)間的個數(shù)h=diff(x);lemda=h(2:n)./(h(1:n-1)+h(2:n));mu=1-lemda;g=3*(lemda.*diff(y(1:n))./h(1:n-1)+mu.*diff(y(2:n+1))./h(2:n));g(1)=g(1)-lemda(1)*dy0;g(n-1)=g(n-1)-mu(n-1)*dyn;%求解三對角方程dy=nachase(lemda,2*ones(1:n-1),mu,g);%假設(shè)給插值節(jié)點，計算插值m=[dy0;dy;dyn];ifnargin>=5s=zeros(size(xx));fori=1:nifi==1kk=find(xx<=x(2));elseifi==nkk=find(xx>=x(n));elsekk=find(xx>x(i)&xx<=x(i+1));endxbar=(xx(kk)-x(i))/h(i);s(kk)=alpha0(xbar)*y(i)+alpha1(xbar)*y(i+1)+...h(i)*beta0(xbar)*m(i)+h(i)*beta1(xbar)*m(i+1);endm=s;end%追趕法functionx=nachase(a,b,c,d)n=length(a);fork=2:nb(k)=b(k)-a(k)/b(k-1)*c(k-1);d(k)=d(k)-a(k)/b(k-1)*d(k-1);endx(n)=d(n)/b(n);fork=n-1:-1:1x(k)=(d(k)-c(k)*x(k+1)/b(k));endx=x(:);%基函數(shù)functiony=alpha0(x)y=2*x.^3-3*x.^2+1;functiony=alpha1(x)y=-2*x.^3+3*x.^2;functiony=beta0(x)y=x.^3-2*x.^2+x;functiony=beta1(x)y=x.^3-x.^2;4.7多項式曲線擬合對給定的試驗數(shù)據(jù)點(xi,yi)(i=1,2,…,N)，可以構(gòu)造m次多項式：由曲線擬合的定義，應(yīng)該使得下式取極小值通過簡單運算可以得出系數(shù)是下面線性方程組的解