Устойчивость явных схем: условие Куранта (CFL)

Аккуратно записать разностную схему — половина дела. Если шаг по времени слишком велик, та же самая схема не приблизит решение, а взорвёт его растущими осцилляциями.

Условие Куранта-Фридрихса-Леви (CFL) — ограничение на соотношение шага по времени и шага по пространству, при нарушении которого явная разностная схема становится неустойчивой и решение разрушается.

В уроке про теплопроводность мы как бы между делом взяли r = 0.4 и предупредили, что число r нельзя брать слишком большим. Настал момент разобраться, почему. Оказывается, у явных схем есть жёсткий порог устойчивости. Перешагнёшь его — и численное решение, вместо того чтобы плавно сглаживаться, начинает раскачиваться со всё растущей амплитудой, пока не превратится в бессмысленную мешанину из огромных чисел противоположных знаков. Это явление называют численной неустойчивостью, а условие, удерживающее от него, — условием CFL.

Порог для теплопроводности: r ≤ 1/2

Для явной схемы теплопроводности условие устойчивости выглядит так: r = α·dt/dx^2 <= 1/2. Если r не превышает 0.5, схема устойчива: ошибки затухают, профиль плавно сглаживается. Если r > 0.5, схема взрывается: малейшие возмущения начинают расти от шага к шагу, появляются знакопеременные осцилляции, и уже через десяток шагов значения вылетают за всякие разумные пределы.

Откуда берётся именно 1/2? Посмотрите на формулу обновления центрального узла: u_i^{n+1} = u_i^n + r·(u_{i-1} − 2u_i + u_{i+1}) = (1 − 2r)·u_i + r·u_{i-1} + r·u_{i+1}. Новое значение — взвешенное среднее старых. Чтобы это среднее было физичным (новое значение не выскакивало за пределы старых), все веса должны быть неотрицательны. Веса у соседей равны r > 0, а вес у самого узла равен (1 − 2r). Он неотрицателен ровно при r <= 1/2. Стоит r превысить 1/2, как центральный вес становится отрицательным — узел начинает "отталкиваться" от своего прежнего значения и перелетать на противоположную сторону, порождая осцилляции.

Демонстрация: устойчивость против взрыва

Сравним два прогона на одном и том же начальном пике [0, 0, 0, 1, 0, 0, 0]: с r = 0.4 (ниже порога) и r = 0.6 (выше порога). Будем печатать max|u| — максимальное по модулю значение — на каждом шаге. У устойчивой схемы оно должно монотонно затухать, у неустойчивой — расти.

def run(r, steps):
    u = [0, 0, 0, 1.0, 0, 0, 0]
    out = []
    def step(u):
        new = u[:]
        for i in range(1, len(u) - 1):
            new[i] = u[i] + r * (u[i-1] - 2*u[i] + u[i+1])
        new[0] = 0.0
        new[-1] = 0.0
        return new
    for s in range(steps + 1):
        out.append(round(max(abs(x) for x in u), 4))
        u = step(u)
    return out

print('r=0.4:', run(0.4, 8))
print('r=0.6:', run(0.6, 8))

Вывод:

r=0.4: [1.0, 0.4, 0.36, 0.24, 0.232, 0.1795, 0.1736, 0.1484, 0.1357]
r=0.6: [1.0, 0.6, 0.76, 0.72, 0.952, 0.9638, 1.3237, 1.3983, 1.936]

Картина наглядна. При r = 0.4 максимум уверенно ползёт вниз: 1.0 → 0.4 → 0.36 → ... → 0.1357. Тепло размывается, амплитуда падает — ровно так, как требует физика остывания. При r = 0.6 максимум сначала проседает до 0.6, но затем начинает раскачиваться и неумолимо расти: 0.76, 0.72, 0.952, и к восьмому шагу уже 1.936 — почти вдвое выше исходного пика! Профиль при этом покрывается зубцами разных знаков. Это и есть численный взрыв: схема не моделирует физику, а генерирует артефакт, который через несколько десятков шагов уйдёт в бесконечность.

Физический смысл условия CFL

За формулой стоит очень наглядная идея, относящаяся ко всем явным схемам. На каждом шаге по времени узел "общается" только со своими непосредственными соседями по сетке — информация за один шаг распространяется ровно на одну клетку. Но в реальном уравнении возмущение тоже распространяется с какой-то скоростью. Условие CFL требует, чтобы за один шаг по времени физическое возмущение не успевало перепрыгнуть больше одной клетки сетки. Иначе схема "не успевает" передать информацию туда, куда она должна была дойти физически, и расчёт разваливается.

Эту же мысль удобно сформулировать на языке областей зависимости, как и сделали сами Курант, Фридрихс и Леви. У точного решения есть физическая область зависимости — множество точек начальных данных, которые реально влияют на значение в данной точке (для волны это конус, в который укладывается всё, до чего сигнал успел дойти). У разностной схемы есть своя численная область зависимости — множество узлов сетки, через которые формула фактически дотягивается до начального слоя за нужное число шагов. Условие CFL — это требование, чтобы численная область зависимости содержала в себе физическую. Если шаг по времени слишком велик, численная область оказывается уже физической: схема пытается вычислить значение, не «дотянувшись» до тех начальных данных, которые на самом деле его определяют. Информации физически не хватает, и расчёт обречён — он не может дать верный ответ в принципе, как бы аккуратно мы ни считали.

Для волнового уравнения это правило записывается совсем прозрачно: число Куранта C = c·dt/dx <= 1. Здесь c — физическая скорость волны, а dx/dt — "скорость сетки" (одна клетка за шаг). Условие C <= 1 буквально означает: волна за шаг проходит не больше одной клетки. Для теплопроводности диффузия не имеет конечной скорости, и анализ даёт квадратичную по dx границу r <= 1/2 — но дух тот же: шаг по времени ограничен шагом по пространству.

И ещё одно важное наблюдение о характере наказания за нарушение CFL. Это не плавная деградация, при которой ответ просто становится чуть менее точным — это резкий обрыв. Пока r <= 1/2, ошибки на каждом шаге умножаются на множитель, не превышающий единицы, и потому затухают. Стоит перешагнуть порог — и множитель становится больше единицы по модулю, а значит, возмущение домножается само на себя шаг за шагом, то есть растёт экспоненциально. Каждый шаг увеличивает амплитуду в одно и то же число раз, и уже через десяток-другой шагов скромная ошибка округления раздувается до астрономических величин — ровно это мы и видели при r = 0.6, где max|u| за восемь шагов почти удвоился и продолжал расти. Поэтому CFL — не рекомендация «для точности», а жёсткая граница работоспособности: чуть ниже неё схема живёт, чуть выше — гарантированно взрывается.

Как работает под капотом

Важнейшее практическое следствие CFL: измельчая сетку, нельзя бесплатно уточнять решение. Для теплопроводности r = α·dt/dx^2 <= 1/2 означает dt <= dx^2/(2α). Уменьшили dx вдвое ради точности — будьте добры уменьшить dt вчетверо, иначе схема взорвётся. А значит, шагов по времени станет вчетверо больше, и каждый — по вдвое большему числу узлов: вычисления растут стремительно. Эта "тирания" квадратичной связи dt и dx — главная причина, по которой для жёстких параболических задач изобрели неявные схемы, не имеющие ограничения на шаг (но требующие на каждом шаге решать систему уравнений).

Стоит подчеркнуть: неустойчивость — это свойство именно схемы, а не уравнения. Точное решение уравнения теплопроводности всегда плавно затухает, никаких осцилляций в нём нет. Растущие зубцы при r = 0.6 — чистый артефакт дискретизации, ложь, которую породил неудачно выбранный шаг. Поэтому, увидев в численном эксперименте растущие знакопеременные колебания там, где физика обещает затухание, первым делом проверяйте число Куранта.

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

Главная ошибка — гнаться за точностью, измельчая только пространственную сетку и забыв пропорционально уменьшить шаг по времени. Результат — внезапный взрыв ранее работавшего расчёта.

Вторая ошибка — принять численные осцилляции за реальную физику. Знакопеременная рябь, растущая по амплитуде в задаче диффузии, — почти всегда признак нарушенного CFL, а не интересный эффект.

Третья — думать, что неявные схемы "всегда лучше". Они снимают ограничение на устойчивость, но платят за это решением системы уравнений на каждом шаге и могут размывать резкие детали. Выбор схемы — всегда компромисс.

Итоги

• Явные схемы устойчивы лишь при выполнении условия CFL; иначе решение взрывается осцилляциями.
• Для теплопроводности порог r = α·dt/dx^2 <= 1/2; при r=0.4 max|u| затухает, при r=0.6 растёт до 1.936 за 8 шагов.
• Порог 1/2 берётся из требования неотрицательности веса (1 − 2r) у самого узла в формуле обновления.
• Физический смысл: за шаг по времени информация не должна перескакивать больше одной клетки сетки.
• Для волнового уравнения условие — число Куранта C = c·dt/dx <= 1; квадратичная связь dt и dx делает явные схемы дорогими и мотивирует неявные.