曲面擬合原理與實(shí)例說課材料

1、曲面擬合原理與實(shí)例精品文檔問題:給定一組坐標(biāo)（,yg,zg） , g1,2,n ,表示有n個(gè)點(diǎn)。要求用以下二元收集于網(wǎng)絡(luò)，如有侵權(quán)請聯(lián)系管理員刪除多項(xiàng)式函數(shù)對所給的坐標(biāo)進(jìn)行擬合:f(x,y)f(x,y)p,qaijx11,1ana21xMi1ai1xMp1apxx2xMxp,yy2yMyq則函數(shù)又可表示為i1j1ajxya2ya22xy2a13y2a23xyq1aqyq1a2qxyi1%xyi12ai3xyi1q1aiqxyp1ap2xyap3xp1q1apqxy,Aal2Laqa22La2qMOMap2Lapqa21Map1f(x,y)xTAy,擬合的目標(biāo)就是求出系數(shù)矩陣Ao最小二乘法:構(gòu)造

2、關(guān)于系數(shù)aj的多元函數(shù):p q/i 1 j 1g( ajx yi 1 j 1nZg)22s(On,L,apq)gf(xg,yg)Zgg1點(diǎn)（a1,，apq）是多元函數(shù)s（a11,L,apq）的極小點(diǎn)，其中g(shù)為權(quán)函數(shù)，默認(rèn)為1,所以點(diǎn)（an,apq）必須滿足方程組saij在g1的情況下，有aijajn2f(xg,yg)zgg12f(xg,yg)zgf(xg,yg)aij2f(%,yg)i1j1zgxgygg1n2g1因此可得xg1yg1f(xg,yg)i1j1xgygzgnxg1ygg11Mxg,yg)i1xg1p,qp,q1,1令u(i,j)p,qau1,1上式實(shí)際共有U11(1,1)Mun(

3、p,q)a1,1xg1ygni1j1xgygzgg11xgygi1j1xgygzg(xg1yg1(i,j)yj)nxg1ygg11zg1i1j1、gxgyg)v(i,j)v(i,j)(i,j)pq個(gè)等式，可將這Upq(1,1)MUpq(p,q)anMapq也就是U*a=V的形式，其中Uii(1,1)UMuii(p,q)Upq(1,1)MUpq(P，q)ni1j1xgygzgg1(1,1)，,(p,q)q個(gè)等式寫成矩陣的形式有:v(1,1)Mv(p,q)a11Mapqv(1,1)Mv(p,q)U為pqpq階矩陣，實(shí)現(xiàn)函數(shù)為functionA=leftmatrix（x,p,y,q）；V為長pq的列

4、向量，實(shí)現(xiàn)函數(shù)為functionB=rightmatrix（x,p,y,q,z）。這樣就可以算出列矩陣a,然后轉(zhuǎn)化成A。例子：某地區(qū)有一煤礦，為估計(jì)其儲量以便于開采，先在該地區(qū)進(jìn)行勘探。假設(shè)該地區(qū)是一長方形區(qū)域，長為4公里，寬為5公里。經(jīng)勘探得到如下數(shù)據(jù):煤礦勘探數(shù)據(jù)表編號12345678910橫坐標(biāo)（公里）1111122222縱坐標(biāo)（公里）1234512345煤層厚度（米）13.7225.808.4725.2722.3215.4721.3314.4924.8326.19編號11121314151617181920橫坐標(biāo)（公里）3333344444縱坐標(biāo)（公里）1234512345煤層厚度（米

5、）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29請你估計(jì)出此地區(qū)內(nèi)（2x4,1y5）煤的儲量，單位用立方米表示,并用電腦畫出該煤礦的三維圖象。如果直接畫出三維曲面圖形：clear;x=1:4;y=1:5;X,Y=meshgrid(x,y)Z=13.7225.808.4725.2722.32;15.4721.3314.4924.8326.19;23.2826.4829.1412.0414.58;19.9523.7315.3518.0116.29'surf(X,Y,Z);X=1 2342 2343 2344 2345 234Y=6 11

6、17 2228 3339 44410 555Z=13.720015.470023.280019.950025.800021.330026.480023.73008.470014.490029.140015.350025.270024.830012.040018.010022.320026.190014.580016.2900粗略計(jì)算體積：底面積乘以平均高度。p=sum(Z);q=P(：,2,3,4);h=sum(q')/15v=2000*4000*hh=20.0773v=1.6062e+008進(jìn)行線性插值：xi=linspace(1,4,31);yi=linspace(1,5,41);

7、XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,'linear');surf(XI,YI,ZI);進(jìn)行三次多項(xiàng)式插值：41);XI,丫I=meshgrid(xi,yi);xi=linspace(1,4,31);yi=linspace(1,5,ZI=interp2(X,Y,Z,XI,YI,'cubic');surf(XI,YI,ZI);進(jìn)行插值后計(jì)算體積：底面積乘以平均高度xi=linspace(1,4,61);yi=linspace(1,5,81);XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y

8、,Z,XI,YI,'cubic');surf(XI,YI,ZI);H=0;n=0;forj=21:61fo門=1:81H=H+ZI(i,j);n=n+1;endendnH=H/nS=2000*4000;V=S*Hn=3321H=20.8222V=1.6658e+00830、15.1Q上面是插值的方法解題，下面用擬合的方法解題。為此編寫了幾個(gè)M函數(shù)：leftmatrix.mfunctionU=leftmatrix(x,p,y,q)%U*a=Va為系數(shù)列矩陣，長度為p*q%U為左邊p*q乘p*q矩陣%x,y為長度一致的列矩陣，給定點(diǎn)的坐標(biāo)%p,q為擬合的函數(shù)中x,y的幕的最高次數(shù)m

9、=length(x);if(nargin=4)&(m=length(y)error('errorcheckcheck!');endU_length=p*q;U=zeros(U_length,U_length);for i=1 : p*qfor j= 1 : p*qx_z=quotient(j-1,q)+quotient(i-1,q);求商y_z=mod(j-1,q)+mod(i-1,q);數(shù)U(i,j尸qiuhe(x,x_z,y,y_z);endend% U為p*q階方陣%賦值0,目的是分配內(nèi)存% x 的事的次數(shù)，quotient 為% y 的事的次rightmatri

10、x.mfunctionV=rightmatrix(x,p,y,q,z)%U*a=V%V為一個(gè)列向量長為p*q%xyz為點(diǎn)的坐標(biāo)%pq分別為xy幕的最高次數(shù)ifnargin=5error('errorcheckcheck!rightmatrixendV=zeros(p*q,1);fori=1:p*qx_z=quotient(i-1,q);y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z);endquuotient.mfunctionsh=quotient(x,y)%sh為x/y的商sh=(x-mod(x,y)/y;qiuhe.mfunctionhe=qi

11、uhe(x,p,y,q,z)%hexAp*yAq從1>m的和%x,y向量長度相同%p,q分別為x,y的幕的次數(shù)m=length(x);%輸入量至少為四，x,y行向量長）；%沒有z ,默認(rèn)為元素全部為1% 1->m 求和if(nargin<4)&(m=length(y)度必需一樣error('errorcheckcheck!endifnargin=4的向量z=ones(m,1);endhe=0;fori=1:mhe=he+x(i)Ap*y(i)Aq*z(i);end卜面一段程序先進(jìn)行擬合，然后驗(yàn)證擬合的效果，具體操作:先輸入x=-y=-z=-p=-q=(注意x,

12、y,z是向量)；擬合得到系數(shù)a,也就是得到了擬合的函數(shù)；根據(jù)擬合函數(shù)計(jì)算給定點(diǎn)(xx,yy)的函數(shù)值zz=f(xx,yy)并進(jìn)行畫圖檢驗(yàn)程序保存于M文件fit.m。fit.mclear;X,Y=meshgrid(1:4,1:5);Z=13.7225.808.4715.4721.3314.4923.2826.4829.1425.2722.32;24.8326.19;12.0414.58;19.9523.7315.3518.0116.29'x=reshape(X,20,1);y=reshape(Y,20,1);z=reshape(Z,20,1);p=4；q=5;U=leftmatrix(x

13、,p,y,q);%U*a_n=VV=rightmatrix(x,p,y,q,z);%a_n=inv(U)*V;a_n=U'V;fori=1:length(a_n)%把長為p*q的列向量a_n轉(zhuǎn)換成p*q的矩陣aaii=quotient(i-1,q)+1;%quotient求商jj=mod(i-1,q)+1;aa(ii,jj)=a_n(i,1)；endaam=31;n=41;%m=4;n=5;XI,丫I=meshgrid(linspace(1,4,m),linspace(1,5,n);xx=reshape(XI,m*n,1);yy=reshape(YI,m*n,1);zz=zeros(m*n,1);xy=zeros(m*n,1);xt=zeros(m*n,1);yt=zeros(m*n,1);是xx,yy 代入所擬合的函數(shù) 求出的函數(shù) 值函數(shù)為因a(i,j)*xAi*yAj, (i=1p,j=1q)為pxq的系數(shù)的矩陣%zz=0;%zzfori=1:p%forj=1:q%aaxt=xx.A(i-1);yt=yy.A(j-i);xy=xt.*yt;zz=zz+aa(i,j).*xy;endendZI=reshape(zz,n,m);surf(XI,YI,ZI);%axis(1415030)aa=1.0e+00