Monte Carlo求地图面积【数学建模】

https://blog.csdn.net/yillc/article/details/6740425

基本思想：所求解的问题是某随机事件A出现的概率，通过某种实验方法，得出A的频率，以此估计A事件出现的概率。

蒙特卡洛求面积：随机在矩形里面产生大量的随机点（数量为N），计算有多少点（数量为c）落在区域面积内，S是矩形的面积，那么S*（c/N）就是所求区域的面积。

1.求曲线围成的面积：

clc
clear
close all
P=rand(10000,2);
x=2*P(:,1)-1;
y=2*P(:,2);
II=find(y<=2-x.^2&y.^3>=x.^2);
M=length(II);
S=4*M/10000
plot(x(II),y(II),'g.')

2.求国土面积

clc
close all
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=[];
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=[are long*height*s/10000/18^2*1600]
end
mean(are)

每一轮循环中先用unifrnd函数生成一个点（在最大范围的矩阵中），再用interp1函数插值计算在生成点x处的函数值，判断是否在国土面积内，若是，s=s+1.

 

for循环做了十次面积求解，每次的答案存放在are（）数组中，最后取均值。