Численные методы для физиков
Этот ноутбук — рабочая часть первой лекции. Две физические задачи — случайное блуждание и оценка \(\pi\) методом Монте-Карло — имеют точный ответ, и на них видно, как устроены массивы NumPy, где теряется точность, как измерять время и когда помогает Numba.
Как запустить. Скачайте ноутбук: 01_python_numpy_numba.ipynb. Дальше один из двух путей:
python -m pip install numpy matplotlib numba jupyterlab, затем jupyter lab и открыть файл;Ячейки выполняются сверху вниз: Shift+Enter выполняет текущую ячейку. Время выполнения, напечатанное ниже, зависит от машины; на вашей оно будет другим.
import sys
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import numba
print("Python ", sys.version.split()[0])
print("NumPy ", np.__version__)
print("Matplotlib", matplotlib.__version__)
print("Numba ", numba.__version__)Python 3.12.15
NumPy 2.5.4
Matplotlib 3.11.2
Numba 0.68.0Версии печатаются в начале каждого ноутбука: численный результат воспроизводим, только если известно, чем он получен.
Сумма квадратов миллиона чисел — сначала циклом, потом одной операцией над массивом.
n = 1_000_000
values = np.arange(n, dtype=np.float64)
def sum_squares_loop(a):
total = 0.0
for v in a:
total += v * v
return total
%timeit -n 1 -r 3 sum_squares_loop(values)
%timeit -n 10 -r 3 np.sum(values * values)140 ms ± 174 μs per loop (mean ± std. dev. of 3 runs, 1 loop each)
585 μs ± 116 μs per loop (mean ± std. dev. of 3 runs, 10 loops each)Цикл Python на каждом шаге проверяет типы и создаёт объекты; NumPy выполняет тот же цикл в скомпилированном коде над непрерывным блоком памяти. Разница — два порядка.
a = np.arange(12).reshape(3, 4)
print(a)
print("shape:", a.shape, " dtype:", a.dtype)
print("строка 1: ", a[1])
print("столбец 2: ", a[:, 2])
print("каждый 2-й: ", a[0, ::2])[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
shape: (3, 4) dtype: int64
строка 1: [4 5 6 7]
столбец 2: [ 2 6 10]
каждый 2-й: [0 2]a = np.zeros(6)
b = a[::2] # вид на те же данные
b[0] = 100.0
print("после b[0] = 100: a =", a) # a тоже изменился
c = a[::2].copy() # независимая копия
c[1] = -1.0
print("после c[1] = -1: a =", a) # a не изменилсяпосле b[0] = 100: a = [100. 0. 0. 0. 0. 0.]
после c[1] = -1: a = [100. 0. 0. 0. 0. 0.]Массивы разной формы согласуются по правилу: размерности сравниваются справа, размерность 1 растягивается.
x = np.array([0.0, 1.0, 2.0]) # форма (3,)
y = np.array([[0.0], [10.0]]) # форма (2, 1)
print(x + y) # форма (2, 3)
pts = np.array([[0.0, 0.0], [3.0, 4.0], [6.0, 8.0]])
d = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=-1)
print(d) # все попарные расстояния без цикла[[ 0. 1. 2.]
[10. 11. 12.]]
[[ 0. 5. 10.]
[ 5. 0. 5.]
[10. 5. 0.]]print(0.1 + 0.2 == 0.3, 0.1 + 0.2)
eps = np.finfo(np.float64).eps
print("eps(float64) =", eps, " 1 + eps/2 == 1:", 1.0 + eps / 2 == 1.0)
print("eps(float32) =", np.finfo(np.float32).eps)False 0.30000000000000004
eps(float64) = 2.220446049250313e-16 1 + eps/2 == 1: True
eps(float32) = 1.1920929e-07Накопление ошибки: десять миллионов раз прибавляем \(0.1\) в одинарной точности. Точный ответ — \(10^6\).
tenths = np.full(10_000_000, 0.1, dtype=np.float32)
print("последовательно (cumsum):", np.cumsum(tenths)[-1])
print("np.sum (попарно): ", tenths.sum())
print("float64: ", tenths.astype(np.float64).sum())последовательно (cumsum): 1.087937e+06
np.sum (попарно): 1.0000001e+06
float64: 1000000.0149011612Последовательная сумма ошибается на 9%: когда сумма велика, прибавка \(0.1\) округляется. np.sum складывает попарно, и ошибка растёт медленнее.
Генератор создаётся один раз с явным seed; с тем же seed результат повторяется.
rng = np.random.default_rng(seed=2026)
print(rng.integers(0, 10, size=8))
print(np.random.default_rng(seed=2026).integers(0, 10, size=8)) # то же самое
print(rng.integers(0, 10, size=8)) # продолжение[8 1 0 6 3 4 0 3]
[8 1 0 6 3 4 0 3]
[6 3 8 7 7 9 7 1]Частица делает \(N\) шагов \(s_i=\pm 1\) с вероятностью \(1/2\). Положение \(x_N=\sum_{i=1}^N s_i\). Шаги независимы, \(\langle s_i\rangle=0\), \(\langle s_i^2\rangle=1\), поэтому \[ \langle x_N\rangle = 0, \qquad \langle x_N^2\rangle = N . \] Это эталон: проверим, что моделирование его воспроизводит.
rng = np.random.default_rng(seed=1)
n_walkers, n_steps = 2000, 1000
steps = rng.choice(np.array([-1, 1], dtype=np.int8), size=(n_walkers, n_steps))
x = np.cumsum(steps, axis=1, dtype=np.int64) # x[k, n-1] — положение k-го блуждания после n шагов
print("x.shape =", x.shape)
print("<x_N> =", x[:, -1].mean())
print("<x_N^2> =", (x[:, -1] ** 2).mean(), " (эталон", n_steps, ")")x.shape = (2000, 1000)
<x_N> = -0.347
<x_N^2> = 1032.122 (эталон 1000 )Насколько среднее по \(M\) блужданиям может отличаться от \(N\)? Для \(x_N^2\) дисперсия равна \(\langle x_N^4\rangle-N^2=2N^2-2N\), поэтому относительная статистическая ошибка \(\langle x_N^2\rangle/N\) близка к \(\sqrt{2/M}\).
n = np.arange(1, n_steps + 1)
ratio = (x ** 2).mean(axis=0) / n # broadcasting: (n_steps,) / (n_steps,)
band = np.sqrt(2 / n_walkers)
fig, ax = plt.subplots(figsize=(8, 3.2))
ax.plot(n, ratio, lw=1.2, label=r"моделирование, $M=2000$")
ax.axhline(1.0, color="k", lw=1, label="эталон")
ax.fill_between(n, 1 - band, 1 + band, color="0.8", label=r"$\pm\sqrt{2/M}$")
ax.set_xlabel("число шагов $n$")
ax.set_ylabel(r"$\langle x_n^2\rangle / n$")
ax.set_ylim(0.85, 1.15)
ax.legend(loc="upper right")
plt.show()
При больших \(N\) распределение \(x_N\) близко к гауссову с дисперсией \(N\). Тонкость: \(x_N\) имеет ту же чётность, что \(N\), поэтому возможные значения идут через 2 и вероятность каждого значения вдвое больше гауссовой плотности.
xN = x[:, -1]
values, counts = np.unique(xN, return_counts=True)
grid = np.linspace(-120, 120, 400)
gauss = 2 / np.sqrt(2 * np.pi * n_steps) * np.exp(-grid**2 / (2 * n_steps))
fig, ax = plt.subplots(figsize=(8, 3.2))
ax.bar(values, counts / n_walkers, width=1.6, alpha=0.6, label="моделирование")
ax.plot(grid, gauss, "k", lw=1.5, label=r"$2\,\mathcal{N}(0, N)$")
ax.set_xlabel("$x_N$")
ax.set_ylabel("вероятность")
ax.legend()
plt.show()
Бросаем \(n\) точек равномерно в квадрат \([0,1]^2\). Доля точек внутри четверти круга стремится к \(p=\pi/4\). Оценка \(\hat\pi=4\,k/n\) имеет стандартное отклонение \[ \sigma_{\hat\pi} = 4\sqrt{\frac{p(1-p)}{n}} \approx \frac{1.64}{\sqrt n}. \]
def pi_numpy(n, rng):
xy = rng.random((n, 2))
inside = (xy ** 2).sum(axis=1) < 1.0
return 4.0 * inside.mean()
rng = np.random.default_rng(seed=3)
for n in [10**2, 10**4, 10**6]:
est = pi_numpy(n, rng)
print(f"n = {n:>8}: pi ~ {est:.5f}, ошибка {est - np.pi:+.5f}, "
f"ожидаемая {1.64 / np.sqrt(n):.5f}")n = 100: pi ~ 3.20000, ошибка +0.05841, ожидаемая 0.16400
n = 10000: pi ~ 3.13760, ошибка -0.00399, ожидаемая 0.01640
n = 1000000: pi ~ 3.14054, ошибка -0.00105, ожидаемая 0.00164Проверим закон \(1/\sqrt n\): для каждого \(n\) повторим оценку 50 раз и найдём среднеквадратичное отклонение от точного \(\pi\).
rng = np.random.default_rng(seed=4)
ns = np.logspace(2, 6, 9).astype(int)
rms = [np.sqrt(np.mean([(pi_numpy(n, rng) - np.pi) ** 2 for _ in range(50)]))
for n in ns]
fig, ax = plt.subplots(figsize=(8, 3.2))
ax.loglog(ns, rms, "o", label="моделирование (50 повторов)")
ax.loglog(ns, 1.64 / np.sqrt(ns), "k", lw=1, label=r"$1.64/\sqrt{n}$")
ax.set_xlabel("число точек $n$")
ax.set_ylabel("ошибка")
ax.legend()
plt.show()
Чтобы получить ещё один верный знак, точек нужно в 100 раз больше. Для \(n=10^8\) массив координат занял бы 1.6 ГБ: векторизация NumPy упирается в память.
Numba компилирует функцию Python в машинный код при первом вызове. Цикл, записанный «как на бумаге», работает со скоростью C и не требует памяти под все точки сразу.
from numba import njit
@njit
def pi_numba(n, seed):
np.random.seed(seed) # у Numba свой генератор; seed задаётся внутри
hits = 0
for _ in range(n):
x = np.random.random()
y = np.random.random()
if x * x + y * y < 1.0:
hits += 1
return 4.0 * hits / nПервый вызов включает компиляцию, следующие — нет.
import time
def timed(f, *args):
t0 = time.perf_counter()
result = f(*args)
return result, time.perf_counter() - t0
_, t_first = timed(pi_numba, 10**6, 5)
_, t_second = timed(pi_numba, 10**6, 5)
_, t_numpy = timed(pi_numpy, 10**6, np.random.default_rng(5))
est, t_big = timed(pi_numba, 10**8, 5)
print(f"Numba, первый вызов (с компиляцией), n=1e6: {t_first:8.3f} с")
print(f"Numba, второй вызов, n=1e6: {t_second:8.3f} с")
print(f"NumPy, n=1e6: {t_numpy:8.3f} с")
print(f"Numba, n=1e8: pi ~ {est:.6f}, ошибка {est - np.pi:+.6f}, время {t_big:.2f} с")Numba, первый вызов (с компиляцией), n=1e6: 0.479 с
Numba, второй вызов, n=1e6: 0.020 с
NumPy, n=1e6: 0.029 с
Numba, n=1e8: pi ~ 3.141380, ошибка -0.000213, время 1.99 сКогда Numba не помогает: функция уже состоит из одной-двух операций NumPy (они и так скомпилированы), функция вызывается один раз и компиляция дороже счёта, код использует объекты, которых Numba не знает (словари произвольных типов, pandas, большинство функций SciPy).
float64, \(1.2\cdot10^{-7}\) для float32; порядок суммирования важен.seed.