Численное решение ОДУ: методы Эйлера, Хойна и Рунге–Кутта

Спутник вошёл в верхние слои атмосферы. Известны его координаты и скорость сейчас, известно уравнение движения — а аналитического решения нет и не будет: сопротивление атмосферы зависит от плотности, плотность — от высоты, высота — от решения. Спрашивать у природы формулу бесполезно. Спрашивать нужно по-другому: известна точка и наклон траектории в ней — куда мы попадём через малую долю секунды? Ответив на этот вопрос многократно, получим траекторию по кусочкам. Так работают все численные решатели дифференциальных уравнений — от калькулятора до аэродинамического пакета.

Задача Коши и идея шагов#

Задача Коши: дано уравнение первого порядка $y^{\prime }=f(x,y)$ и начальное условие $y(x_0)=y_0$. Требуется восстановить решение $y(x)$ на отрезке $[x_0,x_K]$. Существование и единственность гарантируются теоремой Пикара (при разумных условиях на $f$), но теорема молчит о том, как решение записать. Численный метод отвечает так: идём с постоянным шагом $h$, на каждом шаге превращаем наклон $f$ в приращение значения $\Delta y\approx h\cdot f$, и пересчитываем: $x_{k+1}=x_k+h$, $y_{k+1}=y_k+\Delta y$.

Методы различаются лишь одним: сколько раз и где спросить наклон $f$ внутри шага. Один запрос — Эйлер. Два — Хойн. Четыре — классический Рунге–Кутта. Больше запросов — дороже шаг, но точнее результат; вся классификация методов умещается в этот торг.

Метод Эйлера: ломаная и её цена#

Геометрию метода Эйлера видно без формул: на каждом отрезке ломаная идёт по касательной, проведённой в его левом конце. Касательная честна рядом с точкой касания — дальше начинает врать, и вранье копится. Точность первого порядка, $O(h)$: уменьшили шаг вдвое — ошибка упала примерно вдвое. У выпуклых решений ломаная Эйлера систематически отстаёт: она всё время следует за старым наклоном, тогда как кривая этот наклон меняет.

Ручной расчёт на этом же уравнении: один шаг $h=0,5$ даёт $y_1=1+0,5\cdot 1=1,5$ против точного $e^{0,5}\approx 1,6487$. Ошибка $0,149$ за один шаг — и так каждый шаг. Прежде чем осуждать метод, посмотрите на его достоинства: формула из одной строки, ноль лишних вычислений функции, и как первый метод в истории он научил человечество считать решения, которых не существует в формулах. Полная техника этого метода с ручными таблицами — в соседнем уроке метод Эйлера; здесь же идём дальше.

Улучшенный Эйлер (Хойн): предиктор-корректор#

Ошибка Эйлера — от устаревшего наклона. Логично сначала спрогнозировать конец шага грубо (методом Эйлера), посчитать наклон там, а затем пересчитать шаг по среднему из двух наклонов. Это схема Хойна (она же — улучшенный или модифицированный Эйлер):

Два вычисления функции на шаг вместо одного — и точность поднимается до $O(h^2)$: шаг вдвое меньше, ошибка примерно вчетверо меньше. На нашем шаге $h=0,5$: прогноз $\tilde{y}=1,5$, наклон на конце $f(0,5;1,5)=1,5$, итог $y_1=1+0,25(1+1,5)=1,625$. Ошибка $0,024$ вместо эйлеровской $0,149$ — в шесть раз лучше за вдвое большую работу.

1. 1
   Предиктор: $\tilde{y}=1+0,5\cdot 1=1,5$ — грубый эйлеровский прогноз на конец шага
2. 2
   Наклон на конце шага: $f(0,5;1,5)=1,5$
3. 3
   Корректор: $y_1=1+0,25(1+1,5)=1,625$ — шаг пересчитан по среднему наклону
4. 4
   Точное $e^{0,5}\approx 1,6487$: ошибка $0,024$ против эйлеровской $0,149$ — в шесть раз лучше

Классический РК4: четыре наклона#

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

Почему веса именно такие? Схема подбирается так, чтобы её разложение по степеням $h$ совпало с разложением точного решения вплоть до $h^4$ — уравнений на коэффициенты хватает как раз на четыре порядка. Для линейного уравнения $y^{\prime }=y$ шаг РК4 умножает значение на $1+h+\frac{h^2}{2}+\frac{h^3}{6}+\frac{h^4}{24}$ — обрезанный ряд экспоненты, и это лучшее, что можно сделать с четырьмя запросами наклона. Тот же вопрос «сколько точек брать», только для интегралов, разбирался в уроке про численное интегрирование: РК4 — ближайший родственник формулы Симпсона.

Кто точнее: сравнение на y′ = y#

Прогоним все три метода до точки $x=2$ при двух шагах. Точный ответ: $e^2=7,3891$. Уравнение $y^{\prime }=y$ растёт геометрически, поэтому за каждым методом закреплён свой множитель на шаг: у Эйлера $1+h$, у Хойна $1+h+\frac{h^2}{2}$, у РК4 — обрезанный ряд до $h^4$.

• Эйлер: ошибка упала в $1,6$ раза при удвоении числа шагов — первый порядок только разгоняется
• Хойн: падение в $3,3$ раза — второй порядок приближается к своему теоретическому $\times 4$
• РК4: падение в $13$ раз — четвёртый порядок идёт к своему пределу $\times 16$
• При малых $h$ соотношения выравниваются ровно до $2$, $4$ и $16$: так на практике проверяют порядок метода

Уточним язык, чтобы не путаться в отчётах: у метода порядка $p$ ошибка на одном шаге масштабируется как $h^{p+1}$, а накопленная на всём отрезке — как $h^p$. Причина проста: шагов $n=(b-a)/h$, и $n$ умножает локальную ошибку. Именно поэтому «порядок метода» в таблицах — это всегда порядок глобальной ошибки: локальная точность РК4 — пятый порядок, глобальная — четвёртый; обе цифры встречаются в литературе, и теперь вы знаете, какая где.

Читается таблица и по-другому: РК4 при грубом шаге $h=0,5$ точнее Эйлера при $h=0,25$ на три порядка — и при этом дешевле в пересчёте на единицу точности. Четыре вычисления функции против одного, но шаг можно брать в разы крупнее. На гладких задачах побеждает тот, у кого выше порядок, а не тот, кто чаще шагает.

Жёсткие уравнения: когда явные методы взрываются#

На гладкой задаче всё было прилично. Возьмём теперь $y^{\prime }=-15y$, $y(0)=1$ — решение $e^{-15t}$ стремительно умирает к нулю. Устойчивость явного Эйлера требует $|1+h\lambda |<1$ при $\lambda =-15$, то есть шаг $h<2/15\approx 0,13$. Возьмём $h=0,2$ — чуть крупнее: множитель шага $1-3=-2$, и значения идут $1,-2,4,-8,\dots$ — растут по модулю вдвое каждое шаг, хотя настоящее решение к моменту $t=1,6$ уже меньше $4\cdot 10^{-11}$.

Такие задачи с резко разными скоростями процессов называются жёсткими: быстрая составляющая (химическая реакция за миллисекунды) диктует крошечный шаг явному методу, хотя интерес к медленной (нагрев реактора за часы). Лечится переходом к неявным схемам: например, неявный Эйлер $y_{k+1}=y_k+hf(x_{k+1},y_{k+1})$ стабилен на любом шаге — ценой решения уравнения на каждом шаге. Для жёстких систем это плата обязательная, и все промышленные решатели её платят. Почему решения вообще «взрываются» по собственным значениям матрицы системы — отдельная красивая история про устойчивость по Ляпунову.

Практический стандарт для нежёстких задач — адаптивные схемы: пара методов разных порядков (обычно РК4 и РК5, схема Дормана–Принса) считается одновременно, и их разность даёт локальную оценку ошибки без обращения к точному решению — тот же дух, что правило Рунге в численном интегрировании. Шаг автоматически сжимается там, где решение крутое, и растёт на спокойных участках. Именно так работает ode45 в MATLAB и почти любой современный солвер по умолчанию.

Краевые задачи и метод стрельбы: краткий обзор#

В задаче Коши все условия даны в одной точке. В краевой задаче условия распределены по концам отрезка: например, $y^{\prime \prime }+y=0$, $y(0)=0$, $y(\pi /2)=1$. Прямые схемы РК4 к такому виду не приспособлены — им нужен полный стартовый набор $y(0)$ и $y^{\prime }(0)$. Метод стрельбы обходит это остроумно: угадаем недостающее начальное условие — пусть $y^{\prime }(0)=s$ — решим задачу Коши любым РК4 и посмотрим, куда попали в правом конце. Промахнёмся — подкорректируем $s$ и «выстрелим» снова.

В линейном примере зависимость от прицела прямая: решение с наклоном $s$ равно $s\sin x$, на правом конце $y(\pi /2)=s$, значит попасть в единицу даёт $s=1$ — и решение $y=\sin x$ найдено за один выстрел. Для нелинейных уравнений зависимость от $s$ кривая, и прицеливание ведут методом секущих или дихотомии. Альтернатива стрельбе — конечные разности: уравнение и краевые условия заменяют системой разностных уравнений, то есть большой СЛАУ — полный круг возвращается к уроку численные методы решения СЛАУ. Словарное имя этой конструкции — разностная схема.

Следующий конкретный шаг: вернитесь к симулятору в начале урока, включите $n=8$ и сверьте показание «Ошибка» с таблицей — сходится? Потом посчитайте руками один шаг Хойна для $y^{\prime }=x+y$ из соседнего урока и потренируйте тему в тренажёре по дифференциальным уравнениям. Оценку реальной точности численного результата стоит сверять с приёмами из урока про источники погрешностей — двойной расчёт с разными шагами есть частный случай правила Рунге.