Уравнение теплопроводности. Продольно-поперечная схема.

[MJiri]

Стоит задача решить задачу дирехле для уравнения теплопроводности (в прямоугольной области). Насколько я знаю, продольно-поперечная схема является безусловно устойчивой, значит должна сходиться при любом шаге по времени или координате. Я попытался реализовать алгоритм, в итоге никакой безусловной устойчивости не получается. Ведёт себя, как явная схема: если уменьшить шаг по координатам, то, без уменьшения шага по времени, решение расходится. С чем это может быть связано? Я думал, может просто ошибка машинном округлении накапливается.

Если что, вот питоний код моей реализации. Функция принимает значения на одном временном слое и находит значения на следующем.

Используется синтаксис Python

#
# T - шаг по времени
# hx = hy шаги по координатам
#
def solveLayer(vv):
    for k in range(1, Nxy + 1): # Первый полушаг по времени
#Заполнение коэффициентов СЛУ для неявной схемы. Крайние элементы VV всегда  = 0 (гран условия)
        a = [0 if i == 1 else - 0.5 * T / (hx ** 2) for i in range(1, Nxy + 1)]
        b = [T / (hx ** 2) + 1] * Nxy
        c = [0 if i == Nxy else -0.5 * T / (hx ** 2) for i in range(1, Nxy + 1)]
        d = [vv[i][k] * (1 - T / (hy ** 2)) + (vv[i][k + 1] + vv[i][k - 1]) * 0.5 * T / (hy ** 2) for i in
             range(1, Nxy + 1)]
        for j in range(1, Nxy):
            t = a[j] / b[j - 1]
            a[j] = 0
            b[j] -= t * c[j - 1]
            d[j] -= t * d[j - 1]
        vv[Nxy][k] = d[Nxy - 1] / b[Nxy - 1]
        for i in range(Nxy - 1, 0, -1):
            vv[i][k] = 1 / b[i - 1] * (d[i - 1] - c[i - 1] * vv[i + 1][k])

    for k in range(1, Nxy + 1): # Второй полушаг
        a = [0 if i == 1 else - 0.5 * T / (hx ** 2) for i in range(1, Nxy + 1)]
        b = [T / (hx ** 2) + 1] * Nxy
        c = [0 if i == Nxy else -0.5 * T / (hx ** 2) for i in range(1, Nxy + 1)]
        d = [vv[k][i] * (1 - T / (hy ** 2)) + (vv[k + 1][i] + vv[k - 1][i]) * 0.5 * T / (hy ** 2) for i in
             range(1, Nxy + 1)]
        for j in range(1, Nxy):
            t = a[j] / b[j - 1]
            a[j] = 0
            b[j] -= t * c[j - 1]
            d[j] -= t * d[j - 1]
        vv[k][Nxy] = d[Nxy - 1] / b[Nxy - 1]
        for i in range(Nxy - 1, 0, -1):
            vv[k][i] = 1 / b[i - 1] * (d[i - 1] - c[i - 1] * vv[k][i + 1])
    return vv
 

Pphantom:

Очень полезно пользоваться подсветкой кода, благо она есть. Поправил.

[ewert]

MJiri в сообщении #1269322 писал(а):

Я думал, может просто ошибка машинном округлении накапливается.

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

Если бы речь шла о задаче Неймана (скажем), то эффект мог бы объясняться неудачной аппроксимакцией граничных условий. Но для задачи Дирихле в прямоугольнике никакие артефакты невозможны. Какая-то ошибка в коде.

[worm2]

В код совершенно неохота вникать.
Но у вас должен выполняться "принцип максимума": значение во внутренней точке на следующем временном слое является средним из значений на предыдущем (известных) и соседних значений (неизвестных). Что-то типа такого:

$$T_{ij}^{k+0.5}=\alpha_1 T_{ij}^k+\beta_1 T_{i-1,j}^{k+0.5}+\gamma_1 T_{i+1,j}^{k+0.5},$$

$$T_{ij}^{k+1}=\alpha_2 T_{ij}^{k+0.5}+\beta_2 T_{i,j-1}^{k+1}+\gamma_2 T_{i,j+1}^{k+1},$$

где

$$\alpha_1+\beta_1+\gamma_1 = \alpha_2+\beta_2+\gamma_2 = 1,$$

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

P.S. Я всё-таки код глянул одним глазком. По-моему, вы напутали с hx и hy в коэффициентах a, b, c и d: в первой системе везде должно быть hx, а во второй — hy. Ну или наоборот.

[MJiri]

worm2 в сообщении #1269389 писал(а):

P.S. Я всё-таки код глянул одним глазком. По-моему, вы напутали с hx и hy в коэффициентах a, b, c и d: в первой системе везде должно быть hx, а во второй — hy. Ну или наоборот.

Есть такое, но hx = hy, так что это не влияет никак.

[MJiri]

Спасибо всем, нашёл ошибку. Алгоритм реализован верно, проблема была в том, что коэффициент d завист от всех соседних значений температуры на предудущем шаге, а я перезаписывал эти значения на первом же полушаге. Запись промежуточных результатов в ещё однин массив всё исправила.