计算国土面积MATLAB课程设计二Word格式.docx

1、在现在工作中，全站仪的广泛使用我们能够容易得到一系列离散数据坐标，降低了对使用者的基础和计算机语言的要求，使计算不在是测算人员的负担。MATLAB提供了非常方便的绘图功能，越来越受测量人员的青睐。观察本题所给出的点，都是离散的点，我们可以根据所学的插值或拟合的方法计算。提出以下假设建模思路：1. 国土的边界是光滑的连续曲线，能用3次多项式表示出来。2. 国土边界曲线有连续二阶导数。3. 题目所给数据是有选择性的，真实地反映了国土的大概轮廓。4. 假设测量的地图和数据准确，由最西边界与最东边界分为n条连续的边界曲线，边界内所有的土地均为该国国土。三、问题分析：1.从最西边界点到最东边界点，变量x

2、a,b,划分a,b为n小段xi-1.xi,并由此将国土分成n小块，设每一块均为X型区域，即做垂直于x轴的直线穿过该区域，直线与边界曲线最多只有两个交点。2.数值积分方法计算国土面积：根据测量数据利用matlab软件对上下边界进行三次多项式插值。数值积分法的基本思想是将上边界点与下边界点分别利用插值函数求出n条曲线，则曲线所围面积即为国土面积，然后根据比例缩放关系求出国土面积近似解。在求国土面积时，利用求平面图形面积的数值积分法-将该面积分成若干个小长方形，分别求出长方形的面积后相加即为该面积的近似解。分割替代近似代替，在每个小区间xi，xi-1上任取一点ei（i=1,2,n）,作以xi，xi-

3、1为底，以f（ei）为高的小矩形，用其面积f（ei）ei近似代替同底小曲边梯形的面积。 设上边界函数为f2（x）,下边界函数为f1（x）,由定积分定义可知曲线所围区域面积为=四、模型建立：三次样条插值值法是一种分段插值法，由于在插值节点处具有二阶导数连续，从而具有更好的光滑性。三次样条插值用分段低次（不高于三次）的多项式去逼近，比较简单方便。我们可以发现在给我们的x值中有重复点，所以采用三次样条插值比较简单点。首先，要将原图形画出来（图1），观察图形的起伏程度，方便为后面求面积分段，例如下边疆，我们要在x中相同的点处将x分段，根据分段的x和对应的y值进行三次样条插值，画出图形与原图形进行比较，

4、观察哪里的轮廓相差比较大，在轮廓相差大的地方将x分段（根据个人观察，并不是分的越小越好），由于在x值相同处分段，用三次样条插值画出的图像与原图的轮廓相差较大，所以要在相差较大的地方在进行x分段或者用其他插值方法，比如三次插值，线性插值等（图2），直到画出的图形与原图的轮廓相差不大。上边疆同理可得。（图3） 图1（原图形） 图2（插值后比较图形） 图3（最后所得的插值后与原图比较）其次，计算面积，由上面的分析易知，求面积是用梯形公式，由于x有重复的点，并且不是依次增大或减小，所以在计算面积时，会有重复面积，我们可以把所有面积算出来，与图形结合起来，将重复的面积减去，也可以将x的范围重新赋值，也可

5、以避免重复面积出现。 图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-1566

6、500）/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 20

7、7;plot（x,y,b）hold onx12=72 69 57 60 71;y12=207 191 175 166 160;x12i=57: y12i=interp1（x12,y12,x12i）; %用线性插值 plot（x12i,y12i） hold on图形如下：六、误差分析 由上面结果可以看出，所得结果与精确值有一定的误差，产生误差的原因：1. 原始数据误差。因为所测量数据毕竟有限，不能完全表示出所测量点以外的边界的情况；2.把数据点进行适当分组，选取合适插值方法，使插值后所画图像与原始图像更接近；误差解决方法：1.测量更多的数据点，使更接近于原始图像；2.插值误差。插值的原理是近似给

8、出所给数据点之间的数据，并不能与原始图像相吻合；参考文献：计算方法数值分析 袁东锦 编著 南京师范大学出版社冯康等，数值计算方法，北京：国防工业出版社 1978白玉山，计算方法，辽宁人名出版社，1984李庆扬，王能超，易大义，数值分析（第三版），华中工学院出版社，1986邓建中等，计算方法，西安交通大学出版社，1985聂铁军，计算方法，国防工业出版社，1988唐珍，金坚明，李志杰，计算方法，高等教育出版社，1992李岳生，黄友谦，数值逼近，高等教育出版社，1978附录：x=723 722 710 687 676 659 647 630 619 623 626 633 608 596 581 5

9、58 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 33

10、7 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;r*r-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 6

11、35 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;g.gy11=299 298 288 273 262 254 234 220;x11i=17:y11i=interp1（x1

12、1,y11,x11i,plot（x11i,y11i）y12i=interp1（x12,y12,x12i）;plot（x12i,y12i）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:377;x133=72:y13i=interp1（x13,y13,x13i,y133=interp1（x13,y13,x133,plot（x13i,y13i）x=377 377;y=2

13、1 16;x14=377 392 428 462 501 524;y14=16 14 34 43 46 60;x14i=377:524;y14i=interp1（x14,y14,x14i,plot（x14i,y14i）x15=524 533 555 542 550;y15=60 75 95 114 138;x15i=524:555;y15i=interp1（x15,y15,x15i）;plot（x15i,y15i）x16=550 561 574 590 599 610 635 644 649 669 671 677 678 696 720 723;y16=138 139 133 133 139

14、 157 162 174 188 200 207 205 206 216 218 225;x16i=550:10:723;y16i=interp1（x16,y16,x16i,plot（x16i,y16i）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:y21i=interp1（x21,y21,x21i）;plot（x21i,y21i）x22=596 581 558 537 511 484

15、 464;y22=315 315 290 281 270 270 271;x22i=464:596;y22i=interp1（x22,y22,x22i,plot（x22i,y22i）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 2

16、97 298 301 303 307 314;x23i=165:464;y23i=interp1（x23,y23,x23i）;plot（x23i,y23i）x=165 165;y=314 325;x24=165 150 138;y24=325 328 332;x24i=138:165;y24i=interp1（x24,y24,x24i, plot（x24i,y24i）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:138;y25i=i

17、nterp1（x25,y25,x25i,plot（x25i,y25i）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+s24f1=s13+s11+s14+s15+s16f=（f2-f1）/9*100p=abs（f-1566500）/1566500