Параболический тип. Общие принципы метода конечных разностей#
Явный эволюционный во времени процесс. См. характеристики для понимания. Нулевое собственное число отвечает времени.
Явная зависимость от времени.
Решаются как при заданных начальных условиях, так и при краевых.
«Негладкости» и «разрывы» возможны, но сглаживаются со временем.
Примеры:
Уравнение теплопроводности (диффузии)
Уравнение анизотропной диффузии (квазилинейное - страшно, да?)
Пример решения параболической задачи#
Определение. Задача уравнения диффузии (теплопроводности) ставится следующим образом:
В линейном случае \(a=a(x, t)>0\). В квазилинейном случае \(a=a(u, x, t)>0\) и \(f=f(u, x, t)\).
Для того, чтобы задача была поставлена корректно, необходимо задать начальное условие по времени:
и граничные условия (в одном из трёх родов) по координате:
Условия первого рода:
и условия второго рода:
Условия третьего рода:
Метод конечных разностей основан на дискретизации функции и всех производных, как мы уже делали это раньше в одномерном случае. Только теперь дискретная функция - это какая-то многомерная матрица \(u_m^n \stackrel{ideal\, world}{=} u(t_n, x_m)\). Здесь и далее индекс сверху всегда обозначает время, а индексы снизу координаты.
Схемы для уравнения теплопроводности с постоянными коэффициентами#
Приведём некоторые разностные схемы (дискретизации), работающие для константного \(a\). Каждой дискретизации можно сопоставить картинку (шаблон), показывающую, какие конкретно узлы \(u_m^n\) связаны уравнением.
Явная схема#
Устойчива при \(\sigma=\frac{a \tau}{h^2} \leq \frac{1}{2}\) и порядок аппроксимации \(O\left(\tau, h^2\right)\).
Неявная схема#
Устойчива при при любых \(\tau, h\), порядок аппроксимации \(O\left(\tau, h^2\right)\).
Смешанная шеститочечная схема#
Где вес \(\xi \in[0,1]\) - заранее задаваемое нами число - отвечает за смесь явной и неявной производных по координате. Определяется эмпирически или ищется наиболее оптимальное значение \(\xi\) можно отыскать из «максимизации» устойчивости (об этом позже).
Схема Кранка-Никольсон - смешанная схема при \(\xi =\frac{1}{2}\)#
Схема устойчива при любых шагах \(\tau, h\) и имеет порядок аппроксимации \(O\left(\tau, h^2\right)\). Эта схема, в отличие от двух предыдущих, не является монотонной, т.е. она может давать осцилляции разностного происхождения на решениях, имеющих большие градиенты.
Трёхслойная схема - вводим «смесь» явной и неявной схемы также для производной по времени#
При \(\eta=0,5, \quad \xi=1 \quad\) ее порядок аппроксимации равен \(O\left(\tau^2, h^2\right)\). Недостатком схемы является более высокий порядок разностного уравнения по времени, чем порядок исходного дифференциального, следовательно, необходимость ставить дополнительное условие на \(u_t^{\prime}(0, x)\).
Примечание. Неявные схемы решаются путём построения СЛАУ на все пространственные точки в новом слое.
Примечание. Эти схемы легко перенести на многомерный случай или на любое другое параболическое уравнение. Главное помнить, что они «хорошо» работают только при постоянных коэффициентах. Если коэффициенты зависят от времени и координат, то их строят исходя из нужных законов сохранения.
Следующий анимация демонстрирует процесс решения разностной схемы на шаблоне Кранка-Николсон.
#!pip install gdown
import gdown
from IPython.display import HTML
HTML(gdown.download('https://drive.google.com/uc?id=1Co4v72cPOzWH6BESds1OSAjq2KByZvCS', 'anim.html', quiet=True))
Черные точки - это то, что мы уже знаем (граничные условия + найденные). Белые - ещё не знаем.
Видим, что в процессе решения шаблон как бы идёт по сетке, касаясь как черных точек, так и белых (которые после этого становятся зелёными). Каждый кадр соответствует записыванию уравнения, в котором есть как известные (чёрные) точки, так неизвестные.
Когда схема доходит до правого края, мы имеем 8 линейных уравнений на 8 неизвестных (все зелёные точки) - эта система имеет решение, а значит мы можем найти все точки на данной слою.
Как-то так и работает метод конечных разностей.
Примечание. Из этой визуализации понятно, как выбирать граничные условия. Если бы, например, не было чёрных точек справа (это соответствует знанию нами \(u(t, 1)\) для любого момента времени), то мы бы никогда не получили \(N\) уравнений на \(N\) неизвестных. Тем не менее, схему всё равно можно было бы запустить, зная 2 слоя слева (это соответствует знанию граничных условий \(u(t, 0)\) и \(\cfrac{\partial u}{\partial x}(t, 0)\)).
Примечание. Заметим, что если бы использовали явную схему, то не пришлось бы возиться с составлением и решением системы. Следующий пример демонстрирует простоту реализации явного метода.
# Пример явной схемы f = 0
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
# Parameters
a = 1.0
L = 1.0
T = 0.2
Nx = 100
Nt = 100000
h = L/(Nx-1)
tau = T/(Nt-1)
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)
# Initial and boundary conditions
def u0(x):
return 1 + np.sin(2*np.pi*x/L)**2 + np.sin(np.pi*x/L)
def ul(t):
return 1
def ur(t):
return 1
# Finite difference method
r = a*tau/h**2
u = np.zeros((Nx, Nt))
u[:, 0] = u0(x)
for n in range(Nt-1):
u[0, n+1] = ul(t[n+1])
u[-1, n+1] = ur(t[n+1])
for m in range(1, Nx-1):
u[m, n+1] = u[m, n] + r*(u[m-1, n] - 2*u[m, n] + u[m+1, n])
# Plot every 100th timestep
fig, ax = plt.subplots()
for n in range(0, Nt, 5000):
color = cm.jet(n/Nt)
ax.plot(x, u[:, n], color=color, label=f"t={t[n]:.2f}")
ax.set_xlabel("x")
ax.set_ylabel("u")
norm = plt.Normalize(0, T)
cmap = cm.jet
sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
#cbar = fig.colorbar(sm)
#bar.ax.set_ylabel("Time")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 3D surface plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev=40, azim=40)
X, Y = np.meshgrid(x, t[::100])
Z = u.T[::100, :]
surf = ax.plot_surface(X, Y, Z, cmap='coolwarm')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
ax.set_title('Solution of the Heat Equation')
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
Как составлять и решать СЛАУ для неявных шаблонов#
У явных схем помимо простоты есть ещё одна отличительная особенность - они почти никогда не работают. Надо научиться решать УРЧП неявными схемами.
Рассмотрим, например, схему Кранка-Никольсон
Хотим её реализовать. Пусть также заданы начальные и граничные условия:
Хотим записать нашу разностную задачу в каноническом виде (для фиксированных \(n\) и \(n+1\) - по ним здесь суммирования не идёт)
где \(A, B\) и \(b\) пишутся исходя из вида уравнения. Перенесём все слагаемые с \(n+1\) слоя влево:
Введём так называемое число Куранта \(r=\frac{a \tau}{2h^2}\). Тогда уравнение:
Откуда видим выражения на искомые операторы (c учётом граничных условий):
Оператор \(A_{s, m}\):
Оператор \(B_{s, m}\) (при более сложных граничных условиях первая и последняя строчки могут быть ненулевыми):
Столбец \(b_s^n\):
На каждом шаге можно использовать np.linalg.solve относительно неизвестного столбца \(u_s^{n+1}\), но это плохо-обусловленная задача, т.е. при большой размерности ошибка решения СЛАУ станет слишком большой.
Эффективнее использовать алгоритм прогонки.
Обозначим всю известную правую часть за столбец \(F_s^n = B_{sm} u_m^n+b_s^n\). Тогда решаем систему на \(u_m^{n+1}\):
с трёхдиагональной матрицей \(A\). Обозначим далее для простоты формул \(x_i = u_i^{n+1}\).
Эта система равносильна следующим формулам (построчные равенства, там где вылезли за границы просто присваиваем нули):
Ищем решение в виде реккуренты (т.е. перепишем задачу в терминах поиска коэффициентов \(\alpha_i, \beta_i\) - как только они найдены, применяя эту же формулу найдём реккурентно все \(x_i\)):
При этом \(\alpha_i, \beta_i\) не должны зависеть искомых \(x_i\). Начальные члены этих последовательностей можно определить нулевыми \(\alpha_0=\beta_0=0\).
Используя это соотношение, выразим \(x_{i-1}\) и \(x_i\) через \(x_{i+1}\) и подставим в изначальное уравнение:
Беря частную производную слева от \(x_{i+1}\) понимаем, что множитель в скобках должен равняться нулю. А тогда и всё остальное в уравнении должно зануляться:
Отсюда следуют явно заданные реккуренты на \(\alpha_i, \beta_i\)
С начальным условием (из \(\alpha_0=\beta_0=0\)):
Теперь, когда мы знаем \(\alpha_i, \beta_i\), легко реккурентно посчитать искомый вектор неизвестных:
Show code cell source
# Пример схемы Кранка-Никольсон с f = 0 методом прогонки
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from numba import njit
from tqdm import tqdm
def TDMA(A,B,C, F): # Прогонка
# A, B, C, F - просто числовые последовательности, НЕ матрицы
A = np.insert(A, 0, 0)
n = len(F)
x = np.zeros(n)
alpha = np.zeros(n)
beta = np.zeros(n)
alpha[1] = -C[0]/B[0]
beta[1] = F[0]/B[0]
for i in range(2, n):
alpha[i] = -C[i-1]/(A[i-1]*alpha[i-1] + B[i-1])
beta[i] = (F[i-1] - A[i-1]*beta[i-1])/(A[i-1]*alpha[i-1] + B[i-1])
x[n-1] = (F[n-1] - A[n-1]*beta[n-1])/(B[n-1] + A[n-1]*alpha[n-1])
for i in range(n-1, 0, -1): # помним, что последнее значение не кушает
x[i - 1] = alpha[i]*x[i] + beta[i]
return x
# Parameters
a = 1.0
L = 1.0
T = 0.2
Nx = 100
Nt = 100 # Как же мало нам понадобилось шагов по времени!
h = L/(Nx-1)
tau = T/(Nt-1)
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)
# Initial and boundary conditions
def u0(x):
return 1 + np.sin(2*np.pi*x/L)**2 + np.sin(np.pi*x/L)
def ul(t):
return 1
def ur(t):
return 1
# Define the source term
def f(t, x):
return 0
# Initialize the solution matrix
u = np.zeros((Nt, Nx))
# Set the initial condition
u[0, :] = u0(x)
# Set the boundary conditions
u[:, 0] = ul(t)
u[:, -1] = ur(t)
# Set up the matrix for the Crank-Nicolson scheme
r = a*tau/(2*h**2)
A = np.zeros((Nx, Nx))
for i in range(1, Nx-1): # ПОмним, что последний индекс в цикле НЕ кушается
A[i, i] = 1 + 2*r
A[i, i-1] = -r
A[i, i+1] = -r
A[0, 0] = 1
A[-1, -1] = 1
B = np.zeros((Nx, Nx))
for i in range(1, Nx-1): # ПОмним, что последний индекс в цикле НЕ кушается
B[i, i] = 1 - 2*r
B[i, i-1] = r
B[i, i+1] = r
# Time stepping with the Crank-Nicolson scheme
for n in tqdm(range(1, Nt)):
b = np.zeros(Nx)
b[0] = ul(t[n])
b[-1] = ur(t[n])
b[1: -1] = tau*f(t[n-1], x[1:-1])
#print(b)
F = B @ u[n-1] + b
#u[n, :] = np.linalg.solve(A, F) # Не сработает для большого Nx
u[n, :] = TDMA(A.diagonal(offset = -1),
A.diagonal(offset = 0),
A.diagonal(offset = 1),
F)
# Plot every 100th timestep
fig, ax = plt.subplots()
for n in range(0, Nt, 1):
color = cm.jet(n/Nt)
ax.plot(x, u[n, :], color=color, label=f"t={t[n]:.2f}")
ax.set_xlabel("x")
ax.set_ylabel("u")
norm = plt.Normalize(0, T)
cmap = cm.jet
sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
#cbar = fig.colorbar(sm)
#cbar.ax.set_ylabel("Time")
plt.show()
100%|█████████████████████████████████████████████████████████████████████████████████| 99/99 [00:00<00:00, 7082.68it/s]
Построение разностных схем для УРЧП с переменными коэффициентами (Консерватизм)#
Определение. Схема называется консервативной, если в дискретизованной задаче выполняются те же законы сохранения (первые интегралы), что и в дифференциальной. Вспомните про количество независимых первых интегралов в системе ДУ.
Как правило, если коэффициенты в линейном УРЧП зависят от координат \(a=a(x, t)\), то для обеспечения консерватизма системы необходимо строить разностную задачу исходя из требуемых законов сохранения
На примере того же однородного уравнения диффузии (теплопроводности)
Его законом сохранения (первым интегралом) является уравнение непрерывности на тепловой поток \(W=-a(t, x) \frac{\partial u}{\partial x}\):
или, в интегральной форме,
Любое из этих уравнений должно быть выполнено всегда и где угодно. Удобно взять самое последнее и обеспечить его для самого маленького контура в системе (для любого другого контура оно будет выполнено автоматически по линейности):
Это же можно получить из дифференциальной формы закона сохранения
Здесь взяли промежуточные узлы на \(W\), т.к. такой приём, как правило, даёт лучшую устойчивость при переменных коэффициентах (напоминаю \(W=-a(t, x) \frac{\partial u}{\partial x}\)):
Подставляя это в дискретизованный закон сохранения, имеем нашу консервативную численную схему:
Она получилось шеститочечной.
Примечание. В случае неоднородности просто добавляем её. Опять же, чем на бОльших узлах зададим, тем лучше, но не очень критично, если честно.
Квазилинейный случай#
Рассмотрим случай зависимости коэффициентов от от искомой функции \(a=a(u)>0\) и \(f=f(u)\) - т.е. простая квазилинейность. При зависимости этих коэффициентов ещё и от \(x, t\) обобщение тривиально на основе предыдущего пункта.
Т.к. коэффициенты переменные, то опять-таки используем промежуточные узлы.
Явная схема с нелинейностью на нижнем слое
где \(a_{m+1 / 2}^n\) может быть вычислена одним из следующим способов:
\(a_{m+1 / 2}^n=\frac{1}{2}\left(a\left(u_m^n\right)+a\left(u_{m+1}^n\right)\right)\),
\(a_{m+1 / 2}^n=a\left(\cfrac{u_m^n+u_{m+1}^n}{2}\right)\),
\(a_{m+1 / 2}^n=a\left(\cfrac{2 u_m^n u_{m+1}^n}{u_m^n+u_{m+1}^n}\right)\),
\(a_{m+1 / 2}^n=\cfrac{2 a\left(u_m^n\right) a\left(u_{m+1}^n\right)}{a\left(u_m^n\right)+a\left(u_{m+1}^n\right)}\).
На следующем слое опять-таки решается СЛАУ на неизвестные коэффициенты методом прогонки (полученная система будет трёхдиагональной).
Недостаток схемы заключается в необходимости выполнения условия, ограничивающего шаг по времени \(\tau\left\|f_u^{\prime}\right\|<<1\) для обеспечения устойчивости.
Неявная схема с нелинейностью на нижнем слое
Нужна когда условие \(\tau\left\|f_u^{\prime}\right\|<<1\) невыполнимо.
Решать большую систему нелинейных уравнений ну уж очень не хочется, так что здесь, как правило, используют линеаризацию
И алгоритм тот же самый - решить СЛАУ на верхнем уровне методом прогонки.
Примечание. В случае общей квазилинейности по \(a=a(u, x, t)\) просто берём предыдущий параграф и всё, что связано с зависимостью от \(u\) считаем по формулам из этого. Шеститочечная схема, например:
где \(a_{m+1 / 2}^n\) может быть вычислена одним из следующим способов:
\(a_{m+1 / 2}^n=\frac{1}{2}\left(a\left(u_m^n, \; x_m, \;t_n\right)+a\left(u_{m+1}^n, \; x_{m+1}, \;t_n\right)\right)\),
\(a_{m+1 / 2}^n=a\left(\cfrac{u_m^n+u_{m+1}^n}{2},\; \cfrac{x_m+x_{m+1}}{2},\; t_n\right)\),
\(a_{m+1 / 2}^n=a\left(\cfrac{2 u_m^n u_{m+1}^n}{u_m^n+u_{m+1}^n},\; \cfrac{2 x_m x_{m+1}}{x_m+x_{m+1}}, \; t_n\right)\),
\(a_{m+1 / 2}^n=\cfrac{2 \cdot a_m^n \cdot a_{m+1}^n}{a_m^n+a_{m+1}^n}\).
Многомерный (по координате) случай#
Численное решение даже простейших уравнений параболического типа сильно усложняется, если в задаче имеется в наличии более одного пространственного измерения. Условие устойчивости для многомерных схем накладывает столь жесткие ограничения на шаги по времени, что расчет по ним практически невозможен (на прошлом семинаре мы как раз использовали явную схему - работало всё очень долго). Необходимо применять неявные схемы.
Пусть решаем
Неоднородностью \(f(x, y, t)\) пока пренебрежём (см. примечание в конце пункта).
Дискретизуем задачу в виде
где
Схема выглядит следующим образом:
Можно на каждом шаге решать систему на весь следующий слой из \(N^2\) уравнений на \(N^2\) неизвестных, но это не очень эффективно, т.к. аналога метода прогонки для такой задачи нет.
Можно предложить схему расщепления по направлениям, или локально одномерную схему (метод дробных шагов), где на каждом шаге по времени мы последовательно решаем две разностные задачи сначала на ~, потом на \(n+1\):
На каждом из промежуточных шагов мы сути решаем одномерную задачу для строки или столбца матрицы \(u_{ml}\) столько раз, чтобы всю её заполнить. Сначала делаем в одном направлении (по одной координате), потом в другом (по другой координате). Таким образом применим хорошо обусловленный алгоритм прогонки.
Порядок аппроксимации этой схемы \(O\left(\tau, h_x^2, h_y^2\right)\).
Примечание. Обобщение на бОльшие размерности тривиально - просто больше промежуточных слоёв.
Примечание. Порядок аппроксимации по времени легко увеличить до второго, придав слою ~ смысл \(n+\frac{1}{2}\) шага по времени:
Примечание. Если уравнение неоднородно, то правую часть легко засунуть в схему, например: