Немного аналитики для интуиции

Прежде чем доверять численному методу, полезно иметь несколько уравнений, которые мы умеем решать точно — они станут эталоном для проверки.

Аналитическое решение — это формула функции, удовлетворяющей уравнению, полученная точными преобразованиями, а не приближённым счётом.

Большинство дифференциальных уравнений из реальной жизни не имеют решения в виде формулы — именно поэтому существуют численные методы. Но есть несколько простых уравнений, которые поддаются ручному решению. Эти редкие случаи бесценны: они дают нам точный ответ-эталон, с которым можно сравнить результат численного метода и убедиться, что он работает правильно. Разберём самые важные из них.

Метод разделяющихся переменных

Самый доступный приём для ОДУ первого порядка — разделение переменных. Идея проста: если уравнение можно записать так, что слева стоит всё, что зависит от y, а справа — всё, что зависит от t, то обе части можно проинтегрировать по отдельности. Возьмём y' = k*y. Перепишем производную как отношение бесконечно малых: dy / dt = k*y. Перенесём y влево, dt вправо: dy / y = k * dt. Теперь левая часть зависит только от y, правая — только от t. Интегрируем обе стороны: слева получается натуральный логарифм ln(y), справа — k*t плюс константа. Распутав логарифм, приходим к экспоненте.

Экспоненциальный рост

Результат разделения переменных для y' = k*y при k > 0 — это экспоненциальный рост: y(t) = y0 * exp(k*t), где y0 — значение в момент t = 0. Эта формула описывает всё, что растёт пропорционально своему размеру: бактерии, проценты по вкладу, цепную реакцию. Характерная черта — кривая взлетает всё круче, ведь чем больше y, тем больше скорость y'. Проверим формулу численно: посчитаем значения для k = 0.5, y0 = 100.

import math
k = 0.5
y0 = 100.0
for t in [0, 1, 2, 3, 4]:
    print('t=%d  y=%.4f' % (t, y0 * math.exp(k * t)))

Вывод:

t=0  y=100.0000
t=1  y=164.8721
t=2  y=271.8282
t=3  y=448.1689
t=4  y=738.9056

За каждую единицу времени значение умножается примерно на 1.6487 (это exp(0.5)) — характерное постоянное отношение соседних значений и есть подпись экспоненты.

Экспоненциальное затухание

Если k < 0 (или пишут y' = -k*y с положительным k), та же формула описывает затухание: y(t) = y0 * exp(-k*t). Значение убывает, приближаясь к нулю, но никогда его не достигая. Так распадаются радиоактивные ядра, остывает тело (относительно температуры среды), расходуется реагент в реакции первого порядка. Полезное понятие здесь — период полураспада: время, за которое величина уменьшается вдвое; для экспоненты оно постоянно, независимо от того, с какого уровня считать. Это контринтуитивное, но очень важное свойство: упадёт ли величина со 1000 до 500 или со 100 до 50 — займёт одинаковое время. Именно постоянство периода полураспада позволяет датировать древние находки по углероду-14: доля оставшихся ядер однозначно связана с прошедшим временем через ту же формулу exp(-k*t). Заметьте также, что затухание — это просто рост с отрицательным k, так что отдельной теории здесь не нужно: один и тот же закон y' = k*y описывает оба сценария, меняется лишь знак коэффициента.

Логистический рост

Экспонента нереалистична на длинной дистанции: бактерии не могут расти бесконечно, ведь еда кончается. Логистическая модель добавляет торможение: y' = r*y*(1 - y/K), где K — ёмкость среды, предел численности. Пока y мал, множитель (1 - y/K) близок к единице, и рост почти экспоненциальный. Когда y приближается к K, множитель стремится к нулю, и рост останавливается. Это нелинейное уравнение тоже решается разделением переменных, а его решение — знаменитая S-образная кривая: y(t) = K / (1 + A*exp(-r*t)), где A = (K - y0)/y0. Посчитаем её для K = 1000, y0 = 10, r = 0.8.

import math
K = 1000.0
y0 = 10.0
r = 0.8
A = (K - y0) / y0
for t in [0, 2, 4, 6, 8, 10]:
    y = K / (1 + A * math.exp(-r * t))
    print('t=%2d  y=%.2f' % (t, y))

Вывод:

t= 0  y=10.00
t= 2  y=47.65
t= 4  y=198.59
t= 6  y=551.04
t= 8  y=858.74
t=10  y=967.86

Вот она, S-образная форма: сначала медленный старт (10 → 47), потом бурный взрывной рост в середине (198 → 551), и наконец насыщение у потолка K = 1000 (858 → 967). Эту кривую мы тоже будем воспроизводить численно.

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

Что значит «функция является решением уравнения»? Это значит, что если подставить функцию и её производную в уравнение, получится тождество. Проверим, что y(t) = y0*exp(k*t) действительно удовлетворяет y' = k*y. Производную посчитаем приближённо (центральной разностью) и сравним с k*y — они должны совпасть.

import math
k = 0.5
y0 = 100.0

def y(t):
    return y0 * math.exp(k * t)

h = 1e-6
for t in [0.0, 1.0, 2.0]:
    deriv = (y(t + h) - y(t - h)) / (2 * h)
    rhs = k * y(t)
    print('t=%.1f  y_prime=%.4f  k*y=%.4f  ok=%s' % (t, deriv, rhs, abs(deriv - rhs) < 1e-3))

Вывод:

t=0.0  y_prime=50.0000  k*y=50.0000  ok=True
t=1.0  y_prime=82.4361  k*y=82.4361  ok=True
t=2.0  y_prime=135.9141  k*y=135.9141  ok=True

Во всех точках численная производсная совпала с правой частью уравнения — формула действительно решение. Именно так мы и будем проверять численные методы: считаем точный эталон по формуле, прогоняем метод, сравниваем. Если численный ответ близок к эталону на задачах с известным решением, методу можно доверять и на задачах без формулы. Аналитическое решение — это линейка, которой мы измеряем точность приближения.

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

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

Разделение переменных решает простые ОДУ первого порядка, разводя y и t по разные стороны.
Экспоненциальный рост y = y0*exp(k*t) — у соседних значений постоянное отношение.
Затухание y = y0*exp(-k*t) убывает к нулю с постоянным периодом полураспада.
Логистика y' = r*y*(1 - y/K) даёт S-кривую с насыщением у ёмкости K.
Функция — решение, если подстановка её и производной обращает уравнение в тождество.
Аналитическое решение служит эталоном для проверки точности численных методов.