Системы ОДУ и жёсткие уравнения

Урок про то, как переносить методы на системы ОДУ, сводить к ним уравнения высших порядков и распознавать жёсткие задачи.

Система ОДУ — несколько связанных уравнений первого порядка y' = f(t, y), где y — вектор. Жёсткое уравнение содержит процессы очень разных скоростей, из-за чего явные методы требуют неоправданно малого шага.

Системы решаются тем же кодом

Замечательное свойство методов Эйлера и Рунге-Кутты: они почти без изменений работают с системами. Если y — не число, а вектор [y₀, y₁, ...], а f возвращает вектор наклонов, то все формулы (шаг по касательной, усреднение наклонов) применяются покомпонентно. Заменяем числовые операции на векторные — и тот же RK4 решает систему любого размера.

Это сразу даёт и решение уравнений высших порядков: их сводят к системе первого порядка введением новых переменных. Уравнение x'' = −x (гармонический осциллятор) превращается в систему x' = v, v' = −x с вектором состояния [x, v].

def rk4_система(f, t0, Y0, h, n):
    t, Y = t0, Y0[:]
    for _ in range(n):
        k1 = f(t, Y)
        k2 = f(t + h/2, [Y[i] + h/2*k1[i] for i in range(len(Y))])
        k3 = f(t + h/2, [Y[i] + h/2*k2[i] for i in range(len(Y))])
        k4 = f(t + h,   [Y[i] + h*k3[i]   for i in range(len(Y))])
        Y = [Y[i] + h/6*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in range(len(Y))]
        t = t + h
    return t, Y

import math
# x'' = -x  ->  система: x' = v, v' = -x. Старт x=1, v=0.
# Решение: x = cos(t). При t=pi должно быть x=-1, v=0.
f = lambda t, Y: [Y[1], -Y[0]]
t, Y = rk4_система(f, 0.0, [1.0, 0.0], math.pi/100, 100)
print(f"при t=π:  x = {Y[0]:.6f}  v = {Y[1]:.6f}")
print(f"точное :  x = {math.cos(math.pi):.6f}  v = {-math.sin(math.pi):.6f}")

Вывод:

при t=π:  x = -1.000000  v = -0.000000
точное :  x = -1.000000  v = -0.000000

Тот же RK4, лишь записанный для вектора, точно проинтегрировал колебания осциллятора: x(π) = −1, как и положено косинусу. Так решают и движение планет (системы Ньютона), и химическую кинетику, и электрические цепи — везде, где состояние описывается несколькими связанными величинами.

Что такое жёсткость

Некоторые системы содержат процессы очень разных временных масштабов: что-то затухает за наносекунды, что-то меняется за часы. Такие задачи называют жёсткими. Беда в том, что явный метод (Эйлер, RK4) вынужден держать шаг настолько мелким, чтобы «успевать» за самым быстрым процессом — даже когда тот давно затух и физически уже не важен. Иначе решение разносит (мы видели это на y'=−15y). В итоге явный метод делает миллионы ненужных шагов.

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

Формально жёсткость связана с собственными значениями матрицы Якоби системы: если их модули различаются на много порядков (большой разброс |λ|), задача жёсткая. Явные методы устойчивы лишь при h·|λ_max| в пределах небольшой области — отсюда требование микро-шага. Неявные методы (обратный Эйлер, BDF, Radau) A-устойчивы: их область устойчивости охватывает всю отрицательную полуплоскость, и шаг можно выбирать по точности, а не по устойчивости. Платой служит решение нелинейной системы на каждом шаге (Ньютоном) — но для жёстких задач это многократно окупается. В библиотеке выбор за вас делает scipy.integrate.solve_ivp: method='RK45' для обычных, method='BDF' или 'Radau' для жёстких.

Это лишь обзор: устойчивость, жёсткость и краевые задачи подробно разбираются в отдельном курсе по дифференциальным уравнениям. Здесь важно унести главное — как численно интегрировать ОДУ и когда явный метод не годится.

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

• Решать жёсткую задачу явным RK4. Он либо ползёт микро-шагом, либо взрывается; распознавайте жёсткость и берите неявный метод.
• Не сводить уравнение высшего порядка к системе. x''=... нельзя подать в RK4 напрямую — введите v=x' и решайте систему [x, v].
• Путать точность и устойчивость на системах. Метод может быть точным, но при жёсткости неустойчивым — это про устойчивость, а не про порядок.

Итоги

• Методы Эйлера и RK4 переносятся на системы покомпонентно — тот же код для векторного y.
• Уравнения высших порядков сводят к системе первого, вводя производные как новые переменные.
• Жёсткие задачи (разные масштабы времени) требуют неявных методов (BDF, Radau); явные тут неэффективны.
• Это обзор; детали устойчивости и краевых задач — в отдельном курсе по дифурам.