Уравнение теплопроводности и метод конечных разностей
Тепло течёт от горячего к холодному, и скорость этого течения определяется кривизной температурного профиля — именно это и зашифровано в уравнении теплопроводности.
Метод конечных разностей — способ численно решать дифференциальные уравнения, заменяя производные разностями значений функции в узлах сетки.
Уравнение теплопроводности — простейшее и самое важное параболическое ДУЧП. В одномерном случае оно выглядит так: du/dt = α·d^2u/dx^2. Здесь u(x, t) — температура в точке x в момент t, а α — коэффициент температуропроводности, описывающий, насколько легко тепло проходит через материал. Левая часть — скорость изменения температуры в данной точке. Правая — кривизна профиля, умноженная на α.
Физический смысл уравнения
Чтобы прочувствовать уравнение, посмотрите на знак второй производной. Если в некоторой точке профиль образует острый горб (выпуклость вверх), то d^2u/dx^2 там отрицательна, значит du/dt отрицательна — точка остывает. И наоборот: на дне впадины кривизна положительна, du/dt положительна — точка нагревается. Природа сглаживает: горбы оседают, ямы заполняются.
Эту же мысль можно выразить через соседей. Вторая разность u(x+h) − 2u(x) + u(x−h) — это, по сути, сравнение значения в точке со средним её соседей. Если точка холоднее среднего соседей, она нагревается; если горячее — остывает. Тепло перетекает от тёплого к холодному, стремясь всё уравнять. Скорость этого выравнивания тем выше, чем сильнее точка отличается от окружения, то есть чем больше кривизна.
Дискретизация: сетка по пространству и времени
Чтобы решать уравнение численно, разобьём стержень на узлы. Пусть узлы расположены через шаг dx, и их значения хранятся в списке. Время тоже разобьём на шаги dt. Обозначим значение в узле i на временном слое n как u_i^n. Теперь заменим производные конечными разностями.
Производную по времени заменяем разностью вперёд: du/dt ≈ (u_i^{n+1} − u_i^n) / dt. Вторую производную по пространству — центральной разностью: d^2u/dx^2 ≈ (u_{i+1}^n − 2·u_i^n + u_{i-1}^n) / dx^2. Подставив это в уравнение и выразив новое значение, получаем явную схему:
u_i^{n+1} = u_i^n + r * (u_{i+1}^n - 2*u_i^n + u_{i-1}^n), где r = α*dt/dx^2Эта формула — рецепт обновления. Новое значение в узле равно старому плюс поправка, пропорциональная отклонению от соседей. Параметр r собирает в себе всю физику и параметры сетки. Схема называется явной, потому что новое значение каждого узла вычисляется напрямую через уже известные старые значения — без решения системы уравнений.
Расплывание пика: запускаем на практике
Возьмём стержень из 7 узлов, концы которого держим при нуле (краевые условия). В начальный момент поместим в центр единичный пик: [0, 0, 0, 1, 0, 0, 0]. Возьмём r = 0.4 (почему именно меньше 0.5 — тема следующего урока) и сделаем несколько шагов. Проследим, как острый пик расползается в пологий холм.
def step(u, r):
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
u = [0, 0, 0, 1.0, 0, 0, 0]
r = 0.4
for s in range(4):
print('шаг', s, [round(x, 3) for x in u])
u = step(u, r)
print('шаг 4', [round(x, 3) for x in u])
Вывод:
шаг 0 [0, 0, 0, 1.0, 0, 0, 0] шаг 1 [0.0, 0.0, 0.4, 0.2, 0.4, 0.0, 0.0] шаг 2 [0.0, 0.16, 0.16, 0.36, 0.16, 0.16, 0.0] шаг 3 [0.0, 0.096, 0.24, 0.2, 0.24, 0.096, 0.0] шаг 4 [0.0, 0.115, 0.166, 0.232, 0.166, 0.115, 0.0]
Проследим первый шаг руками для центрального узла (i=3): new = 1.0 + 0.4·(0 − 2·1.0 + 0) = 1.0 + 0.4·(−2.0) = 1.0 − 0.8 = 0.2. Для соседнего узла i=2: new = 0 + 0.4·(0 − 0 + 1.0) = 0.4. Совпадает с выводом. Видно главное: пик высотой 1.0 за один шаг просел до 0.2, а его энергия растеклась к соседям. К четвёртому шагу вместо острой иглы — широкий пологий бугор. Это и есть диффузия в действии.
Как работает под капотом
Обратите внимание на функцию step: она строит новый список из старого, и это принципиально. Если бы мы обновляли список на месте, то, посчитав узел i, мы бы испортили значение u[i], которое нужно для расчёта узла i+1. Все узлы должны обновляться одновременно, исходя из одного и того же старого слоя. Поэтому мы копируем список и пишем результат в копию.
Концы списка мы каждый шаг принудительно зануляем — это краевое условие Дирихле: на границах температура фиксирована (равна нулю, концы погружены в лёд). Цикл идёт от 1 до len(u)−1, то есть внутренние узлы; крайние не пересчитываются формулой, а удерживаются условием.
Интересно проследить за суммой значений. Первые два шага сумма строго равна 1.0 — энергия сохраняется, тепло лишь перераспределяется. Но как только волна тепла доходит до предпоследних узлов, через занулённые концы начинается утечка: на третьем шаге сумма падает до 0.872. Физически это значит, что тепло уходит в ледяные концы. Если бы концы были теплоизолированы (другое краевое условие), сумма сохранялась бы дольше.
ASCII-визуализация профиля
Чтобы увидеть расплывание глазами, нарисуем профиль столбиками: длина бара пропорциональна значению.
def step(u, r):
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
u = [0, 0, 0, 1.0, 0, 0, 0]
r = 0.4
for s in range(3):
print('шаг', s)
for x in u:
print(' %5.3f | %s' % (x, int(x * 40) * '#'))
u = step(u, r)
Вывод:
шаг 0 0.000 | 0.000 | 0.000 | 1.000 | ######################################## 0.000 | 0.000 | 0.000 | шаг 1 0.000 | 0.000 | 0.400 | ################ 0.200 | ####### 0.400 | ################ 0.000 | 0.000 | шаг 2 0.000 | 0.160 | ###### 0.160 | ###### 0.360 | ############## 0.160 | ###### 0.160 | ###### 0.000 |
Высокий одинокий столбик первого шага к третьему превращается в широкое плато — наглядная картина того, как тепло размывается по стержню.
Частые ошибки
Самая коварная ошибка — обновлять список на месте (new = u, а потом писать в u). Тогда расчёт правого узла использует уже изменённое значение левого, и схема даёт неверный результат. Всегда работайте с копией старого слоя.
Вторая ошибка — забыть переустановить краевые условия после шага. Формула не трогает концы, но если в начале вы их задали нулём, а потом случайно изменили, профиль поплывёт неправильно.
Третья — взять слишком большой r и удивиться, что вместо плавного остывания получаются растущие осцилляции. Это не баг кода, а нарушение условия устойчивости, которому посвящён отдельный урок.
Итоги
- Уравнение теплопроводности du/dt = α·d^2u/dx^2: скорость изменения температуры равна кривизне профиля.
- Метод конечных разностей заменяет производные разностями значений в узлах сетки.
- Явная схема: u_i^{n+1} = u_i^n + r·(u_{i+1} − 2u_i + u_{i-1}), где r = α·dt/dx^2.
- Все узлы обновляются одновременно — обязательно через копию старого слоя.
- Острый пик за несколько шагов расплывается в пологий холм; через занулённые концы тепло утекает.