От уравнений Максвелла к анимации: реализация FDTD на PascalABC.NET / Хабр

В этой статье мы будем моделировать распространение электромагнитных волн. Для этого будем численно решать систему дифференциальных уравнений в частных производных Максвелла.

При этом подойдём к этой задаче более профессионально, чем в предыдущих материалах: я учёл справедливую критику и замечания читателей по поводу выбора и обоснования численного метода, оценил корректность и реалистичность схемы и результатов моделирования, привёл необходимые условия для воспроизводимости.

В системе СИ уравнения Максвелла записываются следующим образом:

$\begin{align} \nabla \times \mathbf{E} &= -\frac{\partial \mathbf{B}}{\partial t}, \\ \nabla \times \mathbf{H} &= \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}, \\ \nabla \cdot \mathbf{D} &= \rho, \\ \nabla \cdot \mathbf{B} &= 0. \end{align}$

Первое уравнение‑ это Закон индукции Фарадея, читается как: Изменение магнитной индукции порождает вихревое электрическое поле.
Второе уравнение‑ это Теорема о циркуляции магнитного поля, читается как: Электрический ток и изменение электрической индукции порождают вихревое магнитное поле.
Третье уравнение системы‑это так называемый Закон Гаусса— Электрический заряд является источником электрической индукции.
Четвёртое уравнение‑ Закон Гаусса для магнитного поля (Магнитные заряды не обнаружены).
Введённые обозначения:
ρ — объёмная плотность стороннего электрического заряда (в единицах СИ — Кл/м³);
j — плотность электрического тока (плотность тока проводимости) (в единицах СИ — А/м²);
c — скорость света в вакууме (299 792 458 м/с);
E — напряжённость электрического поля (в единицах СИ — В/м);
H— напряжённость магнитного поля (в единицах СИ — А/м);
D— электрическая индукция (в единицах СИ — Кл/м²);
B— магнитная индукция (в единицах СИ — Тл = Вб/м² = кг•с⁻²•А⁻¹);
∇— дифференциальный оператор набла
Источник:[1]

Систему уравнений необходимо дополнять материальными уравнениями. Для самого простого случая (изотропная и однородная среда) они выглядят так:

$\mathbf{D} = \varepsilon \mathbf{E}, \quad \mathbf{B} = \mu \mathbf{H}.$

Где коэффициенты

$\varepsilon = \varepsilon_0 \varepsilon_r, \quad \mu = \mu_0 \mu_r.$

(это диэлектрическая и магнитная проницаемости)‑ постоянны.

Частным случаем этих уравнений является электростатика(когда электрическое поле не меняется со временем)

а также магнитостатика( когда магнитное поле не меняется со временем), например в постоянных магнитах:

магнитное поле постоянного цилиндрического магнита. Неточная модель.

Однако мы будем моделировать электромагнитные явления в динамике (сигнал Гаусса).

Будем использовать FDTD‑метод‑ Метод конечных разностей во временно́й области.

Для моделирования необходимо использовать конечно‑разностные дискретные формулы. Они являются классическими в научной литературе (https://elibrary.ru/download/elibrary_20799764_92247445.pdf и https://elibrary.ru/vvfyhz):

Для компоненты Ez электрического поля:

$E_z^{n+1}(i,j,k) = E_z^n(i,j,k) + \frac{\Delta t}{\varepsilon} \left[ \frac{H_y^{n+1/2}(i+1,j,k) - H_y^{n+1/2}(i,j,k)}{\Delta x} - \frac{H_x^{n+1/2}(i,j+1,k) - H_x^{n+1/2}(i,j,k)}{\Delta y} \right].$

Для компоненты Hx магнитного поля:

$H_x^{n+1/2}(i,j,k) = H_x^{n-1/2}(i,j,k) - \frac{\Delta t}{\mu} \left[ \frac{E_z^n(i,j,k+1) - E_z^n(i,j,k)}{\Delta z} - \frac{E_y^n(i,j,k+1) - E_y^n(i,j,k)}{\Delta y} \right].$

Источником электромагнитных возмущений будет классический Гауссов импульс:

$S(t) = \exp\!\left(-\frac{(t - t_0)^2}{\sigma_t^2}\right).$

или в дискретной форме:

$S_n = \exp\!\left(-\frac{(n\,\Delta t - t_0)^2}{\sigma_t^2}\right).$

Оговорим, что источник в коде введён как плотность тока, эквивалентный гауссову сигналу.

Параметры Гауссового импульса: J0=10^6 А/м², t0=30 шагов и σt=10 шагов, действующий в одной ячейке в центре области (Возбуждение подаётся в ячейку Ex с индексами (50, 49, 49), а не ровно в центр — это следствие сдвига компонентов на сетке Йи.).

Условие устойчивости‑ критерий Куранта:

$\Delta t \le \frac{1}{c \sqrt{\left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}\right)}},$

Шаг по времени выбирался из критерия Куранта с запасом 5%, чему соответствует число Куранта ≈0,55.

Дополнительно, для контроля корректности работы схемы, на каждом N‑м шаге вычислялась интегральная энергия электромагнитного поля:

$\begin{equation} W(t) = \sum_{i,j,k} \left( |\mathbf{E}_{i,j,k}(t)|^2 + |\mathbf{H}_{i,j,k}(t)|^2 \right), \end{equation}$

Отсутствие экспоненциального роста W(t) подтверждает, что численная ошибка не накапливается и схема остаётся устойчивой на всём интервале моделирования.

Также подсчитывается max |div E| с помощью функции MaxDivergenceE. Вне области источника дискретная дивергенция E должна быть близка к машинному нулю; вблизи источника она ненулевая.

Граничные условия— 3 грани PEC, 3 грани PMC: E_tang = 0 на границах

то есть тангенциальные компоненты на трёх «нижних» гранях x=0, y=0, z=0.

Верхние грани (x=Nx-1 и так далее) не содержат тангенциальных E‑компонент в массивах, поэтому там реализуется «магнитная» стенка (H=0).

Эти условия не диссипируют энергию при идеальных условиях.

Моделирование проводилось в кубической области 100×100×100 мм с шагом сетки Δx=Δy=Δz=1 мм. Выполнялось Nt=300 шагов, что соответствует ≈550 пс; за это время фронт волны достигает стенок области и отражается от них. Визуализировался срез Ez в плоскости z=L/2 через каждые 5 шагов.

Процедура обезразмеривания в данном случае необязательна — она необходима там, где числа «выпадают» из комфортного диапазона машинной арифметики — либо переполняются, либо тонут под единицей.

Данную вычислительную схему реализуем на языке семейства Pascal в среде PascalABC.NET, (старый‑добрый), поскольку на Python или C++ у меня возникли проблемы с индексацией.

program FDTD3D_Animated_PEC_Energy;

uses
  GraphABC, System;

type
  Real = double;

const
  Nx = 100;
  Ny = 100;
  Nz = 100;

  dx = 1.0e-3;
  dy = 1.0e-3;
  dz = 1.0e-3;
  c  = 2.998e8;
  eps0 = 8.854e-12;
  mu0  = 4.0 * Pi * 1.0e-7;

  steps = 300;
  snap_every = 5;

  t0 = 30.0;
  sigma_t = 10.0;

  J0 = 1.0e6;

  // ---------- Раскладка окна ----------
  WinWidth  = 1240;
  WinHeight = 660;
  AnimX     = 10;
  AnimY     = 40;
  AnimSize  = 600;
  PlotX     = 630;
  PlotY     = 40;
  PlotW     = 590;
  PlotH     = 600;

var
  dt: Real;
  Ex, Ey, Ez, Hx, Hy, Hz: array of Real;
  eps_r, mu_r: array of Real;
  src_i, src_j, src_k, slice_z: integer;
  W, H: integer;
  E_hist: array of Real;   // история полной энергии

// ---------- Индексаторы Yee-сетки ----------
function idxEx(i, j, k: integer): integer;
begin Result := (i * (Ny - 1) + j) * (Nz - 1) + k; end;
function idxEy(i, j, k: integer): integer;
begin Result := (i * Ny + j) * (Nz - 1) + k; end;
function idxEz(i, j, k: integer): integer;
begin Result := (i * (Ny - 1) + j) * Nz + k; end;
function idxHx(i, j, k: integer): integer;
begin Result := (i * Ny + j) * Nz + k; end;
function idxHy(i, j, k: integer): integer;
begin Result := (i * (Ny - 1) + j) * Nz + k; end;
function idxHz(i, j, k: integer): integer;
begin Result := (i * Ny + j) * (Nz - 1) + k; end;
function idxEps(i, j, k: integer): integer;
begin Result := (i * Ny + j) * Nz + k; end;

// ---------- Источник ----------
function source(n: integer): Real;
begin
  Result := Exp(-Sqr((n - t0) / sigma_t));
end;

// ---------- PEC: E_tang = 0 на нижних гранях ----------
procedure ApplyPEC;
var i, j, k: integer;
begin
  for j := 0 to Ny - 1 do
    for k := 0 to Nz - 2 do
      Ey[idxEy(0, j, k)] := 0;
  for j := 0 to Ny - 2 do
    for k := 0 to Nz - 1 do
      Ez[idxEz(0, j, k)] := 0;

  for i := 0 to Nx - 1 do
    for k := 0 to Nz - 2 do
      Ex[idxEx(i, 0, k)] := 0;
  for i := 0 to Nx - 2 do
    for k := 0 to Nz - 1 do
      Ez[idxEz(i, 0, k)] := 0;

  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 2 do
      Ex[idxEx(i, j, 0)] := 0;
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 1 do
      Ey[idxEy(i, j, 0)] := 0;
end;

// ---------- Полная энергия ----------
function TotalEnergy: Real;
var i, j, k: integer; sum: Real;
begin
  sum := 0;
  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 2 do
      for k := 0 to Nz - 2 do
        sum += eps0 * Sqr(Ex[idxEx(i,j,k)]);
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 1 do
      for k := 0 to Nz - 2 do
        sum += eps0 * Sqr(Ey[idxEy(i,j,k)]);
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 2 do
      for k := 0 to Nz - 1 do
        sum += eps0 * Sqr(Ez[idxEz(i,j,k)]);
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 1 do
      for k := 0 to Nz - 1 do
        sum += mu0 * Sqr(Hx[idxHx(i,j,k)]);
  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 2 do
      for k := 0 to Nz - 1 do
        sum += mu0 * Sqr(Hy[idxHy(i,j,k)]);
  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 1 do
      for k := 0 to Nz - 2 do
        sum += mu0 * Sqr(Hz[idxHz(i,j,k)]);
  Result := 0.5 * sum * dx * dy * dz;
end;

// ---------- max |div E| ----------
function MaxDivergenceE: Real;
var i, j, k: integer; d, maxd: Real;
begin
  maxd := 0;
  for i := 1 to Nx - 2 do
    for j := 1 to Ny - 2 do
      for k := 1 to Nz - 2 do
      begin
        d := (Ex[idxEx(i+1,j,k)] - Ex[idxEx(i,j,k)]) / dx
           + (Ey[idxEy(i,j+1,k)] - Ey[idxEy(i,j,k)]) / dy
           + (Ez[idxEz(i,j,k+1)] - Ez[idxEz(i,j,k)]) / dz;
        if Abs(d) > maxd then maxd := Abs(d);
      end;
  Result := maxd;
end;

// ---------- Срез Ez (рисует только свою панель) ----------
procedure DrawSlice(frameArr: array of Real; x0, y0, sz: integer);
var
  i, j, cell: integer;
  val, norm: Real;
  col: Color;
begin
  norm := 0.0;
  for i := 0 to W * H - 1 do
    if Abs(frameArr[i]) > norm then norm := Abs(frameArr[i]);
  if norm = 0 then norm := 1.0;

  cell := sz div W;
  if cell < 1 then cell := 1;

  for i := 0 to W - 1 do
    for j := 0 to H - 1 do
    begin
      val := frameArr[i * H + j] / norm;
      if val < 0 then
        col := RGB(0, 0, 255 + Round(255 * val))
      else
        col := RGB(255, 255 - Round(255 * val), 0);
      Brush.Color := col;
      FillRect(x0 + i * cell, y0 + j * cell,
               x0 + (i + 1) * cell, y0 + (j + 1) * cell);
    end;

  Font.Color := clBlack;
  Font.Size := 11;
  TextOut(x0 + 5, y0 - 22, 'Срез Ez, z = ' + IntToStr(slice_z));
end;

// ---------- График энергии W(t) ----------
procedure DrawEnergyPlot(E_hist: array of Real; nCur, x0, y0, w, h: integer);
var
  i, px, py: integer;
  Emax, frac: Real;
  marginL, marginB, marginT, marginR, plotW, plotH: integer;
begin
  marginL := 60;
  marginB := 40;
  marginT := 25;
  marginR := 20;
  plotW := w - marginL - marginR;
  plotH := h - marginT - marginB;

  Emax := 1e-30;
  for i := 0 to nCur do
    if E_hist[i] > Emax then Emax := E_hist[i];
  if Emax <= 0 then Emax := 1;

  // Оси
  Pen.Width := 1;
  Pen.Color := clBlack;
  Line(x0 + marginL, y0 + marginT, x0 + marginL, y0 + marginT + plotH);
  Line(x0 + marginL, y0 + marginT + plotH,
       x0 + marginL + plotW, y0 + marginT + plotH);

  Font.Color := clBlack;
  Font.Size := 9;

  // Y-метки: 0%, 25%, 50%, 75%, 100%
  for i := 0 to 4 do
  begin
    frac := i / 4;
    py := y0 + marginT + plotH - Round(plotH * frac);
    Line(x0 + marginL - 5, py, x0 + marginL, py);
    TextOut(x0 + 10, py - 7, IntToStr(Round(frac * 100)) + '%');
  end;

  // X-метки
  for i := 0 to 4 do
  begin
    px := x0 + marginL + Round(plotW * i / 4);
    Line(px, y0 + marginT + plotH, px, y0 + marginT + plotH + 5);
    TextOut(px - 12, y0 + marginT + plotH + 8,
            IntToStr(Round((steps - 1) * i / 4)));
  end;

  // Кривая
  Pen.Width := 2;
  Pen.Color := clRed;
  for i := 0 to nCur do
  begin
    px := x0 + marginL + Round(plotW * i / (steps - 1));
    py := y0 + marginT + plotH - Round(plotH * E_hist[i] / Emax);
    if i = 0 then MoveTo(px, py) else LineTo(px, py);
  end;

  // «Бегунок» текущего шага
  if nCur > 0 then
  begin
    Pen.Width := 1;
    Pen.Color := clBlue;
    px := x0 + marginL + Round(plotW * nCur / (steps - 1));
    Line(px, y0 + marginT, px, y0 + marginT + plotH);
  end;

  // Заголовок и абсолютный максимум
  Font.Color := clBlack;
  Font.Size := 11;
  TextOut(x0 + marginL, y0 + 2, 'Полная энергия W(t)');
  Font.Size := 9;
  TextOut(x0 + marginL + 5, y0 + marginT + 5, Format('W_max = {0:E2}', Emax));
  TextOut(x0 + marginL + plotW div 2 - 25, y0 + h - 15, 'номер шага');
end;

// ======================== ОСНОВНАЯ ПРОГРАММА ========================
begin
  SetWindowSize(WinWidth, WinHeight);
  Window.Title := 'FDTD 3D PEC: срез Ez и энергия';

  dt := 0.95 * dx / (Sqrt(3.0) * c);
  WriteLn('CFL: dt = ', dt, '  c*dt/dx = ', c*dt/dx,
          '  (предел 1/sqrt(3) = 0.5774)');

  src_i := Nx div 2;
  src_j := Ny div 2;
  src_k := Nz div 2;
  slice_z := Nz div 2;

  SetLength(Ex, Nx * (Ny - 1) * (Nz - 1));
  SetLength(Ey, (Nx - 1) * Ny * (Nz - 1));
  SetLength(Ez, (Nx - 1) * (Ny - 1) * Nz);
  SetLength(Hx, (Nx - 1) * Ny * Nz);
  SetLength(Hy, Nx * (Ny - 1) * Nz);
  SetLength(Hz, Nx * Ny * (Nz - 1));
  SetLength(eps_r, Nx * Ny * Nz);
  SetLength(mu_r,  Nx * Ny * Nz);
  SetLength(E_hist, steps);

  for var i := 0 to High(eps_r) do
  begin
    eps_r[i] := 1.0;
    mu_r[i]  := 1.0;
  end;

  W := Nx - 1;
  H := Ny - 1;
  var frameArr := new Real[W * H];

  // ================= Временной цикл =================
  for var n := 0 to steps - 1 do
  begin
    // ---------- Обновление H ----------
    for var i := 1 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEz_dy := (Ez[idxEz(i, j, k)] - Ez[idxEz(i, j - 1, k)]) / dy;
          var dEy_dz := (Ey[idxEy(i, j, k)] - Ey[idxEy(i, j, k - 1)]) / dz;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hx[idxHx(i - 1, j, k)] += -(dEz_dy - dEy_dz) / mu * dt;
        end;

    for var i := 1 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEx_dz := (Ex[idxEx(i, j, k)] - Ex[idxEx(i, j, k - 1)]) / dz;
          var dEz_dx := (Ez[idxEz(i, j, k)] - Ez[idxEz(i - 1, j, k)]) / dx;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hy[idxHy(i, j - 1, k)] += -(dEx_dz - dEz_dx) / mu * dt;
        end;

    for var i := 1 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEy_dx := (Ey[idxEy(i, j, k)] - Ey[idxEy(i - 1, j, k)]) / dx;
          var dEx_dy := (Ex[idxEx(i, j, k)] - Ex[idxEx(i, j - 1, k)]) / dy;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hz[idxHz(i, j, k - 1)] += -(dEy_dx - dEx_dy) / mu * dt;
        end;

    // ---------- Обновление E ----------
    for var i := 0 to Nx - 1 do
      for var j := 0 to Ny - 2 do
        for var k := 0 to Nz - 2 do
        begin
          var dHz_dy := (Hz[idxHz(i, j + 1, k)] - Hz[idxHz(i, j, k)]) / dy;
          var dHy_dz := (Hy[idxHy(i, j, k + 1)] - Hy[idxHy(i, j, k)]) / dz;
          var eps := eps0 * eps_r[idxEps(i, j, k)];
          Ex[idxEx(i, j, k)] += (dHz_dy - dHy_dz) / eps * dt;
          if (i = src_i) and (j = src_j - 1) and (k = src_k - 1) then
            Ex[idxEx(i, j, k)] += dt * J0 * source(n) / eps;
        end;

    for var i := 0 to Nx - 2 do
      for var j := 0 to Ny - 1 do
        for var k := 0 to Nz - 2 do
        begin
          var dHx_dz := (Hx[idxHx(i, j, k + 1)] - Hx[idxHx(i, j, k)]) / dz;
          var dHz_dx := (Hz[idxHz(i + 1, j, k)] - Hz[idxHz(i, j, k)]) / dx;
          var eps := eps0 * eps_r[idxEps(i, j, k)];
          Ey[idxEy(i, j, k)] += (dHx_dz - dHz_dx) / eps * dt;
        end;

    for var i := 0 to Nx - 2 do
      for var j := 0 to Ny - 2 do
        for var k := 0 to Nz - 1 do
        begin
          var dHy_dx := (Hy[idxHy(i + 1, j, k)] - Hy[idxHy(i, j, k)]) / dx;
          var dHx_dy := (Hx[idxHx(i, j + 1, k)] - Hx[idxHx(i, j, k)]) / dy;
          var eps := eps0 * eps_r[idxEps(i, j, k)];
          Ez[idxEz(i, j, k)] += (dHy_dx - dHx_dy) / eps * dt;
        end;

    ApplyPEC;

    // ---------- Энергия на каждом шаге ----------
    E_hist[n] := TotalEnergy;

    // ---------- Диагностика и отрисовка ----------
    if (n mod snap_every = 0) or (n = steps - 1) then
    begin
      WriteLn('n=', n:4,
              '  W=', E_hist[n]:0:6,
              '  max|divE|=', MaxDivergenceE:0:6);

      for var i := 0 to W - 1 do
        for var j := 0 to H - 1 do
          frameArr[i * H + j] := Ez[idxEz(i, j, slice_z)];

      LockDrawing;
      ClearWindow(clWhite);
      DrawSlice(frameArr, AnimX, AnimY, AnimSize);
      DrawEnergyPlot(E_hist, n, PlotX, PlotY, PlotW, PlotH);
      Redraw;
      Sleep(10);
    end;
  end;

  WriteLn('Расчёт завершён.');
end.

В данном коде решение реализовано на сетке Йи — это специальная разнесённая (смещённая) сетка для численного решения уравнений Максвелла методом FDTD. Её придумал Кейн Йи (Kane Yee) в 1966 году.

Её суть в следующем: В обычной сетке все компоненты поля (Ex, Ey, Ez, Hx, Hy, Hz) ставили бы в одни и те же узлы. В сетке Йи их разносят по пространству на полшага и ещё сдвигают по времени.

Поля считают «в шахматном порядке» по времени:

Сначала считают En на шаге n
Потом считают Hn+1/2 на промежуточном шаге
Потом снова En+1 и так далее

Это и есть алгоритм «leapfrog» («перескок»).

Следует оговорить, что «схема Йи немонотонна‑ на границах раздела диэлектриков фактический порядок точности падает до первого.» Это показано в статье 2017 года https://elibrary.ru/item.asp?edn=vvfyhz Домбровской Ж. О. и Боголюбова А. Н.

Схема Йи консервативна в дискретном смысле — это её классическое свойство, отмеченное ещё в пионерской работе Йи и подробно разбираемое в монографиях по FDTD (смотрите, например, Taflove A., Hagness S. C. Computational Electrodynamics: the Finite‑Difference Time‑Domain Method. Artech House, 2005. ) Косвенно это устанавливается с помощью контроля энергии на каждом шаге, что я захотел продемонстрировать в своей работе.

Результат работы скрипта‑ довольно красивая анимация распространения ЭМ‑волн (срез Ez):

Срез Ez. Гауссов сигнал. Код до отладки. Не уверен в верности-картинка должна быть симметриной. — Срез Ez. Гауссов сигнал. Код до отладки. Не уверен в верности‑картинка должна быть симметриной.

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

Сравним с описаниями из литературы данную анимацию (анализируем артефакты):

Наблюдение	Что должно быть по литературе	Обоснование с точки зрения физики	Вердикт и источник
Фронт расходится из центра	Гауссов импульс из центра распространяется наружу — классическая картинка (Салливан, гл. 1 и 3)	Главный лепесток гауссова импульса расходится от точечного источника. Корректно	✓ https://lshoshia.science.tsu.ge/fdtd/Dennissullivan‑fixed.pdf
Жёлтое и синее чередуются	У гауссова импульса есть положительная и отрицательная полуволны: передний фронт и хвост противоположного знака	В FDTD‑литературе такие поля визуализируются двухцветной картой: красный — положительные значения Ez, синий — отрицательные.	✓ https://www.researchgate.net/publication/5313575_Amplified_total_internal_reflection_theory_analysis_and_demonstration_of_existence_via_FDTD
Овал растягивается и расплывается к углу	Численная дисперсия: высокочастотные компоненты отстают, фронт из сферы становится анизотропным	Фазовая скорость в Yee‑сетке анизотропна; при S ≈ 0,55 фронт приобретает ромбовидную форму, максимум скорости — вдоль телесной диагонали.	~, https://www.academia.edu/112873112/Reduction_of_numerical_dispersion_in_FDTD_method_through_artificial_anisotropy
Растяжение и уход «за угол»	Передний край пакета движется быстрее номинала вблизи Найквиста; вдоль диагонали пакет приходит раньше всего	Скорость вдоль диагонали максимальна (S ≈ 0,55), пакет достигает угла раньше всего; далее наблюдается мода, бегущая вдоль стенки — авторская интерпретация, требует отдельной проверки	~ https://www.researchgate.net/publication/3307223_FDTD_dispersion_revisited_faster‑than‑light_propagation
Красный овал возле центра — ближняя зона	Дискретный диполь на сетке Йи анизотропен, ближнее поле локализовано в нескольких ячейках	Источник — диполь по x, срез Ez проходит в 1,5 ячейки от него; дальняя Ez на таком срезе подавлена (∝ n·r̂), доминирует ближнее поле. Аналогия: анизотропия точечного диполя на Йи показана в статье 2025 г. (там — резонансные диполи, эффект Пёрселла)	✓https://arxiv.org/html/2504.07444
Нижняя половина жёлтая всё время	До прихода отражений картина обязана быть зеркально симметричной относительно источника	Расхождение с физикой и литературой. Подтверждено репликацией схемы: нулевой треугольник ограничен антидиагональю i + j = 99 через ячейку источника — структурный дефект индексации, а не свойство метода	✗ —
К концу остаётся несколько красных ячеек	Остаточное поле не должно доминировать в кадре	Подтверждено расчётом: обрезанный гауссиан депонирует паразитный заряд (сумма источника по времени ≈ 17,7 ≠ 0), который держит статическое поле; по‑кадровая нормировка растягивает его на всю палитру. Плюс артефакт выбора среза/компоненты	~ — расчёт по коду

После отладки кода я добавил график энергии во времени:

Отладка кода+ график энергии. Сомневаюсь- должно быть плато энергии, а не экспонента. — Отладка кода+ график энергии. Сомневаюсь‑ должно быть плато энергии, а не экспонента.

Анализ графика энергии, а также расхождение с физикой показывает, что где‑то в коде был баг. Мне удалось его найти и 'пофиксить'‑ я переписал блок, отвечающий за расчёт H и визуализацию Ez.

// ---------- Обновление H (ИСПРАВЛЕНО: запись в собственный слот) ----------
    for var i := 0 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEz_dy := (Ez[idxEz(i, j, k)] - Ez[idxEz(i, j - 1, k)]) / dy;
          var dEy_dz := (Ey[idxEy(i, j, k)] - Ey[idxEy(i, j, k - 1)]) / dz;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hx[idxHx(i, j, k)] += -(dEz_dy - dEy_dz) / mu * dt;
        end;

    for var i := 1 to Nx - 2 do
      for var j := 0 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEx_dz := (Ex[idxEx(i, j, k)] - Ex[idxEx(i, j, k - 1)]) / dz;
          var dEz_dx := (Ez[idxEz(i, j, k)] - Ez[idxEz(i - 1, j, k)]) / dx;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hy[idxHy(i, j, k)] += -(dEx_dz - dEz_dx) / mu * dt;
        end;

    for var i := 1 to Nx - 2 do
     for var j := 1 to Ny - 2 do
        for var k := 0 to Nz - 2 do
        begin
          var dEy_dx := (Ey[idxEy(i, j, k)] - Ey[idxEy(i - 1, j, k)]) / dx;
          var dEx_dy := (Ex[idxEx(i, j, k)] - Ex[idxEx(i, j - 1, k)]) / dy;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hz[idxHz(i, j, k)] += -(dEy_dx - dEx_dy) / mu * dt;
        end;

Итоговая анимация после устранения бага. Полная энергия выходит на плато, а не растёт экспоненциально. Видна симметрия относительно вертикали. На анимации- изменение знаковой сигнатуры поля. — Итоговая анимация после устранения бага. Полная энергия выходит на плато, а не растёт экспоненциально. Видна симметрия относительно вертикали. На анимации- изменение **знаковой сигнатуры поля.**

После фикса стало видно, что энергия выходит на плато примерно к 60-му шагу и стабильна, ромбическая волна (Численная анизотропия сетки Йи https://www.academia.edu/117859596/Application_of_numerical_dispersion_compensation_of_the_Yee_FDTD_algorithm_on_elongated_domains ) симметрична относительно вертикали, затем появляется круговой контур (изотропное распространение волны от точечного источника в однородной среде). На 120 шаге кольцо расширяется и отражается от стенок. После этого рисунок выглядит хаотичным, но остаётся когерентным, из-за наложения отражённых фронтов с разными фазами, волна отражается и примерно на 225-м шаге опять получается круг (это ревайвал в закрытом резонаторе- есть в литературе как управляемая TR-фокусировка https://www.researchgate.net/publication/234903018_Time_reversed_acoustics_techniques_for_elastic_imaging_in_reverberant_and_nonreverberant_media_An_experimental_study_of_the_chaotic_cavity_transducer_concept). Затем этот круг опять расширяется и отражается от стенок.Отметим, что из-за немонотонности схемы Йи при максимальных значениях дают отклонение от физики на 17% ( динамическая ошибка- "звон"), но на знаковой сигнатуре это не видно.

Поскольку уравнения Максвелла линейны, то амплитуда источника не влияет на изображение знаковой сигнатуры. Отметим, что знак фильтрует ошибки: гасит амплитудные, обнажает фазовые. Последнее имеет и экспериментальное подтверждение-

"One-bit time reversal" и "Due to the linearity of the TR process, we find that the pulse amplitude increases linearly with the number of antennas" из https://comptes-rendus.academie-sciences.fr/physique/item/10.1016/j.crhy.2009.12.008.pdf
Авторы: Matthieu Davy , Julien de Rosny , Jean-Christophe Joly , Mathias Fink

Выводы:

Реализация FDTD‑модели в PascalABC.NET показала, что даже сравнительно сложные физические задачи можно решать и демонстрировать средствами «устаревшего» языка программирования. Использование модуля GraphABC позволило отказаться от внешней визуализации и свести весь пайплайн к одному файлу: от постановки задачи до наглядного результата.

Особое внимание при реализации было уделено корректности индексации на сетке Йи и соблюдению условия устойчивости CFL. Однако только его оказалось недостаточно, чтобы получить корректную анимацию. Необходимо было проконтролировать энергетику системы во времени. Также была проведена двойная валидация с литературой - сначала неправильного изображения, а затем исправленного.

Код получился компактным и прозрачным: его можно использовать как шаблон для дальнейших экспериментов — например, добавить неоднородности среды (для неоднородной среды требуется доработка: усреднение ε и μ на Yee‑сетке, иначе схема перестаёт быть консервативной), PML‑слои или другие типы источников.

Автор не имеет специального физико‑математического образования и не является профессиональным физиком!