计算国土面积--MATLAB课程设计(二).doc

《数值分析》课程设计论文
根据某国地图求出该国的国土面积
09信科2班 龚林园 20090202201
摘 要：生产实践中常常会出现给你一批离散样点，要求通过这些点的光滑曲面，以满足设计要求或加工要求反映在数学上即已知函数在一些点上的值，寻求它的分析表达式。因为由函数的表格形式不能直接得出表中未列点处得函数，即利用基本的插值法，在查表和求值的计算过程中，可以根据已知数据库构造插值函数来求得被插值函数在某点的精确值。使用Matlab软件来完成我们所需的值。
关键字：插值 梯形公式

一、问题重述
图3.8是某国的地图，为了计算它的国土面积，首先对地图作如下测量：以由西向东方向为轴，由南到北方向为轴，选择方便的原点，得到了表3.6、表3.7的地图测量数据，比例尺为30（数据单位）：100公里（实际单位）。试由测量数据采用插值的方法产生一张需要的地图，计算该国国土的近似面积，与它的精确值156.6500万平方公里比较。
表3.6、表3.7见附件。

二、模型假设：
在工程建设和地籍管理中，会经常遇到面积的测算和计算工作，传统的方法是在纸上，利用求积仪进行计算，存在绘图等操做的误差，精度较低。在现在工作中，全站仪的广泛使用我们能够容易得到一系列离散数据坐标，降低了对使用者的基础和计算机语言的要求，使计算不在是测算人员的负担。MATLAB提供了非常方便的绘图功能，越来越受测量人员的青睐。观察本题所给出的点，都是离散的点，我们可以根据所学的插值或拟合的方法计算。
提出以下假设建模思路：
1. 国土的边界是光滑的连续曲线，能用3次多项式表示出来。
2. 国土边界曲线有连续二阶导数。
3. 题目所给数据是有选择性的，真实地反映了国土的大概轮廓。
4. 假设测量的地图和数据准确，由最西边界与最东边界分为n条连续的边界曲线，边界内所有的土地均为该国国土。
三、问题分析：
1.从最西边界点到最东边界点，变量x∈[a,b],划分[a,b]为n小段[xi-1.xi],并由此将国土分成n小块，设每一块均为X型区域，即做垂直于x轴的直线穿过该区域，直线与边界曲线最多只有两个交点。
2.数值积分方法计算国土面积：
根据测量数据利用matlab软件对上下边界进行三次多项式插值。数值积分法的基本思想是将上边界点与下边界点分别利用插值函数求出n条曲线，则曲线所围面积即为国土面积，然后根据比例缩放关系求出国土面积近似解。在求国土面积时，利用求平面图形面积的数值积分法--将该面积分成若干个小长方形，分别求出长方形的面积后相加即为该面积的近似解。
分割替代近似代替，在每个小区间[xi，xi-1]上任取一点ei（i=1,2,…n）,作以[xi，xi-1]为底，以f（ei）为高的小矩形，用其面积f（ei）ei近似代替同底小曲边梯形的面积。

设上边界函数为f2(x),下边界函数为f1(x),由定积分定义可知曲线所围区域面积为
=

四、模型建立：
三次样条插值值法是一种分段插值法，由于在插值节点处具有二阶导数连续，从而具有更好的光滑性。三次样条插值用分段低次（不高于三次）的多项式去逼近，比较简单方便。我们可以发现在给我们的x值中有重复点，所以采用三次样条插值比较简单点。
首先，要将原图形画出来（图1），观察图形的起伏程度，方便为后面求面积分段，例如下边疆，我们要在x中相同的点处将x分段，根据分段的x和对应的y值进行三次样条插值，画出图形与原图形进行比较，观察哪里的轮廓相差比较大，在轮廓相差大的地方将x分段（根据个人观察，并不是分的越小越好），由于在x值相同处分段，用三次样条插值画出的图像与原图的轮廓相差较大，所以要在相差较大的地方在进行x分段或者用其他插值方法，比如三次插值，线性插值等（图2），直到画出的图形与原图的轮廓相差不大。上边疆同理可得。（图3）
图1(原图形)

图2（插值后比较图形）

图3（最后所得的插值后与原图比较）

其次，计算面积，由上面的分析易知，求面积是用梯形公式，由于x有重复的点，并且不是依次增大或减小，所以在计算面积时，会有重复面积，我们可以把所有面积算出来，与图形结合起来，将重复的面积减去，也可以将x的范围重新赋值，也可以避免重复面积出现。
图4（表示重复面积）

根据上述画图分段的x求面积（图5），然后根据梯形公式求出对应线段的积分(与坐标轴所围的面积)，将上边疆所求的总面积f2=s21+s22+s23+s24+s25= 2.1288e+005 (平方单位）减去下边疆所求的总面积f1=s11+s12+s13+s14+s15+s16-s12=74928（平方单位）（s12与s11面积重复，图6），即得到所求面积f=f2-f1。因为每30（单位）就是100千米，所以还要将得到的面积化成平方千米f=f/9*100= 1.5708e+006（平方千米），将这个面积与实际准确面积比较，得出相对误差p=abs(f-1566500)/1566500= 0.021538。
图5（上下边疆分段图）

图6（下边疆的分段面积）

五、程序分析：

1.根据插值函数，画出下边疆的图像（上边疆同理）
程序如下：
x11=[17 18 20 31 41 58 66 72]; %对x11赋值
y11=[299 298 288 273 262 254 234 220]; %对y11赋值
x11i=17:1:72; %确定范围以及步长
y11i=interp1(x11,y11,x11i,'spline'); %用三次样条插值
plot(x11i,y11i) %画图命令
hold on %等一下的意思，以便所有图形都画在一个图上
x=[72 72];
y=[220 207];
plot(x,y,'b')
hold on
x12=[72 69 57 60 71];
y12=[207 191 175 166 160];
x12i=57:1:72;
y12i=interp1(x12,y12,x12i); %用线性插值
plot(x12i,y12i)
hold on
图形如下：

六、 误差分析
由上面结果可以看出，所得结果与精确值有一定的误差，产生误差的原因：
1. 原始数据误差。因为所测量数据毕竟有限，不能完全表示出所测量点以外的边界的情况；
2.把数据点进行适当分组，选取合适插值方法，使插值后所画图像与原始图像更接近；
误差解决方法：
1.测量更多的数据点，使更接近于原始图像；
2.插值误差。插值的原理是近似给出所给数据点之间的数据，并不能与原始图像相吻合；
参考文献：
《计算方法—数值分析》 袁东锦 编著 南京师范大学出版社
冯康等，《数值计算方法》，北京：国防工业出版社 1978
白玉山，《计算方法》，辽宁人名出版社，1984
李庆扬，王能超，易大义，《数值分析》（第三版），华中工学院出版社，1986
邓建中等，《计算方法》，西安交通大学出版社，1985
聂铁军，《计算方法》，国防工业出版社，1988
唐珍，金坚明，李志杰，《计算方法》，高等教育出版社，1992
李岳生，黄友谦，《数值逼近》，高等教育出版社，1978
附录：
x=[723 722 710 687 676 659 647 630 619 623 626 633 608 596 581 558 537 511 484 464 456 449 434 425 411 394 368 351 332 329 312 284 281 263 251 249 244 240 247 233 222 217 209 189 180 169 165 165 150 138 138 132 127 122 102 86 65 64 54 32 28 17];
y=[225 220 240 256 256 241 245 237 245 254 273 309 308 315 315 290 281 270 270 272 278 290 293 301 303 308 297 303 311 337 342 353 358 365 356 347 346 332 314 297 290 297 298 301 303 307 314 325 328 332 337 336 341 338 332 328 322 316 314 314 307 299];
plot(x,y,'r*')
plot(x,y,'r-')
hold on
x=[17 18 20 31 41 58 66 72 72 69 57 60 71 104 130 146 160 163 168 179 196 223 258 282 307 315 330 352 377 377 392 428 462 501 524 533 555 542 550 561 574 590 599 610 635 644 649 669 671 677 678 696 720 723];
y=[299 298 288 273 262 254 234 220 207 191 175 166 160 150 137 121 117 106 83 64 63 56 50 52 46 38 32 21 21 16 14 34 43 46 60 75 95 114 138 139 133 133 139 157 162 174 188 200 207 205 206 216 218 225];
plot(x,y,'g.')
plot(x,y,'g')
hold on
x11=[17 18 20 31 41 58 66 72];
y11=[299 298 288 273 262 254 234 220];
x11i=17:1:72;
y11i=interp1(x11,y11,x11i,'spline');
plot(x11i,y11i)
hold on
x=[72 72];
y=[220 207];
plot(x,y,'b')
hold on
x12=[72 69 57 60 71];
y12=[207 191 175 166 160];
x12i=57:1:72;
y12i=interp1(x12,y12,x12i);
plot(x12i,y12i)
hold on
x13=[71 104 130 146 160 163 168 179 196 223 258 282 307 315 330 352 377];
y13=[160 150 137 121 117 106 83 64 63 56 50 52 46 38 32 21 21];
x13i=71:1:377;
x133=72:1:377;
y13i=interp1(x13,y13,x13i,'spline');
y133=interp1(x13,y13,x133,'spline');
plot(x13i,y13i)
hold on
x=[377 377];
y=[21 16];
plot(x,y,'b')
x14=[377 392 428 462 501 524];
y14=[16 14 34 43 46 60];
x14i=377:1:524;
y14i=interp1(x14,y14,x14i,'spline');
plot(x14i,y14i)
hold on
x15=[524 533 555 542 550];
y15=[60 75 95 114 138];
x15i=524:1:555;
y15i=interp1(x15,y15,x15i);
plot(x15i,y15i)
hold on
x16=[550 561 574 590 599 610 635 644 649 669 671 677 678 696 720 723];
y16=[138 139 133 133 139 157 162 174 188 200 207 205 206 216 218 225];
x16i=550:10:723;
y16i=interp1(x16,y16,x16i,'spline');
plot(x16i,y16i)
hold on
x21=[723 722 710 687 676 659 647 630 619 623 626 633 608 596];
y21=[225 220 240 256 256 241 245 237 245 254 273 309 308 315];
x21i=596:1:723;
y21i=interp1(x21,y21,x21i);
plot(x21i,y21i)
hold on
x22=[596 581 558 537 511 484 464];
y22=[315 315 290 281 270 270 271];
x22i=464:1:596;
y22i=interp1(x22,y22,x22i,'spline');
plot(x22i,y22i)
hold on
x23=[464 456 449 434 425 411 394 368 351 332 329 312 284 281 263 251 249 244 240 247 233 222 217 209 189 180 169 165];
y23=[271 278 290 293 301 303 308 297 303 311 337 342 353 358 365 356 347 346 335 314 297 290 297 298 301 303 307 314];
x23i=165:1:464;
y23i=interp1(x23,y23,x23i);
plot(x23i,y23i)
hold on
x=[165 165];
y=[314 325];
plot(x,y,'b')
x24=[165 150 138];
y24=[325 328 332];
x24i=138:1:165;
y24i=interp1(x24,y24,x24i, 'spline');
plot(x24i,y24i)
hold on
x25=[138 132 127 122 102 86 65 64 54 32 28 17];
y25=[337 336 341 338 332 328 322 316 314 314 307 299];
x25i=17:1:138;
y25i=interp1(x25,y25,x25i,'spline');
plot(x25i,y25i)
hold on
s11=trapz(x11i,y11i)
s13=trapz(x133,y133)
s14=trapz(x14i,y14i)
s16=trapz(x16i,y16i)
s15=trapz(x15i,y15i)
s21=trapz(x21i,y21i)
s22=trapz(x22i,y22i)
s23=trapz(x23i,y23i)
s24=trapz(x24i,y24i)
s25=trapz(x25i,y25i)
f2=s22+s21+s25+s23+s24
f1=s13+s11+s14+s15+s16
f=(f2-f1)/9*100
p=abs(f-1566500)/1566500
第11页