Buffon's needle in MATLAB

问题:

I am currently working on a project for my Chemical Engineering class called Buffon's needle. The purpose of this project is to use MATLAB to get an estimate for pi and then to make a "cartoon" which will show the needles on a 10x10 graph with lines every 1 unit apart, with needles crossing the line being one color, and needles not crossing being another. I have found the pi estimate and i have created the graph, but my lines are not the one unit in length like they should be, instead the needles are all different lengths. if anyone could help me with this problem it would be much appreciated. my two scripts are below

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');

This below is my function 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];

回答1:

You made an indexing mistake. Your function returns endpoints:

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

Simplified,

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

These are collected in a 3D matrix, which you then plot:

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

Thus, here you are plotting a line with x-coordinates [start_x,end_x], and y-coordinates [start_x,start_y]! This latter should have been [start_y,end_y].

This should have been:

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

The same mistake happens when plotting endpointsNotHitList.