第8章數據插值與擬合11/26/20221第8章數據插值與擬合11/22/202218.1數據插值一維數據的插值問題二維網格數據的插值問題二維一般分布數據的插值問題高維插值問題11/26/202228.1數據插值一維數據的插值問題11/22/202228.1.1一維數據的插值問題interp1函數spline函數lagrange插值hermite插值11/26/202238.1.1一維數據的插值問題interp1函數11/22/1、interp1函數interp1函數的調用格式：y=interp1(x0,y0,x,method)y=interp1(x0,y0,x,method,'extrap')*臨近點插值：method=‘nearest’*線性插值：method=‘linear’*三次樣條插值：method=‘spline’*立方插值：method=‘pchip’or‘cubic’11/26/202241、interp1函數interp1函數的調用格式：*臨近點1、interp1函數選擇插值方法時主要考慮因素：運算時間、占用計算機內存和插值的光滑程度。運算時間占用計算機內存光滑程度臨近點插值快少差線性插值稍長較多稍好三次樣條插值最長較多最好立方插值較長多較好插值方法比較11/26/202251、interp1函數選擇插值方法時主要考慮因素：運算時間、例:已知的數據點來自函數

根據生成的數據進行插值處理，得出較平滑的曲線

直接生成數據。數據點：>>x=0:0.12:1;>>y=(x.^2-3*x+5).*exp(-5*x).*sin(x);>>plot(x,y,x,y,'o')11/26/20226例:已知的數據點來自函數

根據生成的數據進行插值處理，得調用interp1()函數：>>x0=0:0.02:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);>>y1=interp1(x,y,x0);y2=interp1(x,y,x0,'cubic');>>y3=interp1(x,y,x0,'spline');y4=interp1(x,y,x0,'nearest');>>plot(x0,[y1',y2',y3',y4'],':',x,y,'o',x0,y0)誤差分析>>[max(abs(y0(1:49)-y2(1:49))),max(abs(y0-y3)),max(abs(y0-y4))]ans=0.01770.00860.15981、interp1函數11/26/20227調用interp1()函數：1、interp1函數11/2功能三次樣條數據插值格式

yy=spline(x,y,xx)

例：對離散分布在y=exp(x)sin(x)函數曲線上的數據點進行樣條插值計算：>>x=[024581212.817.219.920];>>y=exp(x).*sin(x);>>xx=0:.25:20;>>yy=spline(x,y,xx);>>plot(x,y,'o',xx,yy)2、spline函數11/26/20228功能三次樣條數據插值2、spline函數11/22/203、lagrange插值方法介紹對給定的n個插值點及對應的函數值，利用構造的n-1次lagrange插值多項式，則對插值區間內任意x的函數值y可通過下式求的：11/26/202293、lagrange插值方法介紹11/22/20229lagrange插值函數functiony=lagrange(x0,y0,x)ii=1:length(x0);y=zeros(size(x));fori=iiij=find(ii~=i);y1=1;forj=1:length(ij),y1=y1.*(x-x0(ij(j)));endy=y+y1*y0(i)/prod(x0(i)-x0(ij));end11/26/202210lagrange插值函數functiony=lagrangLagrange插值>>x0=-1+2*[0:10]/10;y0=1./(1+25*x0.^2);>>x=-1:0.01:1;y=lagrange(x0,y0,x);>>ya=1./(1+25*x.^2);>>plot(x,ya,x,y,':')調用interp1函數>>y1=interp1(x0,y0,x,'cubic');y2=interp1(x0,y0,x,'spline');>>plot(x,ya,x,y1,':',x,y2,'--')3、lagrange插值例：對進行Lagrange插值11/26/202211Lagrange插值3、lagrange插值例：對許多實際插值問題中，為使插值函數能更好地和原來的函數重合，不但要求二者在節點上函數值相等，而且還要求相切，對應的導數值也相等，甚至要求高階導數也相等。這類插值稱作切觸插值，或埃爾米特(Hermite)插值。滿足這種要求的插值多項式就是埃爾米特插值多項式4、埃爾米特插值(hermite插值)11/26/202212許多實際插值問題中，為使插值函數能更好地和原4、埃爾米特插值(hermite插值)方法介紹已知n個插值點及對應的函數值和一階導數值。則對插值區間內任意x的函數值y的Hermite插值公式：11/26/2022134、埃爾米特插值(hermite插值)方法介紹11/22/2hermite插值函數functiony=hermite(x0,y0,y1,x)n=length(x0);m=length(x);fork=1:myy=0.0;fori=1:nh=1.0;a=0.0;forj=1:nifj~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;end11/26/202214hermite插值函數functiony=hermite(例：利用Hermite插值法求sin0.34的近似值。>>x0=[0.3,0.32,0.35];y0=[0.29552,0.31457,0.34290];>>y1=[0.95534,0.94924,0.93937];>>y=hermite(x0,y0,y1,0.34)y=0.3335>>sin(0.34)％與精確值比較ans=0.33354、hermite插值11/26/202215例：利用Hermite插值法求sin0.34的近似值。4、h8.1.2二維網格數據的插值問題二維插值函數interp2的調用格式：zi=interp2(x0,y0,z0,xi,yi)zi=interp2(x0,y0,z0,xi,yi,method)*臨近點插值：method=‘nearest’*線性插值：method=‘linear’（缺省算法）*三次樣條插值：method=‘spline’*立方插值：method=‘pchip’or‘cubic’11/26/2022168.1.2二維網格數據的插值問題二維插值函數interp2例由二元函數獲得一些較稀疏的網格數據，對整個函數曲面進行各種插值，并比較插值結果繪制已知數據的網格圖>>[x,y]=meshgrid(-3:0.6:3,-2:0.4:2);>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);>>surf(x,y,z);axis([-3,3,-2,2,-0.7,1.5])8.1.2二維網格數據的插值問題11/26/202217例由二元函數8.1.2二維網格數據的插值問題11/22/默認線性插值算法進行插值>>[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])立方和樣條插值：>>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])>>figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])11/26/202218默認線性插值算法進行插值11/22/202218算法誤差的比較>>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])11/26/202219算法誤差的比較11/22/2022198.1.3二維一般分布數據的插值問題griddata函數的調用格式：z=griddata(x0,y0,z0,x,y,method)method=‘v4’：插值算法，公認效果較好*臨近點插值：method=‘nearest’*線性插值：method=‘linear’（缺省算法）*三次樣條插值：method=‘spline’*立方插值：method=‘cubic’11/26/2022208.1.3二維一般分布數據的插值問題griddata函數的例在x為［-3，3］，y為［-2，2］矩形區域隨機選擇一組數據點，用’v4’與’cubic’插值法進行處理，并對誤差進行比較。已知數據點>>x=-3+6*rand(200,1);y=-2+4*rand(200,1);>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);>>plot(x,y,'x')%樣本點的二維分布>>figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5])，grid11/26/202221例cubic和v4算法>>[x1,y1]=meshgrid(-3:0.2:3,-2:0.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])誤差分析>>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.15])>>figure;surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15])11/26/202222cubic和v4算法11/22/2022228.1.4高維插值問題三維插值interp3函數的調用格式：三維網格[x,y,z]=meshgrid(x1,y1,z1)griddata3()三維非網格數據插值n維插值interpn函數n維網格[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]griddatan()n維非網格數據插值interp3()、interpn()調用格式同interp2()一致griddata3()、griddatan()調用格式同griddata()一致11/26/2022238.1.4高維插值問題三維插值interp3函數的調用格式例：通過函數生成一些網格型樣本點，據此進行插值并給出插值誤差。>>[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+z0.^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.041911/26/202224例：8.2數據擬合多項式擬合函數線性組合的曲線擬合方法最小二乘曲線擬合11/26/2022258.2數據擬合多項式擬合11/22/202225多項式標準形式：多項式擬合函數polyfit，其調用格式p=polyfit(x,y,n)其中x和y為原始樣本點構成的向量n為多項式的階次p為多項式系數按降冪排列得到的向量8.2.1多項式擬合11/26/202226多項式標準形式：8.2.1多項式擬合11/22/20222三次多項式擬合>>x0=0:.1:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);>>p3=polyfit(x0,y0,3);vpa(poly2sym(p3),10)ans=2.839962923*x^3-4.789842696*x^2+1.943211631*x例:已知的數據點來自函數

用多項式擬合的方法在不同階次下進行擬合11/26/202227三次多項式擬合例:已知的數據點來自函數

用多項式擬合的方繪制擬合曲線：>>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')11/26/202228繪制擬合曲線：11/22/202228不同次數進行擬合：>>p4=polyfit(x0,y0,4);y2=polyval(p4,x);>>p5=polyfit(x0,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)11/26/202229不同次數進行擬合：11/22/202229擬合次數為8的多項式：>>vpa(poly2sym(p8),5)ans=-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-125.29*x^4+74.450*x^3-27.672*x^2+4.9869*x+.42037e-6Taylor冪級數展開：>>symsx;y=(x^2-3*x+5)*exp(-5*x)*sin(x);>>vpa(taylor(y,9),5)ans=5.*x-28.*x^2+77.667*x^3-142.*x^4+192.17*x^5-204.96*x^6+179.13*x^7-131.67*x^8多項式表示數據模型是不唯一的，即是兩個多項式函數完全不同。在某一區域內其曲線可能特別近似。11/26/202230擬合次數為8的多項式：11/22/202230編制程序>>x0=-1+2*[0:10]/10;y0=1./(1+25*x0.^2);>>x=-1:.01:1;ya=1./(1+25*x.^2);>>p3=polyfit(x0,y0,3);y1=polyval(p3,x);>>p5=polyfit(x0,y0,5);y2=polyval(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,':')例：對多項式擬合，多項式擬合的效果并不一定總是很精確的。11/26/202231編制程序例：對用Taylor冪級數展開效果將更差。>>symsx;y=1/(1+25*x^2);p=taylor(y,x,10)p=1-25*x^2+625*x^4-15625*x^6+390625*x^8多項式擬合效果>>x1=-1:0.01:1;ya=1./(1+25*x1.^2);>>y1=subs(p,x,x1);>>plot(x1,ya,'--‘,x1,y1)11/26/202232用Taylor冪級數展開效果將更差。11/22/2022328.2.2函數線性組合的曲線擬合方法已知某函數的線性組合為：g(x)=c1f1(x)+c2f2(x)+…+cnfn(x)其中f1(x),f2(x),…,cnfn(x)為已知函數，c1,c2,…,cn為待定系數假設已經測出數據(x1,y1),(x2,y2),…,(xm,ym)，則可建立如下的線性方程：Ac=y11/26/2022338.2.2函數線性組合的曲線擬合方法已知某函數的線性組合為11/26/20223411/22/202234例假設測出一組實驗數據(xi,yi)，其函數原型為y(x)=c1+c2e-3x+c3cos(-2x)e-4x+c4x2，用已知數據求出待定系數ci的值x00.20.40.70.90.92y2.882.25761.96831.92582.08622.109x0.991.21.41.481.5y2.19792.54092.96273.1553.205211/26/202235例假設測出一組實驗數據(xi,yi)，其函數原型為y(x擬合ci參數>>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=A\y;c1=c'c1=1.22002.3397-0.67970.870011/26/202236擬合ci參數11/22/202236圖形顯示>>x0=[0:0.01:1.5]';>>A1=[ones(size(x0)),exp(-3*x0),cos(-2*x0).*exp(-4*x0)x0.^2];>>y1=A1*c;>>plot(x0,y1,x,y,'x')11/26/202237圖形顯示11/22/202237例對下列實驗數據進行函數擬合已知數據點>>x=[1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,2.2255,2.4596,2.7183,3.6693];>>y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,0.2864,0.2532,0.2238,0.1546];>>plot(x,y,x,y,'*')x1.10521.22141.34991.49181.64871.8221y0.67950.60060.53090.46930.41480.3666x2.01382.22552.45962.71833.6693y0.324140.28650.25320.22380.154611/26/202238例對下列實驗數據進行函數擬合x1.10521.22141對x,y進行對數變換：>>x1=log(x);y1=log(y);plot(x1,y1)線性函數擬合法>>A=[x1',ones(size(x1'))];c=[A\y1']‘c=-1.2339-0.2630>>exp(c(2))ans=0.768711/26/202239對x,y進行對數變換：11/22/202239例對f(x)=(x2-3x+5)e-5xsinx進行多項式擬合，可以選擇各個函數為fi(x)=xn+1-i，i=1,2,…,n，并觀察多項式擬合的結果編制程序：x=[0:0.1:1]';y=(x.^2-3*x+5).*exp(-5*x).*sin(x);n=8;A=[];fori=1:n+1,A(:,i)=x.^(n+1-i);endc=A\y;vpa(poly2sym(c),5)ans=-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-125.29*x^4+74.450*x^3-27.672*x^2+4.9869*x+.42037e-611/26/202240例對f(x)=(x2-3x+5)e-5xsinx進行多項式8.2.3最小二乘曲線擬合有一組數據xi，yi（i=1,2,…,N）滿足某一函數原型，其中a為待定系數向量，求出一組待定系數的值使得目標函數最?。?1/26/2022418.2.3最小二乘曲線擬合有一組數據xi，yi（i=1,最小二乘曲線擬合函數lsqcurvefit的調用格式：[a,Jm]=lsqcurvefit(Fun,a0,x,y)Fun為原型函數的matlab表示，可以是M-函數或inline()函數a0為最優化初值x和y為原始輸入輸出數據向量a為返回的待定系數向量Jm為在待定系數下目標函數的值11/26/202242最小二乘曲線擬合函數lsqcurvefit的調用格式：11/例由下面語句生成一組數據，其中ai為待定系數a=0:0.1:10;y=0.12*exp(-0.213*x)+0.54exp(-0.17*x).*sin(1.23*x);并且該數據滿足采用最小二乘曲線擬合獲得這些待定系數，使目標函數的值為最小。編寫函數：f=inline(‘a1*exp(-a2*x)+a3exp(-a4*x).*sin(a5*x)’,’a’,’x’);11/26/202243例由下面語句生成一組數據，其中ai為待定系數11/22/求解待定系數向量[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y);xx',resans=0.11970.21250.54040.17021.2300res=7.1637e-00711/26/202244求解待定系數向量11/22/202244修改最優化選項：ff=optimset;ff.TolFun=1e-20;ff.TolX=1e-15;[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff);xx,resans=0.12000.21300.54000.17001.2300res=9.5035e-02111/26/202245修改最優化選項：11/22/202245繪制曲線：>>x1=0:0.01:10;y1=f(xx,x1);plot(x1,y1,x,y,'o')11/26/202246繪制曲線：11/22/202246例已知數據可能滿足求滿足數據的最小二乘解a、b、c和d的值輸入已知數據點：x=0.1:0.1:1;y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];x0.10.20.30.40.5y2.32012.64702.97073.28853.6008x0.60.70.80.91.0y3.90904.21474.51914.82325.127511/26/202247例已知數據可能滿足x0.10.20.30.40.5y2.編寫函數functiony=c8f3(a,x)y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);待定系數求解a=lsqcurvefit('c8f3',[1;2;2;3],x,y);a'ans=2.45752.45571.44372.072011/26/202248編寫函數11/22/202248繪制曲線：>>y1=c8f3(a,x);plot(x,y,x,y1,’o’)11/26/202249繪制曲線：11/22/202249第8章數據插值與擬合11/26/202250第8章數據插值與擬合11/22/202218.1數據插值一維數據的插值問題二維網格數據的插值問題二維一般分布數據的插值問題高維插值問題11/26/2022518.1數據插值一維數據的插值問題11/22/202228.1.1一維數據的插值問題interp1函數spline函數lagrange插值hermite插值11/26/2022528.1.1一維數據的插值問題interp1函數11/22/1、interp1函數interp1函數的調用格式：y=interp1(x0,y0,x,method)y=interp1(x0,y0,x,method,'extrap')*臨近點插值：method=‘nearest’*線性插值：method=‘linear’*三次樣條插值：method=‘spline’*立方插值：method=‘pchip’or‘cubic’11/26/2022531、interp1函數interp1函數的調用格式：*臨近點1、interp1函數選擇插值方法時主要考慮因素：運算時間、占用計算機內存和插值的光滑程度。運算時間占用計算機內存光滑程度臨近點插值快少差線性插值稍長較多稍好三次樣條插值最長較多最好立方插值較長多較好插值方法比較11/26/2022541、interp1函數選擇插值方法時主要考慮因素：運算時間、例:已知的數據點來自函數

根據生成的數據進行插值處理，得出較平滑的曲線

直接生成數據。數據點：>>x=0:0.12:1;>>y=(x.^2-3*x+5).*exp(-5*x).*sin(x);>>plot(x,y,x,y,'o')11/26/202255例:已知的數據點來自函數

根據生成的數據進行插值處理，得調用interp1()函數：>>x0=0:0.02:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);>>y1=interp1(x,y,x0);y2=interp1(x,y,x0,'cubic');>>y3=interp1(x,y,x0,'spline');y4=interp1(x,y,x0,'nearest');>>plot(x0,[y1',y2',y3',y4'],':',x,y,'o',x0,y0)誤差分析>>[max(abs(y0(1:49)-y2(1:49))),max(abs(y0-y3)),max(abs(y0-y4))]ans=0.01770.00860.15981、interp1函數11/26/202256調用interp1()函數：1、interp1函數11/2功能三次樣條數據插值格式

yy=spline(x,y,xx)

例：對離散分布在y=exp(x)sin(x)函數曲線上的數據點進行樣條插值計算：>>x=[024581212.817.219.920];>>y=exp(x).*sin(x);>>xx=0:.25:20;>>yy=spline(x,y,xx);>>plot(x,y,'o',xx,yy)2、spline函數11/26/202257功能三次樣條數據插值2、spline函數11/22/203、lagrange插值方法介紹對給定的n個插值點及對應的函數值，利用構造的n-1次lagrange插值多項式，則對插值區間內任意x的函數值y可通過下式求的：11/26/2022583、lagrange插值方法介紹11/22/20229lagrange插值函數functiony=lagrange(x0,y0,x)ii=1:length(x0);y=zeros(size(x));fori=iiij=find(ii~=i);y1=1;forj=1:length(ij),y1=y1.*(x-x0(ij(j)));endy=y+y1*y0(i)/prod(x0(i)-x0(ij));end11/26/202259lagrange插值函數functiony=lagrangLagrange插值>>x0=-1+2*[0:10]/10;y0=1./(1+25*x0.^2);>>x=-1:0.01:1;y=lagrange(x0,y0,x);>>ya=1./(1+25*x.^2);>>plot(x,ya,x,y,':')調用interp1函數>>y1=interp1(x0,y0,x,'cubic');y2=interp1(x0,y0,x,'spline');>>plot(x,ya,x,y1,':',x,y2,'--')3、lagrange插值例：對進行Lagrange插值11/26/202260Lagrange插值3、lagrange插值例：對許多實際插值問題中，為使插值函數能更好地和原來的函數重合，不但要求二者在節點上函數值相等，而且還要求相切，對應的導數值也相等，甚至要求高階導數也相等。這類插值稱作切觸插值，或埃爾米特(Hermite)插值。滿足這種要求的插值多項式就是埃爾米特插值多項式4、埃爾米特插值(hermite插值)11/26/202261許多實際插值問題中，為使插值函數能更好地和原4、埃爾米特插值(hermite插值)方法介紹已知n個插值點及對應的函數值和一階導數值。則對插值區間內任意x的函數值y的Hermite插值公式：11/26/2022624、埃爾米特插值(hermite插值)方法介紹11/22/2hermite插值函數functiony=hermite(x0,y0,y1,x)n=length(x0);m=length(x);fork=1:myy=0.0;fori=1:nh=1.0;a=0.0;forj=1:nifj~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;end11/26/202263hermite插值函數functiony=hermite(例：利用Hermite插值法求sin0.34的近似值。>>x0=[0.3,0.32,0.35];y0=[0.29552,0.31457,0.34290];>>y1=[0.95534,0.94924,0.93937];>>y=hermite(x0,y0,y1,0.34)y=0.3335>>sin(0.34)％與精確值比較ans=0.33354、hermite插值11/26/202264例：利用Hermite插值法求sin0.34的近似值。4、h8.1.2二維網格數據的插值問題二維插值函數interp2的調用格式：zi=interp2(x0,y0,z0,xi,yi)zi=interp2(x0,y0,z0,xi,yi,method)*臨近點插值：method=‘nearest’*線性插值：method=‘linear’（缺省算法）*三次樣條插值：method=‘spline’*立方插值：method=‘pchip’or‘cubic’11/26/2022658.1.2二維網格數據的插值問題二維插值函數interp2例由二元函數獲得一些較稀疏的網格數據，對整個函數曲面進行各種插值，并比較插值結果繪制已知數據的網格圖>>[x,y]=meshgrid(-3:0.6:3,-2:0.4:2);>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);>>surf(x,y,z);axis([-3,3,-2,2,-0.7,1.5])8.1.2二維網格數據的插值問題11/26/202266例由二元函數8.1.2二維網格數據的插值問題11/22/默認線性插值算法進行插值>>[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])立方和樣條插值：>>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])>>figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])11/26/202267默認線性插值算法進行插值11/22/202218算法誤差的比較>>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])11/26/202268算法誤差的比較11/22/2022198.1.3二維一般分布數據的插值問題griddata函數的調用格式：z=griddata(x0,y0,z0,x,y,method)method=‘v4’：插值算法，公認效果較好*臨近點插值：method=‘nearest’*線性插值：method=‘linear’（缺省算法）*三次樣條插值：method=‘spline’*立方插值：method=‘cubic’11/26/2022698.1.3二維一般分布數據的插值問題griddata函數的例在x為［-3，3］，y為［-2，2］矩形區域隨機選擇一組數據點，用’v4’與’cubic’插值法進行處理，并對誤差進行比較。已知數據點>>x=-3+6*rand(200,1);y=-2+4*rand(200,1);>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);>>plot(x,y,'x')%樣本點的二維分布>>figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5])，grid11/26/202270例cubic和v4算法>>[x1,y1]=meshgrid(-3:0.2:3,-2:0.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])誤差分析>>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.15])>>figure;surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15])11/26/202271cubic和v4算法11/22/2022228.1.4高維插值問題三維插值interp3函數的調用格式：三維網格[x,y,z]=meshgrid(x1,y1,z1)griddata3()三維非網格數據插值n維插值interpn函數n維網格[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]griddatan()n維非網格數據插值interp3()、interpn()調用格式同interp2()一致griddata3()、griddatan()調用格式同griddata()一致11/26/2022728.1.4高維插值問題三維插值interp3函數的調用格式例：通過函數生成一些網格型樣本點，據此進行插值并給出插值誤差。>>[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+z0.^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.041911/26/202273例：8.2數據擬合多項式擬合函數線性組合的曲線擬合方法最小二乘曲線擬合11/26/2022748.2數據擬合多項式擬合11/22/202225多項式標準形式：多項式擬合函數polyfit，其調用格式p=polyfit(x,y,n)其中x和y為原始樣本點構成的向量n為多項式的階次p為多項式系數按降冪排列得到的向量8.2.1多項式擬合11/26/202275多項式標準形式：8.2.1多項式擬合11/22/20222三次多項式擬合>>x0=0:.1:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);>>p3=polyfit(x0,y0,3);vpa(poly2sym(p3),10)ans=2.839962923*x^3-4.789842696*x^2+1.943211631*x例:已知的數據點來自函數

用多項式擬合的方法在不同階次下進行擬合11/26/202276三次多項式擬合例:已知的數據點來自函數

用多項式擬合的方繪制擬合曲線：>>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')11/26/202277繪制擬合曲線：11/22/202228不同次數進行擬合：>>p4=polyfit(x0,y0,4);y2=polyval(p4,x);>>p5=polyfit(x0,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)11/26/202278不同次數進行擬合：11/22/202229擬合次數為8的多項式：>>vpa(poly2sym(p8),5)ans=-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-125.29*x^4+74.450*x^3-27.672*x^2+4.9869*x+.42037e-6Taylor冪級數展開：>>symsx;y=(x^2-3*x+5)*exp(-5*x)*sin(x);>>vpa(taylor(y,9),5)ans=5.*x-28.*x^2+77.667*x^3-142.*x^4+192.17*x^5-204.96*x^6+179.13*x^7-131.67*x^8多項式表示數據模型是不唯一的，即是兩個多項式函數完全不同。在某一區域內其曲線可能特別近似。11/26/202279擬合次數為8的多項式：11/22/202230編制程序>>x0=-1+2*[0:10]/10;y0=1./(1+25*x0.^2);>>x=-1:.01:1;ya=1./(1+25*x.^2);>>p3=polyfit(x0,y0,3);y1=polyval(p3,x);>>p5=polyfit(x0,y0,5);y2=polyval(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,':')例：對多項式擬合，多項式擬合的效果并不一定總是很精確的。11/26/202280編制程序例：對用Taylor冪級數展開效果將更差。>>symsx;y=1/(1+25*x^2);p=taylor(y,x,10)p=1-25*x^2+625*x^4-15625*x^6+390625*x^8多項式擬合效果>>x1=-1:0.01:1;ya=1./(1+25*x1.^2);>>y1=subs(p,x,x1);>>plot(x1,ya,'--‘,x1,y1)11/26/202281用Taylor冪級數展開效果將更差。11/22/2022328.2.2函數線性組合的曲線擬合方法已知某函數的線性組合為：g(x)=c1f1(x)+c2f2(x)+…+cnfn(x)其中f1(x),f2(x),…,cnfn(x)為已知函數，c1,c2,…,cn為待定系數假設已經測出數據(x1,y1),(x2,y2),…,(xm,ym)，則可建立如下的線性方程：Ac=y11/26/2022828.2.2函數線性組合的曲線擬合方法已知某函數的線性組合為11/26/20228311/22/202234例假設測出一組實驗數據(xi,yi)，其函數原型為y(x)=c1+c2e-3x+c3cos(-2x)e-4x+c4x2，用已知數據求出待定系數ci的值x00.20.40.70.90.92y2.882.25761.96831.92582.08622.109x0.991.21.41.481.5y2.19792.54092.96273.1553.205211/26/202284例假設測出一組實驗數據(xi,yi)，其函數原型為y(x擬合ci參數>>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=A\y;c1=c'c1=1.22002.3397-0.67970.870011/26/202285擬合ci參數11/22/202236圖形顯示>>x0=[0:0.01:1.5]';>>A1=[ones(size(x0)),exp(-3*x0),cos(-2*x0).*exp(-4*x0)x0.^2];>>y1=A1*c;>>plot(x0,y1,x,y,'x')11/26/202286圖形顯示11/22/202237例對下列實驗數據進行函數擬合已知數據點>>x=[1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,2.2255,2.4596,2.7183,3.6693];>>y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,0.2864,0.2532,0.2238,0.1546];>>plot(x,y,x,y,'*')x1.10521.22141.34991.49181.64871.8221y0.67950.60060.53090.46930.41480.3666x2.01382.22552.45962.71833.6693y0.324140.28650.25320.22380.154611/26/20