Математика / Численные методы ОДУ
Метод TR-BDF2 для ОДУ: трапеция и обратные разности
TR-BDF2 делает один шаг двумя неявными подшагами. Сначала схема трапеций находит внутреннюю точку, затем формула обратных разностей второго порядка использует её вместе с начальным значением.
Формула
Внутренняя точка делит внешний шаг в пропорции γ к 1−γ. Трапеция использует первые две точки, BDF2 — обе уже известные точки и новую правую границу.
Обозначения
- tₙ, h
- начальное время и полный внешний шаг, с, если время в секундах
- $γ$
- доля внешнего шага до внутренней точки, безразмерная
- $Y$
- состояние в момент tₙ+γh после трапеции, единица состояния
- yₙ, yₙ₊₁
- состояния на границах внешнего шага, единица состояния
- $f(t,y)$
- правая часть уравнения y′ = f(t,y), единица состояния в секунду
Условия применения
Начальное состояние известно, правая часть достаточно гладкая, а доля внутреннего подшага удовлетворяет 0 < γ < 1. Для свойств стандартной схемы берут γ = 2−√2. Каждый подшаг решают с контролем невязки. Переиспользование факторизации требует одинаковых матриц, например постоянного якобиана или согласованной заморозки якобиана в итерациях.
Ограничения
Один шаг требует двух неявных решений. На нелинейной задаче переиспользование матрицы между подшагами зависит от выбранного решателя и изменения якобиана. L-устойчивость характеризует подавление линейного теста и сама по себе не сохраняет положительность. Слишком большой шаг может пропустить переходный процесс или дать отрицательное значение положительной величины.
Подробное объяснение
Сначала трапеция проходит только долю γh и находит Y в момент tₙ+γh. После этого известны две точки текущего внешнего шага: yₙ на левой границе и Y внутри. Формула BDF2 использует обе и ищет значение на правой границе, где правая часть уравнения снова вычисляется неявно.
Во второй строке коэффициенты получаются из производной квадратичного полинома через точки с относительными временами 0, γ и 1. При γ=1/2 расстояния равны и восстанавливается обычная BDF2 с шагом h/2. Стандартное γ≈0,586 делает расстояния разными, поэтому подставлять коэффициенты равномерной BDF2 нельзя.
При γ=2−√2 выполняется γ/2=(1−γ)/(2−γ). В обеих неявных строках правая часть в новом состоянии умножается на один и тот же коэффициент γh/2. Если f(t,y)=Ay+b(t) с постоянной A, обе системы имеют матрицу I−γhA/2. Её можно разложить один раз и использовать для двух правых частей.
На тесте y′ = λy при z=hλ обозначим a=γ/2. Внутреннее значение умножает yₙ на (1+az)/(1−az), а итоговый множитель равен R(z)=[1+(1−2a)z]/(1−az)². Стандартная схема A-устойчива, и R(z) стремится к нулю при z→−∞. За подавление быстрой моды отвечает полный шаг. Трапеция внутри него может сменить знак.
Для гладкого решения с точно найденными подшагами ошибка на фиксированном интервале имеет второй порядок. Смена внешнего шага не требует пересчитывать историю: обе точки, нужные BDF2, создаются заново из текущего состояния. После резкого изменения правой части размер шага всё равно выбирают по точности перехода.
Как пользоваться формулой
- Задайте tₙ, yₙ и внешний шаг h. Для стандартной схемы вычислите γ=2−√2 и внутреннее время tₙ+γh.
- Найдите внутреннее значение Y по неявной трапеции, используя правую часть в начале и в конце первого подшага.
- Подставьте Y и yₙ в формулу обратных разностей с долей γ и решите вторую систему относительно yₙ₊₁.
- Проверьте невязки обеих строк. Для линейной задачи сравните матрицы до повторного использования факторизации.
- Оцените точность полного перехода и физически допустимый диапазон, прежде чем принять вычисленное значение.
Историческая справка
В 1985 году Рэндольф Бэнк и соавторы описали вычисление переходных процессов в кремниевых приборах и схемах. Их составной шаг сочетал неявную трапецию и BDF2. В расчёте приборов это позволяло учитывать быстрые переходы, сохраняя преимущество одношаговой процедуры: новый шаг строится из текущего состояния.
В 1996 году М. Хосеа и Л. Шампайн исследовали устойчивость TR-BDF2 и её реализацию как общего метода для ОДУ. Они предложили способ корректировать оценку локальной погрешности для очень жёстких задач и построили непрерывный интерполянт с непрерывной первой производной.
В современной записи обе части часто рассматривают как диагонально-неявную схему Рунге — Кутты с явной стадией в начале шага. Запись через трапецию и обратные разности сохраняет исходную конструкцию и прямо показывает, почему BDF2 хватает точек внутри одного внешнего шага.
Историческая линия формулы
Составной метод TR-BDF2 связан с работой Рэндольфа Бэнка и соавторов 1985 года по моделированию полупроводников и схем. Анализ устойчивости, интерполяции и оценки погрешности продолжили М. Хосеа и Л. Шампайн в 1996 году. Название отражает две составные формулы, трапецию и BDF2.
Пример
Конденсатор заряжается до напряжения 2 В с постоянной времени 1 с. Найдите напряжение на первом шаге TR-BDF2 длиной 1 с. Дано y′=(2−y)/τ, τ=1 с, y₀=0 В. Положим a=γ/2≈0,292893, b=1/[γ(2−γ)]≈1,207107. При h=1 с начальный наклон равен 2 В/с, наклон в Y численно равен 2−Y. Первая строка даёт (1+a)Y=4a, откуда Y≈0,906164 В. Во второй строке член с y₀ равен нулю. Вычисляем (1+a)y₁=bY+2a. Подстановка даёт y₁≈(1,207107·0,906164+2·0,292893)/1,292893≈1,299119 В. Проверка. Оба подшага дают вольты, так как коэффициенты безразмерны, а h·f имеет размерность напряжения. Точное значение 2(1−e⁻¹) В ≈ 1,264241 В.
Частая ошибка
Провести трапецию на всём h. Она должна закончиться в tₙ+γh, иначе внутренняя точка окажется на месте правой границы. Взять коэффициенты BDF2 для равных интервалов. Они подходят при γ=1/2. Для стандартной схемы интервалы отличаются. Заменить yₙ значением с предыдущего внешнего шага. Вторая строка использует левую границу текущего шага. Отождествить совпадающие коэффициенты перед f с одинаковым якобианом. На нелинейной задаче якобианы в двух состояниях могут различаться.
Практика
Задачи с решением
Равные интервалы как частный случай
Условие. Возьмите для сравнения γ=1/2, h=2 и решите один шаг для y′=t, y(0)=0. Все величины безразмерны.
Решение. Внутреннее время равно 1. Трапеция даёт Y=0+(1/2)(0+1)=1/2. Коэффициенты второй строки равны 4/3 перед Y, −1/3 перед y₀ и 2/3 перед f(2,y₁). Поэтому y₁=(4/3)(1/2)−(1/3)·0+(2/3)·2=2.
Ответ. Y=1/2 и y₁=2. Оба значения совпадают с точным решением t²/2.
Проверка положительности на крупном шаге
Условие. Примените стандартный TR-BDF2 к y′=−4y, y(0)=1 на шаге h=1 с. Проверьте знак.
Решение. При z=−4 получаем Y=(1−4a)/(1+4a)≈−0,079009. Полный множитель R(−4)=[1−4(1−2a)]/(1+4a)²≈−0,139290, где a≈0,292893. Точное значение e⁻⁴≈0,018316 положительно.
Ответ. y₁≈−0,139290. Численная амплитуда снизилась, но такой шаг нарушил положительность.
Дополнительные источники
- https://doi.org/10.1109/TCAD.1985.1270142
- https://pages.mtu.edu/~struther/Courses/OLD/Other/LessOld/5628_11/TRBDF2_Hosea_AppNumMath_20_pp21-37_1996.pdf
- https://www.mate.polimi.it/biblioteca/add/qmox/56-2015.pdf
Связанные формулы
Математика
Неявная трапециевидная схема Кранка — Николсон для ОДУ
Неявная трапециевидная схема усредняет производные на концах шага, включая неизвестный новый уровень. Она имеет второй порядок и A-устойчивость, но при крупном шаге затухание может сопровождаться сменой знака.
Математика
BDF2: неявная формула второго порядка для жёсткого ОДУ
BDF2 аппроксимирует производную в новом узле комбинацией двух предыдущих значений и неизвестного y следующего шага. Это неявный метод второго порядка, который применяют для жёстких задач, где явные схемы требуют слишком малого шага.
Математика
Неявный метод Эйлера: шаг и устойчивость затухания
Неявный метод Эйлера использует производную в конце шага. Новое значение входит в обе части равенства, поэтому каждый переход требует решения уравнения, зато схема хорошо подавляет быстро затухающие компоненты.