1、曲面擬合原理與實例精品文檔問題:給定一組坐標(biāo)(xg,yg,zg) , g1,2,n ,表示有n個點。要求用以下二元0收集于網(wǎng)絡(luò)，如有侵權(quán)請聯(lián)系管理員刪除多項式函數(shù)對所給的坐標(biāo)進(jìn)行擬合:aijXi 1yjp,qi 1 j 1f (x, y)aijXyij 1,1f (x, y)a11a12y2a13yLq 1%ya21xa22xy2a23XyLq 1a2qXyMi 1i 1i 12Li 1 q 1ai1xa2x yai3X yaiqxyqMP 1P 1P 12LP 1 qap1xap2xyap3x yapqX y,y,aa11a2La1qa21a22La2qMMOMap1ap2Lapq則函數(shù)又可

2、表示為f(x,y)XTAy ,擬合的目標(biāo)就是求出系數(shù)矩陣 A最小二乘法:構(gòu)造關(guān)于系數(shù)aij的多元函數(shù):P qZi 1 j 1g(aij X yi 1 j 1nZg)22s( ai1 ,L , aPq)gf(Xg，yg)Zg g 1點(a11 ,，apq)是多元函數(shù)s(a11,L ,apq)的極小點，其中g(shù)為權(quán)函數(shù)，默認(rèn)為1,所以點(a11 ,，apq)必須滿足方程組Saij在g 1的情況下，有精品文檔收集于網(wǎng)絡(luò)，如有侵權(quán)請聯(lián)系管理員刪除aijaijn2f (Xg,yg) Zgg 12f(Xg,yg)ZgFf(Xg切2f (Sg)i 1 j 1Zg Xg ygg 1n2g 1因此可得XgIyg

3、If(Xg,yg)i 1 j 1Xg yg Zgp,q1,1U (i,j)p,qa U1,1上式實際共有1f(Xg,yg)i 1Xg1yg 1Zgp,qU11(1,1)MU11(p,q)a1,1Xg1ygni 1 j 1Xg yg Zgg 11Xgygi 1 j 1Xg yg Zg1(Xg 1yg1(i,j)i 1Xgygj1)ng1ygg 11ZgIiIjg Xg yg )，v(i,j)ni 1 jXgygZgg 1v(i, j)(i,j)(1,1),(P,q)P q個等式，可將這Upq(1,1)MUpq(p,q)a11Mapqq個等式寫成矩陣的形式有:v(1,1)Mv(p,q)也就是U*a

4、=V的形式，其中Un(1,1)UMU11( P,q)UPq (1,1)MUpq(p,q)a11Mapqv(1,1)Mv(p,q)U為Pq Pq階矩陣，實現(xiàn)函數(shù)為funCtiOn A=IeftmatriX（X,p,y,q） ； V為長Pq的列 向量，實現(xiàn)函數(shù)為function B=rightmatrix（x,p,y,q,z）。這樣就可以算出列矩陣 a, 然后轉(zhuǎn)化成A。例子：某地區(qū)有一煤礦，為估計其儲量以便于開采，先在該地區(qū)進(jìn)行勘探。假設(shè)該地區(qū)是一長方形區(qū)域，長為 4公里，寬為5公里。經(jīng)勘探得到如下數(shù)據(jù): 煤礦勘探數(shù)據(jù)表編號12345678910橫坐標(biāo)（公里）1111122222縱坐標(biāo)（公里）12

5、34512345煤層厚度（米）13.7225.808.4725.2722.3215.4721.3314.4924.8326.19編號11121314151617181920橫坐標(biāo)（公里）3333344444縱坐標(biāo)（公里）1234512345煤層厚度（米）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29請你估計出此地區(qū)內(nèi)（2x4,1 y 5 ）煤的儲量，單位用立方米表示,并用電腦畫出該煤礦的三維圖象。如果直接畫出三維曲面圖形：CIear;x=1:4;y=1:5;X,Y=meshgrid(x,y)Z=13.7225.808.4725.2722.3

6、2;15.4721.3314.4924.8326.19;23.2826.4829.1412.0414.58;19.9523.7315.3518.0116.29'SUrf(X ,Y,Z);X =12341234123412341234Y =11112222333344445555Z =13.7200 15.4700 23.2800 19.950025.800021.330026.480023.73008.470014.490029.140015.350025.270024.830012.040018.010022.320026.190014.580016.290030粗略計算體積：底面積

7、乘以平均高度。P=SUm(Z);q=P(：,2,3,4);h=sum(q')15v=2000*4000*hh =20.0773V =1.6062e+008進(jìn)行線性插值：xi=li nspace(1,4,31)yi=li nspace(1,5,41);XI,YI=meshgrid(xi,yi);ZI=i nterp2(X,Y,Z,XI,YI,'li near');SUrf(XI,YI,ZI);30. 1 1進(jìn)行三次多項式插值：41);XI,YI=meshgrid(xi,yi);xi=li nspace(1,4,31)yi=li nspace(1,5,Zl=interp2(

8、X,Y,Z,XI,YI,'CUbiC ');SUrf(Xl,YI,ZI);30、進(jìn)行插值后計算體積：底面積乘以平均高度xi=li nspace(1,4,61)yi=li nspace(1,5,81);XI,YI=meshgrid(xi,yi);Zl=i nterp2(X,Y,Z,XI,YI,'cubic');SUrf(XI,YI,ZI);H=O; n=0;for j=21:61for i=1:81H=H+ZI(i,j); n=n+1;endendnH=H/nS=2000*4000;V=S*H n =3321H =20.8222V =1.6658e+00830.上

9、面是插值的方法解題，下面用擬合的方法解題。 為此編寫了幾個M函數(shù)：Ieftmatrix.mfun Cti OnU=IeftmatriX(X,p,y,q)% U*a=V a為系數(shù)列矩陣，長度為p*q% U為左邊p*q乘p*q矩陣% x,y為長度一致的列矩陣，給定點的坐標(biāo)% p,q 為擬合的函數(shù)中x,y的幕的最高次數(shù)m=le ngth(x);if (nargin=4) & (m=le ngth(y)error( 'error CheCk check!');endU_len gth=p*q;U=ZerOS(U_le ngth,U_le ngth);for i=1 : p*qf

10、or j= 1 : p*qx_z=quotie nt(j-1,q)+quotie nt(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ù),% y 的幕的次quotie nt為rightmatrix.mfun Cti OnV=rightmatrix(x,p,y,q,z)% U*a=V% V為一個列向量長為p*q% X y Z 為點的坐標(biāo)%p q分別為Xy幕的最高次數(shù) if nargin=5error( 'error CheCk check!

11、 rightmatrix')endV=ZerOS(P*q,1);for i=1 : p*qx_z=quotie nt(i-1,q); y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z); endquuotie nt.mfun Cti onsh=quotie nt(x,y) % Sh為XZy 的商Sh=(X-mod(x,y)/y;qiuhe.mfun Cti onhe=qiuhe(x,p,y,q,z)% he xp*yq從 1 >m 的和% x,y 向量長度相同% p,q分別為x,y的幕的次數(shù)m=le ngth(x);if (nargin<4

12、 )&(m =le ngth(y)度必需一樣error( 'error CheCk check!' end%輸入量至少為四，x,y行向量長);if nargin=4的向量z=on es(m,1); end%沒有Z ,默認(rèn)為元素全部為1he=0;for i=1:mhe=he+x(i)p * y(i)q*z(i);end% 1->m 求和F面一段程序先進(jìn)行擬合，然后驗證擬合的效果，具體操作:先輸入 x=y=z= p=q=(注意 x,y,z 是向量);擬合得到系數(shù)a,也就是得到了擬合的函數(shù)；根據(jù)擬合函數(shù)計算給定點(xx, yy)的函數(shù)值zz=f(xx, yy)并進(jìn)行畫圖

13、檢驗 程序保存于M文件fit.m。fit.mclear;X,Y=meshgrid(1:4,1:5);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'x=reshape(X,20,1);y=reshape(Y,20,1);z=reshape(Z,20,1);p=4;q=5;U=IeftmatriX(X,p,y,q);% U*a_ n=VV=rightmatrix(x,p,y,q,z);%a_n=i nv(U)*V;a_n=UV;

14、把長為p*q 的列向量a_n轉(zhuǎn)換成p*q的矩陣aa 求商for i=1 : Ien gth(a_ n)%ii=quotient(i-1,q)+1; % quotient jj=mod(i-1,q)+1;aa(ii,jj)=a_ n( i,1);end aa m=31; n=41;%m=4; n=5;XI,YI=meshgrid(I in SPaCe(1,4,m),li nspace(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)*xi*yj, (i=1.p,j=1.q)為PXq的系數(shù)的矩陣%zz=0;% ZZfor i=1 : P%for j=1 : q % aaxt=xx.(i-1);yt=yy.(j-i);xy=xt.*yt; zz=zz+aa(i,j).*xy;endendZl=reshape(z z,n ,m);SUrf(Xl,YI,ZI); %axis(1 4 1 5 0 30)