matlab在科學計算中的應用7多項式、插值與數(shù)據(jù)擬合

1、第七章 多項式、插值與數(shù)據(jù)擬合 實驗數(shù)據(jù)的研究：數(shù)據(jù)差值和函數(shù)逼近 數(shù)據(jù)差值：在樣本點的基礎上求出不在 樣本點上的其他點處的函數(shù)值 函數(shù)逼近：由已知的樣本點數(shù)據(jù)求取能對其有較好擬合效果的函數(shù)表達式 多項式主要內(nèi)容： 多項式MATLAB命令 插值 Lagrange插值 Hermite插值 Runge現(xiàn)象和分段插值 分段插值 樣條插值的MATLAB表示 數(shù)據(jù)擬合 多項式擬合 函數(shù)線性組合的曲線擬合方法 最小二乘曲線擬合7.1 關于多項式MATLAB命令 一個多項式的冪級數(shù)形式可表示為：1121nnnnyc xc xc xc1231()nnyc xc xc xc xc112()()()nyc xrx

2、rxr也可表為嵌套形式或因子形式N階多項式n個根，其中包含重根和復根。若多項式所有系數(shù)均為實數(shù)，則全部復根都將以共軛對的形式出現(xiàn) Matlab實現(xiàn)： 在MATLAB里，多項式用行向量表示， 其元素為多項式的系數(shù)， 并從左至右按降冪排列（ 冪系數(shù) ）其matlab表示為 p=2 1 4 5轉換為符號表達式： poly2sym(p)32245yxxx例：roots: 求解多項式的零點 例： r=roots(p) 得到 r = 0.2500 + 1.5612i 0.2500 - 1.5612i -1.0000 所有零點由一個列向量給出。與多項式相關的Matlab函數(shù)：poly: 由零點可得原始多項式

3、的各系數(shù) （但可能相差一個常數(shù)倍） 例： poly(r)ans = 1.0000 0.5000 2.0000 2.5000注意：若存在重根，這種轉換可能會降低精度。665432(1)615201561yxxxxxxx 例： r=roots(1 -6 15 -20 15 -6 1) r = 1.0042 + 0.0025i 1.0042 - 0.0025i 1.0000 + 0.0049i 1.0000 - 0.0049i 0.9958 + 0.0024i 0.9958 - 0.0024i舍入誤差的影響，與計算精度有關。polyval: 計算多項式的值。例： c=3,-7,2,1,1; xi=2

4、.5,3; yi=polyval(c,xi) yi = 23.8125 76.00004323721yxxxx 例：已知，計算y(2.5)。求解： c=3,-7,2,1,1； xi=2.5; yi=polyval(c,xi) yi = 23.8125注：如果xi是含有多個橫坐標值的數(shù)組， 則yi也為與xi長度相同的向量。polyfit: ppolyfit(x,y,n) 求得滿足p(x)=y的n階多項式 其中x,y為n1向量320.20151.43852.74775.4370yxxx 例： x=1.1,2.3,3.9,5.1; y=3.887,4.276,4.651,2.117; a=polyf

5、it(x,y,length(x)-1)a = -0.2015 1.4385 -2.7477 5.4370 poly2sym(a) ans = -403/2000*x3+2877/2000*x2-27477/10000*x+5437/1000 Polyder: 求多項式一階導數(shù)的系數(shù)。 調(diào)用格式為: b=polyder(c ) c為多項式y(tǒng)的系數(shù)，b是微分后的系數(shù)，1121nnnnyc xc xc xc1212(1)nnnync xnc xc12,(1),nncnccb的值為：例如對多項式：conv(a,b)： c=conv(a,b)1121ddcabddyy yc xc xc xc1121mm

6、ammya xa xa xa1121nnbnnyb xb xb xb計算多項式的乘積，其中a，b分別為多項式 的系數(shù)向量abyy和ccy為乘積多項式 的系數(shù)向量deconv(a,b)： q,r=deconv(a,b)1121mmammya xa xa xa1121nnbnnyb xb xb xb計算多項式的除法運算，其中a，b分別為多項式 的系數(shù)向量abyy和aqbryy yybyay多項式除多項式其中qyry是商,是除法的余數(shù)qyq和r分別為和的系數(shù)向量ry 例 a=2,-5,6,-1,9; b=3,-90,-18; c=conv(a,b)c = 6 -195 432 -453 9 -792

7、 -162 q,r=deconv(c,b)q = 2 -5 6 -1 9r = 0 0 0 0 0 0 0 poly2sym(c) ans = 6*x6-195*x5+432*x4-453*x3+9*x2-792*x-162 7.2 插值7.2.1 Lagrange插值 方法介紹1111( )()nnjkkjkjjkxxy xyxx12,ny yy其對應的函數(shù)值為：則插值區(qū)間內(nèi)任意x的函數(shù)值y為：12,nx xx給定的n個插值點：n次Lagrange差值多項式function y=lagrange(x0,y0,x)ii=1:length(x0); y=zeros(size(x);for i=i

8、i ij=find(ii=i);y1=1; for j=1:length(ij) y1=y1.*(x-x0(ij(j); %, x-x_j ,ji的乘積 end y=y+y1*y0(i)/prod(x0(i)-x0(ij);endMatlab實現(xiàn)(n-1次差值)：11( )()nnjkkjkjj kxxy xyxx 算例：給出f(x)=ln(x)的數(shù)值表， 用Lagrange計算ln(0.54)的近似值。 x=0.4:0.1:0.8； % 取樣本點 y=log(x)ans = 0.916291,-0.693147,-0.510826,-0.356675,-0.223144 y0=lagrang

9、e(x,y,0.54,0.55,0.78)ans = -0.6161 -0.5978 -0.2484 y0(1)-log(0.54); %計算誤差ans = 4.352891839787265e-0057.2.2 Hermite插值 方法介紹 不少實際問題不但要求在節(jié)點上函數(shù)值相等，而且要求一階、二階甚至更高階導數(shù)值也相等，滿足這一要求的插值多項式就是Hermite插值多項式。 只討論函數(shù)值與一階導數(shù)值均相等的Hermite插值。1211( )()(2)1 () ;niiiiiiinnjiijjijijjijiy xhxxa yyyxxhaxxxx其 中12,nx xx12,ny yy12,n

10、y yy已知n個插值點 ：及對應的函數(shù)值：和一階導數(shù)值：則對插值區(qū)間內(nèi)任意x的函數(shù)值y的Hermite插值公式 MATLAB實現(xiàn)% hermite.mfunction y=hermite(x0,y0,y1,x)n=length(x0); m=length(x);for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a; end end yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i); end y(k

11、)=yy;end1211( )()(2)1 () ;niiiiiiinnjiijjijijj ij iy xh xxayyyx xhaxxxx其中x0,y0,y1分別為自變量、函數(shù)值及函數(shù)導數(shù)值向量 算例:對給定數(shù)據(jù),試構造Hermite多項式求出sin0.34的近似值。 x0=0.3,0.32,0.35; y0=sin(x0); y1=0.95534,0.94924,0.93937; format long; y=hermite(x0,y0,y1,0.34)y = 0.33348889007407 sin(0.34) 與精確值比較ans = 0.33348709214081 y-ansans

12、 = 5.406212721936754e-008 x=0.3:0.005:0.35;y=hermite(x0,y0,y1,x); plot(x,y) y2=sin(x); hold on plot(x,y2,-r)7.2.3 Runge現(xiàn)象 問題的提出：在應用插值多項式逼近函數(shù)時，是否插值多項式的次數(shù)越高近似效果就越好？ 答案：不一定21()1fxx反例： 該函數(shù)在區(qū)間-5,5上的各階導數(shù)存在，但在此區(qū)間上取n個節(jié)點所構成的Lagrange插值多項式在全區(qū)間內(nèi)并非都收斂。 xl=-4:1:4; yl=1./(1+xl.2); xh=-4:.5:4; yh=1./(1+xh.2); x0=-5

13、:0.1:5; y0=1./(1+x0.2); y01=lagrange(xl,yl,x0); y0h=lagrange(xh,yh,x0);繪制圖形 plot(x0,y0,r)原曲線 hold on plot(x0,y01,-b,x0,y0h,-g)插值曲線Matlab觀察結果：為解決Rung問題，引入分段插值。 算法分析：所謂分段插值就是通過插值點用折線或低次曲線連接起來逼近原曲線。 MATLAB實現(xiàn) 可調(diào)用內(nèi)部函數(shù)。 命令1 interp1 功能 ： 一維數(shù)據(jù)插值（表格查找）。該命令對數(shù)據(jù)點之間計算內(nèi)插值。它找出一元函數(shù)f(x)在中間點的數(shù)值。其中函數(shù)f(x)由所給數(shù)據(jù)決定。 格式1 y

14、i = interp1(x,Y,xi) %返回插值向量yi，每一元素對應于參量xi，同時由向量x與Y的內(nèi)插值決定。參量x指定數(shù)據(jù)Y的點。若Y為一矩陣，則按Y的每列計算。7.2.4 分段插值 temp=300,400,500,600; beta=1000*3.33,2.50,2.00,1.67; alpha=10000*0.2128,0.3605,0.5324,0.7190; ti=321,400,571; propty=interp1(temp,beta,alpha,ti); ti,proptyans = 1.0e+003 * 0.3210 3.1557 2.4382 0.4000 2.500

15、0 3.6050 0.5710 1.7657 6.6489例 : 對于t，beta 、alpha分別有兩組數(shù)據(jù)與之對應，用分段線性插值法計算當t=321, 440, 571時beta 、alpha的值。格式2 yi = interp1(Y,xi) %假定x=1:N，其中N為向量Y的長度，或矩陣Y的行數(shù)。格式3 yi = interp1(x,Y,xi,method) %用指定的算法計算插值： nearest：最近鄰點插值，直接完成計算； linear：線性插值（缺省方式），直接完成計算； spline：三次樣條函數(shù)插值。 cubic or pchip： 分段三次Hermite多項式插值。 其它，

16、如v5cubic ：Matlab5中的三次插值。 注：對于超出x范圍的xi的分量，插值算法nearest、 linear、v5cubic、將返回NaN。 其他的方法，interp1對超出的分量執(zhí)行外插值算法。 yi = interp1(x,Y,xi,method,extrap) 外插值算法 yi = interp1(x,Y,xi,method,extrapval) %確定超出x范圍的xi中的分量的外插值 extrapval，其值通常取NaN或0。 算例 year=1900:10:2010; product=75.995,91.972,105.711,123.203,131.669,. 150.

17、697,179.323,203.212,226.505,249.633,256.344,267.893; p1995 = interp1(year,product,1995)p1995 = 252.9885 x = 1900:1:2010; y = interp1(year,product,x,cubic); plot(year,product,o,x,y) 例:已知的數(shù)據(jù)點來自函數(shù)根據(jù)生成的數(shù)據(jù)進行插值處理，得出較平滑的曲線直接生成數(shù)據(jù)。求解： x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); plot(x,y,x,y,o)Step1：根據(jù)函數(shù)直接生成

18、數(shù)據(jù)，作圖默認為線性差值Step2：Step3：誤差檢測 legend(linear,cubic,spline,nearest,y,y0) max(abs(y0-y1),max(abs(y0-y2), max(abs(y0-y3),max(abs(y0-y4)ans = 0.0614 0.0177 0.0086 0.1598 x1=0:.02:1; y0=(x1.2-3*x1+5).*exp(-5*x1).*sin(x1); y1=interp1(x,y,x1); y2=interp1(x,y,x1,cubic); y3=interp1(x,y,x1,spline); y4=interp1(x

19、,y,x1,nearest); plot(x1,y1,y2,y3,y4,:,x,y,o,x1,y0) x0=-1+2*0:10/10; y0=1./(1+25*x0.2); %取樣本點 x=-1:.01:1; 取插值點 y=lagrange(x0,y0,x); % Lagrange 插值 ya=1./(1+25*x.2); plot(x,ya,x,y,:) 例: y1=interp1(x0,y0,x,cubic); y2=interp1(x0,y0,x,spline); plot(x,ya,x,y1,:,x,y2,-)命令2 interp2 功能 二維數(shù)據(jù)內(nèi)插值（表格查找） 格式1 ZI =

20、interp2(X,Y,Z,XI,YI) %返回矩陣ZI，其元素包含對應于參量XI與YI的元素（可以是向量、或 同型矩陣）。參量X與Y必須是單調(diào)的，且相同的劃分格式，就像由命令meshgrid生成的一樣。若Xi與Yi中有在X與Y范圍之外的點，則相應地返回NaN。 格式2 ZI = interp2(Z,XI,YI) %默認X=1:n,Y=1:m,其中m,n=size(Z) 格式3 ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method計算二維插值： linear：雙線性插值算法（缺省算法）； nearest：最臨近插值； spline：三次樣條插值； cub

21、ic：雙三次插值。 算例： years=1950:10:1990; service=10:10:30; wage = 150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243; w = interp2(service,years,wage,15,1975)w = 190.6288 例Step1：生成樣本點 x,y=meshgrid(-3:.6:3,-2:.4:2); z=(x.2-2*x).*exp(-x.2-y.

22、2-x.*y); surf(x,y,z), axis(-3,3,-2,2,-0.7,1.5)求解：Step2：選較密的插值點，用默認的線性插值算法進行插值 x1,y1=meshgrid(-3:.2:3,-2:.2:2); z1=interp2(x,y,z,x1,y1); surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)Step3：用立方和樣條插值方法： z1=interp2(x,y,z,x1,y1,cubic); z2=interp2(x,y,z,x1,y1,spline); surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5) figu

23、re;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5)Step4：算法誤差的比較 z=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1); surf(x1,y1,abs(z-z1),axis(-3,3,-2,2,0,0.08) figure;surf(x1,y1,abs(z-z2),axis(-3,3,-2,2,0,0.025)二維一般分布數(shù)據(jù)的插值格式:z=griddata(x0,y0,z0,x,y,method）功能:可對非網(wǎng)格數(shù)據(jù)進行插值 Matlab函數(shù)：griddata（）x0,y0,z0:已知樣本點坐標(不一定是網(wǎng)格型的)x,y:期

24、望的插值位置(單個點、向量或網(wǎng)格型矩陣)method選項:表示插值的算法,其選項有: v4:MATLAB4.0提供的插值算法，公認效果較好; linear：雙線性插值算法(缺省算法); nearest：最臨近插值； spline：三次樣條插值； cubic：雙三次插值。222( , )(2 )xyxyzf x yxx ex為-3,3,y為2,2區(qū)域內(nèi)隨機選擇的一組坐標.分別用v4與cubic插值法進行處理，并對誤差進行比較。 例：考慮函數(shù),其中解法:Step1: 取樣本點 x=-3+6*rand(200,1);y=-2+4*rand(200,1); z=(x.2-2*x).*exp(-x.2-

25、y.2-x.*y);().*( , )a,bmnabarand m n生成區(qū)間內(nèi)的階隨機矩陣Step2: 取期望插值的點并求解 x1,y1=meshgrid(-3:.2:3,-2:.2:2); z1=griddata(x,y,z,x1,y1,cubic); surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5) z2=griddata(x,y,z,x1,y1,v4); figure;surf(x1,y1,z2), axis(-3,3,-2,2,-0.7,1.5)Step3:誤差分析 z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1); surf

26、(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.15) figure;surf(x1,y1,abs(z0-z2),axis(-3,3,-2,2,0,0.15) 考察隨機樣本點的分布狀況考察隨機樣本點的分布狀況: plot(x,y,x) % 樣本點的二維分布 figure, plot3(x,y,z,x), axis(-3,3,-2,2,-0.7,1.5)，gridx=-3+6*rand(200,1); y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-x.*y); 分布比較均勻分布比較均勻 例：222( , )(2 )xyxyzf

27、 x yxx e在x為3,3, y為2,2矩形區(qū)域隨機選擇一組分布不均勻數(shù)據(jù)，進行插值分析。求解:Step1:重新生成分布均勻的隨機樣本點 x=-3+6*rand(200,1); y=-2+4*rand(200,1); z=(x.2-2*x).*exp(-x.2-y.2-x.*y);Step2:去除在(-1,-1/2)點為圓心,以0.5為半徑的圓內(nèi)的點, 即得出分布不均勻的樣本點 ii=find(x+1).2+(y+0.5).20.52); x=x(ii); y=y(ii); z=z(ii); plot(x,y,x) t=0:.1:2*pi,2*pi; x0=-1+0.5*cos(t); y0

28、=-0.5+0.5*sin(t); line(x0,y0) % 在圖形上疊印該圓，可見，圓內(nèi)樣本點均已剔除Step3:用新樣本點擬合出曲面 x1,y1=meshgrid(-3:.2:3, -2:.2:2); z1=griddata(x,y,z,x1,y1,v4); surf(x1,y1,z1), axis(-3,3,-2,2,-0.7,1.5)Step4:誤差分析 z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1); surf(x1,y1,abs(z0-z1), axis(-3,3,-2,2,0,0.1) contour(x1,y1,abs(z0-z1),30); h

29、old on, plot(x,y,x); line(x0,y0) 誤差的二維等高線圖 擬合效果仍很好擬合效果仍很好,但樣本點稀疏的地方效果較差但樣本點稀疏的地方效果較差 命令3 interp3 調(diào)用格式同interp2( )函數(shù)一致函數(shù)一致;三維網(wǎng)格生成用meshgrid( )函數(shù),調(diào)用格式: x,y,z=meshgrid(x1,y1,z1) 其中x1,y1,z1為這三維所需要的分割形式，應以向量形式給出，返回x,y,z為網(wǎng)格的數(shù)據(jù)生成，均為三維數(shù)組。 griddata3( ) 三維非網(wǎng)格形式的插值擬合; 調(diào)用格式同griddata( )函數(shù)一致函數(shù)一致 命令4 interpn調(diào)用格式同int

30、erp2( )函數(shù)一致函數(shù)一致;n維網(wǎng)格生成用ndgrid( )函數(shù),調(diào)用格式: x1,x2,xn=ndgridv1,v2,vn griddatan( ) n維非網(wǎng)格形式的插值擬合 調(diào)用格式同griddata( )函數(shù)一致函數(shù)一致 例：22222( , , )cos()x z y x z yV x y zex yzz yx x,y,z=meshgrid(-1:0.2:1); x0,y0,z0=meshgrid(-1:0.05:1); V=exp(x.2.*z+y.2.*x+z.2.*y).*cos(x.2.*y.*z+z.2.*y.*x); V0=exp(x0.2.*z0+y0.2.*x0+z

31、0.2.*y0).*cos(x0.2.*y0.*z0+z0.2.*y0.*x0); V1=interp3(x,y,z,V,x0,y0,z0,spline); err=V1-V0; max(err(:)ans = 0.0419通過函數(shù)生成一些網(wǎng)格型樣本點，試根據(jù)樣本點進行擬合，并給出擬合誤差。7.2.5樣條插值的MATLAB表示011:nnaxxxxb樣條函數(shù)的概念 樣條(spline)：原指工程設計中使用的一種繪圖工具(如富有彈性的細木條或金屬條),應用其將已知點連接成一條光滑曲線(稱為樣條曲線),并使連接點處有連續(xù)的曲率.三次樣條插值即由此抽象出來. 樣條函數(shù)：數(shù)學上將具有一定光滑性的分段多

32、項式函數(shù)稱為樣條函數(shù). 具體的,給定區(qū)間a,b的一個劃分:1S( ): (1). ,(0,1, -1)S( )(2).S( ) , -1iixx xinxmxa bmxmm如果函數(shù)滿足每個小區(qū)間 上是 次多項式; 在 上具有階連續(xù)導數(shù).則稱S( )為關于劃分 的 次樣條函數(shù),其圖像為 次樣條曲線. 三次樣條插值 - 樣條插值: 利用樣條函數(shù)進行插值 (例如分段線性插值是一次樣條插值)- 三次樣條插值: 已知函數(shù)f(x)在區(qū)間a,b上的n個節(jié)點的函數(shù)值( )(1,2,)iiyf xin, 求插值函數(shù)S(x)滿足:iS- 定義三次樣條函數(shù)類：S=csapi(x,y) 其中x=x1,x2,.,xn,

33、 y=y1,y2,yn為樣本點。 S返回樣條函數(shù)對象的插值結果，包括子區(qū)間點、各區(qū)間點三次多項式系數(shù)等。- 可用 fnplt()繪制出插值結果，其調(diào)用格式：fnplt(S)- 對給定的向量xp, 可用fnval()函數(shù)計算, 其調(diào)用格式： yp=fnval(S,xp) 其中得出的yp是xp上各點的插值結果。 Matlab 樣條插值工具箱 例：已知y0=sin(x0),x0= =0,0.4,1,2,pi, 求該函數(shù) 三次樣條插值結果解法: x0=0,0.4,1,2,pi; y0=sin(x0); sp=csapi(x0,y0), fnplt(sp, r:); hold on,sp = form:

34、 pp breaks: 0 0.4000 1 2 3.1416 coefs: 4x4 double pieces: 4 order: 4 dim: 1 ezplot(sin(t),0,pi); plot(x0,y0,o)- 在(0.4000, 1)區(qū)間內(nèi)，插值多項式可以表示為：322( )0.1627(0.4)0.1876(0.4)0.9245(0.4)0.3894S xxxx- 查看插值多項式 sp.coefsans = -0.1627 0.0076 0.9965 0 -0.1627 -0.1876 0.9245 0.3894 0.0244 -0.4804 0.5238 0.8415 0.0

35、244 -0.4071 -0.3637 0.9093例:點，用三次樣條插值的方法對這些數(shù)據(jù)進行擬合解法: x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); sp=csapi(x,y); fnplt(sp)Columns 7 through 12 0.4800 0.6000 -0.2404 0.7652 -0.5776 0.1588 0.6000 0.7200 -0.4774 0.6787 -0.4043 0.1001 0.7200 0.8400 -0.4559 0.5068 -0.2621 0.0605 0.8400 0.9600 -0.4559 0.3

36、427 -0.1601 0.0356 c=sp.breaks(1:4) sp.breaks(2:5) sp.coefs(1:4,:),. sp.breaks(5:8) sp.breaks(6:9) sp.coefs(5:8,:) c = Columns 1 through 6 0 0.1200 24.7396 -19.3588 4.5151 0 0.1200 0.2400 24.7396 -10.4526 0.9377 0.3058 0.2400 0.3600 4.5071 -1.5463 -0.5022 0.3105 0.3600 0.4800 1.9139 0.0762 -0.6786 0

37、.2358- 格式 S=csapi(x1,x2,xn,z)注: csapi() 可處理多個自變量的網(wǎng)格數(shù)據(jù)三次樣條插值類 x0=-3:.6:3; y0=-2:.4:2; x,y=ndgrid(x0,y0); % 注意這里只能用 ndgrid,否則生成的 z 矩陣順序有問題 z=(x.2-2*x).*exp(-x.2-y.2-x.*y); sp=csapi(x0,y0,z); fnplt(sp);例:對函數(shù)解法: 三次樣條數(shù)據(jù)插值-Matlab實現(xiàn): 函數(shù)spline格式 yy = spline(x,y,xx) 例:對離散分布在y=exp(x)sin(x)函數(shù)曲線上的數(shù)據(jù)點進行樣條插值計算： x

38、 = 0 2 4 5 8 12 12.8 17.2 19.9 20; y = exp(x).*sin(x); xx = 0:.25:20; yy = spline(x,y,xx); plot(x,y,o,xx,yy)B樣條函數(shù)及其MATLAB表示 S=spapi(k,x,y)- Matlab實現(xiàn): x0=0,0.4,1,2,pi; y0=sin(x0); ezplot(sin(t),0,pi); hold on sp1=csapi(x0,y0); fnplt(sp1,-); % 三次分段多項式樣條插值例例: :y=sin(t)和解法: sp2=spapi(5,x0,y0); fnplt(sp2

39、,r:) % 5 次 B 樣條插值 x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); ezplot(x2-3*x+5)*exp(-5*x)*sin(x),0,1), hold on sp1=csapi(x,y); fnplt(sp1,-); sp2=spapi(5,x,y); fnplt(sp2,r:)B樣條的插值效果好于三次分段多項式7.2.6 基于樣條插值的數(shù)值微積分運算 基于樣條插值的數(shù)值微分運算格式： Sd=fnder(S,k)該函數(shù)可以求取S的k階導數(shù)。格式： Sd=fnder(S,k1,kn)可以求取多變量函數(shù)的偏導數(shù) 例：解法： syms

40、x; f=(x2-3*x+5)*exp(-5*x)*sin(x); ezplot(diff(f),0,1), hold on 解析解微分圖 x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x);取樣本點 sp1=csapi(x,y);建立三次樣條函數(shù) dsp1=fnder(sp1,1);數(shù)值微分 fnplt(dsp1, r-)繪制樣條圖 sp2=spapi(5,x,y);5階次B樣條 dsp2=fnder(sp2,1); 數(shù)值微分 fnplt(dsp2,g:); axis(0,1,-0.8,5) 例：解法：Step1：擬合曲面 x0=-3:.3:3; y0=-

41、2:.2:2; x,y=ndgrid(x0,y0); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); sp=spapi(5,5,x0,y0,z); B樣條dspxy=fnder(sp,1,1); fnplt(dspxy)%生成樣條圖理論方法： syms x y; z=(x2-2*x)*exp(-x2-y2-x*y); ezsurf(diff(diff(z,x),y),-3 3,-2 2)對符號變量表達式做三維表面圖 基于樣條插值的數(shù)值積分運算格式： f=fnint(S)其中S為樣條函數(shù)。解法： x=0,0.4,1 2,pi; y=sin(x); sp1=csapi(x,y);

42、 a=fnint(sp1); 建立三次樣條函數(shù)并積分 xx=fnval(a,0,pi); xx(2)-xx(1)ans = 2.01910sin xdx 例：考慮中較稀疏的樣本點,用樣條積分的方式求出定積分及積分函數(shù)。Step2: sp2=spapi(5,x,y); b=fnint(sp2); 建立B樣條函數(shù)并積分 xx=fnval(b,0,pi); xx(2)-xx(1)ans = 1.9999Step3: 繪制曲線 ezplot(-cos(t)+2,0,pi); hold on fnplt(a,r-); fnplt(b,g:)7.3 數(shù)據(jù)擬合 用用插值的方法插值的方法對一函數(shù)進行近似對一函

43、數(shù)進行近似,要求所得到的要求所得到的插值多項式經(jīng)過已知插值節(jié)點插值多項式經(jīng)過已知插值節(jié)點;在在n比較大的情況比較大的情況下下,插值多項式往往是高次多項式插值多項式往往是高次多項式,這也就容易出這也就容易出現(xiàn)振蕩現(xiàn)象（現(xiàn)振蕩現(xiàn)象（Runge現(xiàn)象），即雖然在插值節(jié)點現(xiàn)象），即雖然在插值節(jié)點上沒有誤差上沒有誤差,但在插值節(jié)點之外插值誤差變得很但在插值節(jié)點之外插值誤差變得很大大,從從“整體整體”上看上看,插值逼近效果將變得插值逼近效果將變得“很很差差”。數(shù)據(jù)擬合數(shù)據(jù)擬合, 又稱函數(shù)逼近又稱函數(shù)逼近, 求近似函數(shù)的一類數(shù)值方法。求近似函數(shù)的一類數(shù)值方法。 所謂所謂數(shù)據(jù)擬合數(shù)據(jù)擬合是求一個簡單的函數(shù)是求一

44、個簡單的函數(shù),例如是一個例如是一個低次多項式低次多項式,不要求通過已知的這些點不要求通過已知的這些點,而是要求而是要求在整體上在整體上“盡量好盡量好”的逼近原函數(shù)的逼近原函數(shù)。這時。這時,在每在每個已知點上就會有誤差個已知點上就會有誤差,數(shù)據(jù)擬合就是從整體上數(shù)據(jù)擬合就是從整體上使誤差使誤差,盡量的小一些。盡量的小一些。 例：給定一組實驗數(shù)據(jù),如下表ixiyi143226481.12.84.97.2086422468求x，y的函數(shù)關系。解：做草圖(見左圖), 由此可設10ya xa困難: 無論a1,a0取何值,直線都不能經(jīng)過所有樣本點.問題: 選取適當a1,a0的值,使直線能最好的反映數(shù)據(jù)點的基

45、本趨勢. 殘差：近似值與樣本點函數(shù)值（已知數(shù)據(jù)）之差*iiiyy2minii|miniimax | minii 常用準則：1). 殘差絕對值之和最?。?). 殘差最大絕對值最?。?). 殘差平方和最?。?)使用不方便；2) 準則求近似函數(shù)的方法稱為函數(shù)的最佳一致逼近3) 準則求近似函數(shù)的方法稱為最佳平方逼近，或數(shù)據(jù)擬合的最小二乘法7.3.1 多項式擬合1121( )nnng xc xc xc( ,), 1,2,iixyiL n次多項式： 已知數(shù)據(jù)點：( ),1,2,iiiryg xiL2121111LLnnkiikiiikRryc x 0,1,2,1jRjnc 曲線與數(shù)據(jù)點的殘差為：jc 最小

46、二乘法函數(shù)逼近，即取適當 使得最小.由極值原理，得方程組： 數(shù)學算法 或 1221111(),1,2,1nLLnj knkijiijiixcxy kn 221111112111211110111.LLLLnnnniiiiiiiiiLLLnnniiiiiiinLLLniiiiiixxxx ycxxxyccxxy 或矩陣形式：多項式擬合MATLAB命令：polyfit格式：p=polyfit(x,y,n) Matlab實現(xiàn)：Step1:多項式擬合 x0=0:.1:1; y0=(x0.2-3*x0+5).*exp(- 5*x0).*sin(x0); p3=polyfit(x0,y0,3); vpa(

47、poly2sym(p3),10) % 可以如下顯示多項式ans =2.839962923*x3-4.789842696*x2+1.943211631*x+.5975248921e-1例：例：解法：Step2：繪制擬合曲線并與精確解曲線比較： x=0:.01:1; ya=(x.2-3*x+5).*exp(-5*x).*sin(x); y1=polyval(p3,x); plot(x,y1,x,ya,x0,y0,o) legend(擬合曲線,精確解曲線,樣本點)Step3：就不同的次數(shù)進行擬合： p4=polyfit(x0,y0,4); y2=polyval(p4,x); p5=polyfit(x

48、0,y0,5); y3=polyval(p5,x); p8=polyfit(x0,y0,8); y4=polyval(p8,x); plot(x,ya,x0,y0,o,x,y2,x,y3,x,y4) 最小二乘8次多項式擬合效果較好注：察看擬合最高次數(shù)為8的多項式： vpa(poly2sym(p8),5)ans =-8.2586*x8+43.566*x7-101.98*x6+140.22*x5-125.29*x4+74.450*x3-27.672*x2+4.9869*x +.42037e-6 Taylor冪級數(shù)展開： syms x; y=(x2-3*x+5)*exp(-5*x)*sin(x);

49、vpa(taylor(y,9),5)ans =5.*x-28.*x2+77.667*x3-142.*x4+192.17*x5-204.96*x6+179.13*x7-131.67*x8 多項式表示數(shù)據(jù)模型是不唯一的，即是兩個多項式函數(shù)完全不同。在某一區(qū)域內(nèi)其曲線可能特別近似解法：Step1：取樣本點和精確數(shù)值解 x0=-1+2*0:10/10; y0=1./(1+25*x0.2); x=-1:.01:1; ya=1./(1+25*x.2);Step2：各次多項式擬合 p3=polyfit(x0,y0,3); y1=polyval(p3,x); p5=polyfit(x0,y0,5); y2=p

50、olyval(p5,x); p8=polyfit(x0,y0,8); y3=polyval(p8,x); p10=polyfit(x0,y0,10); 例：例： y4=polyval(p10,x); plot(x,ya,x,y1,x,y2,-.,x,y3,-,x,y4,:) legend(精確解曲線,3次多項式擬合,5次多項式擬合,8次多項式擬合,10次多項式擬合)多項式擬合的效果并不一定總是很精確的。 考察Taylor冪級數(shù)展開擬合效果 syms x; y=1/(1+25*x2); p=taylor(y,x,10)p =1-25*x2+625*x4-15625*x6+390625*x8 x1

51、=-1:0.01:1; ya=1./(1+25*x1.2); y1=subs(p,x,x1); plot(x1,ya,-,x1,y1)Taylor多項式擬合效果更差7.3.2 函數(shù)線性組合的曲線擬合方法且線性無關該方程的最小二乘解為:其中例例 x=0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5; y=2.88;2.2576;1.9683;1.9258;2.0862;2.109; 2.1979;2.5409;2.9627;3.155;3.2052; A=ones(size(x)，exp(-3*x), cos(-2*x).*exp(-4*x) ，x.2; c=Ay; c1=cc1 = 1.2200 2.3397 -0.6797 0.8700 圖形顯示 x0=0:0.01:1.5; A1=ones(size(x0) exp(-3*x0), cos(-2*x0)