曲面拟合原理与实例说课材料.docx
- 文档编号:5571482
- 上传时间:2022-12-27
- 格式:DOCX
- 页数:14
- 大小:210.19KB
曲面拟合原理与实例说课材料.docx
《曲面拟合原理与实例说课材料.docx》由会员分享,可在线阅读,更多相关《曲面拟合原理与实例说课材料.docx(14页珍藏版)》请在冰豆网上搜索。
曲面拟合原理与实例说课材料
曲面拟合原理与实例
问题:
给定一组坐标
,
,表示有n个点。
要求用以下二元多项式函数对所给的坐标进行拟合:
即
设
则函数又可表示为
,拟合的目标就是求出系数矩阵A。
最小二乘法:
构造关于系数
的多元函数:
点(
,…,
)是多元函数
的极小点,其中
为权函数,默认为1,所以点(
,…,
)必须满足方程组
在
的情况下,有
因此可得
令
,
则
上式实际共有
个等式,可将这
个等式写成矩阵的形式有:
也就是U*a=V的形式,其中
,
,
U为
阶矩阵,实现函数为functionA=leftmatrix(x,p,y,q);V为长
的列向量,实现函数为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
请你估计出此地区内(
)煤的储量,单位用立方米表示,并用电脑画出该煤矿的三维图象。
如果直接画出三维曲面图形:
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=
1234
1234
1234
1234
1234
Y=
1111
2222
3333
4444
5555
Z=
13.720015.470023.280019.9500
25.800021.330026.480023.7300
8.470014.490029.140015.3500
25.270024.830012.040018.0100
22.320026.190014.580016.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);
进行三次多项式插值:
xi=linspace(1,4,31);yi=linspace(1,5,41);[XI,YI]=meshgrid(xi,yi);
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*H
n=
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('errorcheckcheck!
');
end
U_length=p*q;%U为p*q阶方阵
U=zeros(U_length,U_length);%赋值0,目的是分配内存
fori=1:
p*q
forj=1:
p*q
x_z=quotient(j-1,q)+quotient(i-1,q);%x的幂的次数,quotient为求商
y_z=mod(j-1,q)+mod(i-1,q);%y的幂的次数
U(i,j)=qiuhe(x,x_z,y,y_z);
end
end
rightmatrix.m
functionV=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.m
functionsh=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);
if(nargin<4)&(m~=length(y))%输入量至少为四,x,y行向量长度必需一样
error('errorcheckcheck!
');
end
ifnargin==4%没有z,默认为元素全部为1的向量
z=ones(m,1);
end
he=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.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=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的矩阵aa
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(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);
%zz=0;%zz是xx,yy代入所拟合的函数求出的函数值
fori=1:
p%函数为Σaa(i,j)*x^i*y^j,(i=1...p,j=1...q)
forj=1:
q%aa为pxq的系数的矩阵
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.26780.2132-0.06240.0058
-0.72871.3972-0.92750.2412-0.0210
0.4416-0.84150.5487-0.14070.0122
-0.06800.1295-0.08390.0214-0.0018
注意:
权函数在拟合的函数非常重要,不过她只能按你遇到的具体问题来取,我这里为1;当p,q越大时,拟合的函数与原数据的方差越小,但是有可能函数本身抖动非常厉害,可以画图看出来。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 曲面 拟合 原理 实例 材料