matlab实践（二）：使用蒙卡罗特方法计算面积

1.题目

2.蒙卡罗特方法

法国数学家Ｂｕｆｆｏｎ于１８世纪首次提出使用随机投针的方式来估算圆周率，虽然受限于当时的实验条件，其结果精度并不高，但却充分体现了蒙特卡洛方法的思想。进入２０世纪中后期随着计算机的进步和核物理领域对随机实验需求的增长，蒙特卡洛方法得到了充分发展，目前已广泛运用到了不确定性分析、粒子输运、数理统计、可靠性分析、经济学、稀薄气体动力学等领域。
使用解析方法或数值方法可以有效的、确定的解决各种问题，但在实际应用领域中的各种条件与关系往往十分复杂，数学模型只能在其大体正确的方向上进行充分简化才有可能求解，而这样的结果却可能与实际情况已相去甚远。而蒙特卡洛方法却能够在条件关系相对复杂或高维性突出的情况下对复杂系统直接进行模拟，决定了其具有强大的发展生命力。随着计算能力和算法的发展，随机数的产生和随机抽样方法已经得以可靠解决，在基本蒙特卡洛方法的基础上，马尔科夫链蒙特卡洛方法、拟蒙塔卡洛方法、序贯蒙塔卡洛方法等方法得到了广泛应用。
同时，蒙特卡洛方法又被比喻为“最后的方法”，意为在能够使用解析方法或数值方法时尽量不要使用蒙特卡洛方法，也意为在其他方法都不能解决问题的情况下，也一定能通过蒙特卡洛方法进行模拟，它是解决问题的最后的方法。

出自张　骁，刘丙杰的《基于蒙特卡洛方法的火箭残骸落点范围预测》。

3.origin图像数字化工具

我们对福建省进行操作。

得出数据并保存在“福建.xlsx”中。

用matlab作图

4.判定点在区域内

我们编写函数inPloy,若点的上下左右都有区域上的点，那么其在区域内，返回1，反之，则在外面，返回0。

function flags = inPoly(p,poly)
flags=0;
a=0;
b=0;
c=0;
d=0;
for i=1:134
    if poly(i,1)>=p(1)&&poly(i,2)>=p(2)
        a=1;
    elseif poly(i,1)<=p(1)&&poly(i,2)>=p(2)
        b=1;
    elseif poly(i,1)>=p(1)&&poly(i,2)<=p(2)
        c=1;
    elseif poly(i,1)<=p(1)&&poly(i,2)<=p(2)
        d=1;
    end
end
if a+b+c+d==4
    flags=1;
end

end

结果如图所示：

 

5.判定在直线两侧

我们用两向量的叉积来评定，若其大于0，这说明在顺时针方向；若其小于于0，这说明在逆时针方向；计算大于0的数目，和区域内总的数目比较即可。调整经过福州的直线，若其比值在0.5左右即可。

​
P1=x1-Y(end,1);
Q1=y1-Y(end,2);
P2=m(i)-Y(end,1);
Q2=n(i)-Y(end,2);
if P1*Q2-P2*Q1>0
      o=o+1;
end

​

6.总代码

clc;clear;
Y=xlsread('福建.xlsx');
figure
plot(Y(1:end-1,1),Y(1:end-1,2),'.');
l=10000;
P=rand(l,2);
x=10.*P(:,1);
y=10.*P(:,2);
pionts=zeros(1,l);
k=1;
for i=1:l
    p=[x(i),y(i)];
    pionts(i)=inPoly(p,Y(1:end-1,1:2));
    if pionts(i)==1
        m(k)=x(i);
        n(k)=y(i);
        k=k+1;
    end
end
x1=0;
o=0;
for b=-pi/2:pi/10:pi/2
    for i=1:k-1
       y1=tan(b)*(x1-Y(end,1))+Y(end,2);
       P1=x1-Y(end,1);
       Q1=y1-Y(end,2);
       P2=m(i)-Y(end,1);
       Q2=n(i)-Y(end,2);
       if P1*Q2-P2*Q1>0
           o=o+1;
       end
    end
    if o/k>0.45&&o/k<0.55
        b
        break;
    else o=0;
    end
end
figure(1);
plot(x,y,'b.');
hold on;
plot(m,n,'rX');
hold on;
x2=0:0.1:10;
y2=tan(b)*(x2-Y(end,1))+Y(end,2);
plot(x2,y2,'k-','Linewidth',2);
hold on;
plot(Y(end,1),Y(end,2),'y.','MarkerSize',30);
hold off

结果：