曲面拟合原理与实例文档格式.docx
- 文档编号:6868321
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:14
- 大小:195KB
曲面拟合原理与实例文档格式.docx
《曲面拟合原理与实例文档格式.docx》由会员分享,可在线阅读,更多相关《曲面拟合原理与实例文档格式.docx(14页珍藏版)》请在冰点文库上搜索。
ap2xy
ap3xy
apqxy
p,q
p
q
即
最小二乘法:
构造关于系数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公里。
经勘探得到如下数据:
煤矿勘探数据表
编号
3
4
5
6
7
8
9
10
横坐标(公里)
纵坐标(公里)
煤层厚度(米)
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
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=
Y=
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=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);
ZI=interp2(X,Y,Z,XI,YI,'
cubic'
H=0;
n=0;
forj=21:
61
fori=1:
81
H=H+ZI(i,j);
n=n+1;
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:
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('
rightmatrix'
V=zeros(p*q,1);
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的幂的次数
%输入量至少为四,x,y行向量长度必需一样
默认为元素全部为1的向量
if(nargin<
4)&
(m~=length(y))
);
ifnargin==4%没有z,
z=ones(m,1);
endhe=0;
m
he=he+x(i)^p*y(i)^q*z(i);
%1-->
m求和
下面一段程序先进行拟合,然后验证拟合的效果,具体操作:
先输入x=⋯y=⋯z=⋯p=⋯q=(注意x,y,z是向量);
拟合得到系数a,也就是得到了拟合的函数;
根据拟合函数计算给定点(xx,yy)的函数值zz=f(xx,yy)并进行画图检验程序保存于M文件fit.m。
fit.m
[X,Y]=meshgrid(1:
4,1:
5);
Z=[13.72
22.32;
26.19;
14.58;
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;
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);
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
p%
forj=1:
q%aa
xt=xx.^(i-1);
yt=yy.^(j-1);
xy=xt.*yt;
zz=zz+aa(i,j).*xy;
ZI=reshape(zz,n,m);
%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文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 曲面 拟合 原理 实例