Matlab曲面面积估计 (2).doc

2013-2014第一学期数学软件与数学实验课程设计2013年11月4日11月8日实验题目 曲面面积估计 组员1组员2组员3姓名学号专业信息与计算科学信息与计算科学成绩数学实验报告实验名称曲面面积估计问题背景描述：在一次对山区种植的农作物进行产量评估时，需要对某一区域的地表面积进行估计，由于山区地表并不平整，所以常规的面积计算公式不能直接应用，而需要求地表曲面的面积。现在通过遥感技术在海拔为300m的高空探测到这一矩形区域1, 22, 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) colorbar5、利用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.001;y3=y;y4=y+0.001;将x1、x2、x3、x4、y1、y2、y3、y4带入（1）方程，求出相应的zx=x+0.001;y=y+0.001;（转下页）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);主要内容（要点）：（模型的求解原理、公式、推导、基本求解步骤、算法的流程图等）（接上页）： 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=2-0.001?Y N x=3-0.001?YN sum end1、将曲面分成若干个小三角形，由每个点算出三角形的边长，由海伦公式计算三角形的面积。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=0x11=x1+i*lengthx;x22=x11+lengthx*0.001;y11=y1+i*lengthy;y22=y11+lengthy*0.001;带入方程算出相应的zleng=(lengthx*0.001)2+(lengthy*0.001)2+(z1-z2)2;length=sqrt(leng)+length;i=1Y输出lengthend（转下页）主要内容（要点）：（模型的求解原理、公式、推导、基本求解步骤、算法的流程图等）（接上页）：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.855.39 5.49 5.51 5.46 5.32 5.14 4.94 4.74 4.59 4.49 4.485.61 5.77 5.81 5.71 5.51 5.23 4.90 4.59 4.36 4.21 4.195.73 5.92 5.97 5.86 5.62 5.27 4.88 4.51 4.23 4.05 4.035.74 5.93 5.98 5.86 5.62 5.28 4.88 4.51 4.21 4.04 4.025.63 5.79 5.84 5.74 5.53 5.23 4.91 4.59 4.33 4.18 4.165.42 5.53 5.56 5.49 5.35 5.16 4.93 4.73 4.55 4.45 4.445.14 5.18 5.19 5.17 5.12 5.05 4.97 4.90 4.84 4.81 4.804.84 4.80 4.79 4.82 4.87 4.94 5.02 5.10 5.16 5.19 5.204.56 4.45 4.43 4.49 4.64 4.84 5.06 5.28 5.45 5.55 5.564.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 offrotate3d onfigurecontour(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*x2 + p11*x*y + p02*y2 + p30*x3 + p21*x2*y + p12*x*y2 + p03*y3 + p40*x4 + p31*x3*y + p22*x2*y2+ p13*x*y3 + p04*y4 + p50*x5 + p41*x4*y + p32*x3*y2+ p23*x2*y3 + p14*x*y4;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=-(8568512622238087*x15)/1152921504606846976-(179477452349969*x14*y1)/18014398509481984+(1144995169262675*x14)/144115188075855872-(1211*x13*y12)/10000-(519*x13*y1)/5000+(1109*x13)/5000-(113*x12*y13)/5000+(7681*x12*y12)/100000+(1307*x12*y1)/10000-(1469*x12)/10000-(5521*x1*y14)/100000-(1053*x1*y13)/10000+(6997*x1*y12)/10000+(6711*x1*y1)/10000-(9761*x1)/10000 + (1148886279340723*y14)/144115188075855872+(3689*y13)/100000 - (291*y12)/2000 - (1809*y1)/10000 + 1043/200;z2=-(8568512622238087*x25)/1152921504606846976-(179477452349969*x24*y2)/18014398509481984+(1144995169262675*x24)/144115188075855872-(1211*x23*y22)/10000 - (519*x23*y2)/5000 +(1109*x23)/5000 - (113*x22*y23)/5000 + (7681*x22*y22)/100000+(1307*x22*y2)/10000-(1469*x22)/10000-(5521*x2*y24)/100000-(1053*x2*y23)/10000+(6997*x2*y22)/10000+(6711*x2*y2)/10000-(9761*x2)/10000+(1148886279340723*y24)/144115188075855872+(3689*y23)/100000 - (291*y22)/2000 - (1809*y2)/10000 + 1043/200; z3=-(8568512622238087*x35)/1152921504606846976-(179477452349969*x34*y3)/1801439850948198+(1144995169262675*x34)/144115188075855872-(1211*x33*y32)/10000-(519*x33*y3)/5000+(1109*x33)/5000-(113*x32*y33)/5000+ (7681*x32*y32)/100000+(1307*x32*y3)/10000-(1469*x32)/10000-(5521*x3*y34)/100000 - (1053*x3*y33)/10000 + (6997*x3*y32)/10000 +1043/200; （转下页）实验过程记录（含：主要程序清单以及异常情况及程序调试过程记录等）（接上页）：z4=-(8568512622238087*x45)/1152921504606846976-(179477452349969*x44*y4)/18014398509481984+(1144995169262675*x44)/144115188075855872-(1211*x43*y42)/10000 - (519*x43*y4)/5000 +(1109*x43)/5000 - (113*x42*y43)/5000 + (7681*x42*y42)/100000+(1307*x42*y4)/10000-(1469*x42)/10000-(5521*x4*y44)/100000 - (1053*x4*y43)/10000 + (6997*x4*y42)/10000 + (6711*x4*y4)/10000 - (9761*x4)/10000 + (1148886279340723*y44)/144115188075855872 + (3689*y43)/100000 - (291*y42)/2000 - (1809*y4)/10000 + 1043/200; 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); endendsum（转下页）实验过程记录（含：主要程序清单以及异常情况及程序调试过程记录等）（接上页）：（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=-(8568512622238087*x115)/1152921504606846976 - (179477452349969*x114*y11)/18014398509481984 .+ (1144995169262675*x114)/144115188075855872 - (1211*x113*y112)/10000 - (519*x113*y11)/5000 + .(1109*x113)/5000 - (113*x112*y113)/5000 + (7681*x112*y112)/100000 + (1307*x112*y11)/10000 - .(1469*x112)/10000 - (5521*x11*y114)/100000 - (1053*x11*y113)/10000 + (6997*x11*y112)/10000 + .(6711*x11*y11)/10000 - (9761*x11)/10000 + (1148886279340723*y114)/144115188075855872 + .(3689*y113)/100000 - (291*y112)/2000 - (1809*y11)/10000 + 1043/200;z2=-(8568512622238087*x225)/1152921504606846976 - (179477452349969*x224*y22)/18014398509481984 .+ (1144995169262675*x224)/144115188075855872 - (1211*x223*y222)/10000 - (519*x223*y22)/5000 + .(1109*x223)/5000 - (113*x222*y223)/5000 + (7681*x222*y222)/100000 + (1307*x222*y22)/10000 - .(1469*x222)/10000 - (5521*x22*y224)/100000 - (1053*x22*y223)/10000 + (6997*x22*y222)/10000 + .(6711*x22*y22)/10000 - (9761*x22)/10000 + (1148886279340723*y224)/144115188075855872 + .(3689*y223)/100000 - (291*y222)/2000 - (1809*y22)/10000 + 1043/200;leng=(lengthx*0.001)2+(lengthy*0.001)2+(z1-z2)2;length=sqrt(leng)+length;end实验结果报告与实验总结：曲面图与等高线拟合工具1.曲面方程：f=(8568512622238087*x5)/115292150460684697-(179477452349969*x4*y)/18014398509481984+(1144995169262675*x4)/144115188075855872-(1211*x3*y2)/10000-(519*x3*y)/5000 + (1