用几种数值积分的方法计算西藏聂荣县面积

摘要：利用测量地图所得数据，根据所学插值和数值积分的方法，在MATLAB上画出地图，并计算地图边界曲线积分，将所求积分值代入地图面积推导公式，算出地图所表示的实际面积。

关键词：面积，插值，积分，MATLAB

1.原始数据的测量

图1是中国西藏聂荣县的地图，为了算出它的土地面积，首先对地图作如下测量：以由西向东方向为x轴，由南向北方向为y轴，选择方便的原点，并将从最西边界点到最东边界点在x轴上的区间适当地划分为若干段，在每个分点的y方向测出南边界点和北边界点的y轴坐标y1,y2,这样就得到了表1的数据。

表1 地图边界点数据

X0 0.9 1.3 1.7 2.1 2.6 3.0 3.2 3.3 3.5 4.0 4.6 5.5 5.5

Y1 6.5 6.0 5.6 5.6 5.0 4.5 4.2 4.1 4.1 3.9 3.6 3.0 3.4

Y2 6.5 7.3 7.2 7.3 7.7 7.3

7.6 9.9 10.1 10.5 10.2 9.8 9.7

X0 6.0 6.3 6.6 7.0 7.3 7.8 8.1 8.4 8.8 9.3 9.6 10.0 10.2

Y1 3.5 3.5 3.5 3.4 3.1 3.2 3.4 3.3 3.3 3.8 3.2 2.8 2.5

Y2 8.8 8.6 8.2 7.7 7.5 7.7 8.2 8.9 8.2 9.7 9.7 9.7 10.0

X0 10.6 10.8 11.1 11.5 11.7 12.0 12.5 12.9 13.6 14.0 14.5

Y1 2.2 2.4 2.6 2.0 2.6 2.9 2.9 3.4 1.5 1.1 1.0

Y2 10.0 10.4 10.3 10.3 10.5 10.5 10.4 10.2 9.6 9.7 8.9

X0 15.1 15.6 16.1 16.7 17.0 17.5 17.6 17.8

Y1 0.7 1.9 1.8 2.2 2.6 2.6 2.7 3.0

Y2 7.9 7.6 5.9 5.8 3.4 3.3 3.2 3.0

2.面积计算公式的推导

可把地图看作两条交于最西边界点a和最东边界点b的曲线所围成的图形，记靠南的一条曲线为L1,记靠北的一条曲线为L2。

首先根据插值法在MATLAB上画出L1、L2的图形，具体做法是：利用表1的数据，把(x0,y1)看作L1上的已知节点，用分段线性插值或者三次样条插值求出a、b之间步长为h=0.1的未知插值点x的插值y01,再用plot函数画出点(x,y01)所连成的曲线即L1；利用表1的数据，把(x0,y2)看作L2上的已知节点，用分段线性插值或者三次样条插值求出a、b之间步长为0.1的未知插值点x的插值y02，再用plot函数画出点(x,y02)所连成的曲线即L2。然后计算曲线L1、L2的积分，并求其差值即为所求地图面积，具体做法是：利用h、y01，根据梯形公式或者辛普森公式求出L1的曲线积分z1，利用h、y02，根据梯形公式或者辛普森公式求出L2的曲线积分z2,两曲线积分之差z2-z1即为所求地图面积。最后根据地图比例尺1：1000000换算出地图所表示的土地实际面积S：

S=(L2的曲线积分－L1的曲线积分)×100

＝(z2－z1)×100 (平方公里)

3.MATLAB编程

3.1用三次样条插值和复化辛普森公式计算：

>>clear %清除当前工作区的所有变量

>> x0=[0.9 1.3 1.7 2.1 2.6 3 3.2

3.3 3.5 4.0 4.6 5.0 5.5 6.0 6.3 6.6 7.0 7.3 7.8 8.1 8.4 8.8 9.3 9.6

10.0 10.2 10.6 10.8 11.1 11.5 11.7 12.0 12.5 12.9 13.6 14.0 14.5

15.1 15.6 16.1 16.7 17.0 17.5 17.6 17.8];

>> y1=[6.5 6.0 5.6 5.6 5.0 4.5 4.2

4.1 4.1 3.9 3.6 3.0 3.4 3.5 3.5 3.5 3.4 3.1 3.2 3.4 3.3 3.3 3.8 3.2

2.8 2.5 2.2 2.4 2.6 2.0 2.6 2.9 2.9 3.4 1.5 1.1 1.0 0.7 1.9 1.8 2.2

2.6 2.6 2.7 3.0]; %已知节点(x0,y1)

>> y2=[6.5 7.3 7.2 7.3 7.7 7.3 7.6

9.9 10.1 10.5 10.2 9.8 9.7 8.8 8.6 8.2 7.7 7.5 7.7 8.2 8.9 8.2 9.7

9.7 9.7 10.0 10.1 10.4 10.3 10.3 10.5 10.5 10.4 10.2 9.6 9.7 8.9

7.9 7.6 5.9 5.8 3.4 3.3 3.2 3.0]; %已知节点(x0,y2)

>> h=0.1; %产生插值点的步长

>>

x=0.9:h:17.8; %产生插值点x

>>

y01=spline(x0,y1,x); %计算L1的三次样条插值

>> y02=spline(x0,y2,x);

%计算L2的三次样条插值

>>

plot(x,y01,'k',x,y02,'r') %三次样条插值作图

>> k=length(x);

>> y011=[y01(2:2:k-1)];

>> s011=sum(y011);

>> y012=[y01(3:2:k-1)];

>> s012=sum(y012);

>>

z1=(y01(1)+y01(k)+4*s011+2*s012)*h/3 %用辛普森公式计算L1的积分

z1 =

52.3696

>> y021=[y02(2:2:k-1)];

>> s021=sum(y021);

>> y022=[y02(3:2:k-1)];

>> s022=sum(y022);

>>

z2=(y02(1)+y02(k)+4*s021+2*s022)*h/3 %用辛普森公式计算L2的积分

z2 =

143.5365

>>

S=(z2-z1)*100 %代入面积推导公式

S =

9116.6

※图略※

根据三次样条插值方法作图即3.1中图

3.2用分段线性插值和梯形公式计算：

>>clear %清除当前工作区的所有变量

>> x0=[0.9 1.3 1.7 2.1 2.6 3 3.2

3.3 3.5 4.0 4.6 5.0 5.5 6.0 6.3 6.6 7.0 7.3 7.8 8.1 8.4 8.8 9.3 9.6

10.0 10.2 10.6 10.8 11.1 11.5 11.7 12.0 12.5 12.9 13.6 14.0 14.5

15.1 15.6 16.1 16.7 17.0 17.5 17.6 17.8];

>> y1=[6.5 6.0 5.6 5.6 5.0 4.5 4.2

4.1 4.1 3.9 3.6 3.0 3.4 3.5 3.5 3.5 3.4 3.1 3.2 3.4 3.3 3.3 3.8 3.2

2.8 2.5 2.2 2.4 2.6 2.0 2.6 2.9 2.9 3.4 1.5 1.1 1.0 0.7 1.9 1.8 2.2

2.6 2.6 2.7 3.0]; %已知节点(x0,y1)

>> y2=[6.5 7.3 7.2 7.3 7.7 7.3 7.6

9.9 10.1 10.5 10.2 9.8 9.7 8.8 8.6 8.2 7.7 7.5 7.7 8.2 8.9 8.2 9.7

9.7 9.7 10.0 10.1 10.4 10.3 10.3 10.5 10.5 10.4 10.2 9.6 9.7 8.9

7.9 7.6 5.9 5.8 3.4 3.3 3.2 3.0]; %已知节点(x0,y2)

>>

h=0.1; %产生插值点的步长

>>

x=0.9:h:17.8; %产生插值点x

>>

y01=interp1(x0,y1,x); %计算L1的分段线性插值

>>

y02=interp1(x0,y2,x); %计算L2的分段线性插值

>>

plot(x,y01,'k',x,y02,'r') %分段线性插值作图

>>

z1=trapz(y01)*h %用梯形公式计算L1 的积分

z1 =

52.4750

>>

z2=trapz(y02)*h %用梯形公式计算L2的积分

z2 =

143.6400

>>

S=(z2-z1)*100 %代入面积推导公式

S=

9116.50

※图略※

根据分段线性插值方法作图即3.2中图

4.结果分析

西藏聂荣县的实际面积是14540平方公里，两种计算方法的结果近似，其中3.1中计算所得面积占实际面积的62.7%，误差还是相当大的，可以从增加已知节点的数目、更准确绘取地图、选取更恰当比例尺方面进行改进。