蒙特卡罗法的投针法模拟试验

最新推荐文章于 2023-01-13 17:24:12 发布

easy_R

最新推荐文章于 2023-01-13 17:24:12 发布

阅读量697

点赞数 2

分类专栏： optimization

本文链接：https://blog.csdn.net/easy_R/article/details/114137664

版权

optimization 专栏收录该内容

6 篇文章 0 订阅

订阅专栏

    突然关注到蒙特卡罗法，觉得挺有意思，在CSDN搜了一下相关原理和程序，组织了一下并运行模拟结果，为加深理解此处做一下记录。

直接上程序：

clc;
clear all;
tic;
l = 1.32;     % 针的长度（任意给的）
a = 2.63;    % 平行线的宽度(大于针的长度l即可)
n = 10000;    % 做n次投针试验
m = 0;    % 记录针与平行线相交的次数
x = rand(1, n) * a / 2 ;   % 在[0, a/2]内服从均匀分布随机产生n个数， x中每一个元素表示针的中点和最近的一条平行线的距离
phi= rand(1, n) * pi;    % 在[0, pi]内服从均匀分布随机产生n个数，phi中的每一个元素表示针和最近的一条平行线的夹角
% axis([0,pi, 0,a/2]);   box on;  % 画一个坐标轴的框架，x轴位于0-pi，y轴位于0-a/2， 并打开图形的边框
for i=1:n  % 开始循环，依次看每根针是否和直线相交
    if x(i) <= l / 2 * sin(phi (i))     % 如果针和平行线相交
        m = m + 1;    % 那么m就要加1
        plot(phi(i), x(i), 'r.')   % 模仿书上的那个图，横坐标为phi，纵坐标为x , 用红色的小点进行标记
        hold on  % 在原来的图形上继续绘制
    end
end
p = m / n;    % 针和平行线相交出现的频率
mypi = (2 * l) / (a * p);  % 我们根据公式计算得到的pi
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])%输出结果
toc;
hold on;
Datas=[phi;l / 2 * sin(phi)]';  %注意此处的纵坐标边缘线是计算的边界值。其实这里可以直接画出其边界不用后续的凹凸判定

t=Datas;
t(t(:,1)==0|t(:,2)==0,:)=[];
x=t(:,1);
y=t(:,2);
figure(1)
ang=[0 60 120 180];
for i=1:4
t2=ScatterHull(t,ang(i));
t2=[t2;t2(1,:)];%画图使封闭
figure (2)
subplot(2,2,i)
hold on
plot(t(:,1),t(:,2),'ro')
plot(t2(:,1),t2(:,2),'g-')
legend('scatter','HullCurve','Location','NorthWest')
title(num2str(ang(i)))
end


function y2=ScatterHull(Data,ang_2)
% y1(:,1)=x;y1(:,2)=y;ang_2:concave min angle
% ang_2=180:conhull;ang_2=0:all scatters
y1=unique(Data,'rows');% Remove duplicate rows 
%%
% sort of scatters
    cen=mean(y1);
    ang=atan2(y1(:,1)-cen(1),y1(:,2)-cen(2)); %每个点到坐标中心极角
    y1=[y1,ang];
    y1=sortrows(y1,3);    %按极角排序
    y1=y1(:,1:2);
%%
y2=y1;
im=1;
[n m]=size(y1);
%%
while length(im)>0
 
im=[];
%向量计算，判断凹点
for i=1:n    
    if i==1  %处理第一个点
        v1=y1(n,:)-y1(1,:);
        v2=y1(2,:)-y1(1,:); 
    elseif i==n %最后一个点
        v1=y1(n-1,:)-y1(n,:);    
        v2=y1(1,:)-y1(n,:);        
    else 
        v1=y1(i-1,:)-y1(i,:);     
        v2=y1(i+1,:)-y1(i,:);
    end
    r=det([v1;v2]);  %叉乘后向量方向
    c=dot(v1,v2)/(norm(v1)*norm(v2));
    ang=rad2deg(acos(c));
        
    if r<=0&(ang<ang_2|isnan(ang)==1)%concave
            ang(ang<ang_2)=[];
            im=[im;i];
    end
end
y2(im,:)=[];
y1=y2;
[n m]=size(y1);
end
%%
y2=[y2;y2(1,:)];%result's curve is closed
end