Схема Кранка-Николсона: устойчивость и точность

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

Уравнение теплопроводности и сеточная дискретизация

Рассмотрим одномерное уравнение теплопроводности (диффузии):

$$\frac{\partial u}{\partial t} = a \frac{\partial^2 u}{\partial x^2}, \quad x \in [0, L],\; t \geq 0,$$

где $a > 0$ - коэффициент температуропроводности (диффузии). Введём равномерную сетку с шагом по пространству $h = L/N$ и шагом по времени $\tau$. Узлы сетки: $x_i = i h$, $t_n = n\tau$. Значение решения в узле обозначим $u_i^n \approx u(x_i, t_n)$.

Центральная разностная аппроксимация второй производной по пространству в узле $x_i$ имеет второй порядок точности:

$$\frac{\partial^2 u}{\partial x^2}\bigg|_{x_i} \approx \frac{u_{i-1} - 2u_i + u_{i+1}}{h^2} + O(h^2).$$

Ключевой вопрос - как аппроксимировать правую часть: на новом временном слое $n+1$ (неявная схема), на старом $n$ (явная схема), или в полусумме?

Вывод схемы Кранка-Николсона

Схема Кранка-Николсона берёт среднее арифметическое пространственного лапласиана на двух соседних временных слоях. Производная по времени аппроксимируется первым порядком (что даёт второй порядок по времени за счёт симметрии):

$$\frac{u_i^{n+1} - u_i^n}{\tau} = \frac{a}{2}\left(\frac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{h^2} + \frac{u_{i-1}^{n+1} - 2u_i^{n+1} + u_{i+1}^{n+1}}{h^2}\right).$$

Введём сеточное число $r = a\tau / h^2$ (коэффициент Куранта-Фридрихса-Леви для задачи теплопроводности). Тогда уравнение в узле $i$ принимает вид:

$$-\frac{r}{2} u_{i-1}^{n+1} + (1 + r) u_i^{n+1} - \frac{r}{2} u_{i+1}^{n+1} = \frac{r}{2} u_{i-1}^n + (1 - r) u_i^n + \frac{r}{2} u_{i+1}^n.$$

Правая часть полностью известна из предыдущего слоя. Левая часть порождает трёхдиагональную систему уравнений на $N-1$ неизвестных $u_1^{n+1}, \ldots, u_{N-1}^{n+1}$ (граничные условия $u_0^{n+1}$ и $u_N^{n+1}$ заданы). Такая система решается за $O(N)$ операций методом прогонки (алгоритм Томаса).

Рассчитать для своей задачи

Порядок аппроксимации

Локальная погрешность аппроксимации схемы Кранка-Николсона получается разложением Тейлора. По времени: производная аппроксимируется разностным отношением с погрешностью $O(\tau^2)$ (симметрия полусуммы убирает члены с нечётными степенями $\tau$). По пространству: центральная разностная формула даёт $O(h^2)$. Итоговый порядок:

$$\psi_i^n = O(\tau^2 + h^2).$$

Это второй порядок по обеим переменным - на порядок лучше, чем у явной схемы по времени ($O(\tau + h^2)$) и полностью неявной ($O(\tau + h^2)$). На практике это означает: при вдвое меньшем шаге погрешность падает вчетверо, а не вдвое.

Рассчитать для своей задачи

Устойчивость по фон Нейману

Анализ устойчивости методом фон Неймана подставляет в схему пробное решение $u_i^n = G^n e^{ik\xi h}$, где $G$ - коэффициент роста, $\xi$ - волновое число. Для схемы Кранка-Николсона:

$$G = \frac{1 - 2r\sin^2(\xi h / 2)}{1 + 2r\sin^2(\xi h / 2)}.$$

Обозначим $z = 2r\sin^2(\xi h / 2) \geq 0$. Тогда $G = (1 - z)/(1 + z)$. Поскольку $z \geq 0$, имеем $-1 \leq G \leq 1$, то есть $|G| \leq 1$ при любых $r > 0$. Схема безусловно устойчива: шаг по времени $\tau$ не ограничен условием устойчивости (в отличие от явной схемы, где требуется $r \leq 1/2$).

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

Матричная запись и метод прогонки

На каждом временном шаге система с трёхдиагональной матрицей $A$ имеет вид:

$$A\, \mathbf{u}^{n+1} = B\, \mathbf{u}^n + \mathbf{b},$$

где $A$ и $B$ - трёхдиагональные матрицы, $\mathbf{b}$ - вектор граничных условий. Матрица $A$ имеет единичный вид на главной диагонали $(1 + r)$ и $-r/2$ на побочных. Метод прогонки (алгоритм Томаса) решает такую систему за $O(N)$ арифметических операций, что делает каждый временной шаг схемы Кранка-Николсона очень быстрым даже при большом $N$.

Рассчитать для своей задачи

Обобщения: параметрический метод тета

Схема Кранка-Николсона является частным случаем параметрического $\theta$-метода:

$$\frac{u_i^{n+1} - u_i^n}{\tau} = a\left[\theta\,\Lambda u^{n+1} + (1-\theta)\,\Lambda u^n\right],$$

где $\Lambda u_i = (u_{i-1} - 2u_i + u_{i+1})/h^2$. При $\theta = 0$ - явная схема (устойчива при $r \leq 1/2$), при $\theta = 1$ - полностью неявная (безусловно устойчива, порядок $O(\tau + h^2)$), при $\theta = 1/2$ - схема Кранка-Николсона (безусловно устойчива, порядок $O(\tau^2 + h^2)$). Таким образом, $\theta = 1/2$ даёт наилучшую точность при сохранении безусловной устойчивости.

Частые ошибки

• Не учитывают граничные условия в правой части. При составлении системы $A\mathbf{u}^{n+1} = \mathbf{b}$ значения $u_0^{n+1}$ и $u_N^{n+1}$ переносятся из левой части в правую со своим коэффициентом $r/2$. Забытые слагаемые дают неверное решение в пристеночных узлах.
• Путают $r$ и условие устойчивости явной схемы. Для схемы Кранка-Николсона нет ограничения $r \leq 1/2$. Брать $r \gg 1$ при явной схеме - катастрофа; при Кранке-Николсоне - просто потеря точности по времени.
• Забывают делить на 2 в усреднении лапласиана. Коэффициент при пространственной разности должен быть $r/2$, а не $r$.
• Применяют схему к нелинейным уравнениям без линеаризации. Прямое применение схемы Кранка-Николсона к нелинейному уравнению даёт нелинейную систему на каждом шаге; нужна линеаризация Кранка-Николсона (Crank-Nicolson-Picard или аналоги).
• Не проверяют сходимость при измельчении сетки. Второй порядок следует верифицировать сравнением результатов при $h$ и $h/2$: ошибка должна упасть примерно в 4 раза.

FAQ

Почему схема Кранка-Николсона лучше явной при одинаковом числе операций? Явная схема требует $\tau \leq h^2/(2a)$ для устойчивости, то есть при $h = 0.01$ и $a = 1$ нужно $\tau \leq 5 \cdot 10^{-5}$. Схема Кранка-Николсона с тем же $h$ позволяет взять $\tau = 0.01$ (в 200 раз крупнее), при этом погрешность по времени $O(\tau^2)$ остаётся контролируемой.

Можно ли применять схему Кранка-Николсона к уравнению переноса (гиперболическому)? Нет, схема разработана для параболических уравнений. Для гиперболических уравнений применяют схемы Лакса-Вендроффа, Бима-Уорминга или схемы с расщеплением потока.

Как выбрать соотношение $h$ и $\tau$ для максимальной точности? Поскольку погрешность $O(\tau^2 + h^2)$, оба члена должны быть одного порядка: $\tau \sim h$. Тогда при вдвое меньшем $h$ и $\tau$ погрешность упадёт в 4 раза. Если взять $\tau \sim h^2$, точность по пространству станет ограничивающей; если $\tau \gg h$, наоборот.

Коротко

Схема Кранка-Николсона аппроксимирует уравнение теплопроводности полусуммой пространственных лапласиан на двух временных слоях. Это даёт второй порядок точности $O(\tau^2 + h^2)$ и безусловную устойчивость: коэффициент роста $|G| = |(1-z)/(1+z)| \leq 1$ при любом $r = a\tau/h^2 > 0$. На каждом шаге решается трёхдиагональная система методом прогонки за $O(N)$ операций. Схема является золотым стандартом для одномерной задачи теплопроводности и основой для обобщений на многомерный случай.