Системы ОДУ и вектор состояния
Один объект, много чисел: как описать состояние, которое меняется по нескольким осям сразу.
Вектор состояния — упорядоченный набор чисел Y = [y1, y2, ..., yn], которого в каждый момент времени достаточно, чтобы полностью описать систему и однозначно предсказать её будущее.
Весь предыдущий курс мы решали уравнения вида y' = f(t, y), где y — одно число. Такая модель описывает мир, у которого ровно одна «степень свободы»: остывающий кофе (одна температура), счёт в банке (одна сумма), концентрация одного вещества. Но стоит сделать шаг в сторону реальности — и одной переменной перестаёт хватать. Чтобы описать движение точки на плоскости, нужны две координаты и две скорости. Чтобы описать эпидемию по модели SIR — три числа: доли здоровых, больных и переболевших. Чтобы описать качание маятника — угол и угловая скорость. Во всех этих случаях состояние системы — это не число, а вектор.
От одного уравнения к нескольким связанным
Рассмотрим простейший содержательный пример — гармонический осциллятор, то есть грузик на пружинке. Его положение y подчиняется уравнению y'' = -y (ускорение направлено к положению равновесия и пропорционально отклонению). Здесь стоит вторая производная, и напрямую наш аппарат для первого порядка к нему не применим. Но давайте введём вторую переменную — скорость v = y'. Тогда исходное уравнение превращается в пару связанных уравнений первого порядка:
y' = v v' = -y
Прочитаем их вслух: «скорость есть производная положения» (это просто определение v) и «производная скорости, то есть ускорение, равна минус положению» (это физика осциллятора). Два уравнения нельзя решать по отдельности: чтобы посчитать, как меняется y, нужно знать v, а чтобы посчитать, как меняется v, нужно знать y. Они сцеплены. Именно сцепленность — суть слова «система».
Вектор состояния и правая часть F(t, Y)
Соберём наши переменные в один список: Y = [y, v]. Это и есть вектор состояния. А правую часть системы оформим как функцию F(t, Y), которая по текущему состоянию возвращает вектор производных — список того, с какой скоростью меняется каждая компонента:
F(t, [y, v]) = [v, -y]
Обратите внимание на красоту этой записи. Внешне она устроена точно так же, как скалярное y' = f(t, y): слева производная состояния, справа функция от времени и состояния. Разница лишь в том, что и слева, и справа теперь стоят векторы одинаковой длины. Это значит, что любой наш численный метод, написанный в терминах «возьми f, посчитай наклон, шагни», автоматически обобщается на системы — достаточно научиться складывать векторы и умножать их на число. Мы займёмся этим в следующем уроке.
Фазовое пространство и фазовый портрет
У векторного взгляда есть мощный геометрический образ. Пространство, по осям которого отложены компоненты состояния, называют фазовым пространством. Для осциллятора это плоскость с осями (y, v): положение по горизонтали, скорость по вертикали. В каждый момент времени система — это одна точка на этой плоскости. Со временем точка движется, и след, который она оставляет, называется фазовой траекторией, а вся картина траекторий — фазовым портретом.
v
^
| . - .
| / \
----+-----------> y
| \ /
| ' - 'Для осциллятора фазовая траектория — это замкнутая кривая (окружность или эллипс): система вечно ходит по кругу, потому что энергия сохраняется. Замкнутость кривой на фазовом портрете — визуальный признак периодического движения. Если бы был трение, спираль закручивалась бы внутрь к точке равновесия. Фазовый портрет — это способ увидеть качественное поведение системы целиком, не привязываясь к конкретному моменту времени.
Как работает под капотом
Пока мы не интегрируем систему — мы лишь учимся её «щупать». Самая базовая операция — взять состояние Y в точке и вычислить вектор наклонов F(t, Y). Этот вектор говорит, в какую сторону и как быстро система собирается двигаться прямо сейчас. На чистом Python вектор — это обычный список, а F — обычная функция, возвращающая список той же длины. Никакого numpy: распаковываем список в переменные, считаем компоненты, собираем обратно в список.
import math
# Правая часть системы: гармонический осциллятор y'' = -y
# Состояние Y = [y, v], где v = y'
# Возвращаем вектор производных [y', v'] = [v, -y]
def F(t, Y):
y, v = Y # распаковываем вектор состояния
dy = v # y' = v
dv = -y # v' = -y
return [dy, dv] # вектор наклонов той же длины
# Пощупаем систему в нескольких точках фазовой плоскости
probes = [[1.0, 0.0], [0.0, 1.0], [0.7, 0.7], [-1.0, 0.0]]
for Y in probes:
slope = F(0.0, Y)
print(f"Y = {Y} -> F(t,Y) = {slope}")
Вывод:
Y = [1.0, 0.0] -> F(t,Y) = [0.0, -1.0] Y = [0.0, 1.0] -> F(t,Y) = [1.0, 0.0] Y = [0.7, 0.7] -> F(t,Y) = [0.7, -0.7] Y = [-1.0, 0.0] -> F(t,Y) = [0.0, 1.0]
Прочитаем результат геометрически. В точке [1, 0] (грузик максимально оттянут вправо, скорость нулевая) наклон [0, -1]: положение не меняется (мгновенно скорость ноль), а скорость начинает уменьшаться — грузик вот-вот пойдёт влево. В точке [0, 1] (грузик в центре, летит вверх с единичной скоростью) наклон [1, 0]: положение растёт, скорость пока неизменна. Если в каждой точке плоскости нарисовать такую стрелочку, получится векторное поле — «течение», вдоль которого движутся все траектории. RK4 в следующем уроке будет просто плыть по этому течению маленькими шагами.
Частые ошибки
- Перепутан порядок компонент. Если вы решили, что Y = [y, v], то и F обязана возвращать [y', v'] именно в этом порядке. Случайно вернуть [v', y'] — значит решать совсем другую систему. Зафиксируйте порядок один раз и держитесь его везде.
- F меняет входной список. Если внутри F писать в Y[0] = ..., вы испортите состояние, которое снаружи ещё понадобится методу. Правило: F только читает Y и возвращает новый список.
- Длина вектора «плавает». F всегда должна возвращать список ровно той же длины, что и Y. Вернёте на одну компоненту меньше — следующий шаг метода развалится.
- Путают фазовое пространство с графиком от времени. На фазовом портрете нет оси времени: по осям — переменные состояния. Замкнутая петля там означает периодичность, а не «вернулись в ту же секунду».
Итоги:
- Состояние системы — вектор Y = [y1, ..., yn]; его достаточно, чтобы предсказать будущее.
- Правая часть F(t, Y) возвращает вектор производных той же длины — это прямое обобщение скалярного y' = f(t, y).
- Уравнение второго порядка y'' = -y превращается в систему [y' = v, v' = -y] вводом скорости как отдельной переменной.
- Фазовое пространство — плоскость переменных; фазовый портрет рисует траектории; замкнутая кривая означает периодичность.
- F(t, Y) задаёт векторное поле — «течение», по которому численный метод поведёт решение.