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.
scipy.constants — физические константы
В scipy.constants собраны фундаментальные физические константы в системе СИ по данным CODATA, что избавляет от необходимости перепечатывать значения из справочника и исключает опечатки в восьмом знаке.
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) # масса электрона, кг
Помимо «именованных» констант имеется полный словарь physical_constants, возвращающий по названию кортеж (значение, единицы, погрешность). Поиск нужного ключа удобно выполнять функцией find.
print(constants.find('Bohr')) # все константы, в названии которых есть 'Bohr'
value, unit, uncertainty = constants.physical_constants['proton mass']
print(value, unit, uncertainty) # 1.672...e-27 kg и её погрешность
Там же находятся множители, переводящие одни единицы в другие: constants.eV хранит электронвольт в джоулях, constants.angstrom — ангстрем в метрах, рядом префиксы constants.milli, constants.femto и т. п. Формулы записываются в СИ, а входные данные задаются в привычных единицах.
E = 13.6 * constants.eV # энергия ионизации водорода в джоулях
lam = constants.h * constants.c / E # длина волны соответствующего фотона, м
print(lam / constants.nano) # ~91 нм — граница серии Лаймана
print(constants.convert_temperature(300, 'Kelvin', 'Celsius')) # 26.85
scipy.integrate — интегрирование и ОДУ
Подмодуль scipy.integrate вычисляет определённые интегралы и интегрирует обыкновенные дифференциальные уравнения.
Определённые интегралы вычисляет функция quad (адаптивная квадратура из QUADPACK), принимающая функцию и пределы (в том числе бесконечные) и возвращающая пару из значения интеграла и оценки погрешности. В качестве примера рассмотрим интеграл, возникающий при выводе закона Стефана—Больцмана из формулы Планка.
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.4939..., совпадает с точным значением
print(error) # оценка погрешности ~1e-9
Экспериментальные данные обычно представлены таблицей отсчётов; для них предусмотрены 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] и получим систему первого порядка.
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, 500)) # где выдать решение
t = sol.t
x, v = sol.y # sol.y имеет форму (2, 500): строки — компоненты решения
print(sol.success) # True, если интегрирование дошло до конца
Параметры модели передаются в правую часть через аргумент args, а не читаются из глобальных переменных: solve_ivp вызывает rhs(t, y, gamma, omega0), и ту же функцию можно проинтегрировать с другими параметрами без изменения её кода. Способы передачи параметров в интеграторы рассмотрены в главе «Функции».
По умолчанию используется явный метод Рунге—Кутты RK45, автоматически подбирающий шаг; точность регулируется аргументами rtol и atol. Для жёстких систем, где временные масштабы сильно различаются (химическая кинетика, цепочки распадов), явные методы медленны; в этом случае следует перейти на method='Radau' или method='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.1224...]
print(2**(1/6)) # аналитический ответ: r = 2^(1/6)
print(res.fun) # значение в минимуме: -1.0 (глубина ямы)
minimize работает и с функциями многих переменных, и тогда x0 становится массивом, по умолчанию применяется метод BFGS; для негладких функций имеет смысл попробовать method='Nelder-Mead'. Необходимо помнить, что найденный минимум является локальным и лежит вблизи x0.
Корни уравнения f(x) = 0 на отрезке, где функция меняет знак, надёжно находит brentq (метод Брента). Выведем постоянную закона смещения Вина. Положение максимума спектра Планка сводится к трансцендентному уравнению (x − 5)eˣ + 5 = 0:
from scipy import constants
f = lambda x: (x - 5) * np.exp(x) + 5
x_max = optimize.brentq(f, 1, 10) # корень на отрезке [1, 10], f меняет знак
print(x_max) # 4.9651...
# Постоянная Вина b из λ_max·T = hc/(x_max·k)
b = constants.h * constants.c / (x_max * constants.k)
print(b) # 2.898e-3 м·К — как в справочнике
Тот же результат даёт более общий интерфейс 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) # ≈ 10 и ≈ 1.5
print(perr) # погрешности подгонки
curve_fit возвращает оптимальные параметры popt и ковариационную матрицу pcov, а корни её диагональных элементов дают погрешности параметров для отчёта. Если для точек известны ошибки измерений, они передаются через sigma= (и absolute_sigma=True), тогда подгонка становится взвешенной. Для нелинейных моделей необходимо задавать начальное приближение p0, иначе метод может сойтись к другому минимуму или не сойтись вовсе.
Подгонка настоящего эксперимента: измерение эмиттанса
На ускорителе измеряется эмиттанс пучка, площадь, занимаемая им в фазовой плоскости «координата — угол». Эмиттанс определяет, до какого размера пучок может быть сфокусирован, и потому является основной мерой его качества.
Углы частиц напрямую не измеряются, измерить можно только поперечный размер пучка на люминофорном экране. Решение даёт метод квадрупольного сканирования: изменяя силу квадрупольной линзы перед экраном, снимают зависимость размера пучка от этой силы. В тонколинзовом приближении линза с нормализованным градиентом \(k\), а за ней промежуток длиной \(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.
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)
k = np.linspace(-1.5, 1.5, 25) # сканируем силу линзы
meas = ... # измеренные размеры пучка, м
sigma_meas = 20e-6 # погрешность измерения размера, 20 мкм
Модель записана для \(\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(lambda k, e, b, a: sigma2(k, e, b, a) / w,
k, meas**2 / w, p0=(1e-6, 5.0, -1.0))
Результат на модельных данных с истинными значениями \(\varepsilon = 0.12\) мкм, \(\beta_0 = 6.94\) м, \(\alpha_0 = -1.69\):
по sigma^2 без весов eps = 0.126 ± 0.045 мкм beta = 6.58 ± 2.35 м alpha = -1.61 ± 0.64
по sigma^2 с весами eps = 0.119 ± 0.007 мкм beta = 6.97 ± 0.44 м alpha = -1.71 ± 0.12
Оба ответа верны в пределах своих погрешностей, однако погрешность различается в шесть раз: 36 % против 6 %. Возведение в квадрат нелинейно, и ошибки оно растягивает неодинаково. Точка с большим размером пучка получает в \(\sigma^2\)-пространстве огромный вес, а точки в перетяжке, где пучок мал, получают почти нулевой. Между тем информацию об эмиттансе несёт как раз перетяжка. Наивная подгонка отбрасывает самые ценные измерения.
Правило: преобразование данных, «выпрямляющее» модель, изменяет и веса точек. Тот же эффект даёт подгонка экспоненты через логарифмирование. log(y) превращает модель в прямую, однако делает равномерный шум неравномерным, и результат смещается. Если преобразование неизбежно, погрешности необходимо пересчитывать вместе с данными и передавать их в sigma=.
На линейном ускорителе инжектора ЦКП «СКИФ» эмиттанс измеряется так же, только вместо тонколинзового приближения расчёт размера ведётся полной моделью ускорителя, а подгонка выполняется методом 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 строит по таблице непрерывную функцию, например калибровочную кривую (термопара: напряжение → температура).
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])
lin = interpolate.interp1d(voltage, temp) # кусочно-линейная интерполяция
spl = interpolate.CubicSpline(voltage, temp) # кубический сплайн
print(lin(2.5), spl(2.5)) # температура при 2.5 мВ по двум методам
print(spl(np.linspace(0, 5, 11))) # интерполянт принимает и массивы
Кусочно-линейная интерполяция не добавляет к данным ничего лишнего; кубический сплайн даёт гладкую кривую с непрерывными первой и второй производными. По умолчанию за пределами таблицы interp1d возбуждает исключение. Экстраполяция опасна, и включать её следует осознанно (fill_value='extrapolate').
Сплайн является полиномом на каждом участке, поэтому его можно аналитически дифференцировать и интегрировать. Таким образом из таблицы координат получают скорость.
# Координата свободно падающего тела, измеренная в 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 почти наверняка найдётся готовая функция нужной точности и скорости.
Полезные ссылки
- Документация SciPy — справочник по всем подмодулям, у каждой функции есть примеры;
- Scientific Python Lectures (бывшие SciPy Lecture Notes) — подробный курс по научному Python: NumPy, SciPy, Matplotlib и продвинутые темы;
- SciPy Cookbook — рецепты для типовых задач: от подгонки данных до решения уравнений в частных производных.