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

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

← Все примеры

Уравнение Пуассона 3D (Гаусс--Зейдель)

Решение трёхмерного уравнения Пуассона методом Гаусса--Зейделя. Это улучшенная версия метода Якоби: обновление происходит «на месте» (in-place), то есть новые значения сразу используются для вычисления соседних точек. Благодаря этому метод сходится примерно вдвое быстрее и требует только один массив вместо двух.

Задача

То же уравнение Пуассона, что и в примере с методом Якоби:

$-\Delta u = f(x,y,z) \quad \text{в } [0,1]^3, \quad u = 0 \text{ на границе}$

Правая часть: $f(x,y,z) = 3\pi^2 \sin(\pi x) \sin(\pi y) \sin(\pi z)$.

Точное решение: $u^*(x,y,z) = \sin(\pi x) \sin(\pi y) \sin(\pi z)$.

Дискретизация по семиточечному шаблону (7-point stencil) даёт формулу:

$u_{i,j,k} = \frac{1}{6}\left(u_{i+1,j,k} + u_{i-1,j,k} + u_{i,j+1,k} + u_{i,j-1,k} + u_{i,j,k+1} + u_{i,j,k-1} + h^2 f_{i,j,k}\right)$

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

  1. Инициализируем $u = 0$ во всех точках.
  2. На каждой итерации обходим все внутренние точки сетки и обновляем $u_{i,j,k}$ на месте (записываем сразу в тот же массив).
  3. Порядок обхода важен: при вычислении $u_{i,j,k}$ соседи с меньшими индексами ($u_{i-1,j,k}$, $u_{i,j-1,k}$, $u_{i,j,k-1}$) уже обновлены на текущей итерации, а соседи с большими ($u_{i+1,j,k}$, $u_{i,j+1,k}$, $u_{i,j,k+1}$) -- ещё нет. Это «смешивание» старых и новых значений и ускоряет сходимость.
  4. После каждой итерации проверяем точность. Если максимальное отклонение от точного решения меньше $10^{-3}$ -- выходим досрочно.

Разбор

Инициализация

В отличие от метода Якоби, здесь нужен только один массив u -- второй не требуется, потому что обновление происходит на месте.

вещ таб u[0:Nx-1, 0:Ny-1, 0:Nz-1]
вещ таб f[0:Nx-1, 0:Ny-1, 0:Nz-1]

| Инициализация: точное решение u* = sin(pi x) sin(pi y) sin(pi z)
| Правая часть f = 3*pi^2 * u*
нц для 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[i,j,k] := 0.0
      f[i,j,k] := 3.0 * pi * pi * sin(pi * x) * sin(pi * y) * sin(pi * z)
    кц
  кц
кц

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

Формула та же, что и в методе Якоби, но запись идёт прямо в u (а не во вспомогательный un). Именно это отличие -- использование уже обновлённых соседей -- даёт ускорение сходимости.

цел it, iters
iters := 1000
вещ tol
tol := 1e-3

| Gauss-Seidel (in-place) для -Δu = f: u = (соседи + h^2 f)/6
нц для it от 0 до iters-1 шаг 1
  нц для i от 1 до Nx-2 шаг 1
    нц для j от 1 до Ny-2 шаг 1
      нц для k от 1 до Nz-2 шаг 1
        u[i,j,k] := (u[i+1,j,k] + u[i-1,j,k] + u[i,j+1,k] + u[i,j-1,k] + u[i,j,k+1] + u[i,j,k-1] + h2 * f[i,j,k]) / 6.0
      кц
    кц
  кц

Проверка сходимости и ранний выход

  | Оценка ошибки и возможный ранний выход
  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 := sin(pi * x) * sin(pi * y) * sin(pi * z)
        diff := abs(u[i,j,k] - uex)
        err := max(err, diff)
      кц
    кц
  кц
  если err < tol то
    выход
  все
кц

Финальный вывод

вывод "Poisson3D-GS: N=", Nx, " max|u-uex|=", err, "\n"
если err < tol то
  вывод "OK\n"
иначе
  вывод "FAIL\n"
все

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

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

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

  вещ таб u[0:Nx-1, 0:Ny-1, 0:Nz-1]
  вещ таб f[0:Nx-1, 0:Ny-1, 0:Nz-1]

  цел i, j, k
  вещ x, y, z
  вещ err, uex, diff

  | Инициализация: точное решение u* = sin(pi x) sin(pi y) sin(pi z)
  | Правая часть f = 3*pi^2 * u*
  нц для 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[i,j,k] := 0.0
        f[i,j,k] := 3.0 * pi * pi * sin(pi * x) * sin(pi * y) * sin(pi * z)
      кц
    кц
  кц

  | Граничные условия Дирихле: u = 0 на границе уже инициализировано нулями

  цел it, iters
  iters := 1000
  вещ tol
  tol := 1e-3

  | Gauss-Seidel (in-place) для -Δu = f: u = (соседи + h^2 f)/6
  нц для it от 0 до iters-1 шаг 1
    нц для i от 1 до Nx-2 шаг 1
      нц для j от 1 до Ny-2 шаг 1
        нц для k от 1 до Nz-2 шаг 1
          u[i,j,k] := (u[i+1,j,k] + u[i-1,j,k] + u[i,j+1,k] + u[i,j-1,k] + u[i,j,k+1] + u[i,j,k-1] + h2 * f[i,j,k]) / 6.0
        кц
      кц
    кц

    | Оценка ошибки и возможный ранний выход
    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 := sin(pi * x) * sin(pi * y) * sin(pi * z)
          diff := abs(u[i,j,k] - uex)
          err := max(err, diff)
        кц
      кц
    кц
    если err < tol то
      выход
    все
  кц

  | Финальная проверка
  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 := sin(pi * x) * sin(pi * y) * sin(pi * z)
        diff := abs(u[i,j,k] - uex)
        err := max(err, diff)
      кц
    кц
  кц

  вывод "Poisson3D-GS: N=", Nx, " max|u-uex|=", err, "\n"
  если err < tol то
    вывод "OK\n"
  иначе
    вывод "FAIL\n"
  все
кон

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