曲面拟合原理与实例.docx
- 文档编号:10892999
- 上传时间:2023-02-23
- 格式:DOCX
- 页数:14
- 大小:195KB
曲面拟合原理与实例.docx
《曲面拟合原理与实例.docx》由会员分享,可在线阅读,更多相关《曲面拟合原理与实例.docx(14页珍藏版)》请在冰豆网上搜索。
曲面拟合原理与实例
问题:
多项式函数对所给的坐标进行拟合:
f(x,y)
i1jaijxyij1,1
1
ai
i1j1
i1j1jxy
f(x,y)
a11
a12y
2a13y
L
q1
a1qy
a21x
a22xy
2a23xy
L
q1a2qxy
M
i1
i1
i12
L
i1q1
ai1xi1
ai2xi1y
ai3xi1y2
aiqxi1yq1
M
p1
p1
p12
L
p1q
ap1x
ap2xy
ap3xy
apqxy
p,q
p
q
即
1
最小二乘法:
构造关于系数aij的多元函数:
点(a11,⋯,apq)是多元函数s(a11,L,apq)的极小点,其中g为权函数,默认为1,所以点(a11,⋯,apq)必须满足方程组
saij
在g1的情况下,有
2
[f(xg,yg)zg]g1
U为pqpq阶矩阵,实现函数为functionA=leftmatrix(x,p,y,q);V为长pq的列向量,实现函数为functionB=rightmatrix(x,p,y,q,z)。
这样就可以算出列矩阵a,然后转化成A。
例子:
某地区有一煤矿,为估计其储量以便于开采,先在该地区进行勘探。
假设
该地区是一长方形区域,长为4公里,宽为5公里。
经勘探得到如下数据:
煤矿勘探数据表
编号
1
2
3
4
5
6
7
8
9
10
横坐标(公里)
1
1
1
1
1
2
2
2
2
2
纵坐标(公里)
1
2
3
4
5
1
2
3
4
5
煤层厚度(米)
13.72
25.80
8.47
25.27
22.32
15.47
21.33
14.49
24.83
26.19
编号
11
12
13
14
15
16
17
18
19
20
横坐标(公里)
3
3
3
3
3
4
4
4
4
4
纵坐标(公里)
1
2
3
4
5
1
2
3
4
5
煤层厚度(米)
23.28
26.48
29.14
12.04
14.58
19.95
23.73
15.35
18.01
16.29
请你估计出此地区内(2x4,1y5)煤的储量,单位用立方米表示,并用电脑画出该煤矿的三维图象。
如果直接画出三维曲面图形:
clear;x=1:
4;y=1:
5;[X,Y]=meshgrid(x,y)
Z=[13.7225.808.4725.2722.32;
surf(X,Y,Z);
X=
1
2
3
4
1
2
3
4
1
2
3
4
1
2
3
4
1
2
3
4
Y=
1
1
1
1
2
2
2
2
3
3
3
3
4
4
4
4
5
5
5
5
Z=
13.720015.470023.280019.9500
25.8000
21.3300
26.4800
23.7300
8.4700
14.4900
29.1400
15.3500
25.2700
24.8300
12.0400
18.0100
22.3200
26.1900
14.5800
16.2900
粗略计算体积:
底面积乘以平均高度。
p=sum(Z);
q=p(:
[2,3,4]);
h=sum(q')/15
v=2000*4000*h
h=
20.0773
v=
1.6062e+008
进行线性插值:
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);
41);[XI,YI]=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,Z,XI,YI,'cubic');
surf(XI,YI,ZI);
H=0;n=0;
forj=21:
61
fori=1:
81
H=H+ZI(i,j);n=n+1;
end
end
n
H=H/n
S=2000*4000;
V=S*Hn=
3321
H=
20.8222
V=
1.6658e+008
上面是插值的方法解题,下面用拟合的方法解题。
为此编写了几个M函数:
leftmatrix.m
functionU=leftmatrix(x,p,y,q)
%U*a=Va为系数列矩阵,长度为p*q
%U为左边p*q乘p*q矩阵
%x,y为长度一致的列矩阵,给定点的坐标
%p,q为拟合的函数中x,y的幂的最高次数
m=length(x);
if(nargin~=4)&(m~=length(y))
error(end
'errorcheckcheck!
'
);
U_length=p*q;
U=zeros(U_length,U_length);
fori=1:
p*q
forj=1:
p*q
x_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的幂的次数
为求商
rightmatrix.mfunctionV=rightmatrix(x,p,y,q,z)%U*a=V
%V为一个列向量长为p*q
%xyz为点的坐标
%pq分别为xy幂的最高次数
ifnargin~=5
error('errorcheckcheck!
rightmatrix'
end
V=zeros(p*q,1);
fori=1:
p*q
x_z=quotient(i-1,q);y_z=mod(i-1,q);
V(i,1)=qiuhe(x,x_z,y,y_z,z);end
quuotient.mfunctionsh=quotient(x,y)%sh为x/y的商
sh=(x-mod(x,y))/y;
qiuhe.m
functionhe=qiuhe(x,p,y,q,z)
%hex^p*y^q从1->m的和
%x,y向量长度相同
%p,q分别为x,y的幂的次数
m=length(x);
%输入量至少为四,x,y行向量长度必需一样
默认为元素全部为1的向量
if(nargin<4)&(m~=length(y))
error('errorcheckcheck!
');end
ifnargin==4%没有z,
z=ones(m,1);
endhe=0;
fori=1:
m
he=he+x(i)^p*y(i)^q*z(i);%1-->m求和
end
下面一段程序先进行拟合,然后验证拟合的效果,具体操作:
先输入x=⋯y=⋯z=⋯p=⋯q=(注意x,y,z是向量);拟合得到系数a,也就是得到了拟合的函数;根据拟合函数计算给定点(xx,yy)的函数值zz=f(xx,yy)并进行画图检验程序保存于M文件fit.m。
fit.m
clear;[X,Y]=meshgrid(1:
4,1:
5);
Z=[13.72
25.80
8.47
25.27
22.32;
15.47
21.33
14.49
24.83
26.19;
23.28
26.48
29.14
12.04
14.58;
19.95
23.73
15.35
18.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=V
V=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);
end
aam=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代入所拟合的函数求出的函数值
函数为Σaa(i,j)*x^i*y^j,(i=1...p,j=1...q)
为pxq的系数的矩阵
%zz=0;%zz
fori=1:
p%
forj=1:
q%aa
xt=xx.^(i-1);
yt=yy.^(j-1);
xy=xt.*yt;zz=zz+aa(i,j).*xy;
end
end
ZI=reshape(zz,n,m);surf(XI,YI,ZI);%axis([1415030])aa=
1.0e+003*
0.1465
-0.2678
0.2132
-0.0624
0.0058
-0.7287
1.3972
-0.9275
0.2412
-0.0210
0.4416
-0.8415
0.5487
-0.1407
0.0122
-0.0680
0.1295
-0.0839
0.0214
-0.0018
注意:
权函数在拟合的函数非常重要,不过她只能按你遇到的具体问题来取,我这里为1;当p,q越大时,拟合的函数与原数据的方差越小,但是有可能函数本身抖动非常厉害,可以画图看出来。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 曲面 拟合 原理 实例