突然关注到蒙特卡罗法,觉得挺有意思,在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
【参考文献】
1.Matlab程序学习(散点边界线/包络线)
https://blog.csdn.net/weixin_30467087/article/details/95560323