国土面积计算

为了算出瑞士的国土面积，首先对瑞士地图作如下测量：以由西向东方向为 x轴，由南到北方向为y 轴，选择方便的原点，并将从最西边界点到最东边界点在 轴上的区间适当地划分为若干段，在每个分点的y方向测出南边界点和北边界点的 y1坐标 和y2 ，这样就得到了表中的数据（单位mm）。
根据地图的比例我们知道18mm相当于40km，试由测量数据计算瑞士国土的近似面积，与它的精确值41288km 比较。

x  7.0  10.5  13.0  17.5  34.0  40.5  44.5   48.0  56.0  61.0  68.5  76.5  80.5  91.0
 
y  44    45   47    50    50    3 8   30     30     34    36   34    41   45     46
 
y  44    59   70    72    93    100   110    110    110  117   118   116  118    118
 
x  96   101  104  106.5  111.5  118   123.5  136.5  142  146   150   157   158 
 
y  43    37    33    28    32    65    55    54     52   50    66    66    68
 
y  121  124   121   121    121   122   116   83     81   82    86    85    68
                   表  瑞士地图测量数据 

一 采取线性差值计算

1.模型的假设:

(1) 假设测量的地图和数据准确，由最西边界与最东边界分为两条连续的边界线，边界内所有的土地均为该国国土。

(2) 假设从最西边界点到最东边界点，变量x∈[a,b],划分[a,b]为n小段[xi-1.xi],并由此将国土分成n小块，设每一块均为矩形区域，长为xi-xi-1。

2.计算国土面积：

如果分成n小块后，近似的矩形高度为左端两点的距离，则用matlab编程如下：
x  =[7.0  10.5  13.0  17.5  34.0  40.5  44.5   48.0  56.0  61.0  68.5  76.5  80.5  91.0  96   101  104  106.5  111.5  118   123.5  136.5  142  146   150   157   158 ]*40/18;
y1= [44    45   47    50    50    38   30     30     34    36   34    41   45     46    43    37    33    28    32    65    55    54     52   50    66    66    68]*40/18;
y2= [44    59   70    72    93    100   110    110    110  117   118   116  118    118  121   124   121   121    121   122   116   83     81   82    86    85    68]*40/18;
plot(x,y1,'*');
hold on
plot(x,y2,'*');
sum=0;
for k=1:26
    sum=sum+(x(k+1)-x(k))*(y2(k)-y1(k))
end

所计算出的面积值为42560km2,与标准值41288km2相差3.08%

x  =[7.0  10.5  13.0  17.5  34.0  40.5  44.5   48.0  56.0  61.0  68.5  76.5  80.5  91.0  96   101  104  106.5  111.5  118   123.5  136.5  142  146   150   157   158 ]*40/18;
y1= [44    45   47    50    50    38   30     30     34    36   34    41   45     46    43    37    33    28    32    65    55    54     52   50    66    66    68]*40/18;
y2= [44    59   70    72    93    100   110    110    110  117   118   116  118    118  121   124   121   121    121   122   116   83     81   82    86    85    68]*40/18;
%subplot(1,2,1),
plot(x,y1,'*');
%subplot(1,2,2),
hold on
plot(x,y2,'*');
sum=0;
for k=1:26
    sum=sum+(x(k+1)-x(k))*(y2(k+1)-y1(k+1))
end

如果分成n小块后，近似的矩形高度为右端两点的距离，则用matlab编程如下：   

所计算出的面积值为42267km2,与标准值41288km2相差2.37%

如果两段当y1=y2时，不采用矩形计算，而采用三角形计算，会更加精确：

x  =[7.0  10.5  13.0  17.5  34.0  40.5  44.5   48.0  56.0  61.0  68.5  76.5  80.5  91.0  96   101  104  106.5  111.5  118   123.5  136.5  142  146   150   157   158 ]*40/18;
y1= [44    45   47    50    50    38   30     30     34    36   34    41   45     46    43    37    33    28    32    65    55    54     52   50    66    66    68]*40/18;
y2= [44    59   70    72    93    100   110    110    110  117   118   116  118    118  121   124   121   121    121   122   116   83     81   82    86    85    68]*40/18;
plot(x,y1,'*');
hold on
plot(x,y2,'*');
sum=0;
sum=(x(2)-x(1))*(y2(2)-y1(2))/2;%最左边的三角形
sum=sum+(x(27)-x(26))*(y2(26)-y1(26))/2;%最右边的三角形
for k=2:25
    sum=sum+(x(k+1)-x(k))*(y2(k)-y1(k));
end
sum1=0;
sum1=(x(2)-x(1))*(y2(2)-y1(2))/2;%最左边的三角形
sum1=sum1+(x(27)-x(26))*(y2(26)-y1(26))/2;%最右边的三角形
for k=2:25
    sum1=sum1+(x(k+1)-x(k))*(y2(k+1)-y1(k+1));
end   

 
所计算出的面积值为42635km2 , 42193km2
 
总体来说，采用线性差值是此类问题最简单的计算方法，由于本题给出的领土较为规则，所以用线性差值计算的领土面积较为精确。

二 采取三次多项式差值后数值积分法计算

1.模型的假设:

(1) 假设测量的地图和数据准确，由最西边界与最东边界分为两条连续的边界曲线，边界内所有的土地均为该国国土。

(2) 假设从最西边界点到最东边界点，变量x∈[a,b],划分[a,b]为n小段[xi-1.xi],并由此将国土分成n小块，设每一块均为X型区域，即做垂直于x轴的直线穿过该区域，直线与边界曲线最多只有两个交点。

2.数值积分方法计算国土面积：

根据测量数据利用matlab软件对上下边界进行三次多项式差值。数值积分法的基本思想是将上边界点与下边界点分别利用差值函数求出两条曲线，则曲线所围面积即为国土面积，然后根据比例缩放关系求出国土面积近似解。在求国土面积时，利用求平面图形面积的数值积分法--将该面积分成若干个小长方形，分别求出长方形的面积后相加即为该面积的近似解。

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

程序如下：

clear
x  =[7.0  10.5  13.0  17.5  34.0  40.5  44.5   48.0  56.0  61.0  68.5  76.5  80.5  91.0  96   101  104  106.5  111.5  118   123.5  136.5  142  146   150   157   158 ]
y1= [44    45   47    50    50    38   30     30     34    36   34    41   45     46    43    37    33    28    32    65    55    54     52   50    66    66    68]
y2= [44    59   70    72    93    100   110    110    110  117   118   116  118    118  121   124   121   121    121   122   116   83     81   82    86    85    68]
plot(x,y1,'*');
hold on
plot(x,y2,'*');
newx=7:0.1:158;
newy1=interp1(x,y1,newx,'spline');
newy2=interp1(x,y2,newx,'spline');
s=sum(newy2-newy1)*0.1/18^2*1600

（图2）

得到答案：42468KM2

三 蒙特卡洛方法

蒙特卡洛方法基本思想：当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。

 蒙特卡洛方法求积分的规则：任何一个积分，都可看做是某个随机变量的期望值，因此，可以用这个随机变量的平均值来近似它。

 

求解国土面积的思路：由插值求出国土边界的大致形状，如图二，然后通过计算机模拟，产生随机点落入矩形区域 (矩形面积为S）内，进行N次试验，设落入曲线范围内m点，则面积公式：Area  ≈ S*m/N，N越大，越趋向真实值。

 

1.模型的假设:

(1) 假设测量的地图和数据准确，由最西边界与最东边界分为两条连续的边界曲线，边界内所有的土地均为该国国土。

(2) 假设从最西边界点到最东边界点，变量x∈[a,b],划分[a,b]为n小段[xi-1.xi],并由此将国土分成n小块，设每一块均为X型区域，即做垂直于x轴的直线穿过该区域，直线与边界曲线最多只有两个交点。

 

clear
x  =[7.0  10.5  13.0  17.5  34.0  40.5  44.5   48.0  56.0  61.0  68.5  76.5  80.5  91.0  96   101  104  106.5  111.5  118   123.5  136.5  142  146   150   157   158 ];
y1= [44    45   47    50    50    38   30     30     34    36   34    41   45     46    43    37    33    28    32    65    55    54     52   50    66    66    68];
y2= [44    59   70    72    93    100   110    110    110  117   118   116  118    118  121   124   121   121    121   122   116   83     81   82    86    85    68];
long=max(x)-min(x);
height=max(y2)-min(y1);
are=zeros(10);
for k=1:10
    s=0;
    for i=1:1000
        ranx(i)=unifrnd(min(x),max(x));
        rany(i)=unifrnd(min(y1),max(y2));
        newy1=interp1(x,y1,ranx(i),'spline');
        newy2=interp1(x,y2,ranx(i),'spline');
        if(rany(i)>=newy1&rany(i)<=newy2) s=s+1;
        end
    end
    are(k)=long*height*s/10000/18^2*1600
end

此代码模拟了10次试验，将所求面积保存在了are()数组里，

一次结果为：

4295.1
4187.7
4245
4130.5
4245
4359.5
4338.1
4209.2
4151.9
4101.8
求其平均值：4226，此结果比较精确，如果想提高精确程度，可以增加模拟次数，但对计算机的运算水平有更高的要求。