中插值與擬合

1、計(jì)算可視化1 插值與數(shù)據(jù)擬合1.1 一維數(shù)據(jù)的插值問題1.1.1一維插值問題的求解求解方法一維插值interp1()函數(shù)的調(diào)用格式為：說明函數(shù)根據(jù)x,y的值，計(jì)算函數(shù)在x1處的值。x,y是兩個(gè)等長(zhǎng)的已知向量，分別描述采樣點(diǎn)和樣本值，x1是一個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)，y1是一個(gè)與x1等長(zhǎng)的插值結(jié)果。 插值方法linear，nearest，cubic，spline 注：X1的取值范圍不能超出X的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。例1-1：已知的數(shù)據(jù)點(diǎn)來自函數(shù)，根據(jù)生成的數(shù)據(jù)進(jìn)行插值處理，得出較平滑的曲線直接生成數(shù)據(jù)。調(diào)用函數(shù)：Step1：x=0:.12:1;y=(x.2-3*x+5).*e

2、xp(-5*x).*sin(x);plot(x,y,x,y,'o')% 繪制原函數(shù)的圖形%'linear'x1=0:.02:1; y0=(x1.2-3*x1+5).*exp(-5*x1).*sin(x1);y1=interp1(x,y,x1); %默認(rèn)為'linear'，在樣本點(diǎn)上斜率變化很大plot(x1, y1',':',x,y,'o',x1,y0)%'cubic'y2=interp1(x,y,x1, 'cubic'); %最占內(nèi)存，輸出結(jié)果與spline差不多plot(

3、x1,y2',':',x,y,'o',x1,y0)%'spline'y3=interp1(x,y,x1,'spline'); %最花時(shí)間，但輸出結(jié)果也最平滑plot(x1,y3',':',x,y,'o',x1,y0)% 'nearest'y4=interp1(x,y,x1, 'nearest');%執(zhí)行速度最快，輸出結(jié)果為直角轉(zhuǎn)折plot(x1,y4',':',x,y,'o',x1,y0)%'linear

4、', 'cubic', 'spline', 'nearest'plot(x1,y1',y2',y3',y4',':',x,y,'o',x1,y0);Step2：max(abs(y0(1:49)-y2(1:49),max(abs(y0-y3),max(abs(y0-y4)ans = 0.0177 0.0086 0.1598例1-2：編寫一段程序，允許利用插值方法手工繪制一條光滑的曲線。function sketcher(vis) x=; y=; i=1; h=; axis(0,

5、1 0 1)while 1 x0,y0,but=ginput(1);if but=1, x=x,x0; y=y,y0; h(i)=line(x0,y0); set(h(i),'Marker','o'); i=i+1;else, break; endendif nargin=1, delete(h); end% 若需要，可以刪除樣本點(diǎn)標(biāo)識(shí) xx=x(1):(x(end)-x(1)/100: x(end); yy=interp1(x,y,xx,'spline'); line(xx,yy)1.1.2Lagrange插值算法及應(yīng)用求解方法：已知xi，yi

6、點(diǎn)，可求出x向量上各點(diǎn)處的插值為：function y=lagrange(x0,y0,x)i1=1:length(x0);y=zeros(size(x);for i=ii ij=find(ii=i); y1=1; for j=1:length(ij),y1=y1.*(x-x0(ij(j); end y=y+y1*y0(i)/prod(x0(i)-x0(ij);end例1-3：對(duì)進(jìn)行Lagrange插值。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)

7、; plot(x,ya,x,y,':')y1=interp1(x0,y0,x,'cubic'); y2=interp1(x0,y0,x,'spline');plot(x,ya,x,y1,':',x,y2,'-')1.2已知樣本點(diǎn)的定積分計(jì)算求解方法編寫函數(shù)：function y=quadspln(x0,y0,a,b)f=inline('interp1(x0,y0,x,"spline")',.'x','x0','y0');y=quad

8、l(f,a,b,1e-8,x0,y0);%先用樣條函數(shù)進(jìn)行插值后再積分其中為樣本點(diǎn)構(gòu)成的橫縱坐標(biāo)向量，為積分區(qū)間例1-4：利用樣條插值算法求解。求解：x0=0:pi/30:pi; y0=sin(x0);I=quadspln(x0,y0,0,pi)結(jié)果：I= 2.0000比較梯形法和插值法：梯形法：x0=0:pi/10:pi;y0=sin(x0);I1=trapz(x0,y0)結(jié)果：I1= 1.9835插值法：I=quadspln(x0,y0,0,pi)結(jié)果：I= 2.0000已知5個(gè)不均勻分布的樣本點(diǎn)：x0=0,0.4,1 2,pi; y0=sin(x0); % 生成樣本點(diǎn)plot(x0,y0

9、,x0,y0,'o') % 繪制出的樣本點(diǎn)折線I=quadspln(x0,y0,0,pi)I=2.0191，大約有1%的相對(duì)誤差，應(yīng)該說是相當(dāng)精確I1=trapz(x0,y0)I1=1.8416，用trapz() 函數(shù)將得出很大的誤差(7.9%)樣條插值的結(jié)果與理論之間的比較:x=0:0.01:pi,pi;y0a=sin(x); y=interp1(x0,y0,x,'spline');plot(x0,y0,x,y,':',x,y0a,x0,y0,'o')例1-5：已知其中的150個(gè)數(shù)據(jù)點(diǎn)，用quadspln()計(jì)算出該定積分的值。

10、求數(shù)值解：x=0:3*pi/2/200:3*pi/2; y=cos(15*x); I=quadspln(x,y,0,3*pi/2)結(jié)果：I= 0.0667繪制曲線：x0=0:3*pi/2/1000:3*pi/2; y0=cos(15*x0);y1=interp1(x,y,x0,'spline'); plot(x,y,x0,y1,':')1.3二維網(wǎng)格數(shù)據(jù)的插值問題求解方法二維插值interp2()函數(shù)的調(diào)用格式為：說明x0,y0,z0為已知數(shù)據(jù)，x1,y1為插值點(diǎn)構(gòu)成的新的網(wǎng)格參數(shù)，z1矩陣為在所選插值網(wǎng)格點(diǎn)出的函數(shù)的近似值 插值方法linear，nearest

11、，cubic，spline 注：X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。例1-6：由可計(jì)算出一些較稀疏的網(wǎng)格數(shù)據(jù)，對(duì)整個(gè)函數(shù)曲面進(jìn)行各種插值擬合，并比較擬合結(jié)果。繪制已知數(shù)據(jù)的網(wǎng)格圖:x,y=meshgrid(-3:.6:3,-2:.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)默認(rèn)插值算法進(jìn)行插值: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)立方插值z(mì)1=interp2(x,y,z,x1,y1,'cubic'); surf(x1,y1,z1), axis(-3,3,-2,2,-0.7,1.5)樣條插值：z2=interp2(x,y,z,x1,y1,'spline');surf(x1,y1,z2), axis(-3,3,-2,2,