Математика / Численные методы ОДУ

Двухстадийный метод Гаусса четвёртого порядка для ОДУ

Двухстадийный метод Гаусса решает задачу Коши для ОДУ с глобальным порядком четыре. На каждом шаге два наклона находят из связанной системы в узлах Гаусса — Лежандра и усредняют с равными весами.

Опубликовано: 9 ноября 2026 г.Обновлено: 9 ноября 2026 г.

Формула

\begin{gathered}c_{1,2}=\tfrac12\mp\tfrac{\sqrt3}{6},\quad A=\begin{pmatrix}\tfrac14&\tfrac14-\tfrac{\sqrt3}{6}\\[2pt]\tfrac14+\tfrac{\sqrt3}{6}&\tfrac14\end{pmatrix},\\ k_i=f\!\left(t_n+c_i h,\,y_n+h\sum_{j=1}^2 a_{ij}k_j\right),\quad i=1,2,\\ y_{n+1}=y_n+\tfrac h2(k_1+k_2) \end{gathered}

Схема Множитель затухания на отрицательной вещественной оси

Ось s = −hλ. Сплошная линия показывает R(−s), штриховая — e⁻ˢ. При s = 10 численная амплитуда остаётся заметной, хотя точная почти исчезла.

Обозначения

tₙ, h: текущее время и положительный шаг до tₙ₊₁, с, если t задано в секундах
yₙ, yₙ₊₁: текущее и следующее приближения решения, единица состояния y
$f(t,y)$: правая часть уравнения y′ = f(t,y), единица y в секунду
c₁, c₂: доли шага, задающие времена двух стадий, безразмерные
$A = (aᵢⱼ)$: матрица коэффициентов связи между стадиями, безразмерная
k₁, k₂: наклоны полинома коллокации в двух внутренних узлах, единица y в секунду
i, j: индексы строк и столбцов матрицы стадий, целые числа 1 и 2

Условия применения

Заданы начальные данные задачи Коши y(t₀) = y₀ и правая часть f. Для стандартной оценки четвёртого порядка f должна быть достаточно гладкой в окрестности решения, а шаги достаточно малы, чтобы выбирать ветвь стадий, близкую к текущему состоянию. При численном решении системы погрешность стадий контролируют отдельно, чтобы она не перекрыла погрешность схемы.

Ограничения

Для системы из d уравнений на шаге приходится находить 2d компонент стадий. Нелинейные стадии обычно требуют итераций. На линейном тесте y′ = λy схема A-устойчива, но не L-устойчива. При большом отрицательном hλ быстрое затухание передаётся плохо, даже когда численная последовательность не растёт. При разрывах правой части стандартная оценка четвёртого порядка через разрыв неприменима.

Подробное объяснение

Два числа c₁ и c₂ задают моменты внутри шага, симметричные относительно его середины. В этих моментах нужны наклоны k₁ и k₂. Каждая формула для наклона содержит оба неизвестных, поэтому стадии решают совместно. После этого новый уровень получается из среднего наклонов, умноженного на длину шага.

За формулой стоит квадратичный полином u(t), который начинается в известной точке yₙ. Его производная линейна, так что её можно задать двумя значениями k₁ и k₂. Условие коллокации приравнивает эти наклоны к f в двух внутренних точках полинома. Интегрирование линейного интерполянта производной до каждой из этих точек даёт строки A, а интегрирование до конца шага — два веса по 1/2.

Узлы Гаусса — Лежандра выбраны так, что двухузловая квадратура точна для полиномов до третьей степени. Для общей задачи зависимость f от y связывает стадии, и порядок четыре следует из условий порядка метода Гаусса. При точном входе и точно найденных стадиях ошибка одного шага имеет порядок h⁵. На фиксированном интервале при стандартных условиях глобальная ошибка имеет порядок h⁴.

Для y′ = λy решение системы стадий даёт множитель R(z) = (12 + 6z + z²)/(12 − 6z + z²), где z = hλ. При Re(z) ≤ 0 его модуль не больше единицы. Если z вещественно и уходит в минус бесконечность, R(z) стремится к единице, тогда как eᶻ стремится к нулю. Поэтому допустимый по устойчивости шаг может плохо передать затухание.

Если f зависит только от времени, наклоны известны сразу, и шаг сводится к двухузловой квадратуре для интеграла f. Если f зависит от состояния, обе стадии обычно приходится искать итерационно. Это главное отличие от классического RK4, где четыре наклона вычисляют последовательно.

Как пользоваться формулой

Запишите уравнение в виде y′ = f(t,y), зафиксируйте текущую пару (tₙ,yₙ) и выберите шаг h.
Найдите времена двух стадий tₙ + c₁h и tₙ + c₂h по симметричным узлам на отрезке [0,1].
Подставьте обе строки A в уравнения для наклонов. Для линейной правой части решите систему линейных уравнений.
Для нелинейной правой части начните с приближения к стадиям и решите связанную систему с контролем невязки.
Прибавьте к yₙ величину h(k₁+k₂)/2, затем сравните один шаг h с двумя шагами h/2 на том же интервале.

Историческая справка

В 1964 году Джон Бутчер опубликовал статью Implicit Runge–Kutta Processes в Mathematics of Computation. В ней он связал квадратурные формулы Гаусса — Лежандра с неявными схемами для задачи Коши и получил семейство порядка 2s с s стадиями. Двухстадийная схема стала конкретным примером того, как внутренние наклоны можно связать системой и достичь четвёртого порядка.

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

В 1969 году Байрон Эле исследовал множители устойчивости схем Гаусса и получил диагональные аппроксимации Паде для экспоненты. Этот результат объясняет, почему высокий порядок и A-устойчивость могут сочетаться со слабым затуханием при крупном шаге. Варианты на узлах Радау изучались в том числе ради устранения этого недостатка.

Историческая линия формулы

Имя Гаусса — Лежандра относится к узлам и весам квадратуры. Построение семейства неявных методов Рунге — Кутты по этим узлам связано со статьёй Джона Бутчера 1964 года. Множитель устойчивости удобно записывать как аппроксимацию Паде [2/2], в соответствии с исследованием Байрона Эле 1969 года.

Пример

В идеализированном растворе концентрация убывает по закону y′ = −y/τ. Найдите концентрацию через секунду одним шагом двухстадийного метода Гаусса. Дано y₀ = 1 моль/л, τ = 1 с, h = 1 с. Для вычисления числовых значений наклонов в моль/(л·с) обозначим d = √3/6. После переноса членов со стадиями влево получаем систему (5/4)k₁ + (1/4 − d)k₂ = −1 и (1/4 + d)k₁ + (5/4)k₂ = −1. Определитель системы равен (5/4)² − (1/4 − d)(1/4 + d) = 19/12. Отсюда k₁ = −12(1 + d)/19 и k₂ = −12(1 − d)/19. Обновление даёт y₁ = 1 + (k₁ + k₂)/2 = 1 − 12/19 = 7/19 ≈ 0,368421 моль/л. Проверка. Произведение шага в секундах на наклон даёт моль/л, а точное значение e⁻¹ ≈ 0,367879 моль/л. Абсолютная ошибка одного шага составляет около 0,000542 моль/л.

Частая ошибка

Считать k₂ после k₁ как в явном RK. Коэффициенты вне диагонали здесь ненулевые, поэтому обе строки нужно решать совместно. Взять узлы ±√3/3 из квадратуры на [−1,1] без преобразования. Для времени шага используют узлы 1/2 ± √3/6 на [0,1]. Подставить в обновление значения y на стадиях вместо наклонов. Множитель h должен умножаться на величину с единицей y/с. Выбирать шаг только по A-устойчивости. Для быстро затухающей компоненты отдельно проверьте отношение численной амплитуды к точной.

Практика

Задачи с решением

Полиномиальная правая часть

Условие. Решите за один шаг h = 1 задачу y′ = t³, y(0) = 2. Переменные безразмерны.

Решение. Стадии не зависят от состояния, поэтому k₁ = c₁³ и k₂ = c₂³. Узлы удовлетворяют c₁ + c₂ = 1 и c₁c₂ = 1/6. Отсюда c₁³ + c₂³ = 1³ − 3·(1/6)·1 = 1/2. Тогда y₁ = 2 + (1/2)(1/2) = 9/4. Точный интеграл t³ от 0 до 1 равен 1/4.

Ответ. y(1) = 9/4 = 2,25. Один шаг совпадает с точным решением.

Крупный шаг для быстрого затухания

Условие. Для y′ = −10y, y(0) = 1 найдите остаток амплитуды через секунду при h = 1 с. Сравните с точным решением.

Решение. Здесь λ = −10 с⁻¹, поэтому z = hλ = −10. Из множителя Гаусса получаем R(−10) = (12 − 60 + 100)/(12 + 60 + 100) = 52/172 = 13/43. Точно y(1) = e⁻¹⁰ ≈ 0,00004540. Численное значение остаётся меньше начального, но значительно превосходит точное.

Ответ. Метод оставляет около 30,23% начальной амплитуды. Точное решение оставляет около 0,00454%.

Дополнительные источники

https://www.jstor.org/stable/2003405
https://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/IRK.pdf
https://www.unige.ch/~hairer/preprints/coimbra.pdf
https://cs.uwaterloo.ca/research/tr/1969/CS-RR-2010.pdf

Связанные формулы

Математика

Неявная средняя точка: шаг и сохранение структуры

$y_{n+1}=y_n+h f\!\left(t_n+\frac h2,\frac{y_n+y_{n+1}}2\right)$

Неявная средняя точка вычисляет производную в среднем старого и нового состояний. Новый уровень входит в правую часть, поэтому шаг находят решением уравнения. Схема второго порядка обратима по времени при согласованном решении.

Математика

Классический метод Рунге — Кутты четвёртого порядка

$\begin{aligned}k_1&=f(t_n,y_n),&k_2&=f(t_n+h/2,y_n+hk_1/2),\\k_3&=f(t_n+h/2,y_n+hk_2/2),&k_4&=f(t_n+h,y_n+hk_3),\\y_{n+1}&=y_n+\frac h6(k_1+2k_2+2k_3+k_4).\end{aligned}$

Классический RK4 вычисляет четыре наклона за шаг и объединяет их с весами 1:2:2:1. Для гладкой задачи Коши он даёт четвёртый порядок общей ошибки без явного вычисления высших производных.

Математика

Квадратура Гаусса—Лежандра с двумя узлами

$G_2=\frac{b-a}{2}\left[f\!\left(m-\frac{b-a}{2\sqrt3}\right)+f\!\left(m+\frac{b-a}{2\sqrt3}\right)\right],\quad m=\frac{a+b}{2}$

Двухузловая квадратура Гаусса—Лежандра вычисляет функцию в двух специально выбранных внутренних точках. На конечном отрезке с единичным весом она точно интегрирует любой многочлен степени не выше трёх.

Математика

Неявный метод Эйлера: шаг и устойчивость затухания

$y_{n+1}=y_n+h f(t_{n+1},y_{n+1})$

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