压实系数计算公式Matlab曲面面积估计 (2)

2013-2014第一学期
《数学软件与数学实验》
课程设计


2013年11月4日－11月8日


实验题目 曲面面积估计



姓名
学号
专业
科学
成绩
科学

组员1


组员2


组员3



信息与计算信息与计算



数学实验报告
实验
名称
问题背景描述：
在一次对山区种植的农作物进行产量评估时，需要对某一区域 的地表面积进行估
计，由于山区地表并不平整，所以常规的面积计算公式不能直接应用，而需要求地表< br>曲面的面积。现在通过遥感技术在海拔为300m的高空探测到这一矩形区域[1, 2]×[2,
3] 内一些地面点的垂直距离如表所示（表中数据均减去了200）。假设所观测的矩形
区域 内均是陆地且整个地表面可以看作是一张光滑曲面。

曲面面积估计

实验目的与任务：
1.用matlab计算出所观测区域近似的地表曲面方程(提示，可基于 Matlab的二维插值
命令实现)；并作出其三维曲面图。
2.算出观测区域地表表面积的近似值。
3.编程计算从矩形区域一边上任意一点到其对边上 任意一点沿直线走向的地表曲线的
近似长度。
实验原理与数学模型：
（含模型的假设、符号说明、模型的建立）

将山体的 表面看成是一张光滑的曲面，先将给出的点用matlab画出后用拟合工具
箱将其拟合成方程，用微分 的思想将曲面分成若干个小三角形，算出每个小面积Si相
加从而算出整个曲面面积近似值。
S=∑Si,(i=1,2,…，n)
曲线长度也用微分的方法计算。






实验所用软件及版本：

Matlab（R2010b）

主要内容（要点）：
（模型的求解原理、公 式、推导、基本求解步骤、算法的流程图等）

(1)


































开始
定义x，y
并给定对应z值
用linspace在
x=[2,3]内产生80
个元素
用linspacw在
y=[1,2]内产生80
个元素
绘制网格线，画
出曲面图
Sftool拟合工具
箱产生拟合函数
结束

主要内容（要点）：
（模型的求解原理、公式、推导、基本求解步骤、算法的 流程图等）
（接
上页）
：
（1）
1、有题目数据表格给定，定义x，y范围以及步长，并将x、y所对应的值赋值给z
x=2:0.1:3;y=1.0:0.1:2;
2、用linspace产生多个元素个数，并用产生的数据绘制网格
xi=linspace(2,3,80);yi=linspace(1,2,80);
[XI,YI]=meshgrid(xi,yi);
3、利用三次样条（cubic）插值，并绘出曲面图像
ZI=interp2(X,Y,Z,XI,YI,'cubic');surfc(XI,YI,ZI)
4、绘制原始网格曲面与数据点以及等高线
plot3(X,Y,Z,'go','markeredgecolor','c')
plot3(X',Y',Z','g')
hold off
rotate3d on
figure
contour(XI,YI,ZI,40)
colorbar
5、利用sftool拟合工具箱产生拟合函数
sftool(X,Y,Z)


（2）
开始



sum=0；x=2；

y=1；



x1=x;x2=x;x3=x+0.001;x4

=x+0.001;y1=y;y2=y+0.00

1;y3=y;y4=y+0.001;



将x1、x2、x3、x4、y1、y2、

y3、y4带入（1）方程，求出相

应的z






主要内容（要点）：
（模型的求解原理、公式、推导、基本求 解步骤、算法的流程图等）
（接
上页）
：

a1=(x1-x2)^2+(y1-y2)^2+(z1-z2)^2;a1

=sqrt(a1);

b1=(x4-x2)^2+(y4-y2)^2+(z4-z2)^2;b1

=sqrt(b1);



c=(x1-x4)^2+(y1-y4)^2+(z1-z4)^2;c=sqrt(c);

a2=(x1-x3)^2+(y1-y3)^2+(z1-z3)^2;a2=sqrt(a2);

b2=(x4-x3)^2+(y4-y3)^2+(z4-z3)^2;b2=sqrt(b2);

p1=(a1+b1+c)2;p2=(a2+b2+c)2;


sum=sum+sqrt(abs(p1*(p1-a1)*(p1-b1)*(p1-c)));

sum=sum+sqrt(abs(p2*(p2-a2)*(p2-b2)*(p2-c)));



x=x+0.001;y=y+0.001;


Y


N
Y

x<=3-0.001?

N

sum


end


1、将 曲面分成若干个小三角形，由每个点算出三角形的边长，由海伦公式计算三角形
的面积。a1=(x1- x2)^2+(y1-y2)^2+(z1-z2)^2;a1=sqrt(a1);
b1=(x4-x2)^2+(y4-y2)^2+(z4-z2)^2;b1=sqrt(b1);
c=(x1-x4)^2+(y1-y4)^2+(z1-z4)^2;c=sqrt(c);
a2=(x1-x3)^2+(y1-y3)^2+(z1-z3)^2;a2=sqrt(a2);
b2=(x4-x3)^2+(y4-y3)^2+(z4-z3)^2;b2=sqrt(b2);
p1=(a1+b1+c)2;
p2=(a2+b2+c)2;

主要内容 （要点）：
（模型的求解原理、公式、推导、基本求解步骤、算法的流程图等）
（接
上 页）
：
2、将每个小面积利用循环相加，得总面积。
sum=sum+sqrt(abs(p1*(p1-a1)*(p1-b1)*(p1-c)));
sum=sum+sqrt(abs(p2*(p2-a2)*(p2-b2)*(p2-c)));

（3）

开始




length(x1,y1,x2,y2)




lengthx=abs(x1-x2);

lengthy=abs(y1-y2);

length=0;



i=0


x11=x1+i*lengthx;x22=x11+lengthx*0.001;

y11=y1+i*lengthy;y22=y11+lengthy*0.001;




带入方程算出相应的

z



leng=(lengthx*0.001)^2+(lengthy*0.001)^2+ (z1-z2)^2;

length=sqrt(leng)+length;

Y


i<=1



输出length


end





主要内容（要点）：
（模型的求解原理、公式、推导、基本求解步骤、算法的流程图 等）
（接
上页）
：
1、定义function length=length(x1,y1,x2,y2)
2、步长取0.001，分别赋值
x11=x1+i*lengthx;x22=x11+lengthx*0.001;
y11=y1+i*lengthy;y22=y11+lengthy*0.001;
3、带入方程，算出相应的z。
4、利用公式算出最终长度
leng= (lengthx*0.001)^2+(lengthy*0.001)^2+(z1-z2)^2;
length=sqrt(leng)+length;



































实验过程记录（含：主要程序清单及异常情况记录等）：
（1）
x=2:0.1:3;
y=1.0:0.1:2;
[X,Y]=meshgrid(x,y);
Z=[
5.11 5.13 5.14 5.13 5.09 5.04 4.98 4.93 4.89 4.85 4.85
5.39 5.49 5.51 5.46 5.32 5.14 4.94 4.74 4.59 4.49 4.48
5.61 5.77 5.81 5.71 5.51 5.23 4.90 4.59 4.36 4.21 4.19
5.73 5.92 5.97 5.86 5.62 5.27 4.88 4.51 4.23 4.05 4.03
5.74 5.93 5.98 5.86 5.62 5.28 4.88 4.51 4.21 4.04 4.02
5.63 5.79 5.84 5.74 5.53 5.23 4.91 4.59 4.33 4.18 4.16
5.42 5.53 5.56 5.49 5.35 5.16 4.93 4.73 4.55 4.45 4.44
5.14 5.18 5.19 5.17 5.12 5.05 4.97 4.90 4.84 4.81 4.80
4.84 4.80 4.79 4.82 4.87 4.94 5.02 5.10 5.16 5.19 5.20
4.56 4.45 4.43 4.49 4.64 4.84 5.06 5.28 5.45 5.55 5.56
4.36 4.19 4.16 4.25 4.47 4.76 5.09 5.41 5.66 5.81 5.83
];

xi=linspace(2,3,80);
yi=linspace(1,2,80);
[XI,YI]=meshgrid(xi,yi);
ZI=interp2(X,Y,Z,XI,YI,'cubic');
surfc(XI,YI,ZI) %光滑曲面+等高线
hold on
%原始网格曲面与数据点
plot3(X,Y,Z,'go','markeredgecolor','c')
plot3(X',Y',Z','g')
hold off
rotate3d on
figure
contour(XI,YI,ZI,40)
colorbar






实验过程记录（含：主要程序清单以及异常情况及程序调试过程
记录等）
（接上页）
：

异常：利用拟合出来的方程画图，与原图像相差甚大
解决：将变量经行预处理
syms x y p00 p01 p10 p20 p11 p02 p30 p21 p12 p03 p40 p31 p22 p13 p04 p50 p41
p32 p23 p14;
f=p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 + p30*x^3 + p21*x^2*y +
p12*x*y^2 + p03*y^3 + p40*x^4 + p31*x^3*y + p22*x^2*y^2+ p13*x*y^3 + p04*y^4
+ p50*x^5 + p41*x^4*y + p32*x^3*y^2+ p23*x^2*y^3 + p14*x*y^4;
t00 =5.215;t10 =-0.9761;t01 =-0.1809;t20 =-0.1469;t11 =0.6711;
t02 =-0.1455;t30 =0.2218;t21 =0.1307;t12 =0.6997;t03 =0.03689;
t40 =0.007945;t31 =-0.1038;t22 =0.07681;t13 =-0.1053;t04 =0.007972;
t50 =-0.007432;t41 =-0.009963;t32 =-0.1211;t23 =-0.0226;t14 =-0.05521;
x=2:0.01:3;
y=1.0:0.01:2;
[X,Y]=meshgrid(x,y);
X=(X- 2.5)0.3175;
Y=(Y- 1.5)0.3175;

f=subs(f,[p00 p01 p10 p20 p11 p02 p30 p21 p12 p03 p40 p31 p22 p13 p04 p50 p41
p32 p23 p14],[t00 t01 t10 t20 t11 t02 t30 t21 t12 t03 t40 t31 t22 t13 t04 t50
t41 t32 t23 t14])
Z=subs(f,{'x','y'},{X,Y})
mesh(X,Y,Z)
ezmesh(f,[2 3 1 2])










（转下页）

实验过程记录（ 含：主要程序清单以及异常情况及程序调试过程
记录等）
（接上页）
：
（2）
sum=0;
for y=1:0.001:2-0.001
for x=2:0.001:3-0.001
x1=x;x2=x;x3=x+0.001;x4=x+0.001;
y1=y;y2=y+0.001;y3=y;y4=y+0.001;
z1=-(85685*x 1^5)11529215-(9969*x1^4*y1)
18984+(62675*x1^4)5 855872-(1211*x1^3*y1
^2)10000-(519*x1^3*y1)5000 +(1109*x1^3)5000-(113*x1^2*y1^3)5000+(7681*x
1^ 2*y1^2)100000+(1307*x1^2*y1)10000-(1469*x1^2)10000 -(5521*x1*y1^4)100
000-(1053*x1*y1^3)10000+(699 7*x1*y1^2)10000+(6711*x1*y1)10000-(9761*x1)
100 00 + (40723*y1^4)5855872+(3689*y1^3)100000 -
(291*y1^2)2000 - (1809*y1)10000 + 1043200;
z2=-(85685*x2^5)11529215-(9969*x2^4*y2)
18984+ (62675*x2^4)5855872-(1211*x2^3*y2
^2)10000 - (519*x2^3*y2)5000 +(1109*x2^3)5000 - (113*x2^2*y2^3)5000 +
(7681*x2^2*y2^2)100000+( 1307*x2^2*y2)10000-(1469*x2^2)10000-(5521*x2*y2
^4)100000-(1053*x2*y2^3)10000+(6997*x2*y2^2)10000 +(6711*x2*y2)10000-(9
761*x2)10000+(40723*y2^4) 5855872+(3689*y2^3)10000
0 - (291*y2^2)2000 - (1809*y2)10000 + 1043200;
z3=-(85685*x3^5)11 529215-(9969*x3^4*y3)
1898+(62675*x3^4)5855872- (1211*x3^3*y3^
2)10000-(519*x3^3*y3)5000+(1109* x3^3)5000-(113*x3^2*y3^3)5000+
(7681*x3^2*y3^2 )100000+(1307*x3^2*y3)10000-(1469*x3^2)10000-(5521 *x3*y3
^4)100000 - (1053*x3*y3^3)10000 + (6997*x3*y3^2)10000 +1043200;









（转下页）

实验过程记录（ 含：主要程序清单以及异常情况及程序调试过程
记录等）
（接上页）
：
z4 =-(85685*x4^5)11529215-(9969*x4^4*y4)
18984+(62 675*x4^4)5855872-(1211*x4^3*y4
^2)10000 - (519*x4^3*y4)5000 +(1109*x4^3)5000 - (113*x4^2*y4^3)5000 +
(7681*x4^2*y4^2)100000+( 1307*x4^2*y4)10000-(1469*x4^2)10000-(5521*x4*y4
^4)100000 - (1053*x4*y4^3)10000 + (6997*x4*y4^2)10000 + (6711*x4*y4)10000
-

a1=(x1-x2)^2+(y1-y2)^2+(z1-z2)^2;a1=sqrt(a1);
b1=(x4-x2)^2+(y4-y2)^2+(z4-z2)^2;b1=sqrt(b1);
c=(x1-x4)^2+(y1-y4)^2+(z1-z4)^2;c=sqrt(c);
a2=(x1-x3)^2+(y1-y3)^2+(z1-z3)^2;a2=sqrt(a2);
b2=(x4-x3)^2+(y4-y3)^2+(z4-z3)^2;b2=sqrt(b2);
p1=(a1+b1+c)2;
p2=(a2+b2+c)2;
sum=sum+sqrt(abs(p1*(p1-a1)*(p1-b1)*(p1-c)));
sum=sum+sqrt(abs(p2*(p2-a2)*(p2-b2)*(p2-c)));
end
end
sum















（转下页）
(9761*x4)10000 + (40723*y4^4)5855872 +
(3689*y4^3)100000 - (291*y4^2)2000 - (1809*y4)10000 + 1043200;

实验过程记录（含：主要程序清 单以及异常情况及程序调试过程
记录等）
（接上页）
：

（3）
function length=length(x1,y1,x2,y2)
lengthx=abs(x1-x2);lengthy=abs(y1-y2);
length=0;
for i=0:0.001:1
x11=x1+i*lengthx;x22=x11+lengthx*0.001;
y11=y1+i*lengthy;y22=y11+lengthy*0.001;
z1=-(85685*x11^5)11529215
(9969*x11^4*y11)18984 ...
+ (62675*x11^4)5855872 - (1211*x11^3*y11^2)10000 -
(519*x11^3*y11)5000 + ...
(1109*x11^3)5000 - (113*x11^2*y11^3)5000 + (7681*x11^2*y11^2)100000 +
(1307*x11^2*y11)10000 - ...
(1469*x11^2)10000 - (5521*x11*y11^4)100000 - (1053*x11*y11^3)10000 +
(6997*x11*y11^2)10000 + ...
(6711*x11*y11)10000 - (9761*x11)10000 +
(40723*y11^4)5855872 + ...
(3689*y11^3)100000 - (291*y11^2)2000 - (1809*y11)10000 + 1043200;
z2=-(85685*x22^5)11529215
(9969*x22^4*y22)18984 ...
+ (62675*x22^4)5855872 - (1211*x22^3*y22^2)10000 -
(519*x22^3*y22)5000 + ...
(1109*x22^3)5000 - (113*x22^2*y22^3)5000 + (7681*x22^2*y22^2)100000 +
(1307*x22^2*y22)10000 - ...
(1469*x22^2)10000 - (5521*x22*y22^4)100000 - (1053*x22*y22^3)10000 +
(6997*x22*y22^2)10000 + ...
(6711*x22*y22)10000 - (9761*x22)10000 +
(40723*y22^4)5855872 + ...
(3689*y22^3)100000 - (291*y22^2)2000 - (1809*y22)10000 + 1043200;
leng=(lengthx*0.001 )^2+(lengthy*0.001)^2+(z1-z2)^2;
length=sqrt(leng)+length;
end

-
-
实验结果报告与实验总结：

曲面图与等高线
拟合工具

1.曲面方程：
f=(85685*x^5)11529215 -(9969*x^4*y)180143
98509481984+(62675*x^4)5855 872-(1211*x^3*y^2)10000
-(519*x^3*y)5000
(5521*x*y^4)100000
(6711*x*y)10000
+
+
-
(1109*x^3)5000
(1053*x*y^3)10000
-
+
-
+
(113*x^2*y^3)5000
- (1469*x^2)10000
(6997*x*y^2)10000
(3689*y^3)100000
+
-
+
+
-
(7681*x^2*y^2)100000 (1307*x^2*y)10000
(9761*x)10000
(40723*y^4)5855872
(291*y^2)2000 - (1809*y)10000 + 1043200
2.面积：sum =6.1627
3.长度：
length(2,1,3,2)
ans =9.9096

思考与深入：
1、做惯了平面的面积计算，山体表面积的计算对我们来说试一次全新的尝试。 通过本
次实验，我们知道算一个物体的表面积可以有多种方法，而且各种各样的表面积都可
以进 行计算。
2、这个题目考察了matlab软件中for语句的灵活运用，以及各种语句之间的配合，
体现了该软件灵活丰富的编程功能。
3、本次实验锻炼了我们分析实际问题、转化问题、并用所学解决问题的能力。
4、此次编程 让我们学会使用了matlab自带的拟合工具，从而可以直接得到方程，但
是刚开始的使用遇到各种我 们没有见过甚至不懂的语句，通过百度以及老师的细心讲
解最终得到了想要的结果。
5、实际 问题转化为数学问题，再将数学问题经过模型建立、简化，运用软件解决问题
的能力，为以后的数模打下 了基础。
6、实验中小组成员之间的相互协作，提高了解决问题的效率，同时增强了我们的团队
意识，在几天的共同努力之下，完成了相应的题目。
7、一次的成功将是下一次最好的开始，我们会继续学习探究matlab这一数学软件。
















参考文献：
1、matlab——help sftool
2、百度百科——contour
http:k?url=YFPwiOFCqkjZvICG782d1W7ABrrnw-ggo3o nZDQdNn
Byt7HYbVA93rXKu_EswHjPLIvt74Hgyo3wT1Ahn d8RIa
3、网络资源：基于双线性插值拟合的山形曲面面积计算_王胜兵