Инструмент на практике: scipy solve_ivp

Зачем писать свой RK4, если в SciPy есть промышленный решатель с адаптивным шагом, выбором метода и детектором событий — и как им правильно пользоваться.

solve_ivp — универсальная функция из scipy.integrate для решения задач Коши (систем ОДУ с начальными условиями), поддерживающая адаптивный шаг, несколько методов, поиск событий и плотный вывод.

В учебных целях мы писали RK4 вручную — это бесценно для понимания. Но в реальной работе никто не пишет интеграторы с нуля. Стандарт де-факто в Python — функция solve_ivp из библиотеки SciPy. Она сама подбирает шаг под нужную точность, выбирает подходящий метод для жёстких задач и умеет ловить моменты особых событий. В этом уроке разберём её устройство и научимся выбирать параметры.

Сигнатура и основные параметры

Минимальный вызов выглядит так. Обратите внимание: SciPy не исполняется в браузере, поэтому весь scipy-код помечен как text и приведён для чтения, а не для запуска.

from scipy.integrate import solve_ivp
import numpy as np

# fun(t, y) возвращает производные; здесь маятник без трения
def fun(t, y):
    theta, omega = y
    return [omega, -9.81 * np.sin(theta)]

sol = solve_ivp(
    fun,                      # правая часть системы
    t_span=(0, 10),           # интервал интегрирования (t0, tf)
    y0=[1.0, 0.0],            # начальные условия
    method='RK45',            # метод интегрирования
    t_eval=np.linspace(0, 10, 100),  # в каких точках вернуть решение
    rtol=1e-6, atol=1e-9      # относительная и абсолютная точность
)

print(sol.t)   # массив моментов времени
print(sol.y)   # массив решений, форма (n_уравнений, n_точек)

Разберём параметры. fun(t, y) — правая часть: принимает время t и вектор состояния y, возвращает список или массив производных. Важно: в solve_ivp порядок аргументов именно (t, y), а не (y, t) как было в старой odeint. t_span — кортеж (начало, конец). y0 — начальный вектор. t_eval — точки, в которых вы хотите получить ответ (решатель всё равно сам выбирает внутренние шаги, а в эти точки интерполирует). rtol/atol — допуски точности: чем меньше, тем точнее и медленнее.

Выбор метода

Параметр method — самое важное решение. Доступны:

RK45 (по умолчанию) — Рунге-Кутта 4(5) Дорманда-Принса. Лучший выбор для большинства нежёстких задач.
RK23 — более низкого порядка, для грубой точности и быстрых прикидок.
DOP853 — Рунге-Кутта 8-го порядка, для очень высокой точности на гладких задачах.
Radau — неявный метод (Радо 5-го порядка) для жёстких систем.
BDF — метод обратного дифференцирования, классика для жёстких задач.
LSODA — автоматически переключается между нежёстким и жёстким режимами; удобен, когда жёсткость заранее неизвестна.

Правило выбора простое: нежёсткая задача — RK45; жёсткая (разномасштабные процессы, явные методы взрываются) — BDF или LSODA. Если не уверены, начните с RK45; если решатель делает миллионы крошечных шагов и тормозит — это признак жёсткости, переключайтесь на LSODA.

Как работает под капотом

В отличие от нашего RK4 с фиксированным шагом, solve_ivp использует адаптивный шаг. Он вычисляет решение двумя методами разного порядка (например, 4-го и 5-го), сравнивает их и оценивает локальную ошибку. Если ошибка укладывается в rtol/atol — шаг принимается и при запасе увеличивается; если нет — шаг уменьшается и пересчитывается. Так решатель сам тратит мелкие шаги там, где решение быстро меняется, и крупные — где оно гладкое. Вручную такое реализовать трудоёмко.

События (events)

Параметр events — мощная фишка. Вы передаёте функцию g(t, y), и решатель точно находит моменты, когда она пересекает ноль. Классический пример — падение тела: событие 'высота равна нулю' означает удар об землю.

def hit_ground(t, y):
    return y[0]          # y[0] = высота; ноль = земля
hit_ground.terminal = True   # остановить интегрирование при ударе
hit_ground.direction = -1    # ловить только когда высота убывает

sol = solve_ivp(fall, (0, 100), [100.0, 0.0], events=hit_ground)
print('Упал в момент t =', sol.t_events[0][0])

Флаг terminal=True останавливает счёт в момент события, а direction=-1 ловит только пересечения сверху вниз. SciPy не просто замечает событие на шаге — он интерполяцией находит его точное время. Реализовать это в ручном RK4 крайне сложно.

Плотный вывод (dense_output)

Если задать dense_output=True, решатель вернёт объект sol.sol — непрерывную функцию, которую можно вызывать в любой точке интервала, даже между шагами: y = sol.sol(3.7). Это удобно для гладких графиков и для случаев, когда нужные точки заранее неизвестны.

Свой мини-аналог на stdlib для контраста

Чтобы почувствовать, что solve_ivp делает 'внутри', напишем крошечный аналог с фиксированным шагом на чистом Python — он реально исполнится. Возьмём ту же задачу маятника.

import math

def fun(t, y):
    theta, omega = y
    return [omega, -9.81 * math.sin(theta)]

def rk4_step(f, t, state, h):
    k1 = f(t, state)
    s2 = [state[i] + 0.5*h*k1[i] for i in range(2)]
    k2 = f(t + 0.5*h, s2)
    s3 = [state[i] + 0.5*h*k2[i] for i in range(2)]
    k3 = f(t + 0.5*h, s3)
    s4 = [state[i] + h*k3[i] for i in range(2)]
    k4 = f(t + h, s4)
    return [state[i] + (h/6.0)*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in range(2)]

state = [1.0, 0.0]
h = 0.01
t = 0.0
print("  t      theta")
for step in range(1, 1001):
    state = rk4_step(fun, t, state, h)
    t += h
    if step % 200 == 0:
        print(str(round(t,1)).rjust(4), "  ", round(state[0],4))

Вывод:

  t      theta
 2.0    -0.5571
 4.0    -0.4239
 6.0    0.9966
 8.0    -0.1396
10.0    -0.8479

Маятник колеблется, угол theta меняет знак — это видно по чередованию положительных и отрицательных значений. Наш аналог делает то же, что RK45 внутри solve_ivp, но без адаптивного шага, без событий и без плотного вывода. Именно эти удобства и делают solve_ivp предпочтительным в реальной работе.

Частые ошибки

Порядок аргументов fun(y, t). В solve_ivp правильно fun(t, y). Перепутанный порядок — частая ошибка при переходе со старой odeint.
RK45 на жёсткой задаче. Если решатель невыносимо медленный и делает миллионы шагов, это жёсткость — переключайтесь на BDF/LSODA.
Путаница t_eval и внутреннего шага. t_eval лишь задаёт точки вывода; на точность и устойчивость влияют rtol/atol, а не t_eval.
Запуск scipy в браузерном раннере. SciPy здесь не исполняется — поэтому такой код помечается как text и приводится для чтения.

Зачем это нужно

solve_ivp — инструмент, которым реально пользуются инженеры и учёные. Адаптивный шаг экономит вычисления и гарантирует заданную точность, выбор метода покрывает и жёсткие задачи, а события и плотный вывод решают типичные практические потребности (поймать удар, нарисовать гладкую кривую). Понимая, как устроен RK4 изнутри, вы осознанно настраиваете готовый решатель, а не используете его как чёрный ящик. Это и есть зрелое владение численными методами.

solve_ivp — стандартный решатель ОДУ в SciPy с адаптивным шагом и выбором метода.
RK45 — для нежёстких задач, BDF/LSODA — для жёстких; LSODA переключается сам.
events точно находят моменты особых условий (например, удар об землю) с флагами terminal и direction.
dense_output даёт непрерывную функцию решения; rtol/atol управляют точностью, а не t_eval.