Конечные разности и метод прогонки

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

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

Сетка и разностная замена производной

Разобьём отрезок [a, b] на N равных частей точками t_0=a, t_1, ..., t_N=b с шагом h=(b-a)/N. Значения искомой функции в узлах обозначим y_0, y_1, ..., y_N. Краевые условия дают нам крайние значения сразу: y_0=A и y_N=B известны. Неизвестны только внутренние узлы y_1, ..., y_{N-1}.

Теперь главный трюк. Вторую производную в узле i заменим симметричной разностью через трёх соседей:

y''(t_i) ≈ ( y_{i-1} - 2*y_i + y_{i+1} ) / h^2

Эта формула имеет второй порядок точности: ошибка убывает как h^2. Геометрически она измеряет «кривизну» — насколько средний узел провисает относительно своих соседей. Подставим её в уравнение. Возьмём линейную задачу y''=f(t):

( y_{i-1} - 2*y_i + y_{i+1} ) / h^2  =  f(t_i)

умножим на h^2:
y_{i-1} - 2*y_i + y_{i+1}  =  h^2 * f(t_i)

И вот ключевой момент: это уравнение записывается для каждого внутреннего узла i от 1 до N-1. Получается система линейных уравнений, где неизвестные — внутренние узлы, а каждое уравнение связывает только трёх соседей. Такая матрица называется трёхдиагональной: ненулевые элементы стоят только на главной диагонали и двух соседних с ней.

Структура матрицы (× — ненулевые элементы):

  | ×  ×  0  0  0 |
  | ×  ×  ×  0  0 |
  | 0  ×  ×  ×  0 |
  | 0  0  ×  ×  × |
  | 0  0  0  ×  × |

Метод прогонки (алгоритм Томаса)

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

Система: a_i*x_{i-1} + b_i*x_i + c_i*x_{i+1} = d_i
  a_i — поддиагональ, b_i — диагональ,
  c_i — наддиагональ, d_i — правая часть.

Прямой ход (модифицируем коэффициенты):
  c'_0 = c_0 / b_0
  d'_0 = d_0 / b_0
  для i = 1..n-1:
     m   = b_i - a_i * c'_{i-1}
     c'_i = c_i / m
     d'_i = (d_i - a_i * d'_{i-1}) / m

Обратный ход:
  x_{n-1} = d'_{n-1}
  для i = n-2..0:
     x_i = d'_i - c'_i * x_{i+1}

Реализация на чистом Python

Решим ту же задачу, что и стрельбой, чтобы сравнить: y''=2, y(0)=0, y(1)=0, точное y=t^2-t. Возьмём N=4 (шаг h=0.25), внутренних узлов три. Для y''=const разностная формула точна в узлах, поэтому численное решение должно совпасть с аналитическим до машинной точности. Прогонку напишем вручную на списках, без numpy.

def thomas(a, b, c, d):
    """Решение трёхдиагональной СЛАУ методом прогонки.
    a — поддиагональ (a[0] не используется),
    b — диагональ, c — наддиагональ (c[-1] не используется),
    d — правая часть. Возвращает список решений."""
    n = len(d)
    cp = [0.0] * n   # модифицированная наддиагональ
    dp = [0.0] * n   # модифицированная правая часть
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m if i < n - 1 else 0.0
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = [0.0] * n
    x[n - 1] = dp[n - 1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x

# Сначала проверим прогонку на маленькой системе 3x3 с известным ответом:
#   2x0 -  x1        = 1
#  -x0 + 2x1 -  x2   = 0
#        -x1 + 2x2   = 1
# Ответ (по симметрии): x0 = x1 = x2 = 1
test = thomas([0.0, -1.0, -1.0],   # поддиагональ
              [2.0,  2.0,  2.0],   # диагональ
              [-1.0, -1.0, 0.0],   # наддиагональ
              [1.0,  0.0,  1.0])   # правая часть
print("проверка 3x3:", [round(v, 4) for v in test])

# Теперь сама краевая задача y'' = 2, y(0)=0, y(1)=0
N = 4
h = 1.0 / N
rhs = 2.0 * h * h            # h^2 * f, где f = 2

# Уравнение узла: y_{i-1} - 2 y_i + y_{i+1} = h^2 * 2
# Крайние y_0=0 и y_4=0 переносим в правую часть.
A = [0.0,  1.0,  1.0]        # поддиагональ
Bd = [-2.0, -2.0, -2.0]      # диагональ
C = [1.0,  1.0,  0.0]        # наддиагональ
D = [rhs - 0.0, rhs, rhs - 0.0]   # y_0=0 и y_4=0 -> вычитаем 0

sol = thomas(A, Bd, C, D)
nodes = [0.25, 0.5, 0.75]
print("h =", h)
for t, y_num in zip(nodes, sol):
    y_exact = t * t - t
    print("t =", t, " численно =", round(y_num, 6),
          " точно =", round(y_exact, 6))

Вывод:

проверка 3x3: [1.0, 1.0, 1.0]
h = 0.25
t = 0.25  численно = -0.1875  точно = -0.1875
t = 0.5  численно = -0.25  точно = -0.25
t = 0.75  численно = -0.1875  точно = -0.1875

Тест 3x3 дал ровно [1, 1, 1] — прогонка работает. А краевая задача решена в точку: численные значения во всех внутренних узлах совпали с аналитическими y=t^2-t. Это не случайность: для уравнения с постоянной правой частью разностная аппроксимация второй производной не вносит ошибки в узлах.

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

Почему конечные разности приводят именно к трёхдиагональной матрице? Потому что разностная формула второй производной связывает узел только с двумя ближайшими соседями — y_{i-1}, y_i, y_{i+1}. Никакой узел не «дотягивается» дальше соседа, поэтому в каждой строке матрицы ровно три ненулевых элемента, выстроенных в полосу. Это прямое следствие локальности разностного шаблона.

Метод прогонки — это попросту метод Гаусса, аккуратно адаптированный под эту полосу. Обычный Гаусс на матрице размера n тратит порядка n^3 операций; прогонка — порядка n, потому что не трогает заведомо нулевые элементы вне трёх диагоналей. Прямой ход исключает поддиагональ, превращая систему в двухдиагональную, обратный ход выражает неизвестные одну за другой снизу вверх. Алгоритм устойчив, если матрица обладает диагональным преобладанием — а для разностной аппроксимации второй производной (диагональ -2, соседи по +1) это как раз выполняется, поэтому прогонка не накапливает ошибок.

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

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

Забыть перенести краевые значения в правую часть. В уравнениях у границы (i=1 и i=N-1) фигурируют известные y_0 и y_N. Их нужно вычесть из правой части, а не оставлять как неизвестные. Иначе система собрана неверно.
Перепутать диагонали при сборке. Поддиагональ a, диагональ b и наддиагональ c легко перепутать местами, особенно индексацию a[0] и c[-1], которые не используются. Всегда тестируйте прогонку на маленькой системе 3x3 с известным ответом.
Делить на ноль в прямом ходе. Если m=b_i - a_i*cp[i-1] обращается в ноль (нет диагонального преобладания), прогонка ломается. Для классических краевых задач этого не происходит, но для произвольных коэффициентов проверяйте.
Ожидать точного совпадения для произвольного f. Совпадение «в точку» здесь — особенность постоянной правой части. Для общего f(t) разностная схема даёт ошибку порядка h^2: измельчайте сетку, чтобы её уменьшить.

Итоги

Конечные разности покрывают отрезок сеткой и заменяют y'' симметричной разностью (y_{i-1} - 2*y_i + y_{i+1})/h^2.
Уравнение, записанное для каждого внутреннего узла, даёт трёхдиагональную СЛАУ.
Её решает метод прогонки (алгоритм Томаса) за линейное время — прямой и обратный ход.
Краевые значения переносятся в правую часть уравнений у границы.
В отличие от стрельбы, все узлы находятся одновременно; метод устойчив там, где прямое интегрирование расходится.