📖 Документация Qumir

← Вернуться в Playground

← Все примеры

Уравнение теплопроводности 3D (неявная схема)

Неявная схема (Backward Euler) для трёхмерного уравнения теплопроводности. В отличие от явной схемы, неявная безусловно устойчива -- можно брать сколь угодно большие шаги по времени без риска расходимости. За это приходится платить решением системы линейных уравнений на каждом шаге, что реализовано здесь методом Гаусса--Зейделя.

Задача

Решается то же уравнение теплопроводности в кубе $[0,1]^3$:

$\frac{\partial u}{\partial t} = \Delta u$

с начальным условием $u_0(x,y,z) = \sin(\pi x) \sin(\pi y) \sin(\pi z)$ и нулевыми граничными условиями $u = 0$ на границе.

Неявная схема записывается так: вместо того чтобы вычислять $u^{n+1}$ напрямую, нужно решить систему уравнений:

$(I - dt \cdot \Delta_h) , u^{n+1} = u^n$

Здесь $I$ -- единичная матрица, $\Delta_h$ -- дискретный оператор Лапласа (семиточечный шаблон). Это большая разреженная система, которую мы решаем итерационно.

Точное решение для проверки: $u(x,y,z,t) = e^{-3\pi^2 t} \sin(\pi x) \sin(\pi y) \sin(\pi z)$.

Идея алгоритма

  1. Разбиваем куб на сетку $N \times N \times N$ с шагом $h$.
  2. Выбираем параметр $\lambda = dt / h^2 = 0.5$ -- в два раза больше, чем допускает условие устойчивости явной схемы ($\leq 1/6$). Это возможно благодаря безусловной устойчивости неявной схемы.
  3. На каждом шаге по времени копируем предыдущее решение как начальное приближение.
  4. Итерационно решаем систему $(I - dt \cdot \Delta_h) u^{n+1} = u^n$ методом Гаусса--Зейделя (метод, который обновляет значения «на месте», используя уже посчитанные соседние значения).
  5. Повторяем внутренние итерации до сходимости (изменения меньше $10^{-6}$) или до 200 итераций.
  6. Сравниваем результат с точным решением.

Разбор

Параметры и начальное условие

Параметр $\lambda = 0.5$ задаёт крупный шаг по времени. Решение хранится в четырёхмерном массиве u[время, x, y, z], что позволяет сохранить всю историю.

цел Nx, Ny, Nz
Nx := 16
Ny := 16
Nz := 16

вещ pi
pi := 3.141592653589793
вещ h, h2
h := 1.0 / (Nx - 1)
h2 := h * h

| Неявная схема (Backward Euler) для u_t = Δu в 3D: (I - dt Δ) u^{n+1} = u^n
| Итерационно решаем линейную систему методом Гаусса–Зейделя на каждом шаге по времени
вещ lam, dt, T
lam := 0.5 | lam = dt / h^2 (умеренно крупный шаг по времени, неявная схема устойчива)
dt := lam * h2
T := 0.10

Четырёхмерный массив и начальное условие

Массив u[0:nt, 0:Nx-1, 0:Ny-1, 0:Nz-1] хранит решение на всех временных слоях. Нулевой слой u[0,...] заполняется начальным условием $\sin(\pi x) \sin(\pi y) \sin(\pi z)$.

| 4D массив: u[время, x, y, z]
вещ таб u[0:nt, 0:Nx-1, 0:Ny-1, 0:Nz-1]

| Начальное условие u(0,x,y,z) = sin(pi x) sin(pi y) sin(pi z)
нц для i от 0 до Nx-1 шаг 1
  x := i * h
  нц для j от 0 до Ny-1 шаг 1
    y := j * h
    нц для k от 0 до Nz-1 шаг 1
      z := k * h
      u[0,i,j,k] := sin(pi * x) * sin(pi * y) * sin(pi * z)
    кц
  кц
кц

Итерации Гаусса--Зейделя

На каждом временном шаге решаем систему итерационно. Формула обновления:

$u_{i,j,k}^{n+1} = \frac{u_{i,j,k}^{n} + \lambda \cdot (\text{сумма шести соседей из } u^{n+1})}{1 + 6\lambda}$

Обновление происходит «на месте» (in-place): часть соседей уже содержит значения текущей итерации, а часть -- предыдущей. Это и есть суть метода Гаусса--Зейделя, в отличие от метода Якоби, где все соседи берутся с предыдущей итерации.

maxIt := 200 | максимум внутренних итераций Гаусса–Зейделя на шаг
denom := 1.0 + 6.0 * lam

нц для n от 0 до nt-1 шаг 1
  | Инициализируем начальное приближение u[n+1] := u[n]
  нц для i от 0 до Nx-1 шаг 1
    нц для j от 0 до Ny-1 шаг 1
      нц для k от 0 до Nz-1 шаг 1
        u[n+1,i,j,k] := u[n,i,j,k]
      кц
    кц
  кц

  | Внутренние итерации решения (I - dt Δ) u^{n+1} = u^n методом Гаусса–Зейделя
  нц для it от 0 до maxIt-1 шаг 1
    change := 0.0
    нц для i от 1 до Nx-2 шаг 1
      нц для j от 1 до Ny-2 шаг 1
        нц для k от 1 до Nz-2 шаг 1
          oldv := u[n+1,i,j,k]
          rhs := u[n,i,j,k]
          newv := (rhs + lam * (u[n+1,i+1,j,k] + u[n+1,i-1,j,k] + u[n+1,i,j+1,k] + u[n+1,i,j-1,k] + u[n+1,i,j,k+1] + u[n+1,i,j,k-1])) / denom
          u[n+1,i,j,k] := newv
          diff := abs(newv - oldv)
          change := max(change, diff)
        кц
      кц
    кц
    | Ранний выход по сходимости внутреннего решателя
    если change < 1e-6 то выход все
  кц
кц

Проверка результата

Сравниваем численное решение на последнем временном слое с точным: $u(x,y,z,t) = e^{-3\pi^2 t} \sin(\pi x) \sin(\pi y) \sin(\pi z)$.

t := nt * dt
decay := exp(-3.0 * pi * pi * t)
err := 0.0
нц для i от 0 до Nx-1 шаг 1
  x := i * h
  нц для j от 0 до Ny-1 шаг 1
    y := j * h
    нц для k от 0 до Nz-1 шаг 1
      z := k * h
      uex := decay * sin(pi * x) * sin(pi * y) * sin(pi * z)
      diff := abs(u[nt,i,j,k] - uex)
      err := max(err, diff)
    кц
  кц
кц

вывод "Heat3D implicit (4D): N=", Nx, ", T=", t, ", dt=", dt, ", max|u-uex|=", err, "\n"
если err < 1e-2 то
  вывод "OK\n"
иначе
  вывод "FAIL\n"
все

Полная программа

алг
нач
  цел Nx, Ny, Nz
  Nx := 16
  Ny := 16
  Nz := 16

  вещ pi
  pi := 3.141592653589793
  вещ h, h2
  h := 1.0 / (Nx - 1)
  h2 := h * h

  | Неявная схема (Backward Euler) для u_t = Δu в 3D: (I - dt Δ) u^{n+1} = u^n
  | Итерационно решаем линейную систему методом Гаусса–Зейделя на каждом шаге по времени
  вещ lam, dt, T
  lam := 0.5 | lam = dt / h^2 (умеренно крупный шаг по времени, неявная схема устойчива)
  dt := lam * h2
  T := 0.10

  цел nt
  nt := int(T / dt)
  если nt < 1 то nt := 1 все

  | 4D массив: u[время, x, y, z]
  вещ таб u[0:nt, 0:Nx-1, 0:Ny-1, 0:Nz-1]

  цел i, j, k, n, it, maxIt
  вещ x, y, z, t
  вещ decay, uex, diff, err, rhs, newv, oldv, change, denom

  | Начальное условие u(0,x,y,z) = sin(pi x) sin(pi y) sin(pi z)
  нц для i от 0 до Nx-1 шаг 1
    x := i * h
    нц для j от 0 до Ny-1 шаг 1
      y := j * h
      нц для k от 0 до Nz-1 шаг 1
        z := k * h
        u[0,i,j,k] := sin(pi * x) * sin(pi * y) * sin(pi * z)
      кц
    кц
  кц

  | Граница нулевая для всех t
  maxIt := 200 | максимум внутренних итераций Гаусса–Зейделя на шаг
  denom := 1.0 + 6.0 * lam

  нц для n от 0 до nt-1 шаг 1
    | Инициализируем начальное приближение u[n+1] := u[n]
    нц для i от 0 до Nx-1 шаг 1
      нц для j от 0 до Ny-1 шаг 1
        нц для k от 0 до Nz-1 шаг 1
          u[n+1,i,j,k] := u[n,i,j,k]
        кц
      кц
    кц

    | Внутренние итерации решения (I - dt Δ) u^{n+1} = u^n методом Гаусса–Зейделя
    нц для it от 0 до maxIt-1 шаг 1
      change := 0.0
      нц для i от 1 до Nx-2 шаг 1
        нц для j от 1 до Ny-2 шаг 1
          нц для k от 1 до Nz-2 шаг 1
            oldv := u[n+1,i,j,k]
            rhs := u[n,i,j,k]
            newv := (rhs + lam * (u[n+1,i+1,j,k] + u[n+1,i-1,j,k] + u[n+1,i,j+1,k] + u[n+1,i,j-1,k] + u[n+1,i,j,k+1] + u[n+1,i,j,k-1])) / denom
            u[n+1,i,j,k] := newv
            diff := abs(newv - oldv)
            change := max(change, diff)
          кц
        кц
      кц
      | Ранний выход по сходимости внутреннего решателя
      если change < 1e-6 то выход все
    кц
  кц

  | Оценка ошибки на времени t = nt * dt относительно точного решения e^{-3 pi^2 t} sin(pi x) sin(pi y) sin(pi z)
  t := nt * dt
  decay := exp(-3.0 * pi * pi * t)
  err := 0.0
  нц для i от 0 до Nx-1 шаг 1
    x := i * h
    нц для j от 0 до Ny-1 шаг 1
      y := j * h
      нц для k от 0 до Nz-1 шаг 1
        z := k * h
        uex := decay * sin(pi * x) * sin(pi * y) * sin(pi * z)
        diff := abs(u[nt,i,j,k] - uex)
        err := max(err, diff)
      кц
    кц
  кц

  вывод "Heat3D implicit (4D): N=", Nx, ", T=", t, ", dt=", dt, ", max|u-uex|=", err, "\n"
  если err < 1e-2 то
    вывод "OK\n"
  иначе
    вывод "FAIL\n"
  все
кон

▶ Запустить пример