Волновое уравнение и уравнение Лапласа

Если теплопроводность всё сглаживает, то волновое уравнение заставляет возмущение бежать, а уравнение Лапласа описывает мир, в котором всё уже пришло в равновесие.

Метод Якоби — итерационный способ решения уравнения Лапласа, в котором каждый внутренний узел на новой итерации становится средним арифметическим четырёх своих соседей.

В этом уроке мы знакомимся с двумя другими типами ДУЧП. Волновое уравнение d^2u/dt^2 = c^2·d^2u/dx^2 — гиперболическое: в нём по времени стоит вторая производная (ускорение), и потому возмущение не размывается, а распространяется как бегущая волна со скоростью c. Уравнение Лапласа d^2u/dx^2 + d^2u/dy^2 = 0 — эллиптическое: в нём нет времени вообще, оно описывает установившееся, равновесное состояние, например стационарную температуру пластины или электрический потенциал.

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

Качественно решение волнового уравнения для бесконечной струны описывается формулой д'Аламбера: любой начальный профиль распадается на две одинаковые по форме волны, разбегающиеся влево и вправо со скоростью c. Не нужно решать уравнение, чтобы это увидеть, — достаточно представить, что половина «горба» уезжает в одну сторону, половина в другую, и в каждый момент полное поле есть сумма этих двух бегущих частей. Эта картина двух встречных волн объясняет и отражения, и интерференцию, и почему для расчёта одного нового слоя схеме нужны сразу два предыдущих: вторая производная по времени хранит память не только о форме, но и о направлении движения.

Волновое уравнение: схема с двумя предыдущими слоями

Поскольку в волновом уравнении по времени стоит вторая производная, для шага вперёд нужны не один, а два предыдущих временных слоя. Заменив обе вторые производные центральными разностями, получаем явную схему:

u_i^{n+1} = 2*u_i^n - u_i^{n-1} + C^2 * (u_{i+1}^n - 2*u_i^n + u_{i-1}^n),  C = c*dt/dx

Здесь C — число Куранта. Новое значение узла зависит от его значения сейчас (слой n), от значения шаг назад (слой n−1) и от кривизны текущего профиля. Чтобы стартовать, нужны два первых слоя: начальный профиль и профиль через шаг. Если струну отпустили из покоя (нулевая начальная скорость), первый слой получают по упрощённой формуле u_i^1 = u_i^0 + 0.5·C^2·(u_{i+1}^0 − 2u_i^0 + u_{i-1}^0).

Уравнение Лапласа: итерации Якоби (исполнимый пример)

Для надёжно воспроизводимого запуска возьмём именно уравнение Лапласа — у него красивое и точное свойство сходимости. Дискретизация d^2u/dx^2 + d^2u/dy^2 = 0 на квадратной сетке даёт простое условие: в каждом внутреннем узле значение равно среднему четырёх соседей (сверху, снизу, слева, справа). Метод Якоби просто многократно применяет это правило: на каждой итерации пересчитываем все внутренние узлы как среднее соседей со старой итерации, пока поле не перестанет меняться.

Возьмём сетку 5×5. Верхнюю границу нагреем до 100, остальные три границы держим при 0. Внутренние 3×3 узла стартуют с нуля. Будем итерировать и следить за центральным узлом, а в конце выведем установившееся поле.

N = 5
u = [[0.0] * N for _ in range(N)]
for j in range(N):
    u[0][j] = 100.0  # верхняя граница горячая

def jacobi(u):
    new = [row[:] for row in u]
    for i in range(1, N - 1):
        for j in range(1, N - 1):
            new[i][j] = 0.25 * (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1])
    return new

print('Сходимость центрального узла:')
for it in range(1, 7):
    u = jacobi(u)
    print('итерация %d: центр = %.4f' % (it, u[2][2]))

for it in range(200):
    u = jacobi(u)
print('Установившееся поле:')
for row in u:
    print('  ' + ' '.join('%6.2f' % x for x in row))

Вывод:

Сходимость центрального узла:
итерация 1: центр = 0.0000
итерация 2: центр = 6.2500
итерация 3: центр = 12.5000
итерация 4: центр = 15.6250
итерация 5: центр = 18.7500
итерация 6: центр = 20.3125
Установившееся поле:
  100.00 100.00 100.00 100.00 100.00
    0.00  42.86  52.68  42.86   0.00
    0.00  18.75  25.00  18.75   0.00
    0.00   7.14   9.82   7.14   0.00
    0.00   0.00   0.00   0.00   0.00

Центральный узел плавно сходится и в пределе застывает на значении 25.00. И это не случайность, а красивая физика: по симметрии центр квадрата, у которого одна сторона нагрета до 100, а три держатся при 0, в равновесии должен оказаться ровно на среднем — (100 + 0 + 0 + 0) / 4 = 25. Численный метод сам, без подсказок, нашёл этот ответ.

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

Почему правило "узел = среднее соседей" решает Лапласа? Раскроем дискретизацию. Условие new[i][j] = 0.25·(сумма соседей) эквивалентно тому, что (u_{i+1} − 2u_i + u_{i-1}) + (u_{j+1} − 2u_j + u_{j-1}) = 0, то есть сумма вторых разностей по обоим направлениям равна нулю. А это и есть дискретный Лаплас d^2u/dx^2 + d^2u/dy^2 = 0. Когда итерации сходятся и поле перестаёт меняться, каждый узел в точности равен среднему соседей — значит, дискретное уравнение Лапласа выполнено всюду.

Как и в теплопроводности, метод Якоби строит новую сетку из копии старой: все узлы обновляются одновременно от одного и того же старого состояния. Существует родственный метод Гаусса-Зейделя, где значения обновляются на месте и сразу используются в этом же проходе, — он сходится почти вдвое быстрее, но теряет свойство "чистого одновременного обновления". Для ясности мы выбрали Якоби.

У итераций есть выразительная физическая аналогия — релаксация к равновесию. Можно представить, что каждый узел потихоньку «подтягивается» к среднему своих соседей, как если бы поле было упругой плёнкой, натянутой на жёсткий каркас границы: где-то плёнка провисает, где-то выгибается, и шаг за шагом она успокаивается в форму с минимальной «энергией натяжения». Метод Гаусса-Зейделя в этой картинке — та же релаксация, но «более жадная»: свежие значения соседей он использует немедленно, и потому информация о границе быстрее проникает вглубь сетки. Гармонические функции (решения Лапласа) — это как раз такие предельно гладкие, ни о чём не «спорящие» поля, к которым релаксация и стремится.

Из правила «узел равен среднему соседей» следует красивое свойство — принцип максимума: в равновесии значение в любой внутренней точке не может быть больше, чем самое большое значение на границе, и не меньше самого маленького. Интуиция простая: если бы какой-то внутренний узел оказался строго выше всех соседей, он не был бы равен их среднему, а значит, поле ещё не пришло в равновесие. Поэтому в нашем примере все внутренние значения честно лежат между 0 и 100 — никакая внутренняя точка не «перегревается» выше нагретой границы. Это же объясняет, почему стационарную задачу так удобно использовать на практике: уравнение Лапласа описывает не только остывшую до равновесия пластину, но и электростатический потенциал в области без зарядов, и течение идеальной несжимаемой жидкости, и установившееся тепло в детали — всюду, где система уже «успокоилась» и du/dt обратилось в ноль.

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

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

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

Для метода Якоби ошибка — обновлять сетку на месте и думать, что это всё ещё Якоби. На месте — это уже Гаусса-Зейделя; результаты сходятся к тому же ответу, но скорость и промежуточные значения отличаются.

Ещё одна ошибка — пересчитывать граничные узлы. Граница в задаче Лапласа фиксирована (это краевое условие), её трогать нельзя; цикл должен идти только по внутренним узлам от 1 до N−1.

Итоги

Волновое уравнение d^2u/dt^2 = c^2·d^2u/dx^2 гиперболическое: возмущение бежит как волна; явная схема требует двух предыдущих слоёв.
Уравнение Лапласа d^2u/dx^2 + d^2u/dy^2 = 0 эллиптическое: описывает равновесие, времени в нём нет.
Дискретный Лаплас: каждый внутренний узел равен среднему четырёх соседей.
Метод Якоби многократно применяет это правило до сходимости; поле застывает в равновесии.
В примере центр квадрата сходится к 25.00 = среднему граничных значений — это подтверждает физика по симметрии.