Математика / Численные методы ОДУ
Двухстадийный метод Radau IIA третьего порядка для ОДУ
Двухстадийный Radau IIA строит приближение решения ОДУ в точках одной трети и конца шага. Схема имеет третий порядок, а её множитель устойчивости стремится к нулю для очень быстрых затухающих компонент.
Формула
Первый узел находится при доле 1/3. Наклон в нём входит в итог с весом 3/4, наклон на правой границе — с весом 1/4. Во второй стадии уже получено yₙ₊₁.
Обозначения
- tₙ, h
- время начала шага и его длина, с, если время в секундах
- yₙ, yₙ₊₁
- приближения в начале и конце шага, единица состояния
- Y₁, Y₂
- значения полинома коллокации при tₙ + h/3 и tₙ + h, единица состояния
- k₁, k₂
- наклоны f в первой и второй стадии, единица состояния в секунду
- $f(t,y)$
- правая часть уравнения y′ = f(t,y), единица состояния в секунду
Условия применения
Нужна задача Коши с достаточно гладкой правой частью в окрестности решения и начальным состоянием в области её определения. Для стандартного третьего порядка обе стадии должны решаться с погрешностью, которая мала относительно ошибки дискретизации. Если выбран метод Ньютона, невязку контролируют для обеих строк системы.
Ограничения
Две стадии связаны и в общем случае требуют совместного решения нелинейной системы. L-устойчивость описывает линейный тест и допускает отрицательный множитель, поэтому знак положительной концентрации при крупном шаге нужно проверять отдельно. Оценка третьего порядка рассчитана на гладкое решение и может терять силу, если шаг пересекает скачок правой части.
Подробное объяснение
Radau IIA берёт первую стадию в точке tₙ+h/3, а вторую прямо в конце шага. Обе стадии входят в уравнения друг друга. Как только система решена, значение Y₂ можно принять за yₙ₊₁. Именно это равенство выражает свойство жёсткой точности для данного метода.
Представьте полином второй степени, который начинается в точке yₙ. Его линейную производную задают двумя наклонами k₁ и k₂. Если проинтегрировать этот линейный интерполянт до доли 1/3, выйдут коэффициенты 5/12 и −1/12. Если проинтегрировать до правой границы, получатся веса 3/4 и 1/4. Поэтому формула обновления совпадает со второй строкой системы.
При точном входном значении и достаточно гладком решении ошибка одного шага имеет порядок h⁴, а на фиксированном интервале глобальная ошибка имеет порядок h³. Для линейного теста y′ = λy подстановка kᵢ = λYᵢ даёт R(z) = (1 + z/3)/(1 − 2z/3 + z²/6) при z = hλ.
Модуль R(z) не превосходит единицы в левой полуплоскости, и R(z) стремится к нулю при z → −∞. Это объясняет подавление быстрой затухающей моды. Для вещественного z < −3 числитель становится отрицательным, поэтому затухание может сопровождаться сменой знака. В задаче с концентрацией такой шаг может быть непригоден.
Если f зависит только от t, стадии вычисляются напрямую. При зависимости от y первая строка использует наклон в правой точке, а вторая — наклон в первой. Схема с двумя последовательными независимыми вычислениями правой части здесь дала бы другой результат.
Как пользоваться формулой
- Разметьте моменты tₙ+h/3 и tₙ+h. Запишите оба выражения для Y и оба уравнения для наклонов.
- Для линейной правой части перенесите неизвестные значения стадий влево и решите систему из двух строк.
- Для нелинейного уравнения решите обе строки совместно и проверьте невязку до принятия результата.
- Примите Y₂ как yₙ₊₁ и проверьте ожидаемые свойства модели, например знак концентрации или диапазон напряжения.
- Сравните результат с вычислением на двух шагах h/2, если точное решение для оценки ошибки недоступно.
Историческая справка
В 1964 году Джон Бутчер опубликовал работу Integration Processes Based on Radau Quadrature Formulas. В ней квадратуры Радау стали основой неявных методов Рунге — Кутты. Наличие крайнего узла в квадратуре сделало естественным вариант схемы с одной стадией на границе шага.
Для жёстких уравнений требовалось ещё подавление быстро затухающих компонент. В диссертации 1969 года Байрон Эле исследовал поддиагональные приближения Паде к экспоненте и построил схемы с множителем, который стремится к нулю. Герхард Ваннер отмечает, что к таким расширениям формул Радау независимо пришёл Ове Аксельссон.
В коллокационной записи узлы Radau IIA выбирают с включением правой границы. Для двух стадий это одна треть шага и его конец. Последняя строка таблицы коэффициентов совпадает с весами обновления, поэтому вычисленный в правом узле результат можно сразу переносить на следующий шаг.
Историческая линия формулы
Название указывает на семейство квадратур Радау. Его применение к неявным схемам связано с работой Джона Бутчера 1964 года, а развитие L-устойчивого семейства Radau IIA в 1969 году — с работами Байрона Эле и Ове Аксельссона. Двухстадийная формула является частным случаем этого семейства.
Пример
В линейной модели датчика сигнал затухает с постоянной времени 1 с. Найдите напряжение через секунду одним шагом Radau IIA. Дано y′ = −y/τ, τ = 1 с, y₀ = 1 В, h = 1 с. После подстановки k₁ = −Y₁/τ и k₂ = −Y₂/τ числовые уравнения для напряжений имеют вид (17/12)Y₁ − (1/12)Y₂ = 1 и (3/4)Y₁ + (5/4)Y₂ = 1. Умножаем первую строку на 12, вторую на 4: 17Y₁ − Y₂ = 12 и 3Y₁ + 5Y₂ = 4. Подстановка Y₂ = 17Y₁ − 12 даёт 88Y₁ = 64, откуда Y₁ = 8/11 и Y₂ = 4/11. Ответ y₁ = Y₂ ≈ 0,363636 В. Проверка. Точное напряжение e⁻¹ В ≈ 0,367879 В, ошибка около 0,004243 В. В формулах стадий шаг в секундах умножается на наклон в вольтах в секунду.
Частая ошибка
Считать первый узел серединой. Для двухстадийного Radau IIA он расположен на одной трети шага, второй совпадает с концом. Терять знак перед k₂ в первой строке. Коэффициент −1/12 относится к интегрированию до первого узла и входит в систему со своим знаком. Ещё раз прибавлять взвешенные наклоны к Y₂. Вторая стадия уже содержит обновление относительно yₙ. Считать L-устойчивость гарантией положительного результата. Для z < −3 численный множитель имеет отрицательный знак.
Практика
Задачи с решением
Правая часть зависит от времени
Условие. Для безразмерной задачи y′ = t², y(0) = 0 сделайте шаг h = 1.
Решение. Стадии вычисляются напрямую в t = 1/3 и t = 1: k₁ = 1/9, k₂ = 1. Подставляем во вторую строку: Y₂ = (3/4)(1/9) + 1/4 = 1/12 + 1/4 = 1/3. Интеграл t² от 0 до 1 равен 1/3.
Ответ. y₁ = 1/3. В этом примере метод точно интегрирует квадратичную правую часть.
Отрицательный множитель затухания
Условие. Для y′ = −6y, y(0) = 1 возьмите один шаг h = 1 с. Проверьте знак численного результата.
Решение. При λ = −6 с⁻¹ имеем z = −6. Числитель R равен 1−2 = −1, знаменатель 1+4+36/6 = 11. Поэтому y₁ = −1/11, тогда как точное значение e⁻⁶ ≈ 0,002479 остаётся положительным.
Ответ. Численный результат −0,090909. Амплитуда уменьшилась, но знак изменился.
Дополнительные источники
- https://doi.org/10.1090/S0025-5718-1964-0165693-1
- https://cs.uwaterloo.ca/research/tr/1969/CS-RR-2010.pdf
- https://www.unige.ch/~hairer/preprints/coimbra.pdf
Связанные формулы
Математика
Неявная средняя точка: шаг и сохранение структуры
Неявная средняя точка вычисляет производную в среднем старого и нового состояний. Новый уровень входит в правую часть, поэтому шаг находят решением уравнения. Схема второго порядка обратима по времени при согласованном решении.
Математика
Неявный метод Эйлера: шаг и устойчивость затухания
Неявный метод Эйлера использует производную в конце шага. Новое значение входит в обе части равенства, поэтому каждый переход требует решения уравнения, зато схема хорошо подавляет быстро затухающие компоненты.
Математика
BDF2: неявная формула второго порядка для жёсткого ОДУ
BDF2 аппроксимирует производную в новом узле комбинацией двух предыдущих значений и неизвестного y следующего шага. Это неявный метод второго порядка, который применяют для жёстких задач, где явные схемы требуют слишком малого шага.
Математика
Квадратура Гаусса—Лежандра с двумя узлами
Двухузловая квадратура Гаусса—Лежандра вычисляет функцию в двух специально выбранных внутренних точках. На конечном отрезке с единичным весом она точно интегрирует любой многочлен степени не выше трёх.