SciPy
NumPy сам по себе почти ничего не вычисляет: он хранит числа, уложенные подряд, и выполняет над ними арифметические операции.
SciPy является библиотекой численных методов, построенной поверх NumPy. NumPy отвечает за структуру данных, то есть за быстрый многомерный массив и базовые операции над ним, а SciPy добавляет готовые алгоритмы: интегрирование, решение дифференциальных уравнений, оптимизацию и подгонку кривых, интерполяцию, линейную алгебру, обработку сигналов, спектральный анализ. Прежде чем реализовывать численный метод самостоятельно, целесообразно поискать его в SciPy: почти наверняка он уже реализован поверх проверенных десятилетиями FORTRAN- и C-библиотек (QUADPACK, LAPACK, MINPACK, ODEPACK). Код при этом остаётся коротким, а выполняется со скоростью компилируемого.
SciPy организована в подмодули по областям: scipy.integrate, scipy.optimize, scipy.linalg и так далее. Импортировать библиотеку целиком не принято, импортируется нужный подмодуль.
import numpy as np # NumPy понадобится везде
from scipy import optimize # импортируем конкретный подмодуль
Листинги главы выполняются подряд, в одном сеансе интерпретатора: numpy импортирован как np один раз, а имена, введённые в одном листинге, доступны в следующих. Числа в комментариях и распечатках получены на виртуальной машине под Python 3.12 с NumPy 2.5 и SciPy 1.18, на которой записаны демонстрации лекции «NumPy, pandas и SciPy».
Слайды к главе. Физические константы, интегралы и задача Коши, поиск минимума и корня, подгонка кривых с погрешностями параметров на примере измерения эмиттанса и интерполяция изложены также в шестой части лекции «NumPy, pandas и SciPy» с демонстрациями в терминале; слайды лекции доступны на сайте книги и в PDF.
scipy.constants — физические константы
В scipy.constants собраны фундаментальные физические константы в системе СИ по данным CODATA, что избавляет от необходимости перепечатывать значения из справочника и исключает опечатки в восьмом знаке. SciPy 1.18 содержит рекомендованные значения CODATA 2022.
from scipy import constants
print(constants.c) # скорость света в вакууме, м/с
print(constants.h) # постоянная Планка, Дж·с
print(constants.hbar) # приведённая постоянная Планка, Дж·с
print(constants.e) # элементарный заряд, Кл
print(constants.k) # постоянная Больцмана, Дж/К
print(constants.epsilon_0) # электрическая постоянная, Ф/м
print(constants.m_e) # масса электрона, кг
Скорость света, а после переопределения единиц СИ в 2019 году также постоянная Планка, постоянная Больцмана и элементарный заряд точны по определению; точна и приведённая постоянная Планка, выражающаяся через постоянную Планка. Масса электрона и электрическая постоянная определяются измерениями и известны с погрешностью.
Помимо «именованных» констант имеется полный словарь physical_constants, возвращающий по названию кортеж (значение, единицы, погрешность); у точных констант погрешность равна нулю. Поиск нужного ключа удобно выполнять функцией find.
print(constants.find('Bohr')) # все константы, в названии которых есть 'Bohr'
value, unit, uncertainty = constants.physical_constants['proton mass']
print(value, unit, uncertainty) # 1.67262192595e-27 kg 5.2e-37
entry = constants.physical_constants['Planck constant']
print(entry) # (6.62607015e-34, 'J Hz^-1', 0.0)
Там же находятся множители, переводящие одни единицы в другие: constants.eV хранит электронвольт в джоулях, constants.angstrom — ангстрем в метрах, рядом префиксы constants.milli, constants.femto и т. п. Формулы записываются в СИ, а входные данные задаются в привычных единицах.
E = 13.6 * constants.eV # энергия ионизации водорода в джоулях
lam = constants.h * constants.c / E # длина волны соответствующего фотона, м
print(lam / constants.nano) # ≈ 91.2 нм — граница серии Лаймана
print(constants.convert_temperature(300, 'Kelvin', 'Celsius')) # 26.850000000000023
Последняя строка печатает 26.850000000000023, а не 26,85: десятичная дробь 273,15 представима в двоичной арифметике лишь приближённо, и разность 300 − 273,15 наследует ошибку её представления.
scipy.integrate — интегрирование и ОДУ
Подмодуль scipy.integrate вычисляет определённые интегралы и интегрирует обыкновенные дифференциальные уравнения.
Определённые интегралы вычисляет функция quad, реализующая адаптивную квадратуру из QUADPACK. Она принимает функцию и пределы, в том числе бесконечные, и возвращает пару из значения интеграла и оценки его абсолютной погрешности. Интервал интегрирования дробится, пока оценка не станет меньше абсолютного допуска epsabs или относительного допуска epsrel, умноженного на модуль интеграла; по умолчанию оба допуска равны \(1{,}49 \cdot 10^{-8}\). В качестве примера рассмотрим интеграл, возникающий при выводе закона Стефана—Больцмана из формулы Планка. Интегрирование спектральной плотности излучения по частотам после замены \(x = h\nu/kT\) приводит к интегралу
\[ I = \int_0^\infty \frac{x^3 dx}{e^x - 1} = \frac{\pi^4}{15}, \]
через который постоянная Стефана—Больцмана выражается как \(\sigma = 2\pi k^4 I/(h^3 c^2)\).
from scipy import integrate
# Интеграл ∫ x³/(eˣ − 1) dx от 0 до ∞; аналитический ответ — π⁴/15.
# Числитель и знаменатель умножены на e⁻ˣ, чтобы eˣ не переполнялся при больших x
f = lambda x: x**3 * np.exp(-x) / (1 - np.exp(-x))
value, error = integrate.quad(f, 0, np.inf)
print(value, np.pi**4 / 15) # 6.493939402266831 6.493939402266828
print(error) # 2.628470539040253e-09 — оценка погрешности
# постоянная Стефана—Больцмана, Вт/(м²·К⁴), и её значение в scipy.constants
sigma_sb = 2 * np.pi * constants.k**4 / (constants.h**3 * constants.c**2) * value
print(sigma_sb, constants.sigma) # 5.670374419184433e-08 5.6703744191844314e-08
Фактическая ошибка интеграла, около \(3 \cdot 10^{-15}\), почти на шесть порядков меньше оценки \(2{,}6 \cdot 10^{-9}\): оценка погрешности quad обычно пессимистична. Постоянная Стефана—Больцмана, вычисленная по интегралу, отличается от constants.sigma лишь в шестнадцатом знаке, на уровне округления float64; в SciPy она вычислена из точных констант по формуле \(2\pi^5 k^4/(15 h^3 c^2)\).
Множитель \(e^{-x}\) в числителе и знаменателе введён не случайно. Наибольшее число float64 близко к \(1{,}8 \cdot 10^{308}\), поэтому \(e^x\) переполняется уже при \(x > 709{,}78\), а на бесконечном интервале quad вычисляет подынтегральную функцию и в точках порядка тысяч. Без этой предосторожности результат зависит от вида подынтегральной функции.
naive = lambda x: x**3 / (np.exp(x) - 1) # исходная запись
print(integrate.quad(naive, 0, np.inf)) # (6.49393940226683, 2.628470028924825e-09)
debye = lambda x: x**4 * np.exp(x) / (np.exp(x) - 1)**2 # интеграл равен 4π⁴/15
print(integrate.quad(debye, 0, np.inf)) # (nan, nan)
В исходной записи NumPy возвращает при больших \(x\) бесконечность и печатает предупреждение RuntimeWarning: overflow encountered in exp, а результат остаётся верным лишь потому, что конечное число, делённое на бесконечность, даёт ноль. В интеграле \(\int_0^\infty x^4 e^x (e^x - 1)^{-2} dx = 4\pi^4/15\), входящем в закон Дебая для теплоёмкости кристалла при низких температурах, бесконечность возникает и в числителе, и в знаменателе. Отношение двух бесконечностей равно nan, и quad возвращает (nan, nan) с предупреждениями о переполнении и об ошибке округления, не позволяющей достичь заданной точности, тогда как та же функция в записи \(x^4 e^{-x}/(1 - e^{-x})^2\) даёт 25,97575760906732 при точном значении 25,975757609067312. Поэтому запись, устойчивая к переполнению, необходима и там, где исходная запись даёт верный ответ: верным он оказывается лишь при частном виде подынтегральной функции.
Экспериментальные данные обычно представлены таблицей отсчётов; для них предусмотрены integrate.simpson(y, x=x) и integrate.trapezoid(y, x=x), работающие непосредственно по узлам сетки.
Основным инструментом для ОДУ является solve_ivp (initial value problem, задача Коши). Он интегрирует систему уравнений первого порядка dy/dt = f(t, y); уравнение более высокого порядка предварительно сводится к системе. Для затухающего гармонического осциллятора x'' + 2γx' + ω₀²x = 0 введём вектор состояния y = [x, v] и получим систему первого порядка: x' = v, v' = −2γv − ω₀²x.
from scipy.integrate import solve_ivp
omega0 = 2 * np.pi # собственная частота, рад/с
gamma = 0.3 # коэффициент затухания, 1/с
def rhs(t, y, gamma, omega0):
x, v = y # y = [координата, скорость]
return [v, -2 * gamma * v - omega0**2 * x]
sol = solve_ivp(rhs, t_span=(0, 10), y0=[1.0, 0.0], # x(0)=1, v(0)=0
args=(gamma, omega0), # параметры модели
t_eval=np.linspace(0, 10, 501)) # где выдать решение
t = sol.t
x, v = sol.y # sol.y имеет форму (2, 501): строки — компоненты решения
print(sol.success, sol.nfev) # True 428: интегрирование дошло до конца
Параметры модели передаются в правую часть через аргумент args, а не читаются из глобальных переменных: solve_ivp вызывает rhs(t, y, gamma, omega0), что позволяет проинтегрировать ту же функцию с другими параметрами без изменения её кода. Способы передачи параметров в интеграторы рассмотрены в главе «Функции». Аргумент t_eval задаёт лишь моменты, в которые сохраняется решение: шаг интегрирования метод выбирает сам, а значения в заданные моменты получает интерполяцией внутри шага. Поле nfev содержит число вычислений правой части и позволяет сравнивать стоимость расчёта.
По умолчанию используется явный метод Рунге—Кутты RK45 пятого порядка, подбирающий шаг так, чтобы оценка локальной ошибки, полученная по вложенной формуле четвёртого порядка, оставалась меньше atol + rtol·|y|. Допуски по умолчанию, rtol=1e-3 и atol=1e-6, обеспечивают лишь около трёх верных знаков. Оценим фактическую ошибку. При \(\gamma < \omega_0\) уравнение осциллятора имеет точное решение \(x(t) = e^{-\gamma t}\left(\cos\omega t + \frac{\gamma}{\omega}\sin\omega t\right)\), где \(\omega = \sqrt{\omega_0^2 - \gamma^2}\), и с ним можно сравнить решения, полученные с разными допусками и методами.
w = np.sqrt(omega0**2 - gamma**2) # частота затухающих колебаний, рад/с
exact = np.exp(-gamma * t) * (np.cos(w * t) + gamma / w * np.sin(w * t))
for method, tol in (('RK45', {}), # допуски по умолчанию
('RK45', {'rtol': 1e-8, 'atol': 1e-10}),
('DOP853', {'rtol': 1e-10, 'atol': 1e-12})):
s = solve_ivp(rhs, (0, 10), [1.0, 0.0], method=method,
args=(gamma, omega0), t_eval=t, **tol)
print(method, tol, s.nfev, f'{np.abs(s.y[0] - exact).max():.1e}')
RK45 {} 428 2.2e-03
RK45 {'rtol': 1e-08, 'atol': 1e-10} 4796 1.0e-08
DOP853 {'rtol': 1e-10, 'atol': 1e-12} 2954 1.1e-10
С допусками по умолчанию метод вычисляет правую часть 428 раз, и наибольшая ошибка координаты составляет \(2{,}2 \cdot 10^{-3}\) при начальном отклонении 1. С rtol=1e-8 и atol=1e-10 ошибка уменьшается до \(1{,}0 \cdot 10^{-8}\) ценой 4796 вычислений. Метод восьмого порядка DOP853, рекомендуемый документацией для расчётов высокой точности, с ещё более строгими допусками достигает \(1{,}1 \cdot 10^{-10}\) за 2954 вычисления. Таким образом, для расчёта, сравниваемого с измерением, допуски необходимо задавать явно, а при высоких требованиях к точности метод высокого порядка оказывается выгоднее.
Для жёстких систем, где временные масштабы различаются на порядки (химическая кинетика, цепочки распадов), шаг явного метода ограничен условием устойчивости по самому быстрому масштабу, и расчёт становится медленным. В этом случае применяются неявные методы: method='Radau', неявный метод Рунге—Кутты пятого порядка, и method='BDF', формулы дифференцирования назад переменного порядка. Метод LSODA из ODEPACK сам определяет жёсткость задачи и переключается между методами Адамса и BDF.
scipy.optimize — минимизация, корни, подгонка кривых
В scipy.optimize решаются три типовые задачи обработки эксперимента: поиск минимума функции, поиск корня уравнения и подгонка модели к данным.
Минимизацией занимается функция minimize. Найдём равновесное межатомное расстояние в потенциале Леннард-Джонса V(r) = 4ε[(σ/r)¹² − (σ/r)⁶] (в безразмерных единицах ε = σ = 1).
from scipy import optimize
def lj(r):
return 4 * (r**-12 - r**-6) # потенциал Леннард-Джонса
res = optimize.minimize(lj, x0=1.5) # x0 — начальное приближение
print(res.x) # положение минимума: [1.12246209]
print(2**(1/6)) # аналитический ответ: r = 2^(1/6)
print(res.fun) # значение в минимуме: ≈ -1.0 (глубина ямы)
minimize работает и с функциями многих переменных, и тогда x0 становится массивом. Без ограничений по умолчанию применяется квазиньютоновский метод BFGS, при заданных границах переменных — L-BFGS-B, при ограничениях общего вида — SLSQP; для негладких функций имеет смысл попробовать method='Nelder-Mead', не использующий производных. Необходимо помнить, что найденный минимум является локальным: метод сходится к одному из минимумов, обычно ближайшему к x0, и глобального минимума не гарантирует.
Корни уравнения f(x) = 0 на отрезке, на концах которого функция имеет разные знаки, надёжно находит brentq. Метод Брента сочетает деление отрезка пополам с интерполяцией, и смена знака на концах обеспечивает его гарантированную сходимость. Выведем постоянную закона смещения Вина. Спектральная плотность излучения по длинам волн пропорциональна \(\lambda^{-5}/(e^{hc/\lambda kT} - 1)\); приравняв нулю её производную и введя \(x = hc/(\lambda kT)\), получим трансцендентное уравнение \((x - 5)e^x + 5 = 0\). У него есть тривиальный корень \(x = 0\), поэтому отрезок поиска выбран так, чтобы исключить его.
from scipy import constants
f = lambda x: (x - 5) * np.exp(x) + 5
x_max = optimize.brentq(f, 1, 10) # корень на отрезке [1, 10]: f(1) < 0 < f(10)
print(x_max) # 4.965114231744276
# Постоянная Вина b из λ_max·T = hc/(x_max·k)
b = constants.h * constants.c / (x_max * constants.k)
print(b, constants.Wien) # 0.0028977719551851727 0.0028977719551851727
Вычисленная постоянная Вина совпадает с constants.Wien во всех знаках: SciPy вычисляет её по той же формуле из точных констант. Тот же корень даёт более общий интерфейс optimize.root_scalar(f, bracket=(1, 10)), а для систем нелинейных уравнений предназначен optimize.root, принимающий вектор невязок.
Подгонка модели и погрешности параметров
Наиболее востребованной в лаборатории функцией является curve_fit, выполняющая подгонку параметров модели к данным методом наименьших квадратов. Сгенерируем точки экспоненциального распада с шумом и восстановим параметры.
rng = np.random.default_rng(42)
# «Эксперимент»: N(t) = N₀·exp(−t/τ) с N₀ = 10, τ = 1.5 плюс шум измерений
t = np.linspace(0, 5, 50)
data = 10 * np.exp(-t / 1.5) + rng.normal(0, 0.3, t.size)
def model(t, n0, tau): # первый аргумент — независимая переменная,
return n0 * np.exp(-t / tau) # остальные — подгоняемые параметры
popt, pcov = optimize.curve_fit(model, t, data, p0=(5, 1)) # p0 — стартовая точка
n0, tau = popt
perr = np.sqrt(np.diag(pcov)) # стандартные погрешности параметров
print(n0, tau) # ≈ 9.876 и ≈ 1.534 при истинных 10 и 1.5
print(perr) # [0.11241755 0.02603728]
curve_fit возвращает оптимальные параметры popt и ковариационную матрицу pcov: корни её диагональных элементов дают погрешности параметров для отчёта, а внедиагональные элементы — ковариации между параметрами. Восстановленные \(N_0 = 9{,}88 \pm 0{,}11\) и \(\tau = 1{,}534 \pm 0{,}026\) согласуются с истинными 10 и 1,5 в пределах полутора погрешностей. Задача без ограничений решается методом Левенберга—Марквардта из MINPACK, а при заданных границах параметров (bounds) — методом доверительной области 'trf'. Сами погрешности основаны на линейном приближении модели вблизи оптимума.
Без начального приближения p0 все параметры начинаются с единицы. В примере выше такое приближение тоже привело бы к верному ответу, однако значение по умолчанию не учитывает ни физического смысла параметров, ни выбранных единиц. Если записать те же моменты времени в наносекундах, время жизни станет равным 1500 нс и превысит стартовое значение на три порядка.
t_ns = t * 1e3 # те же моменты времени в наносекундах
popt_ns, pcov_ns = optimize.curve_fit(model, t_ns, data) # p0 = (1, 1)
print(popt_ns, np.sqrt(np.diag(pcov_ns))) # [10.09141512 1. ] [inf inf]
popt_ns, pcov_ns = optimize.curve_fit(model, t_ns, data, p0=(5, 1000))
print(popt_ns) # [ 9.87635832 1534.20618922]
При τ = 1 нс множитель \(e^{-t/\tau}\) уже во второй точке, t ≈ 102 нс, меньше \(10^{-44}\), а в первой точке, t = 0, от τ не зависит, поэтому производная модели по τ во всех точках практически равна нулю. Метод остаётся в стартовой точке по τ, подогнав лишь \(N_0\) к первому измерению, а SciPy выдаёт предупреждение OptimizeWarning: Covariance of the parameters could not be estimated и бесконечные погрешности. С приближением, выбранным по смыслу задачи, получается прежний ответ в новых единицах. Поэтому для нелинейной модели p0 задаётся явно.
Смысл погрешностей в pcov определяют аргументы sigma и absolute_sigma. Без sigma все точки считаются равноценными, а pcov масштабируется по разбросу остатков вокруг модели. С sigma минимизируется взвешенная сумма квадратов
\[ \chi^2 = \sum_{i=1}^{M} \left(\frac{y_i - f(x_i)}{\sigma_i}\right)^2, \]
однако по умолчанию, при absolute_sigma=False, учитываются лишь относительные значения sigma: ковариационная матрица умножается на приведённый \(\chi^2\), то есть на \(\chi^2/(M - N)\), где \(M\) — число точек, а \(N\) — число параметров. Погрешности тогда описывают разброс точек вокруг модели, а не точность прибора. При absolute_sigma=True значения sigma считаются стандартными ошибками измерений, и pcov выражает именно их, а приведённый \(\chi^2\), близкий к единице, подтверждает согласие модели и оценки ошибок с данными. Если ошибки измерений известны, они передаются через sigma= вместе с absolute_sigma=True. Разницу между режимами наглядно показывает измерение из физики ускорителей.
Подгонка настоящего эксперимента: измерение эмиттанса
На ускорителе измеряется эмиттанс пучка, площадь, занимаемая им в фазовой плоскости «координата — угол», делённая на π. Частицы пучка занимают в этой плоскости эллипс \(\gamma x^2 + 2\alpha x x' + \beta x'^2 = \varepsilon\) площадью \(\pi\varepsilon\); коэффициенты эллипса \(\alpha, \beta, \gamma\) называются параметрами Твисса и связаны соотношением \(\beta\gamma - \alpha^2 = 1\), а поперечный размер пучка равен \(\sigma = \sqrt{\varepsilon\beta}\). Эмиттанс определяет, до какого размера пучок может быть сфокусирован, и потому является основной мерой его качества.
Углы частиц напрямую не измеряются, измерить можно только поперечный размер пучка на люминофорном экране. Решение даёт метод квадрупольного сканирования: изменяя силу квадрупольной линзы перед экраном, снимают зависимость размера пучка от этой силы. В тонколинзовом приближении линза с оптической силой \(k = 1/f\) (нормализованный градиент, умноженный на длину линзы, в 1/м), а за ней промежуток длиной \(l\) до экрана дают матрицу перехода с элементами \(R_{11} = 1 - kl\) и \(R_{12} = l\), а квадрат размера на экране выражается через параметры Твисса \(\alpha_0, \beta_0, \gamma_0\) в месте линзы:
\[ \sigma^2(k) = \varepsilon\left(R_{11}^2\beta_0 - 2R_{11}R_{12}\alpha_0 + R_{12}^2\gamma_0\right), \qquad \gamma_0 = \frac{1+\alpha_0^2}{\beta_0}. \]
Это парабола по \(k\), и три её коэффициента однозначно определяют три искомые величины. Задача сводится к подгонке средствами curve_fit. Модельный скан ниже совпадает со сканом из демонстрации лекции: расстояние до экрана 1,2 м, 25 значений \(k\) от −1,5 до 1,5 1/м, погрешность измерения размера 20 мкм, истинные \(\varepsilon = 0{,}12\) мкм (0,12 мм·мрад), \(\beta_0 = 6{,}94\) м и \(\alpha_0 = -1{,}69\).
l = 1.2 # расстояние «линза → экран», м
def sigma2(k, eps, beta0, alpha0):
gamma0 = (1 + alpha0**2) / beta0
R11, R12 = 1 - k * l, l
return eps * (R11**2 * beta0 - 2 * R11 * R12 * alpha0 + R12**2 * gamma0)
true = (0.12e-6, 6.94, -1.69) # истинные ε (м·рад), β0 (м) и α0
rng = np.random.default_rng(42)
k = np.linspace(-1.5, 1.5, 25) # сила линзы в скане, 1/м
sigma_meas = 20e-6 # погрешность измерения размера, 20 мкм
meas = np.sqrt(sigma2(k, *true)) + rng.normal(0, sigma_meas, k.size) # размеры, м
Модель записана для \(\sigma^2\), поэтому кажется естественным возвести измерения в квадрат и подогнать параболу. Сравним такую подгонку с подгонкой, в которой измерениям возвращён их собственный вес.
# наивно: подгоняем квадраты
popt, pcov = optimize.curve_fit(sigma2, k, meas**2, p0=(1e-6, 5.0, -1.0))
# правильно: та же модель, но точки взвешены своими погрешностями;
# погрешность квадрата: d(sigma^2) = 2*sigma*d(sigma)
w = 2 * meas * sigma_meas
popt_w, pcov_w = optimize.curve_fit(sigma2, k, meas**2, p0=(1e-6, 5.0, -1.0),
sigma=w, absolute_sigma=True)
for name, p, c in (('без весов', popt, pcov), ('с весами', popt_w, pcov_w)):
e = np.sqrt(np.diag(c))
print(f'{name:10} eps = {p[0] * 1e6:.3f} ± {e[0] * 1e6:.3f} мкм '
f'beta0 = {p[1]:.2f} ± {e[1]:.2f} м alpha0 = {p[2]:.2f} ± {e[2]:.2f}')
без весов eps = 0.126 ± 0.045 мкм beta0 = 6.58 ± 2.35 м alpha0 = -1.61 ± 0.64
с весами eps = 0.119 ± 0.009 мкм beta0 = 6.97 ± 0.51 м alpha0 = -1.71 ± 0.13
Оба ответа верны в пределах своих погрешностей, однако погрешность эмиттанса различается в пять раз: 36 % против 7 %. Возведение в квадрат нелинейно, и ошибки оно растягивает неодинаково: погрешность квадрата \(2\sigma \cdot d\sigma\) велика у широкого пучка и мала в перетяжке, где пучок узкий. Подгонка без весов приписывает всем точкам одинаковую погрешность, поэтому результат определяют точки с большим размером пучка, а точки перетяжки почти не влияют на него. Между тем информацию об эмиттансе несёт как раз перетяжка, и наивная подгонка почти не учитывает самые ценные измерения.
Проверим согласие модели с данными и роль absolute_sigma: вычислим приведённый \(\chi^2\) и повторим взвешенную подгонку без absolute_sigma.
M, N = k.size, popt_w.size # число точек и параметров
chi2 = np.sum(((meas**2 - sigma2(k, *popt_w)) / w)**2)
print(f'chi2/(M - N) = {chi2 / (M - N):.2f}')
# тот же вызов без absolute_sigma: pcov умножается на chi2/(M - N)
_, pcov_rel = optimize.curve_fit(sigma2, k, meas**2, p0=(1e-6, 5.0, -1.0), sigma=w)
print(np.allclose(pcov_rel, pcov_w * chi2 / (M - N)))
print(f'погрешность eps без absolute_sigma: {np.sqrt(pcov_rel[0, 0]) * 1e6:.3f} мкм')
corr = pcov_w[0, 1] / np.sqrt(pcov_w[0, 0] * pcov_w[1, 1])
print(f'корреляция eps и beta0: {corr:.2f}')
chi2/(M - N) = 0.74
True
погрешность eps без absolute_sigma: 0.007 мкм
корреляция eps и beta0: -0.99
Приведённый \(\chi^2\) равен 0,74 и для 22 степеней свободы согласуется с единицей, поэтому заданная ошибка 20 мкм не противоречит данным. Без absolute_sigma=True погрешность эмиттанса была бы 0,007 мкм: curve_fit умножил бы ковариационную матрицу на 0,74, чтобы приведённый \(\chi^2\) стал равен единице, и погрешности всех параметров уменьшились бы в \(1/\sqrt{0{,}74} \approx 1{,}16\) раза, описывая разброс точек вокруг модели, а не заданную ошибку измерений. Коэффициент корреляции \(\varepsilon\) и \(\beta_0\), равный −0,99, показывает, что параметры определяются совместно: для величины, вычисляемой из нескольких параметров, необходима вся матрица pcov, а не только её диагональ.
Найдём, наконец, перетяжку, то есть наименьший размер пучка в скане. Минимум параболы по \(R_{11}\) равен \(\sigma_{\min}^2 = \varepsilon l^2/\beta_0\) и достигается при \(R_{11} = l\alpha_0/\beta_0\), то есть при \(k = 1/l - \alpha_0/\beta_0\): для истинных параметров \(\sigma_{\min} = l\sqrt{\varepsilon/\beta_0} = 0{,}158\) мм при \(k = 1{,}08\) 1/м. Сравним с перетяжкой, восстановленной двумя подгонками, на плотной сетке по \(k\).
kk = np.linspace(-1.5, 1.5, 301) # плотная сетка по k, 1/м
for name, p in (('без весов', popt), ('с весами', popt_w), ('модель', true)):
s = np.sqrt(sigma2(kk, *p)) * 1e3 # размер пучка, мм
print(f'{name:10} перетяжка {s.min():.3f} мм при k = {kk[s.argmin()]:.2f} 1/м')
без весов перетяжка 0.166 мм при k = 1.08 1/м
с весами перетяжка 0.157 мм при k = 1.08 1/м
модель перетяжка 0.158 мм при k = 1.08 1/м
Подгонка без весов хуже всего описывает именно точки у перетяжки, несущие информацию об эмиттансе. График такого скана с отрезками погрешностей и кривой модели построен в главе «Визуализация на Python».
Отсюда общее правило: преобразование данных, «выпрямляющее» модель, изменяет и веса точек. Тот же эффект даёт подгонка экспоненты через логарифмирование: log(y) превращает модель в прямую, однако делает равномерный шум неравномерным, и результат смещается. Если преобразование неизбежно, погрешности необходимо пересчитывать вместе с данными и передавать их в sigma=; другой путь состоит в подгонке модели в исходных величинах, например самого размера \(\sigma(k)\) с постоянной погрешностью 20 мкм.
На линейном ускорителе инжектора ЦКП «СКИФ» эмиттанс измеряется так же, только вместо тонколинзового приближения расчёт размера ведётся полной моделью ускорителя, а подгонка выполняется методом COBYQA (Constrained Optimization BY Quadratic Approximations), не требующим производных целевой функции, что удобно, когда каждое её вычисление представляет собой запуск программы моделирования. На выходе первой ускоряющей структуры измерен эмиттанс \(1{,}54 \pm 0{,}08\) мкм при энергии \(36{,}4 \pm 1{,}8\) МэВ, на выходе ускорителя \(0{,}12 \pm 0{,}01\) мкм при \(197{,}9 \pm 9{,}9\) МэВ. Каждый параметр приведён с погрешностью. Подгонка без погрешностей в физике эксперимента считается незаконченной работой, а curve_fit возвращает их в pcov.
scipy.interpolate — интерполяция табличных данных
Экспериментальные данные представлены таблицей, а значение требуется в произвольной точке между узлами. scipy.interpolate строит по таблице непрерывную функцию, например калибровочную кривую (термопара: напряжение → температура). Кусочно-линейную интерполяцию, соединяющую соседние узлы отрезками, выполняет функция NumPy np.interp, а кубический сплайн строит класс CubicSpline.
from scipy import interpolate
# Калибровка термопары: напряжение (мВ) -> температура (°C)
voltage = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
temp = np.array([0.0, 25.3, 50.1, 76.4, 101.8, 128.0])
spl = interpolate.CubicSpline(voltage, temp) # кубический сплайн
print(np.interp(2.5, voltage, temp), spl(2.5)) # 63.25 63.20875
print(spl(np.linspace(0, 5, 11))) # интерполянт принимает и массивы
print(spl.derivative()(2.5)) # 26.518055555555563, °C/мВ
Кусочно-линейная интерполяция не добавляет к данным ничего, кроме отрезков. Кубический сплайн состоит из кубических полиномов, сшитых в узлах с непрерывными первой и второй производными, и обеспечивает гладкую кривую. Метод derivative возвращает производную сплайна как новую функцию: наклон калибровочной кривой при 2,5 мВ равен 26,5 °C/мВ, что соответствует чувствительности термопары около 38 мкВ/°C.
print(np.interp(6.0, voltage, temp), spl(6.0)) # 128.0 158.0866666666666
print(interpolate.CubicSpline(voltage, temp, extrapolate=False)(6.0)) # nan
За пределами таблицы np.interp возвращает крайнее значение, а CubicSpline по умолчанию продолжает крайний полином: при 6 мВ получается 128,0 и 158,1 °C соответственно, и ни одна из функций не сообщает о выходе за таблицу. Экстраполяция опасна и требует осознанного решения; сплайн, построенный с extrapolate=False, возвращает за таблицей nan. Класс interp1d, встречающийся в старом коде, отнесён в документации SciPy к устаревшим (legacy); за пределами таблицы он по умолчанию возбуждает исключение ValueError. Для новых программ документация рекомендует np.interp, CubicSpline и make_interp_spline, причём последняя по умолчанию также экстраполирует.
Поскольку на каждом участке сплайн совпадает с полиномом, его можно аналитически дифференцировать и интегрировать; это позволяет, например, получить из таблицы координат скорость.
# Координата свободно падающего тела, измеренная в 11 моментах времени
t = np.linspace(0, 2, 11)
x = 0.5 * 9.81 * t**2
spl = interpolate.CubicSpline(t, x)
v = spl.derivative() # v(t) = dx/dt — тоже готовая функция
print(v(1.0)) # ≈ 9.81 — скорость через секунду падения
Интерполяция проводит кривую точно через все точки, включая испорченные шумом. Для зашумлённых данных применяется сглаживание (interpolate.make_smoothing_spline) или подгонка модели через curve_fit.
scipy.linalg — линейная алгебра
scipy.linalg покрывает и расширяет numpy.linalg: решение систем линейных уравнений, разложения матриц (LU, QR, SVD, Холецкого), собственные значения, матричные функции наподобие expm, считающей экспоненту от матрицы. Решим систему уравнений Кирхгофа для схемы с тремя контурами.
from scipy import linalg
# Метод контурных токов: R @ i = v
R = np.array([[ 11., -5., 0.],
[ -5., 18., -3.],
[ 0., -3., 9.]]) # матрица сопротивлений, Ом
v = np.array([10., 0., 0.]) # ЭДС в контурах, В
i = linalg.solve(R, v) # контурные токи, А
print(i)
print(R @ i) # проверка: получаем обратно v
linalg.solve(A, b) быстрее и численно устойчивее, чем linalg.inv(A) @ b; явное вычисление обратной матрицы почти никогда не требуется.
Второй задачей являются собственные значения, дающие в механике нормальные моды, а в квантовой механике уровни энергии. Рассмотрим два одинаковых груза массы m на пружинах жёсткости k, связанные пружиной kc; уравнения движения имеют вид ẍ = −(K/m)x, а квадраты частот нормальных мод совпадают с собственными значениями матрицы K/m.
m, k, kc = 1.0, 1.0, 0.5
K = np.array([[k + kc, -kc],
[-kc, k + kc]]) # матрица жёсткости связанных осцилляторов
evals, evecs = linalg.eigh(K / m) # eigh — для симметричных (эрмитовых) матриц
omega = np.sqrt(evals) # частоты нормальных мод
print(omega) # [1.0, 1.414...]: √(k/m) и √((k+2kc)/m)
print(evecs) # столбцы — моды: синфазная (1,1) и противофазная (1,-1)
Для симметричных и эрмитовых матриц лучше использовать eigh, а не общий eig: он быстрее, гарантированно возвращает вещественные собственные значения в порядке возрастания и ортонормированные собственные векторы. Гамильтонианы в матричном представлении относятся к этому случаю.
scipy.signal — обработка сигналов
scipy.signal очищает и анализирует данные с АЦП, осциллографа или детектора. Три типовые операции: фильтрация, оценка спектра, поиск пиков.
Смоделируем измерение, в котором полезный сигнал 5 Гц закрыт широкополосным шумом, и подавим шум низкочастотным фильтром Баттерворта.
from scipy import signal
rng = np.random.default_rng(0)
fs = 1000 # частота дискретизации, Гц
t = np.arange(0, 2, 1 / fs) # 2 секунды записи
clean = np.sin(2 * np.pi * 5 * t) # полезный сигнал 5 Гц
noisy = clean + 0.8 * rng.normal(size=t.size)
# НЧ-фильтр Баттерворта 4-го порядка, частота среза 20 Гц
b, a = signal.butter(4, 20, btype='low', fs=fs)
filtered = signal.filtfilt(b, a, noisy) # прогон вперёд и назад:
# нулевой фазовый сдвиг
filtfilt пропускает сигнал через фильтр дважды (вперёд и назад), поэтому не вносит фазовой задержки, и форма и положение пиков не искажаются.
Спектральную плотность мощности оценивают методом Уэлча, усредняющим спектры по перекрывающимся окнам, что снижает разброс оценки по сравнению с одиночным БПФ.
freqs, psd = signal.welch(noisy, fs=fs, nperseg=512)
print(freqs[np.argmax(psd)]) # ≈ 5.9 Гц — пик вблизи частоты сигнала
# (разрешение по частоте здесь fs/nperseg ≈ 2 Гц)
Пики (спектральные линии, импульсы детектора) находит find_peaks с настраиваемыми критериями (минимальная высота, расстояние между пиками, ширина).
peaks, props = signal.find_peaks(filtered, height=0.5, distance=100)
print(t[peaks]) # моменты максимумов синусоиды: шаг ≈ 0.2 с
print(props['peak_heights'])
scipy.fft — быстрое преобразование Фурье
Спектральный анализ напрямую выполняется через scipy.fft, заменивший старый scipy.fftpack. Для вещественных сигналов (а измеренные сигналы вещественны) предпочтительнее использовать rfft, возвращающий только физически осмысленную половину спектра, от нуля до частоты Найквиста fs/2.
Задача состоит в том, чтобы найти в зашумлённой записи скрытые периодические составляющие.
from scipy.fft import rfft, rfftfreq
rng = np.random.default_rng(1)
fs = 1000 # частота дискретизации, Гц
t = np.arange(0, 1, 1 / fs) # 1 секунда записи, N = 1000 отсчётов
x = (np.sin(2 * np.pi * 50 * t) # компонента 50 Гц, амплитуда 1
+ 0.5 * np.sin(2 * np.pi * 120 * t) # компонента 120 Гц, амплитуда 0.5
+ rng.normal(0, 0.5, t.size)) # белый шум
spectrum = rfft(x) # комплексный спектр
freqs = rfftfreq(t.size, 1 / fs) # сетка частот: 0 ... fs/2 = 500 Гц
amp = 2 / t.size * np.abs(spectrum) # нормировка на амплитудный спектр
print(freqs[np.argmax(amp)]) # 50.0 — главный пик
print(amp[freqs == 50], amp[freqs == 120]) # ≈ 1.0 и ≈ 0.5 — амплитуды нашлись
Разрешение по частоте равно 1/T, где T — длительность записи (здесь 1 Гц). Всё, что выше частоты Найквиста fs/2, «заворачивается» вниз (алиасинг), поэтому сигнал перед оцифровкой должен быть отфильтрован. Множитель 2/N преобразует модуль БПФ в амплитуды гармоник. Обратное преобразование выполняет irfft, для комплексных сигналов служит пара fft/ifft, имеются и многомерные версии (fft2, fftn), например для дифракционных картин.
Прочие возможности SciPy
Остальные подмодули:
scipy.statsсодержит распределения, статистические тесты, генерацию случайных величин;scipy.sparseхранит разреженные матрицы (гамильтонианы больших систем, сеточные задачи);scipy.specialсобрал специальные функции: Бесселя, Лежандра, эрмитовы полиномы, гамма-функцию;scipy.ndimageзанимается обработкой изображений (данные с ПЗС-камер);scipy.spatialпокрывает вычислительную геометрию: kd-деревья, триангуляции, ближайших соседей.
Если задача сформулирована в стандартных терминах (задача Коши, наименьшие квадраты, собственные значения), в SciPy почти наверняка найдётся готовая функция нужной точности и скорости. Таким образом, остаётся задать явно точность и смысл результата: допуски solve_ivp по умолчанию дают около трёх верных знаков, погрешности параметров curve_fit зависят от весов и аргумента absolute_sigma, а преобразование данных меняет веса точек.
Полезные ссылки
- Документация SciPy — справочник по всем подмодулям, у каждой функции есть примеры;
- Scientific Python Lectures (бывшие SciPy Lecture Notes) — подробный курс по научному Python: NumPy, SciPy, Matplotlib и продвинутые темы;
- SciPy Cookbook — рецепты для типовых задач: от подгонки данных до решения уравнений в частных производных.