曲面拟合原理与实例说课材料Word格式文档下载.docx
《曲面拟合原理与实例说课材料Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《曲面拟合原理与实例说课材料Word格式文档下载.docx（14页珍藏版）》请在冰豆网上搜索。
上式实际共有
个等式，可将这
个等式写成矩阵的形式有：
也就是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
横坐标（公里）
纵坐标（公里）
煤层厚度（米）
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
请你估计出此地区内（
）煤的储量，单位用立方米表示，并用电脑画出该煤矿的三维图象。
如果直接画出三维曲面图形：
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
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）;
进行三次多项式插值：
cubic'
进行插值后计算体积：
xi=linspace（1,4,61）;
yi=linspace（1,5,81）;
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*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!
'
U_length=p*q;
%U为p*q阶方阵
U=zeros（U_length,U_length）;
%赋值0，目的是分配内存
fori=1:
p*q
forj=1:
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
rightmatrix.m
functionV=rightmatrix（x,p,y,q,z）
%U*a=V
%V为一个列向量长为p*q
%xyz为点的坐标
%pq分别为xy幂的最高次数
ifnargin~=5
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）;
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的幂的次数
if（nargin<
4）&
（m~=length（y））%输入量至少为四，x,y行向量长度必需一样
ifnargin==4%没有z,默认为元素全部为1的向量
z=ones（m,1）;
he=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）;
;
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的矩阵aa
ii=quotient（i-1,q）+1;
%quotient求商
jj=mod（i-1,q）+1;
aa（ii,jj）=a_n（i,1）;
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代入所拟合的函数求出的函数值
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;
ZI=reshape（zz,n,m）;
%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越大时，拟合的函数与原数据的方差越小，但是有可能函数本身抖动非常厉害，可以画图看出来。
升级会员 