function [PI_Square, error_Square] = MMG(N)
% Monte Carlo Method to Get PI
close all
clc
format long
PI_Square = zeros(1,N);
error_Square = zeros(1,N);
for i = 1:N-1
[PI_Square(i), error_Square(i)] = getPI_Square(i, 0);
end
[PI_Square(N), error_Square(N)] = getPI_Square(N, 1);
plotAll(PI_Square, error_Square, N);
end
function [PI, error] = getPI_Square(nmax, isplot)
% 正方形内切一圆的模型
XY = 2*rand(2, nmax)-1;
XY_in = XY(:,( XY(1,:).^2 + XY(2,:).^2 ) <= 1);
XY_out = XY(:,( XY(1,:).^2 + XY(2,:).^2 ) > 1);
flag_inCircle = size(XY_in, 2);
PI = 4*flag_inCircle/nmax;
error = abs(PI - pi);
if isplot == 1
plotPoint(XY_in, XY_out);
end
end
function plotPoint(XY_in, XY_out)
line([-1, 1], [1, 1], 'LineWidth',2.5);
hold on
grid on
line([-1, 1], [-1, -1], 'LineWidth',2.5);
line([-1, -1], [-1, 1], 'LineWidth',2.5);
line([1, 1], [-1, 1], 'LineWidth',2.5);
alpha = 0:pi/20:2*pi;
R = 1;
x = R*cos(alpha);
y = R*sin(alpha);
plot(x, y, 'LineWidth',2.5)
axis equal
plot(XY_out(1,:), XY_out(2,:), '.k', 'MarkerSize', 10);
plot(XY_in(1,:), XY_in(2,:), '.m', 'MarkerSize', 10);
end
function plotAll(PI_Square, error_Square, N)
figure
PI_vector = pi * ones(1,N);
plot(1:N,PI_vector,1:N,PI_Square,'LineWidth',2);
grid on
legend('PI精确值', '蒙特卡洛计算结果');
xlabel('蒙特卡洛投点个数');
ylabel('PI的计算结果');
title('蒙特卡洛计算PI值图');
figure
plot(1:N,PI_vector,1:N,error_Square,'LineWidth',2);
grid on
legend('PI精确值', '蒙特卡洛计算误差');
xlabel('蒙特卡洛投点个数');
ylabel('PI的计算误差');
title('蒙特卡洛计算误差图');
end
评论0