Уравнения высших порядков как системы первого порядка

Универсальный приём: любую старшую производную можно «спрятать» в новую переменную.

Сведение к первому порядку — приём, при котором уравнение n-го порядка заменяется системой n уравнений первого порядка вводом новых переменных для производных y', y'', ..., y^(n-1).

Наши численные методы — Эйлер, RK4 — умеют решать только уравнения первого порядка: им нужна формула для производной состояния. Но физика полна уравнений второго порядка: второй закон Ньютона F = m*a связывает положение с ускорением, то есть со второй производной. Колебания, орбиты, маятники — всё это второй порядок. Как же их решать? Ответ — один простой и абсолютно общий приём, который превращает любое уравнение высшего порядка в систему первого порядка, после чего работает наш готовый rk4_system.

Приём: вводим переменные для производных

Идея до смешного проста. Пусть дано y'' = что-то. Объявим новую переменную для первой производной: v = y'. Тогда:

Первое уравнение системы — это определение новой переменной: y' = v.
Второе уравнение — это исходное уравнение, переписанное через v: v' = y'' = (правая часть).

Состояние было одним числом y, стало вектором [y, v]. Уравнение второго порядка стало системой из двух уравнений первого порядка — ровно то, что мы умеем интегрировать. Если бы был третий порядок, мы ввели бы ещё и a = v' и получили систему из трёх уравнений. Алгоритм механический и никогда не подводит.

Стоит оценить, насколько приём универсален. Уравнение любого порядка n сводится ровно к n уравнениям первого порядка: вводим по одной новой переменной на каждую производную вплоть до (n минус первой), и старшая производная попадает в правую часть последнего уравнения. Никаких исключений, никаких особых случаев — поэтому все библиотечные решатели ОДУ принимают на вход систему первого порядка, а сведение к ней считается стандартным первым шагом любой задачи. Один раз освоив этот приём, вы получаете доступ ко всей вычислительной машинерии для систем, не написав ни одного нового интегратора.

Особого внимания заслуживает вопрос: почему новой переменной берут именно скорость v = y'? Дело в физическом смысле. Пара (положение, скорость) — это полное мгновенное описание механической системы: зная, где тело и как быстро оно движется, мы однозначно предсказываем его будущее. Эта пара задаёт точку в так называемом фазовом пространстве (по сути, пространстве «положение плюс импульс»), а решение уравнения — траектория точки в нём. Вектор состояния [y, v] не случайная техническая уловка, а естественные координаты, в которых динамика выглядит как движение одной точки по плоскости.

Маятник

Математический маятник длиной L подчиняется уравнению theta'' = -(g/L)*sin(theta), где theta — угол отклонения, g — ускорение свободного падения. Это нелинейное уравнение второго порядка (из-за sin), элементарного решения у него нет. Сведём к системе. Введём угловую скорость omega = theta':

theta' = omega
omega' = -(g/L)*sin(theta)

Состояние Y = [theta, omega]. Для малых углов sin(theta) ≈ theta, и уравнение становится линейным с известным периодом T = 2*pi*sqrt(L/g) — это та самая формула из школьной физики. Проверим её численно: возьмём малый стартовый угол 0.2 рад и посмотрим, через сколько времени маятник вернётся в исходную точку.

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

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]

def rk4_system(F, t0, Y0, h, n):
    t, Y = t0, 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)))
        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

g, L = 9.8, 1.0

# theta'' = -(g/L)*sin(theta)  ->  [theta'=omega, omega'=-(g/L)*sin(theta)]
def F(t, Y):
    theta, omega = Y
    return [omega, -(g/L)*math.sin(theta)]

rows = rk4_system(F, 0.0, [0.2, 0.0], 0.05, 80)
print("  t     угол theta")
for i in range(0, 81, 8):
    t, Y = rows[i]
    theta, omega = Y
    print(f"{t:5.2f}   {theta:8.5f}")

T = 2*math.pi*math.sqrt(L/g)
print("Период малых колебаний 2*pi*sqrt(L/g) =", round(T, 4))

Вывод:

  t     угол theta
 0.00    0.20000
 0.40    0.06329
 0.80   -0.16005
 1.20   -0.16448
 1.60    0.05606
 2.00    0.19986
 2.40    0.07043
 2.80   -0.15540
 3.20   -0.16867
 3.60    0.04874
 4.00    0.19943
Период малых колебаний 2*pi*sqrt(L/g) = 2.0071

Маятник качается: угол идёт от +0.2 вниз через ноль к отрицательным значениям и обратно вверх. Посмотрите на момент t=2.00: theta=0.19986 — почти ровно стартовые 0.2! Маятник совершил полный цикл примерно за две секунды. И это в точности совпадает с теоретическим периодом T = 2*pi*sqrt(1/9.8) ≈ 2.0071 с. Численное решение нелинейного уравнения подтвердило школьную формулу для малых колебаний — отличный признак, что и приём сведения, и движок rk4_system работают верно. (При больших стартовых углах sin(theta) уже не равен theta, и реальный период станет заметно длиннее теоретического — это знаменитый эффект амплитудной зависимости периода маятника.)

Затухающий осциллятор

Тем же приёмом разбирается осциллятор с трением: y'' + 2*gamma*y' + w^2*y = 0, где gamma — коэффициент затухания, w — собственная частота. Выразим старшую производную: y'' = -2*gamma*y' - w^2*y. Введём v = y' и получим систему:

y' = v
v' = -2*gamma*v - w^2*y

В коде это просто другая функция F: def F(t, Y): y, v = Y; return [v, -2*gamma*v - w*w*y]. Член -2*gamma*v и есть сила трения: она всегда направлена против скорости и потому отбирает энергию. Решением будут затухающие колебания — амплитуда уменьшается от цикла к циклу, а на фазовом портрете вместо замкнутой петли получается спираль, закручивающаяся к началу координат. Движок не меняется ни на строку — мы лишь подставляем новое F. В этом и состоит сила векторного подхода: один универсальный rk4_system решает осциллятор, хищника с жертвой, маятник и затухающие колебания, отличаясь только правой частью.

Поведение этого осциллятора качественно меняется в зависимости от того, насколько сильно трение gamma по сравнению с собственной частотой w, и распадается на три классических режима. При слабом трении (gamma меньше w — недодемпфирование) система успевает много раз качнуться, прежде чем остановится: получаются колебания с медленно тающей амплитудой, а на фазовом портрете — спираль, виток за витком стягивающаяся к началу. При сильном трении (gamma больше w — передемпфирование) колебаний нет вовсе: выведенная из равновесия система вяло, без единого «перелёта» через ноль, сползает обратно к покою — так ведёт себя тяжёлая дверь в густой смазке. На границе этих режимов (gamma равно w — критическое демпфирование) система возвращается к равновесию максимально быстро и тоже без колебаний; именно этот режим закладывают в автомобильные амортизаторы и доводчики дверей, чтобы гасить толчок как можно скорее, не раскачивая систему. Поменять режим в нашем коде — значит просто подобрать другие числа gamma и w; вся динамика по-прежнему рождается из одной и той же функции F.

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

Забыли выразить старшую производную. Перед сведением уравнение надо разрешить относительно y'' (перенести всё остальное вправо). Из y'' + 2*gamma*y' + w^2*y = 0 получаем y'' = -2*gamma*y' - w^2*y — только тогда правая часть готова.
Перепутали определение и физику. Первое уравнение системы всегда тривиально (y' = v — это определение v), второе несёт всю физику (v' = правая часть). Если написать наоборот, система будет бессмысленной.
Линеаризацию приняли за точную модель. Период 2*pi*sqrt(L/g) верен лишь для малых углов. При старте с большого угла численное решение даст более длинный период — и это правильно, формула просто перестаёт применяться.
Неверный порядок [theta, omega]. Зафиксируйте, что Y[0]=theta, Y[1]=omega, и возвращайте производные в том же порядке. Перестановка превратит маятник в другую систему.

Итоги:

Любое уравнение n-го порядка сводится к системе n уравнений первого порядка вводом переменных для производных.
Первое уравнение системы — определение новой переменной (y' = v), последнее — исходное уравнение через эти переменные.
Маятник theta'' = -(g/L)*sin(theta) превращается в [theta' = omega, omega' = -(g/L)*sin(theta)].
Для малых углов численный период совпал с теоретическим 2*pi*sqrt(L/g) ≈ 2.007 с — проверка корректности.
Затухающий осциллятор y'' + 2*gamma*y' + w^2*y = 0 решается тем же движком; меняется только функция F.