Математика / Численные методы ОДУ
Метод Розенброка — Эйлера с одним решением системы
Метод Розенброка — Эйлера находит приращение решения автономного ОДУ из одной линейной системы с якобианом. На линейной задаче он совпадает с неявным Эйлером, а на нелинейной даёт шаг первого порядка.
Формула
На графике невязка G(Δ) = Δ + (1+Δ)²/2 и её касательная при Δ = 0. Нуль касательной даёт Δ = −0,25. Нуль полной невязки даёт Δ ≈ −0,268.
Обозначения
- yₙ, yₙ₊₁
- старое и новое состояния, единицы компонент состояния
- $h$
- длина временного шага, с, если время задано в секундах
- $f(y)$
- правая часть автономного уравнения y′=f(y), единицы компонент в секунду
- $Jₙ$
- якобиан f по компонентам y в точке yₙ, для скаляра с⁻¹ при времени в секундах
- $I$
- единичная матрица размерности состояния, безразмерные диагональные элементы
- $Δₙ$
- приращение состояния за шаг, единицы компонент состояния
Условия применения
Формула записана для автономной задачи Коши. Нужны начальное состояние, дифференцируемая правая часть и якобиан, вычисленный в текущем состоянии. Матрица I−hJₙ должна быть обратимой. Для оценки первого порядка предполагают достаточную гладкость в окрестности решения и контролируемую погрешность линейного решателя.
Ограничения
Один линейный решатель может быть дорогим для большой системы, если матрица плотная. Близость I−hJₙ к вырожденной ухудшает точность приращения. Для явной зависимости f от времени нужна отдельно определённая схема, например применение автономной формулы к расширенной системе t′=1. Устойчивость на линейном тесте сама по себе не даёт малую погрешность при крупном шаге.
Подробное объяснение
Формула обновляет состояние приращением Δₙ. В явном Эйлере это приращение равно hf(yₙ). Здесь оно сначала проходит через линейную систему с I−hJₙ. Якобиан учитывает, как меняется скорость движения при изменении каждой компоненты состояния.
Неявный Эйлер требует решить уравнение G(Δ)=Δ−hf(yₙ+Δ)=0. При разложении правой части вокруг yₙ получается f(yₙ+Δ)≈f(yₙ)+JₙΔ. Отсюда выходит линейная система для приращения. Её решение совпадает с первой итерацией Ньютона для G, начатой с Δ=0.
Для аффинной функции f(y)=Ay+b разложение точно, поэтому шаг совпадает с неявным Эйлером. На тесте y′=λy это даёт R(z)=1/(1−z), где z=hλ. Множитель ограничен по модулю в левой полуплоскости и стремится к нулю на отрицательной вещественной оси. Так проявляется L-устойчивость.
На нелинейной задаче остаток разложения зависит от квадрата приращения. Как правило, после одной линейной операции остаётся невязка исходного нелинейного уравнения. При уменьшении h метод сохраняет глобальный первый порядок. Уменьшение невязки за дополнительные итерации решает уже задачу неявного Эйлера.
Для системы Δ лучше находить линейным решателем, например факторизацией, а не строить обратную матрицу. При неизвестном точном ответе вычислите полный шаг и два половинных, каждый раз обновляя якобиан в текущей точке. Разность показывает, насколько результат чувствителен к сетке.
Как пользоваться формулой
- Запишите автономную систему y′=f(y) и выберите шаг h. Вычислите f в известном состоянии yₙ.
- Возьмите частные производные каждой компоненты f по компонентам y и соберите якобиан в состоянии yₙ.
- Составьте I−hJₙ и правую часть h·f(yₙ), используя одну и ту же систему единиц для времени и производных.
- Решите линейную систему на Δₙ и проверьте невязку. Прибавьте найденное приращение к исходному состоянию.
- Проверьте шаг по известному решению или по двум половинным шагам, для второго обновите якобиан.
Историческая справка
В 1963 году Говард Розенброк опубликовал в The Computer Journal статью Some General Implicit Processes for the Numerical Solution of Differential Equations. Он предложил строить аналоги неявных схем Рунге — Кутты, в которых на шаге решаются линейные системы. Это позволяло получить широкую область устойчивости без последовательного решения полных нелинейных уравнений.
Линейно-неявный шаг Эйлера наглядно показывает основной приём: достаточно одной линеаризации уравнения для нового состояния. В современном обозначении это матрица I−hJₙ и правая часть hf(yₙ), где якобиан оценивается до перехода.
В многостадийных методах Розенброка и Розенброка — Ваннера идея развилась в последовательность линейных решений с одной матрицей на шаге. Например, Шампайн и Райхельт в 1997 году описали решатель ode23s на основе формулы Розенброка второго порядка. Одностадийная формула здесь сохраняет первый порядок и показывает принцип без дополнительных стадий.
Историческая линия формулы
Эйлер в названии указывает на исходное одношаговое обновление. Линейно-неявная запись получается из линеаризации уравнения неявного Эйлера. Имя Говарда Розенброка относится к общему подходу на основе линейных систем, изложенному им в 1963 году и развитому в многостадийных схемах.
Пример
Безразмерная концентрация убывает по нелинейному закону y′=−y². Найдите приближение из y(0)=1 за шаг h=1/2. Правая часть f(1)=−1, якобиан скалярный J₀=−2y₀=−2. Тогда (1−hJ₀)Δ₀=hf(y₀) принимает вид [1−(1/2)(−2)]Δ₀=−1/2. Получаем Δ₀=−1/4 и y₁=1−1/4=3/4. Точное решение y(t)=1/(1+t) даёт y(1/2)=2/3, поэтому абсолютная погрешность равна 1/12. Для сравнения, неявный Эйлер требует решить y₁=1−y₁²/2 и даёт положительный корень √3−1≈0,732051. Проверка. Все величины нормированы, а подстановка Δ₀=−1/4 в линейную строку возвращает правую часть −1/2. Нелинейная невязка Δ₀+(1/2)(1+Δ₀)² равна 1/32, как видно из различия двух приближений.
Частая ошибка
Подставить якобиан в новом неизвестном состоянии. В формуле Розенброка — Эйлера производные берутся в известной точке yₙ. Дважды умножить на h. Правая часть уже равна hf(yₙ), поэтому Δₙ прибавляют к состоянию сразу. Ожидать совпадения с неявным Эйлером при любой f. Совпадение следует из точной линеаризации аффинной функции, у нелинейной остаётся остаток. Считать только диагональные производные. При связанных уравнениях внедиагональные элементы передают влияние между компонентами.
Практика
Задачи с решением
Аффинная правая часть
Условие. Для y′=2−y, y(0)=0 вычислите один безразмерный шаг h=1 и проверьте формулу неявного Эйлера.
Решение. f(0)=2, J₀=−1. Поэтому 2Δ₀=2 и y₁=0+1=1. Неявный Эйлер требует y₁=0+2−y₁, то есть 2y₁=2. Точное решение в t=1 равно 2(1−e⁻¹)≈1,264241.
Ответ. Обе схемы дают y₁=1, поскольку правая часть аффинна.
Связанные компоненты
Условие. Для системы y₁′=−y₁+y₂, y₂′=−3y₂ из состояния (1;1) выполните шаг h=1. Величины безразмерны.
Решение. В начальной точке f=(0;−3). Матрица I−J имеет первую строку (2;−1) и вторую (0;4). Из второй строки Δ₂=−3/4, из первой 2Δ₁−Δ₂=0, поэтому Δ₁=−3/8. Прибавляем приращения к (1;1).
Ответ. Новое состояние (5/8; 1/4). Первая компонента изменилась из-за связи со второй, хотя её начальная производная равна нулю.
Дополнительные источники
- https://academic.oup.com/comjnl/article-abstract/5/4/329/316388
- https://www.wias-berlin.de/people/john/LEHRE/NUMERIK_II_24_25/ode_2.pdf
- https://doi.org/10.1137/S1064827594276424
Связанные формулы
Математика
Явный метод Эйлера: решение задачи Коши по шагам
Явный метод Эйлера приближённо решает уравнение y′ = f(t,y), продолжая решение по наклону в начале каждого шага. При гладком решении глобальная ошибка имеет первый порядок по шагу.
Математика
Неявный метод Эйлера: шаг и устойчивость затухания
Неявный метод Эйлера использует производную в конце шага. Новое значение входит в обе части равенства, поэтому каждый переход требует решения уравнения, зато схема хорошо подавляет быстро затухающие компоненты.
Математика
Метод Тейлора второго порядка для задачи Коши
Метод Тейлора второго порядка добавляет к шагу Эйлера член с полной производной правой части. Для скалярного уравнения y′=f(t,y) вторая производная решения равна fₜ+fᵧf, вычисленной в начале шага.