RK4 для систем: обобщение

Тот же RK4, что и для скаляра, — нужно лишь научить его складывать списки.

RK4 для систем — метод Рунге-Кутты 4-го порядка, в котором промежуточные наклоны k1, k2, k3, k4 являются векторами, а арифметика над ними выполняется покомпонентно.

В одном из прошлых разделов мы вывели классический RK4 для скалярного уравнения. Его суть — за один шаг h собрать четыре оценки наклона (в начале, дважды в середине и в конце интервала) и взять их взвешенное среднее с весами 1, 2, 2, 1. Сейчас мы покажем поразительный факт: чтобы применить тот же метод к системе из любого числа уравнений, не нужно менять ни одной формулы. Меняется только тип данных: вместо чисел в формулах живут векторы. А значит, всё, что нам потребуется добавить — две крошечные функции для арифметики над векторами.

k1..k4 как векторы

Напомним скалярные формулы RK4 и поставим рядом их векторный аналог. Слева Y — число, справа Y — список:

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_new = Y + (h/6)*(k1 + 2*k2 + 2*k3 + k4)

Формулы буква в букву те же. Но теперь k1, k2, k3, k4 — это векторы наклонов: F возвращает список, значит каждый k — список. Выражения вроде Y + (h/2)*k1 надо понимать покомпонентно: умножить каждый элемент k1 на h/2 и поэлементно прибавить к Y. Чистый Python не умеет складывать списки знаком +, поэтому мы заведём две вспомогательные функции, которые и составляют весь «секрет» обобщения:

vadd(a, b)   # поэлементная сумма двух векторов: [a0+b0, a1+b1, ...]
smul(s, v)   # умножение вектора на скаляр s: [s*v0, s*v1, ...]

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

Соберём полноценный движок rk4_system(F, t0, Y0, h, n). Он копирует начальный вектор (чтобы не испортить аргумент вызывающего), делает n шагов, на каждом считает четыре векторных наклона, комбинирует их через vadd и smul и накапливает результат. Проверять будем на осцилляторе [y' = v, v' = -y] со стартом Y0 = [1, 0] — это начальное условие соответствует точному решению y = cos(t), v = -sin(t). Заодно посчитаем «энергию» E = y^2 + v^2: для точного решения это cos^2 + sin^2 = 1 в любой момент, и хороший метод обязан держать её почти постоянной.

import math

# --- арифметика над векторами (списками) ---
def vadd(a, b):
    return [ai + bi for ai, bi in zip(a, b)]

def smul(s, v):
    return [s * vi for vi in v]

# --- универсальный RK4 для систем ---
def rk4_system(F, t0, Y0, h, n):
    t = t0
    Y = Y0[:]                       # копия, чтобы не портить аргумент
    rows = [(t, Y[:])]
    for _ in range(n):
        k1 = F(t, Y)
        k2 = F(t + h/2, vadd(Y, smul(h/2, k1)))
        k3 = F(t + h/2, vadd(Y, smul(h/2, k2)))
        k4 = F(t + h,   vadd(Y, smul(h,   k3)))
        # Y += (h/6)*(k1 + 2*k2 + 2*k3 + k4)
        s = vadd(vadd(k1, smul(2, k2)), vadd(smul(2, k3), k4))
        Y = vadd(Y, smul(h/6, s))
        t = t + h
        rows.append((t, Y[:]))
    return rows

# Осциллятор: Y = [y, v], F = [v, -y], решение y=cos t, v=-sin t
def F(t, Y):
    y, v = Y
    return [v, -y]

rows = rk4_system(F, 0.0, [1.0, 0.0], 0.1, 30)
for i in range(0, 31, 5):
    t, Y = rows[i]
    y, v = Y
    E = y*y + v*v
    print(f"t={t:4.1f}  y={y:8.5f}  cos={math.cos(t):8.5f}  E={E:.6f}")

Вывод:

t= 0.0  y= 1.00000  cos= 1.00000  E=1.000000
t= 0.5  y= 0.87758  cos= 0.87758  E=1.000000
t= 1.0  y= 0.54030  cos= 0.54030  E=1.000000
t= 1.5  y= 0.07074  cos= 0.07074  E=1.000000
t= 2.0  y=-0.41615  cos=-0.41615  E=1.000000
t= 2.5  y=-0.80114  cos=-0.80114  E=1.000000
t= 3.0  y=-0.98999  cos=-0.98999  E=1.000000

Два вывода из таблицы. Во-первых, столбцы y и cos совпадают до пятого знака — RK4 воспроизводит точное решение синуса с великолепной точностью при крупном шаге h = 0.1. Во-вторых, энергия E = y^2 + v^2 держится равной 1.000000 на всём интервале. Это не случайность: RK4 4-го порядка вносит настолько малую ошибку за шаг, что величины, которые физически должны сохраняться, почти не дрейфуют на коротких и средних дистанциях. Проверка сохраняющейся величины — лучший способ убедиться, что движок для системы написан верно: если энергия «поплыла», ищите ошибку в порядке компонент или в формуле комбинации.

Почему это работает в любой размерности

Заметьте, что в коде rk4_system нигде не написано, что вектор двумерный. vadd и smul работают со списком любой длины, F вызывается как чёрный ящик. Тот же самый движок без единого изменения проинтегрирует систему из трёх уравнений (как SIR-модель), из четырёх (точка на плоскости) или из шести (орбита в пространстве). Мы написали движок один раз — и дальше меняется только функция F. Это и есть выигрыш от перехода к векторному мышлению.

Почему компоненты обновляются синхронно

Здесь скрыта тонкость, в которую легко упереться. Возникает соблазн «упростить» код: пройтись по компонентам вектора по очереди и обновить каждую своим маленьким RK4, как если бы это были независимые скалярные уравнения. Так делать нельзя. В системе уравнения связаны: производная одной компоненты зависит от значений других (в осцилляторе y' зависит от v, а v' зависит от y). Если обновить сначала y, а потом считать наклон для v уже по новому y, мы смешаем два разных момента времени и разрушим связанность системы. Поэтому все наклоны k1, k2, k3, k4 вычисляются от одного и того же вектора состояния, и лишь в самом конце шага весь вектор Y обновляется одним движением. Синхронность — не стилистическая прихоть, а условие корректности: только так метод «видит» систему как единое целое, а не как набор разрозненных уравнений.

Эта же синхронность объясняет, почему движок так легко масштабируется. Вектор состояния может включать не два, а десятки и сотни величин: в химической кинетике каждая компонента — концентрация одного вещества в реакции, в многотельной задаче (звёзды, спутники) — координаты и скорости каждого тела. Размерность вектора растёт, но логика шага не меняется ни на букву: те же четыре оценки наклона, та же взвешенная комбинация. Именно благодаря тому, что vadd и smul работают над списком любой длины, переход от двух уравнений к двумстам не требует ни строчки новой математики — лишь более длинного вектора и более содержательной функции F.

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

Y = Y0 вместо Y = Y0[:]. Без копии вы будете изменять список вызывающего, и повторный запуск с тем же Y0 даст мусор. Всегда копируйте начальное состояние.
Складывание списков через +. В Python [1,2] + [3,4] даёт [1,2,3,4] (конкатенацию!), а не [4,6]. Поэлементную сумму делает только vadd. Это самая коварная ошибка обобщения.
Забыли множители 2 у k2 и k3. Веса в RK4 — это 1, 2, 2, 1. Потеряете двойку — порядок метода упадёт, и точность резко ухудшится.
Сохраняют ссылку на Y в истории. В rows нужно класть Y[:], копию. Иначе все строки истории будут указывать на один и тот же изменяемый список и покажут лишь финальное состояние.

Итоги:

RK4 для систем — те же формулы, но k1..k4 и Y становятся векторами.
Вся «магия обобщения» — две функции: vadd (поэлементная сумма) и smul (умножение на скаляр).
Движок rk4_system не зависит от размерности: меняется только F.
Проверка сохраняющейся величины (энергия y^2 + v^2 ≈ 1) — надёжный тест корректности.
Главная ловушка Python — + конкатенирует списки, а не складывает их поэлементно.