import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationN = 5000 # 总投点数（可调大，如 10000 或 100000）# ------------------------------# 生成随机点：x, y ∈ [-1, 1]# ------------------------------np.random.seed(42) # 固定随机种子便于复现x = np.random.uniform(-1, 1, N)y = np.random.uniform(-1, 1, N)# 判断是否在单位圆内：x² + y² ≤ 1inside = x**2 + y**2 <= 1pi_est = 4 * np.sum(inside) / N# ===== 动态动画（逐步显示点） =====fig, ax = plt.subplots(figsize=(8, 8))ax.set_xlim(-1.1, 1.1)ax.set_ylim(-1.1, 1.1)ax.set_aspect('equal')ax.grid(True, linestyle='--', alpha=0.5)ax.set_title(f'Monte Carlo π(N={N})', fontsize=14)ax.set_xlabel('x')ax.set_ylabel('y')# 画单位圆与正方形theta = np.linspace(0, 2*np.pi, 300)ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=2)ax.plot([-1, 1, 1, -1, -1], [-1, -1, 1, 1, -1], 'k--', linewidth=1)# 初始化散点scat_inside = ax.scatter([], [], s=10, c='blue', alpha=0.7, label='within circle')scat_outside = ax.scatter([], [], s=10, c='red', alpha=0.7, label='out of circle')ax.legend()# 当前点数统计text_pi = ax.text(-0.9, 0.9, '', fontsize=12, bbox=dict(boxstyle="round", facecolor="wheat"))# 动画更新函数def update(frame): # 每次加一批点（提升速度），这里每帧加 50 个点 step = 50 end = min(frame * step, N) idx = slice(0, end) x_in = x[idx][inside[idx]] y_in = y[idx][inside[idx]] x_out = x[idx][~inside[idx]] y_out = y[idx][~inside[idx]] scat_inside.set_offsets(np.column_stack([x_in, y_in])) scat_outside.set_offsets(np.column_stack([x_out, y_out])) pi_cur = 4 * len(x_in) / end if end > 0 else 0 text_pi.set_text(f'used points: {end}/{N}\nπ ≈ {pi_cur:.5f}\nerror: {abs(pi_cur - np.pi):.5f}') return scat_inside, scat_outside, text_piframes = N // 50 + 1ani = FuncAnimation(fig, update, frames=frames, interval=50, blit=False, repeat=False)plt.show()
