RK4 — рабочая лошадка: реализация
Теория позади — пора написать функцию, которой можно пользоваться каждый день. RK4 заслуженно называют рабочей лошадкой численного интегрирования: он прост, надёжен и точен. Соберём его на чистом Python, без numpy, и убедимся, что ошибка остаётся крошечной даже при крупном шаге.
Рабочая лошадка (workhorse) — в контексте ОДУ так называют классический RK4: метод, который при фиксированном шаге обеспечивает отличный баланс точности и простоты и потому служит выбором по умолчанию для подавляющего большинства гладких задач.
Реализация умещается в десяток строк. Заведём функцию rk4(f, t0, y0, h, n), которая делает n шагов размера h от точки (t0, y0) и возвращает два списка — узлы по времени и значения решения. Внутри — буквальный перевод формулы из прошлого урока в код.
Базовая реализация
def rk4(f, t0, y0, h, n):
"""n шагов метода RK4. Возвращает (списки t, списки y)."""
ts = [t0]
ys = [y0]
t, y = t0, y0
for _ in range(n):
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)
y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
t = t + h
ts.append(t)
ys.append(y)
return ts, ys
if __name__ == "__main__":
import math
# эталон: y' = y, точное решение e^t
ts, ys = rk4(lambda t, y: y, 0.0, 1.0, 0.5, 4)
print(f"{'t':>5} {'y_числ':>12} {'y_точн':>12} {'ошибка':>10}")
for t, y in zip(ts, ys):
exact = math.exp(t)
print(f"{t:5.2f} {y:12.6f} {exact:12.6f} {abs(exact - y):10.2e}")
Вывод:
t y_числ y_точн ошибка 0.00 1.000000 1.000000 0.00e+00 0.50 1.648438 1.648721 2.84e-04 1.00 2.717346 2.718282 9.36e-04 1.50 4.479375 4.481689 2.31e-03 2.00 7.383970 7.389056 5.09e-03
Шаг h = 0.5 — очень крупный, но к t = 2 ошибка всего около 5·10⁻³, то есть пять знаков после запятой ещё верны. Обратите внимание, как ошибка медленно растёт от шага к шагу: это накопление глобальной погрешности, неизбежное для любого метода, но у RK4 оно идёт от крошечной базы.
Нелинейная задача
Эталон y' = y подозрительно прост: правая часть линейна. Проверим RK4 на задаче, где решение ведёт себя интереснее: y' = t − y, y(0) = 1. Её точное решение — y(t) = t − 1 + 2e⁻ᵗ (можно проверить подстановкой). Здесь правая часть зависит и от t, и от y, поэтому метод работает «по-настоящему».
import math
def rk4(f, t0, y0, h, n):
ts, ys = [t0], [y0]
t, y = t0, y0
for _ in range(n):
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)
y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
t += h
ts.append(t)
ys.append(y)
return ts, ys
f = lambda t, y: t - y
exact = lambda t: t - 1 + 2 * math.exp(-t)
ts, ys = rk4(f, 0.0, 1.0, 0.5, 4)
print(f"{'t':>5} {'y_числ':>12} {'y_точн':>12} {'ошибка':>10}")
for t, y in zip(ts, ys):
print(f"{t:5.2f} {y:12.6f} {exact(t):12.6f} {abs(exact(t) - y):10.2e}")
Вывод:
t y_числ y_точн ошибка 0.00 1.000000 1.000000 0.00e+00 0.50 0.713542 0.713061 4.80e-04 1.00 0.736342 0.735759 5.83e-04 1.50 0.946791 0.946260 5.30e-04 2.00 1.271100 1.270671 4.29e-04
И снова четыре-пять верных знаков при шаге 0.5. Здесь ошибка даже не растёт монотонно: решение этой задачи затухает к прямой y = t − 1, и метод следует за ним устойчиво. RK4 одинаково хорош и на линейных, и на нелинейных гладких правых частях — менять в коде нужно только функцию f.
Сравнение порядков: Эйлер против RK4
Чтобы прочувствовать разницу в порядке, сравним Эйлер и RK4 при одном и том же шаге на эталоне y' = y. Разрыв в точности должен быть колоссальным.
import math
def euler(f, t0, y0, h, n):
t, y = t0, y0
for _ in range(n):
y = y + h * f(t, y)
t += h
return y
def rk4_final(f, t0, y0, h, n):
t, y = t0, y0
for _ in range(n):
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)
y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
t += h
return y
f = lambda t, y: y
h, n = 0.5, 4 # доходим до t = 2
ye = euler(f, 0.0, 1.0, h, n)
yr = rk4_final(f, 0.0, 1.0, h, n)
exact = math.exp(2.0)
print(f"точно = {exact:.6f}")
print(f"Эйлер = {ye:.6f} ошибка = {abs(exact - ye):.6f}")
print(f"RK4 = {yr:.6f} ошибка = {abs(exact - yr):.6f}")
print(f"во сколько раз RK4 точнее: {abs(exact - ye) / abs(exact - yr):.0f}")
Вывод:
точно = 7.389056 Эйлер = 5.062500 ошибка = 2.326556 RK4 = 7.383970 ошибка = 0.005086 во сколько раз RK4 точнее: 457
Как работает под капотом
Почему один и тот же шаг даёт разницу в сотни раз? Дело в порядке. Глобальная ошибка Эйлера ведёт себя как C·h, у RK4 — как C·h⁴. При h = 0.5 это h⁴ = 0.0625 против h = 0.5 — уже множитель восемь в пользу RK4 только за счёт степени, а константы C у RK4 к тому же меньше. Итоговая разница в нашем эксперименте — около 457 раз. Чтобы Эйлер дотянулся до точности, которую RK4 показал при h = 0.5, ему пришлось бы взять шаг примерно в 450 раз мельче, то есть сделать в 450 раз больше итераций. Даже с поправкой на то, что один шаг RK4 стоит четырёх вызовов f против одного у Эйлера, RK4 выигрывает по совокупной работе на порядки. Вот почему при фиксированном шаге его выбирают по умолчанию: на гладких задачах он почти всегда дешевле при равной точности.
Частые ошибки
- Обновлять t раньше наклонов. Все четыре наклона считаются от текущих t и y; t увеличивают только после вычисления y_next. Сдвиг t в начале итерации испортит аргументы k1.
- Переиспользовать изменённый y внутри шага. k2, k3, k4 опираются на исходное y_n плюс сдвиг, а не на промежуточно обновлённое y. Перезапись y до конца формулы ломает метод.
- Сравнивать методы при разном числе шагов. Чтобы сравнение порядков было честным, Эйлер и RK4 должны пройти один отрезок с одним h (и, значит, одинаковым n).
- Считать RK4 универсально лучшим. На жёстких (stiff) задачах явный RK4 теряет устойчивость и требует крошечного шага; там нужны неявные методы. RK4 — лошадка для гладких нежёстких задач.
Итоги
- RK4 реализуется десятком строк: четыре наклона и сборка y_next = y + (h/6)(k1 + 2k2 + 2k3 + k4).
- Достаточно менять функцию f — метод одинаково работает на линейных (y'=y) и нелинейных (y'=t−y) задачах.
- При крупном шаге h = 0.5 RK4 держит 4-5 верных знаков, тогда как Эйлер ошибается в первом знаке.
- На одном отрезке RK4 оказался точнее Эйлера примерно в 457 раз — это прямое следствие порядка h⁴ против h.
- RK4 — выбор по умолчанию для гладких нежёстких задач; для жёстких нужны неявные методы.