曲面拟合原理与实例说课材料

1、曲面拟合原理与实例精品文档问题:给定一组坐标（,yg,zg） , g1,2,n ,表示有n个点。要求用以下二元收集于网络，如有侵权请联系管理员删除多项式函数对所给的坐标进行拟合:f(x,y)f(x,y)p,qaijx11,1ana21xMi1ai1xMp1apxx2xMxp,yy2yMyq则函数又可表示为i1j1ajxya2ya22xy2a13y2a23xyq1aqyq1a2qxyi1%xyi12ai3xyi1q1aiqxyp1ap2xyap3xp1q1apqxy,Aal2Laqa22La2qMOMap2Lapqa21Map1f(x,y)xTAy,拟合的目标就是求出系数矩阵Ao最小二乘法:构造

2、关于系数aj的多元函数:p q/i 1 j 1g( ajx yi 1 j 1nZg)22s(On,L,apq)gf(xg,yg)Zgg1点（a1,，apq）是多元函数s（a11,L,apq）的极小点，其中g为权函数，默认为1,所以点（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上式实际共有U11(1,1)Mun(

3、p,q)a1,1xg1ygni1j1xgygzgg11xgygi1j1xgygzg(xg1yg1(i,j)yj)nxg1ygg11zg1i1j1、gxgyg)v(i,j)v(i,j)(i,j)pq个等式，可将这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个等式写成矩阵的形式有:v(1,1)Mv(p,q)a11Mapqv(1,1)Mv(p,q)U为pqpq阶矩阵，实现函数为functionA=leftmatrix（x,p,y,q）；V为长pq的列

4、向量，实现函数为functionB=rightmatrix（x,p,y,q,z）。这样就可以算出列矩阵a,然后转化成A。例子：某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。假设该地区是一长方形区域，长为4公里，宽为5公里。经勘探得到如下数据:煤矿勘探数据表编号12345678910横坐标（公里）1111122222纵坐标（公里）1234512345煤层厚度（米）13.7225.808.4725.2722.3215.4721.3314.4924.8326.19编号11121314151617181920横坐标（公里）3333344444纵坐标（公里）1234512345煤层厚度（米

5、）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29请你估计出此地区内（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粗略计算体积：底面积乘以平均高度。p=sum(Z);q=P(：,2,3,4);h=sum(q')/15v=2000*4000*hh=20.0773v=1.6062e+008进行线性插值：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);进行三次多项式插值：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);进行插值后计算体积：底面积乘以平均高度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上面是插值的方法解题，下面用拟合的方法解题。为此编写了几个M函数：leftmatrix.mfunctionU=leftmatrix(x,p,y,q)%U*a=Va为系数列矩阵，长度为p*q%U为左边p*q乘p*q矩阵%x,y为长度一致的列矩阵，给定点的坐标%p,q为拟合的函数中x,y的幕的最高次数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);数U(i,j尸qiuhe(x,x_z,y,y_z);endend% U为p*q阶方阵%赋值0,目的是分配内存% x 的事的次数，quotient 为% y 的事的次rightmatri

10、x.mfunctionV=rightmatrix(x,p,y,q,z)%U*a=V%V为一个列向量长为p*q%xyz为点的坐标%pq分别为xy幕的最高次数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的幕的次数m=length(x);%输入量至少为四，x,y行向量长）；%没有z ,默认为元素全部为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卜面一段程序先进行拟合，然后验证拟合的效果，具体操作:先输入x=-y=-z=-p=-q=(注意x,

12、y,z是向量)；拟合得到系数a,也就是得到了拟合的函数；根据拟合函数计算给定点(xx,yy)的函数值zz=f(xx,yy)并进行画图检验程序保存于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转换成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 代入所拟合的函数 求出的函数 值函数为因a(i,j)*xAi*yAj, (i=1p,j=1q)为pxq的系数的矩阵%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