Адаптивный шаг: RKF45 и Дорман-Принс
До сих пор шаг h был фиксирован и выбирался вслепую. Но зачем мельчить там, где решение спокойное, и зачем грубить там, где оно резко меняется? Адаптивные методы сами подбирают шаг под локальную сложность задачи — и именно они работают внутри scipy.
Адаптивный шаг — стратегия, при которой метод на каждом шаге оценивает собственную локальную ошибку и автоматически подстраивает размер шага h: уменьшает его там, где ошибка велика, и увеличивает там, где есть запас точности.
Главная трудность очевидна: чтобы подстроить шаг, нужно знать ошибку — а точного решения у нас нет. Значит, ошибку надо оценить, не зная ответа. Все адаптивные методы Рунге-Кутты решают это одним трюком: считают шаг двумя способами разной точности и смотрят на расхождение.
Две оценки — одна ошибка
Представьте, что мы делаем шаг методом четвёртого порядка и получаем значение y⁴, а заодно — методом пятого порядка значение y⁵. Метод пятого порядка точнее, поэтому разность |y⁵ − y⁴| — это хорошая оценка локальной ошибки именно метода четвёртого порядка. Мы получили оценку ошибки бесплатно, не зная истинного решения: достаточно сравнить два приближения между собой.
Наивный способ сделать это — посчитать два совершенно независимых метода, удвоив работу. Гениальность Эрвина Фельберга (1969) в том, что он подобрал коэффициенты так, чтобы методы 4-го и 5-го порядков использовали одни и те же пробные наклоны. Это метод RKF45 (Рунге-Кутта-Фельберг): шесть наклонов, из которых собираются две оценки — порядка 4 и порядка 5. Их разность даёт ошибку почти даром, всего за шесть вызовов f вместо десяти.
Дорман-Принс — стандарт de-facto
В 1980 году Дорманд и Принс улучшили идею. Их метод DOPRI5 (Dormand-Prince) тоже семиступенчатый и даёт оценки 4-го и 5-го порядков, но коэффициенты подобраны так, чтобы минимизировать ошибку именно у решения 5-го порядка (его и берут как итог — приём FSAL, «first same as last», где последний наклон шага становится первым наклоном следующего). Именно Дорман-Принс лежит в основе метода RK45 в scipy и ode45 в MATLAB. Когда вы вызываете solve_ivp(..., method='RK45'), под капотом крутится этот алгоритм.
Как scipy выбирает шаг: atol и rtol
Контроллер шага сравнивает оценку ошибки err с допуском tol и корректирует h. Допуск в scipy задаётся двумя параметрами: atol (абсолютная точность) и rtol (относительная). Для каждой компоненты решения порог считается как tol = atol + rtol·|y|: когда y близко к нулю, главенствует atol; когда y велик, главенствует относительный rtol. Это разумно — требовать абсолютной точности 10⁻₆ от величины порядка миллиона бессмысленно.
Логика контроллера (упрощённо, как в scipy RK45):
err = норма( (y5 - y4) / (atol + rtol*|y|) )
if err <= 1:
шаг принят: y = y5, t = t + h
h_new = h * min(MAXFACTOR, SAFETY * err**(-1/5)) # увеличиваем
else:
шаг отклонён, не двигаемся
h_new = h * max(MINFACTOR, SAFETY * err**(-1/5)) # уменьшаем
SAFETY ~ 0.9 — коэффициент осторожности
степень 1/5 — потому что ведущий порядок ошибки ~ h^5
MAXFACTOR ~ 10, MINFACTOR ~ 0.2 — рост/спад шага ограниченыПоказатель −1/5 не случаен: локальная ошибка ведёт себя как h⁵, поэтому, чтобы из текущей ошибки err получить целевую ошибку 1, шаг масштабируют множителем err в степени −1/5. Множитель SAFETY чуть меньше единицы страхует от слишком оптимистичного увеличения, а ограничители не дают шагу скакать в десятки раз за итерацию.
Как работает под капотом
Реализовать RKF45 целиком в уроке тяжело из-за громоздких таблиц коэффициентов, поэтому построим контроллер на той же идее, но через правило Рунге (удвоения шага). Сделаем один шаг RK4 размера h и получим y_big. Затем пройдём тот же отрезок двумя полушагами h/2 и получим y_two. Поскольку RK4 имеет порядок 4, ошибка большого шага примерно в 2⁴ = 16 раз больше, чем ошибка двух мелких, и оценка локальной ошибки выражается как |y_two − y_big| / 15. Эта оценка играет роль той самой разности двух методов разного порядка — только «второй метод» мы получаем измельчением, а не другими весами.
Правило Рунге для метода порядка p = 4: y_big = один шаг RK4 размера h y_two = два шага RK4 размера h/2 err = |y_two - y_big| / (2^p - 1) = |y_two - y_big| / 15 err <= tol -> шаг принят, можно увеличить h err > tol -> шаг отклонён, уменьшаем h и пробуем снова
Соберём это в рабочий контроллер. Он печатает, как шаг сам подстраивается: на эталоне y' = y решение растёт экспоненциально, поэтому метод вынужден постепенно дробить шаг, чтобы удерживать ошибку под допуском.
def rk4_step(f, t, y, h):
k1 = f(t, y)
k2 = f(t + h/2, y + h/2 * k1)
k3 = f(t + h/2, y + h/2 * k2)
k4 = f(t + h, y + h * k3)
return y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
def adaptive(f, t0, y0, t_end, h, tol):
t, y = t0, y0
print(f"{'t':>8} {'h':>10} {'y':>12} {'оценка ошибки':>14}")
print(f"{t:8.4f} {h:10.6f} {y:12.6f} {'-':>14}")
while t < t_end - 1e-12:
if t + h > t_end:
h = t_end - t
y_big = rk4_step(f, t, y, h) # один шаг h
y_half = rk4_step(f, t, y, h/2) # первый полушаг
y_two = rk4_step(f, t + h/2, y_half, h/2) # второй полушаг
err = abs(y_two - y_big) / 15.0 # правило Рунге
if err <= tol: # шаг принят
t += h
y = y_two + (y_two - y_big) / 15.0 # уточнение по Ричардсону
print(f"{t:8.4f} {h:10.6f} {y:12.6f} {err:14.2e}")
if err < tol / 16: # большой запас -> растём
h *= 2.0
else: # ошибка велика -> дробим
h *= 0.5
adaptive(lambda t, y: y, 0.0, 1.0, 2.0, 0.5, 1e-6)
Вывод:
t h y оценка ошибки 0.0000 0.500000 1.000000 - 0.2500 0.250000 1.284025 5.27e-07 0.5000 0.250000 1.648721 6.76e-07 0.7500 0.250000 2.117000 8.68e-07 0.8750 0.125000 2.398875 3.42e-08 1.0000 0.125000 2.718281 3.88e-08 1.1250 0.125000 3.080216 4.40e-08 1.2500 0.125000 3.490342 4.98e-08 1.3750 0.125000 3.955076 5.65e-08 1.5000 0.125000 4.481688 6.40e-08 1.6250 0.125000 5.078418 7.25e-08 1.7500 0.125000 5.754602 8.21e-08 1.8750 0.125000 6.520818 9.31e-08 2.0000 0.125000 7.389055 1.05e-07
Посмотрите на колонку h. Стартовый шаг 0.5 сразу отклоняется (оценка ошибки выше допуска 10⁻₆) и дробится до 0.25. Ближе к t = 0.875, когда решение разогналось, запас исчерпывается и метод снова дробит шаг до 0.125, после чего держит его до конца. Никто не задавал эти значения руками — контроллер сам нащупал шаг, при котором локальная ошибка стабильно сидит около 10⁻₇-10⁻₈, чуть ниже допуска. Финальное значение 7.389055 совпадает с e² ≈ 7.389056 до шести знаков. Именно так, только с более хитрой оценкой ошибки от пары методов 4/5 порядка, работает scipy.
Частые ошибки
- Делить разность не на тот множитель. Для метода порядка p оценка ошибки — |y_two − y_big| / (2ᵖ − 1). Для RK4 (p = 4) это деление на 15, а не на 16 и не на 2.
- Принимать менее точное значение. Если шаг принят, итогом логично брать более точное y_two (а ещё лучше — уточнение по Ричардсону), а не грубое y_big по большому шагу.
- Путать atol и rtol. atol спасает у нуля, rtol — на больших значениях. Один atol при больших y приведёт к мучительному дроблению шага; один rtol около нуля — к потере точности.
- Не ограничивать рост и спад шага. Без множителей-ограничителей (вроде 0.2 и 10) шаг может скакнуть так резко, что следующий шаг тут же отклонится — метод зациклится на месте. Полезен и предохранитель на минимальный h.
- Считать, что adaptive всегда быстрее. На задачах с почти постоянной сложностью накладные расходы на оценку ошибки могут не окупиться; адаптивность выигрывает там, где решение неоднородно по гладкости.
Итоги
- Адаптивные методы оценивают локальную ошибку, сравнивая два приближения разного порядка, и подстраивают шаг под неё.
- RKF45 и Дорман-Принс делят пробные наклоны между методами 4-го и 5-го порядков, получая оценку ошибки почти даром.
- Дорман-Принс — основа scipy RK45 и MATLAB ode45; допуск задаётся как tol = atol + rtol·|y|.
- Новый шаг масштабируется множителем вида SAFETY·err⁻¹ˣ⁵ — степень 1/5 отражает ведущий порядок ошибки h⁵.
- Правило Рунге (шаг h против двух полушагов h/2) даёт оценку ошибки |y_two − y_big|/15 и позволяет построить адаптивный контроллер вручную, без scipy.