«Школьный» пример. Вычисление числа \(\pi\)

Идея и реализация

Самый классический пример - вычисление числа \(\pi\) как площадь единичного круга.

Идея, лежащая в основе Монте-Карло подхода, на самом деле очень проста. Будем случайным образом выбирать точки из квадрата \((x, y) \in[-1,1] \times[-1,1]\). Тогда доля точек, попавших в круг, будет равна отношению площади круга к площади квадрата:

\[ \frac{\pi}{4}=\frac{S_{\circ}}{S_{\square}} \equiv \lim _{N \rightarrow \infty} \frac{\widetilde{\pi}_N}{4} \]

где введено приближенное значение числа \(\pi\) на N-м шаге.

\[ \widetilde{\pi}_N=4 \times \frac{\sum_{i=1}^N \theta\left(1-x_i^2-y_i^2\right)}{N} \]

Здесь \(\theta(x)=\left\{\begin{array}{ll}1, & x \geq 0 \\ 0, & x<0\end{array}-\right.\) стандартная функция Хевисайда.

### Визуализация на python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# set up the figure and axes
fig, ax = plt.subplots(figsize=(6, 6))
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)

# draw a unit circle
circle = plt.Circle((0, 0), 1, edgecolor='black', fill=False)
ax.add_patch(circle)

# initialize the scatter plot
scatter = ax.scatter([], [], s=10)  # use s=10 for original size

# generate the random points outside of the update function
n_max = 25000
points = np.random.uniform(-1, 1, size=(n_max, 2))
inside = np.linalg.norm(points, axis=1) <= 1

# define the Monte Carlo method
def monte_carlo(n, points, inside):
    num_inside = np.sum(inside[:n])
    if num_inside == 0:
        pi = 0
    else:
        pi = 4 * num_inside / n
    return pi

# define the update function for the animation
def update(frame, points, inside):
    n = frame * 1000
    pi = monte_carlo(n, points, inside)
    scatter.set_offsets(points[:n])
    scatter.set_color(np.where(inside[:n], 'blue', 'red'))
    scatter.set_sizes(np.full(n, 3))  # decrease size by a factor of 3
    ax.set_title(f"n = {n}, pi = {pi:.4f}")
    return scatter,

# create the animation
ani = animation.FuncAnimation(fig, update, fargs=(points, inside), frames=25, interval=200)

# display the animation
plt.close(ani._fig)
from IPython.display import HTML
HTML(ani.to_jshtml())

Once Loop Reflect