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

Метод Розенброка — Эйлера с одним решением системы

Метод Розенброка — Эйлера находит приращение решения автономного ОДУ из одной линейной системы с якобианом. На линейной задаче он совпадает с неявным Эйлером, а на нелинейной даёт шаг первого порядка.

Опубликовано: Обновлено:

Формула

$$\begin{gathered} J_n=f'(y_n),\qquad (I-hJ_n)\Delta_n=h f(y_n),\\ y_{n+1}=y_n+\Delta_n \end{gathered}$$
Схема Одна линеаризация вместо решения нелинейного уравнения
Невязка G(Δ) −0,50Δ −0,268−0,250 Нелинейная невязка Касательная при Δ = 0

На графике невязка 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 метод сохраняет глобальный первый порядок. Уменьшение невязки за дополнительные итерации решает уже задачу неявного Эйлера.

Для системы Δ лучше находить линейным решателем, например факторизацией, а не строить обратную матрицу. При неизвестном точном ответе вычислите полный шаг и два половинных, каждый раз обновляя якобиан в текущей точке. Разность показывает, насколько результат чувствителен к сетке.

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

  1. Запишите автономную систему y′=f(y) и выберите шаг h. Вычислите f в известном состоянии yₙ.
  2. Возьмите частные производные каждой компоненты f по компонентам y и соберите якобиан в состоянии yₙ.
  3. Составьте I−hJₙ и правую часть h·f(yₙ), используя одну и ту же систему единиц для времени и производных.
  4. Решите линейную систему на Δₙ и проверьте невязку. Прибавьте найденное приращение к исходному состоянию.
  5. Проверьте шаг по известному решению или по двум половинным шагам, для второго обновите якобиан.

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

В 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_{n+1}=y_n+h f(t_n,y_n),\qquad t_{n+1}=t_n+h$

Явный метод Эйлера приближённо решает уравнение y′ = f(t,y), продолжая решение по наклону в начале каждого шага. При гладком решении глобальная ошибка имеет первый порядок по шагу.

Математика

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

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

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

Математика

Метод Тейлора второго порядка для задачи Коши

$y_{n+1}=y_n+h f_n+\frac{h^2}{2}\left(f_t+f_y f\right)_n$

Метод Тейлора второго порядка добавляет к шагу Эйлера член с полной производной правой части. Для скалярного уравнения y′=f(t,y) вторая производная решения равна fₜ+fᵧf, вычисленной в начале шага.