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

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

← Все примеры

Уравнение Пуассона 3D (метод Якоби)

Итерационное решение трёхмерного уравнения Пуассона методом Якоби. Уравнение Пуассона -- одно из фундаментальных уравнений математической физики: оно описывает стационарное (не меняющееся со временем) распределение потенциала -- электрического, гравитационного или температурного. Программа находит численное решение на сетке $16 \times 16 \times 16$ и сравнивает его с точным.

Задача

Уравнение Пуассона в кубе $[0,1]^3$:

$-\Delta u = -\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right) = f(x,y,z)$

с граничными условиями Дирихле $u = 0$ на всех гранях куба.

В нашем примере правая часть подобрана так, чтобы точное решение было известно:

$u^*(x,y,z) = \sin(\pi x) \sin(\pi y) \sin(\pi z), \quad f(x,y,z) = 3\pi^2 \sin(\pi x) \sin(\pi y) \sin(\pi z)$

(Подставив $u^*$ в $-\Delta u$, получим ровно $f$.)

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

  1. Разбиваем куб на равномерную сетку с шагом $h = 1/(N-1)$ и записываем дискретный аналог уравнения с помощью семиточечного шаблона (7-point stencil): центральная точка и шесть соседей.
  2. Получаем формулу: $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)$.
  3. Начинаем с нулевого приближения $u = 0$.
  4. На каждой итерации вычисляем новые значения un по формуле, читая из u (метод Якоби использует значения только с предыдущей итерации).
  5. Копируем un обратно в u.
  6. Проверяем точность -- если максимальное отклонение от точного решения меньше допуска $10^{-3}$, завершаем досрочно.

Разбор

Инициализация сетки и правой части

Создаём три массива: u (текущее приближение), un (новое приближение) и f (правая часть уравнения). Правую часть вычисляем аналитически: $f = 3\pi^2 \sin(\pi x) \sin(\pi y) \sin(\pi z)$.

вещ таб u[0:Nx-1, 0:Ny-1, 0:Nz-1]
вещ таб un[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
      un[i,j,k] := 0.0
      f[i,j,k] := 3.0 * pi * pi * sin(pi * x) * sin(pi * y) * sin(pi * z)
    кц
  кц
кц

Итерации Якоби

Главная особенность метода Якоби: новые значения un вычисляются целиком из старых u. Поэтому нужны два массива. После прохода по всем точкам un копируется в u.

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

| Итерации Якоби для -Δu = f (7-точечный трафарет)
нц для it от 0 до iters-1 шаг 1
  нц для i от 1 до Nx-2 шаг 1
    нц для j от 1 до Ny-2 шаг 1
      нц для k от 1 до Nz-2 шаг 1
        | -Δu = f  ⇒  6 u[i,j,k] - (соседи) = h^2 f[i,j,k]
        | u[i,j,k] = (соседи + h^2 f[i,j,k]) / 6
        un[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
      кц
    кц
  кц
  | Копирование обратно (внутренняя область)
  нц для i от 1 до Nx-2 шаг 1
    нц для j от 1 до Ny-2 шаг 1
      нц для k от 1 до Nz-2 шаг 1
        u[i,j,k] := un[i,j,k]
      кц
    кц
  кц

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

После каждой итерации вычисляем максимальное отклонение от точного решения. Если оно меньше $10^{-3}$ -- выходим досрочно, не дожидаясь всех 2000 итераций.

  | Оценка ошибки и возможный ранний выход
  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: N=", Nx, " max|u-uex|=", err, "\n"
| Простая проверка точности: для N=16 и Якоби с ранним выходом добиваемся < 1e-3
если 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]
  вещ таб un[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
        un[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 := 2000
  вещ tol
  tol := 1e-3

  | Итерации Якоби для -Δu = f (7-точечный трафарет)
  нц для it от 0 до iters-1 шаг 1
    нц для i от 1 до Nx-2 шаг 1
      нц для j от 1 до Ny-2 шаг 1
        нц для k от 1 до Nz-2 шаг 1
          | -Δu = f  ⇒  6 u[i,j,k] - (соседи) = h^2 f[i,j,k]
          | u[i,j,k] = (соседи + h^2 f[i,j,k]) / 6
          un[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
        кц
      кц
    кц
    | Копирование обратно (внутренняя область)
    нц для i от 1 до Nx-2 шаг 1
      нц для j от 1 до Ny-2 шаг 1
        нц для k от 1 до Nz-2 шаг 1
          u[i,j,k] := un[i,j,k]
        кц
      кц
    кц
    | Оценка ошибки и возможный ранний выход
    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: N=", Nx, " max|u-uex|=", err, "\n"
  | Простая проверка точности: для N=16 и Якоби с ранним выходом добиваемся < 1e-3
  если err < tol то
    вывод "OK\n"
  иначе
    вывод "FAIL\n"
  все
кон

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