Немного аналитики для интуиции
Прежде чем доверять численному методу, полезно иметь несколько уравнений, которые мы умеем решать точно — они станут эталоном для проверки.
Аналитическое решение — это формула функции, удовлетворяющей уравнению, полученная точными преобразованиями, а не приближённым счётом.
Большинство дифференциальных уравнений из реальной жизни не имеют решения в виде формулы — именно поэтому существуют численные методы. Но есть несколько простых уравнений, которые поддаются ручному решению. Эти редкие случаи бесценны: они дают нам точный ответ-эталон, с которым можно сравнить результат численного метода и убедиться, что он работает правильно. Разберём самые важные из них.
Метод разделяющихся переменных
Самый доступный приём для ОДУ первого порядка — разделение переменных. Идея проста: если уравнение можно записать так, что слева стоит всё, что зависит от $y$, а справа — всё, что зависит от $t$, то обе части можно проинтегрировать по отдельности. Возьмём $y' = k\,y$. Перепишем производную как отношение бесконечно малых: $\frac{dy}{dt} = k\,y$. Перенесём $y$ влево, $\Delta t$ вправо: $\frac{dy}{y} = k\,dt$. Теперь левая часть зависит только от $y$, правая — только от $t$. Интегрируем обе стороны: слева получается натуральный логарифм $\ln(y)$, справа — $k\,t$ плюс константа. Распутав логарифм, приходим к экспоненте.
Экспоненциальный рост
Результат разделения переменных для $y' = k\,y$ при $k \gt 0$ — это экспоненциальный рост: $y(t) = y_0\,e^{k t}$, где $y_0$ — значение в момент $t = 0$. Эта формула описывает всё, что растёт пропорционально своему размеру: бактерии, проценты по вкладу, цепную реакцию. Характерная черта — кривая взлетает всё круче, ведь чем больше $y$, тем больше скорость $y'$. Проверим формулу численно: посчитаем значения для $k = 0.5$, $y_0 = 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$ (это $e^{0.5}$) — характерное постоянное отношение соседних значений и есть подпись экспоненты.
Экспоненциальное затухание
Если $k \lt 0$ (или пишут $y' = -k\,y$ с положительным $k$), та же формула описывает затухание: $y(t) = y_0\,e^{-k t}$. Значение убывает, приближаясь к нулю, но никогда его не достигая. Так распадаются радиоактивные ядра, остывает тело (относительно температуры среды), расходуется реагент в реакции первого порядка. Полезное понятие здесь — период полураспада: время, за которое величина уменьшается вдвое; для экспоненты оно постоянно, независимо от того, с какого уровня считать. Это контринтуитивное, но очень важное свойство: упадёт ли величина со 1000 до 500 или со 100 до 50 — займёт одинаковое время. Именно постоянство периода полураспада позволяет датировать древние находки по углероду-14: доля оставшихся ядер однозначно связана с прошедшим временем через ту же формулу $e^{-k t}$. Заметьте также, что затухание — это просто рост с отрицательным $k$, так что отдельной теории здесь не нужно: один и тот же закон $y' = k\,y$ описывает оба сценария, меняется лишь знак коэффициента.
Логистический рост
Экспонента нереалистична на длинной дистанции: бактерии не могут расти бесконечно, ведь еда кончается. Логистическая модель добавляет торможение: $y' = r\,y\,(1 - y/K)$, где $K$ — ёмкость среды, предел численности. Пока $y$ мал, множитель $(1 - y/K)$ близок к единице, и рост почти экспоненциальный. Когда $y$ приближается к $K$, множитель стремится к нулю, и рост останавливается. Это нелинейное уравнение тоже решается разделением переменных, а его решение — знаменитая S-образная кривая: $y(t) = \frac{K}{1 + A\,e^{-r t}}$, где $A = \frac{K - y_0}{y_0}$. Посчитаем её для $K = 1000$, $y_0 = 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) = y_0\,e^{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
Во всех точках численная производсная совпала с правой частью уравнения — формула действительно решение. Именно так мы и будем проверять численные методы: считаем точный эталон по формуле, прогоняем метод, сравниваем. Если численный ответ близок к эталону на задачах с известным решением, методу можно доверять и на задачах без формулы. Аналитическое решение — это линейка, которой мы измеряем точность приближения.
Частые ошибки
Первая ошибка — забыть начальное условие при разделении переменных и потерять константу интегрирования; именно она превращается в $y_0$ и определяет конкретную кривую из семейства. Вторая — путать экспоненциальный и логистический рост: на коротком интервале они выглядят одинаково (логистика стартует как экспонента), но логистика обязательно выходит на плато, а экспонента уходит в бесконечность. Третья — думать, что раз есть формула, численные методы не нужны. Наоборот: формула нужна именно как тест-эталон, а для подавляющего большинства реальных уравнений формулы не существует вовсе, и без численных методов не обойтись.
- Разделение переменных решает простые ОДУ первого порядка, разводя $y$ и $t$ по разные стороны.
- Экспоненциальный рост $y = y_0\,e^{k t}$ — у соседних значений постоянное отношение.
- Затухание $y = y_0\,e^{-k t}$ убывает к нулю с постоянным периодом полураспада.
- Логистика $y' = r\,y\,(1 - y/K)$ даёт S-кривую с насыщением у ёмкости $K$.
- Функция — решение, если подстановка её и производной обращает уравнение в тождество.
- Аналитическое решение служит эталоном для проверки точности численных методов.