Visualizing the estimation of pai in Monte Carlo simulation

A visualization of the estimation of pai used Monte Carlo simulation

clear; clc; close all;

n = 50000;

N_total = 1000000;

update_rate_sim = 50;

update_rate_conv = 500;

fig = figure('Name', 'Monte Carlo Estimation of π (High-Speed Background Calculation)', 'NumberTitle', 'off', 'Color', 'w', 'Position', [100, 100, 1200, 600]);

ax1 = subplot(1, 2, 1);

axis(ax1, 'square');

axis(ax1, [-1, 1, -1, 1]);

hold(ax1, 'on');

grid(ax1, 'on'); box(ax1, 'on');

ax2 = subplot(1, 2, 2);

hold(ax2, 'on');

grid(ax2, 'on'); box(ax2, 'on');

xlabel(ax2, 'Number of Throws', 'FontSize', 12);

ylabel(ax2, 'π Estimate', 'FontSize', 12);

title(ax2, 'Convergence of π Estimate', 'FontSize', 14);

rectangle(ax1, 'Position', [-1, -1, 2, 2], 'EdgeColor', '#4DBEEE', 'LineWidth', 1.5);

viscircles(ax1, [0 0], 1, 'Color', 'r', 'LineWidth', 2);

title_handle = title(ax1, 'Preparing to start simulation...', 'FontSize', 14);

pi_line = yline(ax2, pi, '--r', 'True π Value', 'LineWidth', 2);

pi_line.LabelHorizontalAlignment = 'left';

pi_line.FontSize = 12;

axis(ax2, [0, n, 2.8, 3.6]);

points_inside_plot = plot(ax1, NaN, NaN, 'b.', 'MarkerSize', 8);

points_outside_plot = plot(ax1, NaN, NaN, 'g.', 'MarkerSize', 8);

legend(ax1, {'Circle', 'Points Inside', 'Points Outside'}, 'Location', 'northeastoutside', 'AutoUpdate', 'off');

convergence_plot = plot(ax2, NaN, NaN, '-b', 'LineWidth', 1.5);

points_in_circle = 0;

history_i = [];

history_pi = [];

pause(1);

for i = 1:n

x = 2 * rand() - 1;

y = 2 * rand() - 1;

distance_sq = x^2 + y^2;

if distance_sq <= 1

points_in_circle = points_in_circle + 1;

set(points_inside_plot, 'XData', [get(points_inside_plot, 'XData'), x], 'YData', [get(points_inside_plot, 'YData'), y]);

else

set(points_outside_plot, 'XData', [get(points_outside_plot, 'XData'), x], 'YData', [get(points_outside_plot, 'YData'), y]);

end

pi_estimate = 4 * points_in_circle / i;

if mod(i, update_rate_sim) == 0 || i == n

set(title_handle, 'String', sprintf('Points Thrown: %d / %d\nπ Estimate ≈ %.6f', i, n, pi_estimate));

drawnow limitrate;

end

if mod(i, update_rate_conv) == 0 || i == n

history_i = [history_i, i];

history_pi = [history_pi, pi_estimate];

set(convergence_plot, 'XData', history_i, 'YData', history_pi);

end

end

set(title_handle, 'String', sprintf('Visualization finished!\nCalculating up to %d points in background...', N_total));

drawnow;

n_remaining = N_total - n;

x_bg = 2 * rand(n_remaining, 1) - 1;

y_bg = 2 * rand(n_remaining, 1) - 1;

points_in_circle_bg = sum(x_bg.^2 + y_bg.^2 <= 1);

points_in_circle = points_in_circle + points_in_circle_bg;

pi_final_estimate = 4 * points_in_circle / N_total;

plot(ax2, N_total, pi_final_estimate, 'r*', 'MarkerSize', 12, 'LineWidth', 2, 'DisplayName', 'Final Estimate');

legend(ax2, {'Convergence Curve', 'True π Value', 'Final Value (Million Throws)'}, 'Location', 'best');

legend(ax2, 'boxoff');

final_title_str = sprintf('Final Result (%d Points Thrown)\nHigh-Precision π Estimate ≈ %.6f', N_total, pi_final_estimate);

set(title_handle, 'String', final_title_str);

hold(ax1, 'off');

hold(ax2, 'off');