Решение ОДУ: метод Эйлера и scipy solve_ivp

Дифференциальное уравнение описывает, как система меняется во времени. Решим его численно — шаг за шагом методом Эйлера.

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

Зачем это нужно

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

Идея метода Эйлера

Производная — это наклон. Если мы знаем значение y сейчас и наклон y' = f(t, y), можно сделать маленький шаг h по времени, считая наклон постоянным:

y(t + h) ≈ y(t) + h * f(t, y)

Повторяем шаг за шагом — получаем траекторию. Это как идти в тумане: смотрим под ноги (текущий наклон), делаем шажок, снова смотрим.

Пример: радиоактивный распад

Уравнение y' = -k·y (скорость распада пропорциональна количеству). Точное решение — y(t) = y0·e^(−k·t). Сравним с Эйлером:

import math

k = 0.5
y = 100.0      # начальное количество
t = 0.0
h = 0.25       # шаг
steps = 8      # дойдём до t = 2

while t < 2.0 - 1e-9:
    y = y + h * (-k * y)   # шаг Эйлера
    t = t + h

exact = 100.0 * math.exp(-k * 2.0)
print("Эйлер при t=2 :", round(y, 4))
print("Точно         :", round(exact, 4))

Вывод:

Эйлер при t=2 : 34.3609
Точно         : 36.7879

Эйлер ошибается — шаг крупноват. Уменьшите h (и увеличьте число шагов) — приближение улучшится. Это фундаментальный компромисс: мелкий шаг = точнее, но дольше.

А в SciPy — одна строка

SciPy использует методы куда точнее Эйлера (Рунге-Кутты, адаптивный шаг):

import numpy as np
from scipy.integrate import solve_ivp

def rhs(t, y):
    return -0.5 * y          # y' = -k y

sol = solve_ivp(rhs, [0, 2], [100.0], t_eval=[2.0])
print(sol.y[0, -1])          # ~36.7879 — почти точно

solve_ivp по умолчанию применяет метод Рунге-Кутты 4(5) порядка с автоматическим подбором шага — отсюда высокая точность «из коробки».

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

Почему Эйлер так грубо ошибается? Он считает наклон постоянным на всём шаге, хотя на самом деле тот меняется. Ошибка за шаг ~ h², а накапливается за всё решение в ~ h (метод первого порядка). Методы Рунге-Кутты хитрее: они вычисляют наклон в нескольких промежуточных точках шага и усредняют. RK4 имеет ошибку ~ h⁴ — на порядки точнее при той же длине шага. Адаптивные версии ещё и сами уменьшают шаг там, где решение быстро меняется, и увеличивают на «спокойных» участках.

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

  • Слишком большой шаг. Эйлер при крупном h может не просто ошибаться, а вообще «разойтись» (для жёстких уравнений).
  • Жёсткие (stiff) уравнения. Когда в системе есть быстрые и медленные процессы, явные методы требуют крошечного шага; нужны специальные неявные решатели (method='Radau').
  • Забыть начальное условие. Без y(0) решение не определено — ОДУ задаёт лишь семейство кривых.

Итог

  • ОДУ задаёт скорость изменения; решение восстанавливает величину во времени.
  • Метод Эйлера: шаг по наклону, прост, но метод 1-го порядка — грубоват.
  • Методы Рунге-Кутты усредняют наклон и куда точнее; их даёт solve_ivp.
  • Компромисс шага: мельче — точнее, но дольше; адаптивные методы выбирают шаг сами.
Проверьте себя
1. Что вычисляет один шаг метода Эйлера y(t+h) ≈ y(t) + h·f(t,y)?
AТочное значение решения
BПриближение, экстраполируя по текущему наклону на шаг h вперёд
CИнтеграл функции f
DПроизводную y
2. Почему методы Рунге-Кутты (RK4) точнее метода Эйлера?
AОни используют меньше памяти
BОни вычисляют наклон в нескольких промежуточных точках шага и усредняют, давая ошибку ~h^4 вместо ~h
CОни не требуют начального условия
DОни работают только для линейных ОДУ
3. Что такое «жёсткое» (stiff) ОДУ?
AУравнение без решения
BСистема с одновременно быстрыми и медленными процессами, требующая неявных решателей или крошечного шага
CУравнение с целыми коэффициентами
DЛинейное уравнение