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

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

← Все примеры

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

Численное решение трёхмерного уравнения теплопроводности методом конечных разностей. Программа моделирует, как тепло распространяется внутри куба, используя явную схему Эйлера (Forward Euler) -- простейший способ шагать по времени. Результат сравнивается с точным аналитическим решением.

Задача

Уравнение теплопроводности описывает, как температура $u(x,y,z,t)$ меняется со временем в трёхмерной области $[0,1]^3$:

$\frac{\partial u}{\partial t} = \Delta u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}$

Начальное условие (распределение температуры в момент $t = 0$):

$u_0(x,y,z) = \sin(\pi x) \sin(\pi y) \sin(\pi z)$

Граничные условия (температура на стенках куба всегда равна нулю):

$u = 0 \quad \text{на границе } [0,1]^3$

Точное решение (для проверки):

$u(x,y,z,t) = e^{-3\pi^2 t} \sin(\pi x) \sin(\pi y) \sin(\pi z)$

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

  1. Разбиваем куб $[0,1]^3$ на равномерную сетку $N \times N \times N$ с шагом $h = 1/(N-1)$.
  2. Выбираем шаг по времени $dt = \alpha \cdot h^2$, где $\alpha \leq 1/6$ -- условие устойчивости явной схемы в 3D.
  3. На каждом шаге по времени обновляем температуру во всех внутренних точках сетки по семиточечному шаблону (7-point stencil): берём значение в текущей точке и шести её соседях (слева, справа, сверху, снизу, спереди, сзади).
  4. Граничные точки не трогаем -- они всегда равны нулю.
  5. После всех временных шагов сравниваем результат с точным решением.

Разбор

Параметры сетки и шага по времени

Задаём размер сетки $16 \times 16 \times 16$ и вычисляем шаг $h$. Параметр $\alpha = dt/h^2 = 0.1$ выбран с запасом меньше предела устойчивости $1/6 \approx 0.167$.

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

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

| Явная схема для уравнения теплопроводности u_t = Δu в 3D
| Выбор шага по времени для устойчивости: dt ≤ h^2 / (2*d), d=3 ⇒ dt ≤ h^2/6
вещ dt, T, alpha
alpha := 0.1      | alpha = dt / h^2 ≤ 1/6, берём с запасом
dt := alpha * h2
T := 0.05         | конечное время моделирования

Начальное условие

Заполняем массив u начальным распределением температуры $u_0 = \sin(\pi x) \sin(\pi y) \sin(\pi z)$. Массив un -- буфер для нового значения на следующем шаге.

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

нц для 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] := sin(pi * x) * sin(pi * y) * sin(pi * z)
      un[i,j,k] := 0.0
    кц
  кц
кц

Основной цикл по времени

Явная схема Эйлера: новое значение в каждой точке вычисляется по семиточечному шаблону. Формула берёт текущее значение $u_{i,j,k}$ и шесть соседей -- по два вдоль каждой из трёх осей. Это дискретный аналог оператора Лапласа $\Delta u$:

$u_{i,j,k}^{n+1} = u_{i,j,k}^{n} + \alpha \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} - 6 \cdot u_{i,j,k} \right)$

Обновляем только внутренние точки (индексы от 1 до $N-2$), затем копируем результат обратно.

нц для n от 0 до nt-1 шаг 1
  | Обновляем только внутреннюю область, границы остаются нулями
  нц для i от 1 до Nx-2 шаг 1
    нц для j от 1 до Ny-2 шаг 1
      нц для k от 1 до Nz-2 шаг 1
        un[i,j,k] := u[i,j,k] + alpha * ((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]) - 6.0 * u[i,j,k])
      кц
    кц
  кц
  | Копируем внутреннюю область обратно
  нц для i от 1 до Nx-2 шаг 1
    нц для j от 1 до Ny-2 шаг 1
      нц для k от 1 до Nz-2 шаг 1
        u[i,j,k] := un[i,j,k]
      кц
    кц
  кц
кц

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

Вычисляем точное решение $u(x,y,z,t) = e^{-3\pi^2 t} \sin(\pi x) \sin(\pi y) \sin(\pi z)$ и находим максимальное отклонение. Если ошибка меньше $10^{-2}$, тест считается пройденным.

вещ decay
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[i,j,k] - uex)
      err := max(err, diff)
    кц
  кц
кц

вывод "Heat3D (explicit): 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
  h := 1.0 / (Nx - 1)
  вещ h2
  h2 := h * h

  | Явная схема для уравнения теплопроводности u_t = Δu в 3D
  | Выбор шага по времени для устойчивости: dt ≤ h^2 / (2*d), d=3 ⇒ dt ≤ h^2/6
  вещ dt, T, alpha
  alpha := 0.1      | alpha = dt / h^2 ≤ 1/6, берём с запасом
  dt := alpha * h2
  T := 0.05         | конечное время моделирования

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

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

  | Начальное условие и точное решение для проверки:
  | u0(x,y,z) = sin(πx) sin(πy) sin(πz)
  | u(x,y,z,t) = e^{-3π^2 t} u0(x,y,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[i,j,k] := sin(pi * x) * sin(pi * y) * sin(pi * z)
        un[i,j,k] := 0.0
      кц
    кц
  кц

  | Нулевые граничные условия (согласованы с u0 на границе)

  | Число шагов по времени и фактическое конечное время (T может скорректироваться)
  nt := int(T / dt)
  если nt < 1 то
    nt := 1
  все
  t := nt * dt

  | Основной цикл по времени: явная схема Эйлера вперёд
  нц для n от 0 до nt-1 шаг 1
    | Обновляем только внутреннюю область, границы остаются нулями
    нц для i от 1 до Nx-2 шаг 1
      нц для j от 1 до Ny-2 шаг 1
        нц для k от 1 до Nz-2 шаг 1
          un[i,j,k] := u[i,j,k] + alpha * ((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]) - 6.0 * u[i,j,k])
        кц
      кц
    кц
    | Копируем внутреннюю область обратно
    нц для i от 1 до Nx-2 шаг 1
      нц для j от 1 до Ny-2 шаг 1
        нц для k от 1 до Nz-2 шаг 1
          u[i,j,k] := un[i,j,k]
        кц
      кц
    кц
  кц

  | Оценка ошибки по максимум-норме относительно точного решения u(x,y,z,t)=e^{-3π^2 t} u0
  вещ decay
  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[i,j,k] - uex)
        err := max(err, diff)
      кц
    кц
  кц

  вывод "Heat3D (explicit): N=", Nx, ", T=", t, ", dt=", dt, ", max|u-uex|=", err, "\n"
  | Базовая проверка: при Nx=16, T≈0.05 и alpha=0.1 ошибка должна быть небольшой (< 1e-2)
  если err < 1e-2 то
    вывод "OK\n"
  иначе
    вывод "FAIL\n"
  все
кон

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