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

1、曲面擬合原理與實(shí)例精品文檔問題:給定一組坐標(biāo)（,yg,zg） , g1,2,n ,表示有n個(gè)點(diǎn)。要求用以下二元收集于網(wǎng)絡(luò)，如有侵權(quán)請(qǐng)聯(lián)系管理員刪除多項(xiàng)式函數(shù)對(duì)所給的坐標(biāo)進(jìn)行擬合:f (x, y)f(x, y)p,qaijx11,1an a21xM i 1 ai1xM p 1 apxx2 xMxp,yy2 yMyq則函數(shù)又可表示為i 1 j 1ajxya2ya22xy2a13y2 a23xyq 1aqyq 1a2qxyi 1%x yi 1 2ai3x yi 1 q 1aiqx yp 1ap2xyap3xp 1 q 1apqx y,Aal2Laqa22La2qMOMap2Lapqa21Map1f

2、(x, y)xTAy,擬合的目標(biāo)就是求出系數(shù)矩陣Ao最小二乘法:構(gòu)造關(guān)于系數(shù)aj的多元函數(shù):p q/i 1 j 1g( ajx yi 1 j 1nZg)22s(On,L ,apq) g f (xg, yg) Zgg 1點(diǎn)（a1 ,，apq）是多元函數(shù)s（a11,L , apq）的極小點(diǎn)，其中 g為權(quán)函數(shù)，默認(rèn)為1,所以點(diǎn)（an ,apq）必須滿足方程組saij在g 1的情況下，有aijajn2f(xg,yg) zgg 12f(xg,yg)zg f(xg,yg) aij2f (%,yg)i 1 j 1 zgxg ygg 1n2g 1因此可得xg1yg 1f(xg,yg)i 1 j 1xg yg

3、zgnxg1ygg 11Mxg,yg)i 1 xg1p,qp,q1,1令u (i, j)p,qa u1,1上式實(shí)際共有U11(1,1)Mun(p,q)a1,1xg1ygni 1 j 1xg yg zgg 11xgygi 1 j 1xg yg zg(xg 1yg1(i,j)yj)nxg1ygg 11zg1 i 1 j 1、 g xg yg )v(i, j)v(i, j)(i,j)p q個(gè)等式，可將這Upq(1,1)MUpq(p,q)anMapq也就是U*a=V的形式，其中Uii(1,1)U M uii( p,q)Upq (1,1)MUpq(P，q)ni 1 j 1xg yg zgg 1(1,1)

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

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

6、721.3314.4924.8326.19;23.2826.4829.1412.0414.58;19.9523.7315.3518.0116.29surf(X,Y,Z);X =1 2342 2343 2344 2345 234Y =6 1117 2228 3339 44410 555Z =13.7200 15.4700 23.2800 19.950025.8000 21.3300 26.4800 23.73008.4700 14.4900 29.1400 15.350025.2700 24.8300 12.0400 18.010022.3200 26.1900 14.5800 16.2900粗

7、略計(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);XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,linear);surf(XI,YI,ZI);進(jìn)行三次多項(xiàng)式插值：41);XI,YI=meshgrid(xi,yi);xi=linspace(1,4, 31);yi=linspace(1,5,ZI=interp2(X,Y,Z,XI,YI,cub

8、ic );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,Z,XI,YI,cubic);surf(XI,YI,ZI);H=0;n=0;for j=21:61 fo門=1:81 H=H+ZI(i,j);n=n+1;end end n H=H/n S=2000*4000;V=S*H n =3321H =20.8222V =1.6658e+00830、15.1Q上面是插值的方法解題，下面用擬合的方法解題。 為此編寫了幾個(gè)M函數(shù)

9、：leftmatrix.mfunction U=leftmatrix(x,p,y,q)% U*a=V a為系數(shù)列矩陣，長(zhǎng)度為p*q% U為左邊p*q乘p*q矩陣% x,y為長(zhǎng)度一致的列矩陣，給定點(diǎn)的坐標(biāo)% p,q 為擬合的函數(shù)中x,y的幕的最高次數(shù)m=length(x);if (nargin=4) & (m=length(y)error( error check check!);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

10、_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 的事的次rightmatrix.m function V=rightmatrix(x,p,y,q,z) % U*a=V% V為一個(gè)列向量長(zhǎng)為p*q% x y z為點(diǎn)的坐標(biāo)%p q分別為xy幕的最高次數(shù) if nargin=5error( error check check! rightmatrix endV=zeros(p*q,1);for i=1 : p*qx_z=quotient(i-1

11、,q);y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z); endquuotient.m function sh=quotient(x,y)% sh為x/y 的商sh=(x-mod(x,y)/y;qiuhe.mfunction he=qiuhe(x,p,y,q,z)% he xAp*yAq從 1 m 的和% x,y 向量長(zhǎng)度相同% p,q分別為x,y的幕的次數(shù) m=length(x);%輸入量至少為四，x,y行向量長(zhǎng)）；%沒有z ,默認(rèn)為元素全部為1% 1-m 求和if (nargin4 )&(m=length(y)度必需一樣error( error che

12、ck check! endif nargin=4的向量z=ones(m,1);endhe=0;for i=1:mhe=he+x(i)Ap * y(i)Aq*z(i);end卜面一段程序先進(jìn)行擬合，然后驗(yàn)證擬合的效果，具體操作:先輸入 x= - y= - z= - p= - q=(注意 x,y,z是向量)；擬合得到系數(shù)a,也就是得到了擬合的函數(shù)；根據(jù)擬合函數(shù)計(jì)算給定點(diǎn)(xx, yy)的函數(shù)值z(mì)z=f(xx, yy)并進(jìn)行畫圖檢驗(yàn) 程序保存于M文件fit.m。fit.mclear;X,Y=meshgrid(1:4,1:5);Z=13.72 25.80 8.4715.47 21.33 14.4923

13、.28 26.48 29.1425.27 22.32;24.83 26.19;12.04 14.58;19.95 23.73 15.3518.01 16.29;x=reshape(X,20,1);y=reshape(Y,20,1);z=reshape(Z,20,1);p=4；q=5;U=leftmatrix(x,p,y,q); % U*a_n=VV=rightmatrix(x,p,y,q,z);%a_n=inv(U)*V;a_n=UV;for i=1 : length(a_n) %把長(zhǎng)為p*q的列向量a_n轉(zhuǎn)換成p*q的矩陣aaii=quotient(i-1,q)+1; % quotient求

14、商jj=mod(i-1,q)+1;aa(ii,jj)=a_n(i,1)；endaam=31;n=41;%m=4;n=5;XI,YI=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;% zzfor i=1 : p %for j=1 : q % aa xt=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(1 4 1 5 0 30) aa =1.