Математика / Начала анализа

Неявная трапециевидная схема Кранка — Николсон для ОДУ

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

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

Формула

$$y_{n+1}=y_n+\frac h2[f(t_n,y_n)+f(t_{n+1},y_{n+1})]$$
Схема Неявная трапеция: затухание
Неявная трапеция: затуханиеНеявная трапеция: затухание10,500,25шаг: 0,6

x — время t от 0 до 0,25; y от 0,5 до 1. Линия exp(−2t), точка — шаг 0,6.

Обозначения

$y_n,y_{n+1}$
старое и новое состояния, единица y
$t_n,t_{n+1}$
концы шага, единица t
$h$
шаг, единица t
$f$
правая часть ОДУ, единица y/единица t

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

Требуется начальное состояние и возможность решить алгебраическое уравнение для y_{n+1}. Для нелинейного f применяют Ньютон или другой внутренний решатель с контролем невязки. Гладкость решения нужна для второго порядка; шаг не должен пересекать неразрешённый скачок.

Ограничения

Безусловная линейная устойчивость не означает монотонность: при больших шагах на жёстком затухании возможны слабые колебания. Схема не является автоматически точной для нелинейностей и требует контроля сходимости внутреннего решения. При разрыве коэффициентов нужно согласовать узел с событием.

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

Трапециевидная схема следует из приближения интеграла производной площадью трапеции: среднее значение f в начале и конце шага умножается на h. Конечная производная зависит от неизвестного y_{n+1}, поэтому формула не даёт число одним присваиванием. Для каждого шага возникает нелинейное или линейное алгебраическое уравнение.

На линейной модели y'=λy схема имеет множитель R(z)=(1+z/2)/(1-z/2), z=hλ. При отрицательном λ модуль множителя не превосходит единицы, но при большом |z| он стремится к -1. Это объясняет одновременно устойчивость и возможное чередование знака у быстро затухающей компоненты.

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

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

Название Кранка—Николсон особенно распространено для параболических уравнений после пространственной дискретизации. В статье об ОДУ корректнее говорить о неявном трапециевидном методе; связь с Кранком—Николсон — это тот же временной оператор, а не новая формула для любой задачи.

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

  1. Запишите остаток G(z)=z-y_n-h[f(t_n,y_n)+f(t_{n+1},z)]/2 и задайте начальное приближение z.
  2. Для линейной f решите полученное равенство напрямую; для нелинейной примените Ньютон или другой проверенный решатель.
  3. Проверьте невязку G(z), а затем примите z как y_{n+1}; не смешивайте ошибку решателя с ошибкой сетки.
  4. Перенесите время на h и повторите шаг, контролируя знак, масштаб и возможные события модели.
  5. Сравните сетки h и h/2; отдельно проверьте устойчивость на тесте y'=λy при отрицательном λ.

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

Трапециевидное правило интегрирования относится к классической квадратуре и задолго предшествует компьютерным решателям ОДУ. Когда его применяют к интегральной форме y(t_{n+1})-y(t_n)=∫fdt, получается неявная трапециевидная схема. В анализе временной дискретизации уравнений диффузии эта же идея получила имя Кранка—Николсон после работ Джона Кранка и Филлиса Николсона в середине XX века. Для ОДУ метод рассматривается как двухточечный линейный многошаговый или одношаговый неявный метод с вторым порядком. Позднейшая теория устойчивости, жёсткости и неявных решателей уточнила, почему он полезен и где даёт колебания. Современные библиотеки используют трапециевидные шаги вместе с контролем невязки, адаптацией и обработкой событий; одна формула без этих механизмов не является готовым промышленным решателем.

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

Неявная трапециевидная формула выросла из общей квадратуры; имя Кранка—Николсон закрепилось прежде всего за временной схемой для диффузии. Атрибуция не означает, что авторы создали всю теорию неявных ОДУ или единственную реализацию метода.

Пример

Решим y'=-2y, y(0)=1 одним шагом h=0,25. Неявная трапеция даёт y1=1+0,125[-2-2y1]. После раскрытия скобок y1=1-0,25-0,25y1, поэтому 1,25y1=0,75 и y1=0,6. Точное значение e^-0,5≈0,606530660, ошибка около 0,006530660. Множитель схемы здесь (1-0,25)/(1+0,25)=0,6. Ответ показывает, что устойчивость не равна точности: при более крупном шаге решение всё ещё затухает, но может заметно отличаться от экспоненты. При h=0,125 два шага дают примерно 0,604938, поэтому сравнение сеток отделяет дискретизацию от ошибки решателя.

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

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

Практика

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

Линейное затухание

Условие. Для y'=-y, y(0)=1 и h=0,5 найдите один шаг.

Решение. y1=1+0,25(-1-y1), поэтому 1,25y1=0,75.

Ответ. y1=0,6.

Постоянная скорость

Условие. Для y'=3, y(0)=2 и h=0,4 найдите y1.

Решение. Оба наклона равны 3, поэтому y1=2+0,4·3.

Ответ. y1=3,2.

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

  • https://fncbook.com/multistep/
  • https://fncbook.com/implicit/

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

Математика

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

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

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

Математика

Метод Хойна: явное усреднение двух наклонов

$k_1=f(t_n,y_n),\quad k_2=f(t_n+h,y_n+hk_1),\quad y_{n+1}=y_n+\frac h2(k_1+k_2)$

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

Математика

Метод Адамса — Башфорта второго порядка для ОДУ

$y_{n+1}=y_n+h\left(\frac32 f_n-\frac12 f_{n-1}\right)$

Двухшаговый метод Адамса — Башфорта использует два уже вычисленных наклона и получает явное приближение второго порядка с одним новым вызовом правой части.