MATLAB的差值與曲線擬合

1、MATLAB的差值與曲線擬合,MATLAB數(shù)據(jù)處理與應用2011-2012學年選修課,第八講,王文健 ,.SOUTHWEST JIAOTONG UNIVERSITY,2,主要內(nèi)容,插值運算 曲線擬合 分段函數(shù) 拉格拉日插值 牛頓插值 Hermite插值,.SOUTHWEST JIAOTONG UNIVERSITY,3,MATLAB插值與曲線擬合,簡介 把試驗數(shù)據(jù)繪制成圖是點狀分布，將這些點練成曲線，然后轉化為有意義的數(shù)學函數(shù)，才能對數(shù)據(jù)做對比和分析 上述過程涉及插值(interpolation)和曲線擬合(curve-fitting) 插值過程中認為數(shù)據(jù)是準確的，求取描述點之間的數(shù)據(jù) 曲線擬合

2、中，假定已知曲線的規(guī)律，做曲線的最佳逼近，不需要經(jīng)過所有數(shù)據(jù)點,.SOUTHWEST JIAOTONG UNIVERSITY,4,MATLAB插值與曲線擬合,插值計算 假設試驗得到一組數(shù)據(jù)形態(tài)為： f(xk)；其中k=1,2,n；且x1=axn=c 如果某些點xi不屬于上述的xi ，但是a=xi=c，要估計這些點的函數(shù)值f(xi)就需要做插值運算 根據(jù)原始數(shù)據(jù)所描述的函數(shù)的復雜程度，存在有： 一維插值 二維插值 Spline插值,.SOUTHWEST JIAOTONG UNIVERSITY,5,MATLAB插值與曲線擬合,一維插值 一維插值可以是線性的，也可以是三次多項式或spline插值 一

3、維線性插值是假設兩個數(shù)據(jù)中的變化為線性關系，因此可由已知點的坐標(f(a),a)和(f(c),c)計算b點的函數(shù)值f(b) 一維線性插值是最簡單的插值，適用范圍很小 可用較為復雜的三次多項式或spline來近似找到原函數(shù)f(x),.SOUTHWEST JIAOTONG UNIVERSITY,6,MATLAB插值與曲線擬合,一維插值 MATLAB的一維插值函數(shù)為interp1(a,c,b)和interp1(a,c,b,method) 其中a,c為已知數(shù)據(jù)，而b為要插值的數(shù)據(jù)點，method為預先設定的插值方法，分別為線性(linear)、三次多項式(cubic)和spline 如果數(shù)據(jù)的變化較大

4、，以spline插值所形成的曲線最為平滑，效果最好 三次多項式所得的曲線平滑度介于linear和spline之間,.SOUTHWEST JIAOTONG UNIVERSITY,7,MATLAB插值與曲線擬合,一維插值舉例 假設汽車引擎在額定轉速下溫度與時間的測量值關系如下00； 120； 260； 368； 477； 5110；用一維插值法估計時間為2.6s和4.9s時的溫度值。 分別用三種方法估計其值 x=0 1 2 3 4 5; y=0 20 60 68 77 110; y1=interp1(x,y,2.6) y2=interp1(x,y,2.6 4.9) y3=interp1(x,y,2

5、.6 4.9,cubic) y3=interp1(x,y,2.6 4.9,spline),.SOUTHWEST JIAOTONG UNIVERSITY,8,MATLAB插值與曲線擬合,一維插值舉例 畫出插值后的時間與溫度變化曲線圖。 t=0 1 2 2.6 3 4 4.9 5; y1=0 20 60 64.8 68 77 106.7 110; y2=0 20 60 65.2 68 77 105.6 110; y3=0 20 60 67.3 68 77 105.2 110; plot(t,y1,t,y2,t,y3),.SOUTHWEST JIAOTONG UNIVERSITY,9,MATLAB插

6、值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,10,MATLAB插值與曲線擬合,二維插值 假設試驗所得到的一組數(shù)據(jù)形態(tài)為z=f(xk,yk)，期中k=1，2，n；如果某些點(xi,yi)不屬于上述點，要估計這些點的函數(shù)值f(xi,yi)就需要進行二維擦繪制運算，二維插值相當于二元函數(shù)運算 MATLAB中二維插值函數(shù)是interp2(x,y,z,xi,yi)，期中xi,yi為要插值的數(shù)據(jù)點，x，y，z為已知數(shù)據(jù) 通過設定interp2(x,y,z,xi,yi,method)也可以設定不同的插值方法，有l(wèi)inear、cubic和spline三種,.SOUTHWEST

7、JIAOTONG UNIVERSITY,11,MATLAB插值與曲線擬合,二維插值舉例 假設汽車的引擎轉速、溫度與時間的測量值關系如下表，估計時間為2.6s和速度為2500r時的溫度。,.SOUTHWEST JIAOTONG UNIVERSITY,12,MATLAB插值與曲線擬合,二維插值舉例 應用函數(shù)interp2函數(shù)分別用linear、cubic、spline來估計時間2.6s和轉速2500r的溫度值。 t=1 2 3 4 5; speed=2000 3000 4000; temp=20 110 176;60 180 220;68 240 349;77 310 450;110 405 50

8、3; temp1=interp2(speed,t,temp,2500,2.6) temp2=interp2(speed,t,temp,2500,2.6,cubic) temp3=interp2(speed,t,temp,2500,2.6,spline),.SOUTHWEST JIAOTONG UNIVERSITY,13,MATLAB插值與曲線擬合,spline插值 函數(shù)spline有兩種應用形式，spline(x,y,xi)和spline(x,y) spline(x,y,xi)與interp1(x,y,xi,spline)效果一樣 舉例 x=0 1 2 3 4 5; y=0 20 60 68

9、77 110; y1=spline(x,y,2.6) y1=spline(x,y,2.6,4.9),.SOUTHWEST JIAOTONG UNIVERSITY,14,MATLAB插值與曲線擬合,spline插值 舉例：在0,10區(qū)間上按正弦規(guī)律y=sin(x)取10個點yi，再在區(qū)間0,10上取41個點，以yi為已知數(shù)據(jù)，對這41個點做spline插值運算，得到函數(shù)值并繪圖，比較其值與正弦波偏離程度。 x=0:10; y=sin(x); xx=0:0.25:10; yy=spline(x,y,xx); plot(x,y,o,xx,yy,xx,sin(xx);,.SOUTHWEST JIAOT

10、ONG UNIVERSITY,15,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,16,MATLAB插值與曲線擬合,spline插值 采用一維線性插值運算 x=0:10; y=sin(x); xx=0:0.25:10; yy=interp1(x,y,xx); plot(x,y,o,xx,yy,xx,sin(xx);,.SOUTHWEST JIAOTONG UNIVERSITY,17,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,18,MATLAB插值與曲線擬合,spline插值 采用一維三次多項式插值運算

11、x=0:10; y=sin(x); xx=0:0.25:10; yy=interp1(x,y,xx,cubic); plot(x,y,o,xx,yy,xx,sin(xx);,.SOUTHWEST JIAOTONG UNIVERSITY,19,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,20,MATLAB插值與曲線擬合,spline插值 向量t=1900:10:1990，p=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633表示從1900年到199

12、0年人口普查的數(shù)字，利用spline插值法預測2000年的人口數(shù)，并繪制人口變化趨勢圖。 t=1900:10:1990; p=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633; y=spline(t,p,2000) plot(t 2000,p y,o),.SOUTHWEST JIAOTONG UNIVERSITY,21,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,22,MATLAB插值與曲線擬合,spline插值 畫一個圓并在數(shù)據(jù)點y(:,2)

13、，y(:,6)處標上字母o。 x=pi*0:0.5:2; y=0 1 0 -1 0 1 0;1 0 1 0 -1 0 1; pp=spline(x,y); yy=ppval(pp,linspace(0,2*pi,101);%計算分算函數(shù)的值 plot(yy(1,:),yy(2,:),-b,y(1,2:5),y(2,2:5),or),.SOUTHWEST JIAOTONG UNIVERSITY,23,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,24,MATLAB插值與曲線擬合,曲線擬合 曲線擬合(curve-fitting)顧名思義就是用近似的曲線方

14、程來代表一組離散的數(shù)據(jù) 曲線擬合與插值有許多相似之處，二者最大區(qū)別再于曲線擬合要找出一個曲線函數(shù)式，而插值僅求出與數(shù)據(jù)點對應的函數(shù)值而已 曲線擬合運用最小二乘法原理 如果擬合的曲線限定為多項式就稱為多項式最小二乘法曲線擬合 MATLAB提供的曲線擬合有多項式最小二乘擬合、普通最小二乘擬合等,.SOUTHWEST JIAOTONG UNIVERSITY,25,MATLAB插值與曲線擬合,線性回歸 假設有一組數(shù)據(jù)形態(tài)為y=y(x)，其中x=0 1 2 3 4 5，y=0 20 60 68 77 110，用一個最簡單的方程式來表示這組數(shù)據(jù)，一階線性方程式最簡單 隨意假設其方程式為y=20 x x=0

15、 1 2 3 4 5; y=0 20 60 68 77 110; y1=20*x sum_sq=sum(y-y1).2) plot(x,y1,x,y,o),title(Linear estimate),grid,.SOUTHWEST JIAOTONG UNIVERSITY,26,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,27,MATLAB插值與曲線擬合,線性回歸 從圖中看出，所假設的方程與實際數(shù)據(jù)的差距很大，從理論上講，誤差平方和最小才是代表數(shù)據(jù)的最佳方程式，以此標準尋找數(shù)據(jù)方程式的方法稱為最小平方誤差法或線性回歸 MATLAB中的polyfit

16、()函數(shù)提供了從一階到高階多項式的回歸運算，語法為coef=polyfit(x,y,n)，期中x、y為輸入數(shù)據(jù)，n為多項式的階數(shù)，coef為輸出數(shù)組，數(shù)值為多項式的系數(shù)項,.SOUTHWEST JIAOTONG UNIVERSITY,28,MATLAB插值與曲線擬合,線性回歸 假設一組數(shù)據(jù)x=0 1 2 3 4 5和y=0 20 60 68 77 110，期中y是x的函數(shù)，使用最小二乘法建立一階線性回歸方程 用polyfit函數(shù)求出系數(shù) x=0 1 2 3 4 5; y=0 20 60 68 77 110; coef=polyfit(x,y,1),.SOUTHWEST JIAOTONG UNI

17、VERSITY,29,MATLAB插值與曲線擬合,線性回歸 畫出方程的直線及x和y的數(shù)據(jù)點 a1=coef(1);a0=coef(2); ybest=a1*x+a0; sum_sq=sum(y-ybest).2) plot(x,ybest,x,y,o);title(Linear regression estimate),grid 通過比較發(fā)現(xiàn)，最小平方誤差法獲得的誤差平方和要小于y=20 x的值,.SOUTHWEST JIAOTONG UNIVERSITY,30,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,31,MATLAB插值與曲線擬合,多項式回

18、歸 函數(shù)polyfit(x,y,n)中當n大于等于2時，稱為多項式回歸 高階多項式回歸有時并不如低階的更能代表原始數(shù)據(jù) 對上述數(shù)據(jù)的x和y分別求二、三、四、五階多項式回歸方程，畫出其擬合曲線 二階： x=0 1 2 3 4 5; y=0 20 60 68 77 110; coef=polyfit(x,y,2); ybest=coef(1)*x.2+coef(2)*x+coef(3); plot(x,ybest,x,y,o),title(Linear regression estimate),grid,.SOUTHWEST JIAOTONG UNIVERSITY,32,MATLAB插值與曲線擬合

19、,多項式回歸 三階： x=0 1 2 3 4 5; y=0 20 60 68 77 110; coef=polyfit(x,y,3); ybest=coef(1)*x.3+coef(2)*x.2+coef(3)*x+coef(4); plot(x,ybest,x,y,o),title(Linear regression estimate),grid 三階的擬合曲線最為光滑,.SOUTHWEST JIAOTONG UNIVERSITY,33,MATLAB插值與曲線擬合,多項式回歸 四階： x=0 1 2 3 4 5; y=0 20 60 68 77 110; coef=polyfit(x,y,4

20、); ybest=coef(1)*x.4+coef(2)*x.3+coef(3)*x.2 +coef(4)*x+coef(5); plot(x,ybest,x,y,o),title(Linear regression estimate),grid,.SOUTHWEST JIAOTONG UNIVERSITY,34,MATLAB插值與曲線擬合,多項式回歸 五階： x=0 1 2 3 4 5; y=0 20 60 68 77 110; coef=polyfit(x,y,5); ybest=coef(1)*x.5+coef(2)*x.4+coef(3)*x.3+coef(4)*x.2+coef(5)

21、*x+coef(6); plot(x,ybest,x,y,o),title(Linear regression estimate),grid 五階以上的多項式回歸方程會通過所有的原始數(shù)據(jù)點，但曲振蕩劇烈,.SOUTHWEST JIAOTONG UNIVERSITY,35,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,36,MATLAB插值與曲線擬合,多次式擬合及函數(shù)的計算 polyval也可用來做多項式函數(shù)計算，由polyfit計算出多項式的系數(shù)后，再由polyval計算出多項式的函數(shù)值 語法polyval(coef,x)，期中coef為多項式的各個

22、系數(shù)所構成的數(shù)組 x=0:0.1:1; y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; coef =polyfit(x,y,2) xi=0:0.01:1; yi=polyval(coef,xi); plot(x,y,o,xi,yi);,.SOUTHWEST JIAOTONG UNIVERSITY,37,MATLAB插值與曲線擬合,多次式擬合及函數(shù)的計算 舉例：對數(shù)據(jù)x=0 1 2 3 4 5和y=0 20 60 68 77 110做二階到九階的多項式回歸分析，然后計算多項式的函數(shù)值，繪制其曲線 x=0 1 2 3 4 5;

23、 y=0 20 60 68 77 110; newx=0:0.05:5; for n=2:9 f(:,n)=polyval(polyfit(x,y,n),newx); plot(newx,f(:,n),x,y,o);hold on; end,.SOUTHWEST JIAOTONG UNIVERSITY,38,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,39,MATLAB插值與曲線擬合,分段多項式函數(shù) 函數(shù)pp=mkpp(breaks,coefs)是區(qū)間breaks上根據(jù)多項式的系數(shù)coefs構建一個分段多項式pp breaks是一個長度為L+1的向

24、量，元素是每一個區(qū)間的起點與終點，共有L個區(qū)間 coefs是一個L(k+1)的矩陣 函數(shù)ppval()可估算分段多項式pp在特殊點xx上的值 v=ppval(pp,xx),.SOUTHWEST JIAOTONG UNIVERSITY,40,MATLAB插值與曲線擬合,分段多項式函數(shù) 舉例在區(qū)間-8 -4 0 4 8上分別畫出首尾相接的y=-0.5x+1；y=0.5x-1；y=-0.5x+1；y=0.5x-1四條線段 breaks=-8 -4 0 4 8; coefs=-0.5 1;0.5 -1;-0.5 1;0.5 -1; pp=mkpp(breaks,coefs) x=-8:0.1:8; y

25、=ppval(pp,x); plot(x,y,-b);,.SOUTHWEST JIAOTONG UNIVERSITY,41,MATLAB插值與曲線擬合,分段多項式函數(shù) 二次多項式y(tǒng)=-x2/4+x，在區(qū)間-8 -4 0 4 8上分別畫出y;-y;y;-y的曲線 breaks=-8 -4 0 4 8; coefs=-1/4 1 0; pp=mkpp(breaks,coefs;-coefs;coefs;-coefs) x=-8:0.1:8; y=ppval(pp,x); plot(x,y,K-);,.SOUTHWEST JIAOTONG UNIVERSITY,42,MATLAB插值與曲線擬合,分段

26、多項式函數(shù) 在區(qū)間-8 -4上應用函數(shù)pp=mkpp(breaks,coefs)構建二次多項式函數(shù)y=-0.25x2+x coefs=-0.25 1 0; pp=mkpp(-8 -4,coefs) x=-8:0.1:-4; y=ppval(pp,x); plot(x,y);,.SOUTHWEST JIAOTONG UNIVERSITY,43,MATLAB插值與曲線擬合,分段多項式函數(shù) 畫出分段函數(shù)在連續(xù)三個周期內(nèi)的圖像 x=0 0.075 0.1 0.175 0.2 0.275 0.3; y=40/3 0;-40 1;40/3 0;-40 1;40/3 0;-40 1; pp=mkpp(x,y

27、); xx=0:0.001:0.3; yy=ppval(pp,xx); plot(xx,yy,b-);,.SOUTHWEST JIAOTONG UNIVERSITY,44,MATLAB插值與曲線擬合,.SOUTHWEST JIAOTONG UNIVERSITY,45,MATLAB插值與曲線擬合,獲得分段多項式的詳細資料 使用函數(shù)breaks,coefs,l,k,d=unmkpp(pp)可以獲得分段多項式的詳細資料 Breaks是斷點數(shù)，coefs是系數(shù)，l是段數(shù)，k是項數(shù)，d是維數(shù)，pp是由于spine或mkpp所創(chuàng)建 為分段二項多項式-0.25x4+x和-(-0.25x4+x)創(chuàng)建一個描述

28、y=-0.25 1 0; pp=mkpp(-8 -4 0 4 8,y;-y;y;-y); breaks,coefs,l,k,d=unmkpp(pp),.SOUTHWEST JIAOTONG UNIVERSITY,46,MATLAB插值與曲線擬合,拉格拉日插值法 通過平面上不同兩點可以去定一條直線經(jīng)過這兩點，這就是拉格拉日插值問題，對于不在同一條直線的三個點得到哦的插值多項式為拋物線 舉例：已知cos30=sqrt(3)/2、cos45=sqrt(2)/2、 cos60=0.5、cos90=0，利用拉格拉日插值法計算cos(-40)、cos47、cos53、cos79、cos174的值，并給出插

29、值多項式 步驟： 1.編寫拉格拉日插值方法函數(shù)文件lag.m 2.編寫主函數(shù)文件lag_main.m 3.命令窗口運行l(wèi)ag_main,.SOUTHWEST JIAOTONG UNIVERSITY,47,MATLAB插值與曲線擬合,拉格拉日插值法 x=pi/4,pi/6,pi/3,pi/2; y=cos(pi/4),cos(pi/6),cos(pi/3),cos(pi/2); t=-40*pi/180,47*pi/180,53*pi/180,79*pi/180,174*pi/180; du=-40 47 53 79 174 yt=lag(x,y,t) yreal=cos(-40*pi/180)

30、 cos(47*pi/180) cos(53*pi/180) cos(79*pi/180) cos(174*pi/180) dy=yt-yreal yt=lag(x,y) ezplot(yt,-pi/4,pi) hold on ezplot(cos(t),-pi/4,pi);,.SOUTHWEST JIAOTONG UNIVERSITY,48,MATLAB插值與曲線擬合,牛頓插值法 牛頓插值法主要是基于插商的概念進行的，當增加插值基點時原來的計算結果對后面的計算仍然有用 舉例：已知零階Bvessel函數(shù)f(x)在若干點處的函數(shù)值如下， X 1.0 1.3 1.6 1.9 2.2 f(x) 0.

31、7651977 0.6200860 0.4554022 0.2818186 0.1103623 x=1.0 1.3 1.6 1.9 2.2; y=0.7651977 0.6200860 0.4554022 0.2818186 0.1103623; yt=niudun(x,y,1.5),.SOUTHWEST JIAOTONG UNIVERSITY,49,MATLAB插值與曲線擬合,牛頓插值法 5階勒讓德多項式為P5=(63X5-70X3+15X)/8，在-1,1上任意取9個點，如x=-1 -0.3 0.4 -0.7 0 -0.4 0.8 0.7 1，不管取點順序 1.根據(jù)5階勒讓德多項式，計算這9點的函數(shù)值； 2.依據(jù)這9點作牛頓插值，畫出插值多項式圖像，與勒讓德多項式圖像相比較； 3.計算x=0.24 -0.46 0.83時的近似結果，與真實結果相比。 步驟： 1.編寫牛頓插值方