《數(shù)值分析》課程設(shè)計論文國土面積計算

《數(shù)值分析》課程設(shè)計論文 國土面積計算 摘 要：數(shù)學(xué)建模方法是處理科學(xué)理論的一種經(jīng)典方法，也是解決各類實際問題的常用方法。本文采用插值、復(fù)化梯形公式的方法，并利用數(shù)學(xué)軟件MATLAB對國土面積進行計算，使之計算結(jié)果與實際記錄基本吻合。 關(guān)鍵詞 ：插值 復(fù)化梯形公式 MATLAB 一、問題提出 圖3.8是某國的地圖，為了計算它的國土面積，首先對地圖作如下測量：以由西向東方向為軸，由南到北方向為軸，選擇方便的原點，得到了表3.6、表3.7的地圖測量數(shù)據(jù)，比例尺為30（數(shù)據(jù)單位）：100公里（實際單位）。試由測量數(shù)據(jù)采用插值的方法產(chǎn)生一張需要的地圖，計算該國國土的近似面積，與它的精確值156.6500萬平方公里進行比較。 二、模型假設(shè) 題中所給數(shù)據(jù)是真實可用的，假設(shè)國土面積在一定時期內(nèi)沒有變化，各種因素對此題也無影響，可以進行實地測量。 三、問題分析與模型建立 在工程建設(shè)和地籍管理中，會經(jīng)常遇到面積的測算和計算工作，傳統(tǒng)的方法是在紙上，利用求積儀進行計算，存在繪圖等操做的誤差，精度較低。在現(xiàn)在工作中，全站儀的廣泛使用我們能夠容易得到一系列離散數(shù)據(jù)坐標，降低了對使用者的基礎(chǔ)和計算機語言的要求，使計算不在是測算人員的負擔(dān)。MATLAB提供了非常方便的繪圖功能，越來越受測量人員的青睞?，F(xiàn)根據(jù)題目給出的數(shù)據(jù)，我們可以使用MATLAB軟件中的plot函數(shù)先大概畫出該圖形輪廓，從該圖可知本題主要是求不規(guī)則平面圖形的面積，題中已經(jīng)把圖形分為上下邊疆線，所以可以采用線性插值與復(fù)化梯形公式分別求出上下邊疆與橫坐標所圍的面積，最后求出圖形的總面積。 根據(jù)題目中給出的已知條件，應(yīng)用機理分析進行建模具體分析如下： 1. 根據(jù)對某地地圖測量的數(shù)據(jù)計算出某地地圖的近似面積。 2. 將測算數(shù)據(jù)與其精確值156.6500萬平方公里進行比較。 四、模型求解 1.插值作圖 在將數(shù)據(jù)進行插值時，由于數(shù)據(jù)中出現(xiàn)重復(fù)點，故不能直接插值，所以考慮分段處理，分別將上下邊疆線合理分段，每段記為xi(i=1,2,……),而y的值也與之對應(yīng)分段記為yi(i=1,2,……)，本文將下邊疆線分為六段，將上邊疆線分為七段，現(xiàn)以下邊疆線為例分段： 第一段：x1=[17 18 20 31 41 58 66 72]；y1=[299 298 288 273 262 254 234 220]；第二段：x2=[72 69 57]；y2=[207 191 175]；第三段：x3=[57 60 71 104 130 146 160 163 168 179 196 223 258 282 307 315 330 352 377]；y3=[175 166 160 150 137 121 117 106 83 64 63 56 50 52 46 38 32 21 21]；第四段：x4=[377 392 428 462 501 524 533 555];y4=[16 14 34 43 46 60 75 95]；第五段：x5=[555 542];y5=[95 114]；第六段：x6=[542 550 561 574 590 599 610 635 644 649 669 671 677 678 696 720 723]；y6=[114 138 139 133 133 139 157 162 174 188 200 207 205 206 216 218 225]； 下面先將每段進行一維線性插值如第一段: x1=[17 18 20 31 41 58 66 72]; y1=[299 298 288 273 262 254 234 220]; x1i=17:1:72; y1i=interp1(x1,y1,x1i); %線性插值 plot(x1i,y1i) hold on; 其中y1i則為被插值點的函數(shù)值，x1，y1為插值節(jié)點，且x1,y1單調(diào)，x1i為被插值點，x1i的步長為1，x1i的值不能超出x1的范圍，plot函數(shù)為輸出這一段的圖形，再用hold on將下面每一段數(shù)據(jù)輸出的圖形承接在一個圖中，從而構(gòu)成一幅整圖如下(其中藍線則為插值后的線)： 再考慮將每一段進行一維三次樣條插值(1)進行比較，同樣以第一段為例：y1i=interp1(x1,y1,x1i,’spline’) %三次樣條插值 plot(x1i,y1i) hold on; 其中x1,y1,x1i,y1i的含義同上，而’spline’為插值的類型是三次樣條插值，再次得出插值后的圖如下：（藍線為插值后的線) 根據(jù)以上兩圖的比較可以看出線性插值與三次樣條插值存在著一定的差別，但是只根據(jù)圖形無法判斷它們的精確性，因此下面再分別用它們插值后的線來計算面積。 2.求面積 對以上每一段所構(gòu)造的函數(shù)值用復(fù)化梯形公式（I=） （h=1）求積分，記為An，n=（1，2，……6），下以A1為例： s=0; for i=1:54 運用for循環(huán)語句 s=s+y1i(i); end A1=0.5*(299+220+2*s) %復(fù)化梯形公式 先對s賦值0，然后采用for循環(huán)語句將i從1到54（54為x1i區(qū)間的范圍）對s疊加，從而求得s的值，最后用復(fù)化梯形公式進行積分，求得A1；同理求出下，上邊疆其余段的積分分別為Ai (i=1,2,……6) ，Bj (j=1,2,……7)。 從而求出下邊疆與橫坐標的面積為=，上邊疆與橫坐標的面積為 =。 下圖分別為線性插值、三次樣條插值后下邊疆線與橫坐標所圍的面積：以下每種顏色就代表下邊疆每一段所求得的積分，將它們合并即為整個下邊疆線與橫坐標所圍面積。 而該題的面積就為上邊疆與橫坐標的面積減去下邊疆與橫坐標的面積，即為，再根據(jù)題中給出的比例尺將換成實際單位，并與精確值=156.6500萬平方公里進行比較，最后算出誤差。 五、運行結(jié)果 1.用線性插值算的的面積為：下邊疆與橫坐標的面積為=8.2899e+004，而上邊 疆與橫坐標的面積為=2.2475e+005； 所以=2.2475e+005-(8.2899e+004)=1.4185e+005，即求得的近似面積為=157.6150萬平方公里； 而題中給出的精確值為=156.6500萬平方公里，所以得出相對誤差為==0.0062。 2.而用三次樣條插值算的的面積為：下邊疆與橫坐標的面積為= 8.2780e+004，而上邊疆與橫坐標的面積為=2.2475e+005； 所以=2.2475e+005-(8.2780e+004)=1.4197e+005，即求得的近似面積為=157.7450萬平方公里； 所以得出相對誤差為== 0.0070。 六、誤差分析與結(jié)論 1.從以上兩種插值方法可以看出結(jié)果都存在著一定的誤差，產(chǎn)生誤差的原因： (1)原始數(shù)據(jù)誤差，由于所測量數(shù)據(jù)有限，不能完全表示出所測量點以外的邊界情況； (2)插值誤差，插值的原理是近似給出所給數(shù)據(jù)點之間的數(shù)據(jù)，并不能與原始圖像相吻合； (3)程序中x1i所取步長的大小不同。 2.減小誤差方法： (1)選取合適的插值方法，使得插值后的圖像更逼近原圖：將以上兩種方法進行比較可以看出用線性插值所作的圖更接近原圖，且求得的結(jié)果誤差更小，與實際記錄更接近，所以采用線性插值； (2)改變程序中x1i的步長：將線性插值程序中x1i的步長1改為0.1，其余都不改變，從而得出的相對誤差為H=0.0196，與步長為1時比較可以看出步長為1的誤差更小。 3.結(jié)論： 由以上分析可知，就這題而言應(yīng)采用線性插值，得出的近似面積為157.6150萬平方公里，求得的相對誤差為0.0062，但是這并不能肯定的說明線性插值就一定比三次樣條插值好，或者更加精確，在不同的情況下會有不同的結(jié)果，并且有些誤差是無法避免的 所以我們只能因題而議。 至此該題解答完成。 七、模型的應(yīng)用與推廣 此模型可以應(yīng)用于工程建設(shè)和地籍管理面積測量和計算工作中,全站儀的廣泛使用使我們能夠容易得到一系列離散點的坐標,利用MATLAB 這種功能強大的計算軟件，能夠比較簡單的對于邊界不規(guī)則地塊面積進行計算。此模型適用于不是十分精確的測量?？梢詫Υ笮凸こ淌┕ぬ峁┐蟾诺拿娣e和邊界圖形。還可以應(yīng)用于對湖泊、山川、河流的測量計算。 參考文獻： [1]聶鐵軍，《計算方法》,國防工業(yè)出版社，1988 [2]王能超，《數(shù)值分析簡明教程》,高等教育出版社，1984 [3]易大義等，《數(shù)值方法》，浙江科學(xué)技術(shù)出版社，1984 [4]張德榮等，《計算方法與算法語言》（第二版），人民教育出版社，1989 [5]鄧建中等，《計算方法》，西安交通大學(xué)出版社，1985 [6]唐珍，金堅明，李志杰，《計算方法》，高等教育出版社，1992 [7]李岳生，黃友謙，《數(shù)值逼近》，高等教育出版社，1978 附 錄： 完整代碼如下： X1=[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]; Y1=[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(X1,Y1,g.) plot(X1,Y1,g) hold on X2=[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]; Y2=[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(X2,Y2,r.) plot(X2,Y2,r) hold on; x1=[17 18 20 31 41 58 66 72]; y1=[299 298 288 273 262 254 234 220]; x1i=17:1:72; y1i=interp1(x1,y1,x1i); %線性插值 plot(x1i,y1i) hold on; s=0; for i=1:54 s=s+y1i(i); end A1=0.5*(299+220+2*s) %復(fù)化梯形公式 x2=[72 69 57]; y2=[207 191 175]; x2i=57:1:72; y2i=interp1(x2,y2,x2i); %線性插值 plot(x2i,y2i) hold on; s=0; for i=1:14 s=s+y2i(i); end A2=0.5*(207+175+2*s) %復(fù)化梯形公式 x3=[57 60 71 104 130 146 160 163 168 179 196 223 258 282 307 315 330 352 377]; y3=[175 166 160 150 137 121 117 106 83 64 63 56 50 52 46 38 32 21 21]; x3i=57:1:377; y3i=interp1(x3,y3,x3i ); %線性插值 plot(x3i,y3i) hold on; s=0; for i=1:319 s=s+y3i(i); end A3=0.5*(175+21+2*s) %復(fù)化梯形公式 x4=[377 392 428 462 501 524 533 555]; y4=[16 14 34 43 46 60 75 95]; x4i=377:1:555; y4i=interp1(x4,y4,x4i); %線性插值 plot(x4i,y4i) hold on; s=0; for i=1:177 s=s+y4i(i); end A4=0.5*(95+16+2*s) %復(fù)化梯形公式 x5=[555 542 ]; y5=[95 114 ]; x5i=542:1:555; y5i=interp1(x5,y5,x5i); %線性插值 plot(x5i,y5i) hold on; s=0; for i=1:12 s=s+y5i(i); end A5=0.5*(95+114+2*s) %復(fù)化梯形公式 x6=[542 550 561 574 590 599 610 635 644 649 669 671 677 678 696 720 723]; y6=[114 138 139 133 133 139 157 162 174 188 200 207 205 206 216 218 225]; x6i=542:1:723; y6i=interp1(x6,y6,x6i); %線性插值 plot(x6i,y6i) hold on; s=0; for i=1:180 s=s+y6i(i); end A6=0.5*(114+225+2*s) %復(fù)化梯形公式 M=A1+A2+A3+A4+A5+A6 x7=[723 722 710 687 676 659 647 630 619]; y7=[225 220 240 256 256 241 245 237 245]; x7i=619:1:723; y7i=interp1(x7,y7,x7i); %線性插值 plot(x7i,y7i) hold on; s=0; for i=1:103 s=s+y7i(i); end B1=0.5*(225+245+2*s) %復(fù)化梯形公式 x8=[619 623 626 633]; y8=[245 254 273 309]; x8i=619:1:633; y8i=interp1(x8,y8,x8i); %線性插值 plot(x8i,y8i) hold on; s=0; for i=1:13 s=s+y8i(i); end B2=0.5*(245+309+2*s) %復(fù)化梯形公式 x9=[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]; y9=[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]; x9i=240:1:633; y9i=interp1(x9,y9,x9i); %線性插值 plot(x9i,y9i) hold on; s=0; for i=1:392 s=s+y9i(i); end B3=0.5*(309+332+2*s) %復(fù)化梯形公式 x10=[240 247 ]; y10=[332 314 ]; x10i=240:1:247; y10i=interp1(x10,y10,x10i); %線性插值 plot(x10i,y10i) hold on; s=0; for i=1:6 s=s+y10i(i); end B4=0.5*(314+332+2*s) %復(fù)化梯形公式 x11=[247 233 222 217 209 189 180 169 165]; y11=[314 297 290 297 298 301 303 307 314]; x11i=165:1:247; y11i=interp1(x11,y11,x11i); %線性插值 plot(x11i,y11i) hold on; s=0; for i=1:81 s=s+y11i(i); end B5=0.5*(314+314+2*s) %復(fù)化梯形公式 x12=[165 150 138]; y12=[325 328 332]; x12i=138:1:165; y12i=interp1(x12,y12,x12i); %線性插值 plot(x12i,y12i); hold on s=0; for i=1:26 s=s+y12i(i); end B6=0.5*(325+332+2*s) %復(fù)化梯形公式 x13=[138 132 127 122 102 86 65 64 54 32 28 17]; y13=[337 336 341 338 332 328 322 316 314 314 307 299]; x13i=17:1:138; y13i=interp1(x13,y13,x13i); %線性插值 plot(x13i,y13i) hold on; s=0; for i=1:120 s=s+y13i(i); end B7=0.5*(337+299+2*s) %復(fù)化梯形公式 N=B1+B2+B3+B4+B5+B6+B7 C=N-M E=C*(100/9)/10000 C0=156.6500 H=(E-C0)/C0 - 15 -