在MATLAB中布冯的针

Question

我目前正在为我的化学工程课程做一个名为布冯针的项目。这个项目的目的是使用MATLAB来估算pi，然后制作一个“卡通”，它将显示10x10图形上的针，每1个单位间隔一行，针穿过该线为一种颜色，针没有穿越另一个。我已经找到了pi估计并且我已经创建了图形，但是我的线条并不是它们应该的长度单位，而是针的长度都不同。如果有人能帮我解决这个问题，我将不胜感激。我的两个脚本在下面

clear all;  
close all;
clc;
format compact 

% Script to illustrate the estimation of pi value by using Buffon's needle
% experiment

% set number of separate experiments
nExperiments = 1000;

% set number of separate trials
nTrials = 3;

% total number of dropped needles is directly based on number of trials and
% number of experiments
ndropped=nTrials.*nExperiments;

% spacing between the parallel lines 
spacing = 1;

%  length of the needle 
L = spacing;

% the lower bound of x coordinate.
a = 10;

totalhits = 0;
for i = 1:nTrials
    %     keeps track of the number of hits
    hits = 0;

    %     keeps track of the number of times the needle doesn't hit the
    %     any of the lines
    nothits = 0;

    for j = 1:nExperiments
        [outcome,endpoints,angle] = needle_drop(spacing,L);

        if outcome
            hits = hits + 1;
            endpointsHitList(:,:,hits) = endpoints;
        else
            nothits = nothits + 1;
            endpointsNotHitList(:,:,nothits) = endpoints;
        end

        angleList(j) = angle;
    end

    scatter(1:nExperiments,angleList);
    xlabel('Experiments');
    ylabel('Angles');
    piestimate(i) = (2*L/spacing)/(hits/nExperiments);
end

fprintf('The average value of pi is %f plus or minus %f after %d trials of %d experiments with %d total number of dropped needle.\n',mean(piestimate),std(piestimate),nTrials,nExperiments,ndropped);

figure
hold on

% plot the vertical separations
for i = 0:spacing:a
    p1 =  plot([i,i],[0 11],'k','LineWidth',2);
end

% plot the needles that hit the vertical separation
for i = 1:hits
    p2 = plot(endpointsHitList(:,1,i),endpointsHitList(1,:,i),['-','b']);
end

% plot the needles that don't hit the vertical separation
for i = 1:nothits
    p3 = plot(endpointsNotHitList(:,1,i),endpointsNotHitList(1,:,i),['-','r']);
end

axis([-2,12 -2 12]);
legend([p1 p2 p3],'Vertical Separations','Hits','Not Hits')
title('Buffon Needle Experiment');
xlabel('x-axis');
ylabel('y-axis');
figure
histogram(piestimate)
title('Histogram of pi estimate');

以下是我的功能needle_drop：

function [outcome,endpoints,theta] = needle_drop(spacing,L)
% spacing = spacing between the parallel lines spacing
% L = length of the needle
% spacing = 1;
% % In the special case where the length of the needle is equal to the grid spacing
% % between the parallel lines
% L = spacing;

% a is the lower bound of x coordinate. b is the upper bound.
% the needle position will be randomly between 0 and 10.
a = 0;
b = 10;

% generate random number r1 in [0,1]
r1 = rand;

% sample a value of the angle uniformly distributed over the interval
% from zero to ninety degrees
theta = (pi/2)*r1;

% the projection of half the length of the needle horizontally: S
S = (L/2)*cos(theta);

% Another random number r2 is generated
% this corresponds to x,y coordinate
r2 = a + (b-a).*rand(1,2);

% we need to take care of the offset.
% if the x coordinate is between 0 and d then offset is 0 if xcord is
% between d and 2d then offset is d and likewise for other values.
offset = floor(r2(1));

% The sampled position T of the center of the needle is next compared to the
% sampled projection of half the length of the needle
if r2(1)-S <=0+offset || r2(1)+S >=spacing+offset
    outcome = 1;
else
    outcome = 0;
end

% the projection of half the length of the needle vertically: V
V = L/2*sin(theta);

endpoints = [r2(1)-S,r2(2)+V;r2(1)+S,r2(2)-V];

Answer 1

你犯了一个索引错误。你的函数返回endpoints：

endpoints = [ r2(1)-S, r2(2)+V; ...
              r2(1)+S, r2(2)-V ];

简化，

endpoints = [ start_x, start_y; ...
              end_x,   end_y   ];

这些是在3D矩阵中收集的，然后您可以绘制：

p2 = plot( endpointsHitList(:,1,i), endpointsHitList(1,:,i), ['-','b'] );
%          ^ x-coordinates          ^ y-coordinates

因此，在这里你用x坐标[start_x,end_x]和y坐标[start_x,start_y]绘制一条线！后者应该是[start_y,end_y]。

这应该是：

p2 = plot( endpointsHitList(:,1,i), endpointsHitList(:,2,i), ['-','b'] );
%                                                    ^^^ get second column

绘制endpointsNotHitList时会发生同样的错误。