Потенціал електричного поля та метод релаксації

1. Рівняння Лапласа та Пуасона

Електричний потенціал φ пов'язаний з полем E = −∇φ. З рівняння Гаусса ∇·E = ρ/ε₀ отримуємо:

Рівняння Лапласа описує потенціал у просторі без зарядів — між провідниками з фіксованими потенціалами. Рівняння Пуасона — загальний випадок (наявні заряди). Обидва відносяться до еліптичного типу: розв'язок "гладкий" і залежить від граничних умов по всій межі.

Фізичний зміст

2. Дискретизація методом скінченних різниць

Розбиваємо область на рівномірну сітку з кроком h у 2D. Оператор Лапласа апроксимуємо п'ятиточковим шаблоном (stencil):

Підставляємо в рівняння Пуасона ∇²φ = f (правій частині = −ρ/ε): отримуємо систему лінійних рівнянь A·φ = b, де A — розріджена матриця (5 ненульових елементів у рядку для 2D). Розмірність N × N при сітці N×N. Пряме розв'язання (LU) дорого — O(N³). Ітераційні методи — O(N²) за ітерацію, і зазвичай конвергують за O(N) – O(N log N) ітерацій.

3. Метод Якобі

Найпростіший ітераційний метод: виражаємо φ_{i,j} через сусідів, використовуючи старі значення:

Збіжність: спектральний радіус ρ_J ≈ 1 − π²/(2N²) для квадратної сітки N×N. Кількість ітерацій для похибки ε: k ≈ N² ln(1/ε) / π² ∝ N². Для N=100 і ε=10⁻⁶ — ~28 000 ітерацій. Повільно!

Чому у WebGL зручно реалізувати Якобі

Кожен піксель текстури = один вузол сітки. Якобі ідеально паралельний — нові значення не залежать одне від одного. Один прохід = один drag-call у ping-pong framebuffer.

// GLSL fragment shader — одна ітерація Якобі
precision highp float;
uniform sampler2D uPhi;   // поточне наближення φ
uniform sampler2D uRho;   // щільність заряду ρ (нормована)
uniform vec2 uTexelSize;  // 1/width, 1/height
uniform float uH2;        // h² (крок сітки у квадраті)
varying vec2 vUV;

void main() {
  float n = texture2D(uPhi, vUV + vec2(0., uTexelSize.y)).r;
  float s = texture2D(uPhi, vUV - vec2(0., uTexelSize.y)).r;
  float e = texture2D(uPhi, vUV + vec2(uTexelSize.x, 0.)).r;
  float w = texture2D(uPhi, vUV - vec2(uTexelSize.x, 0.)).r;
  float rho = texture2D(uRho, vUV).r;
  gl_FragColor = vec4((n + s + e + w + uH2 * rho) / 4.0, 0., 0., 1.);
}

4. Метод Гаусса-Зейделя

Відмінність від Якобі: одразу використовуємо нові значення сусідів, як тільки вони отримані (лівий і нижній сусіди вже оновлені):

Збіжність вдвічі швидша Якобі: ρ_GS ≈ ρ_J². На практиці вдвічі менша кількість ітерацій при тій самій похибці. Складність реалізації на GPU: залежність між пікселями — потрібен шаховий порядок оновлення ("red-black Gauss-Seidel").

Red-Black ordering

Розбиваємо сітку на "червоні" (i+j парне) і "чорні" (i+j непарне) вузли. Червоні залежать тільки від чорних і навпаки. Два проходи за ітерацію:

function redBlackGS(phi, rho, N, h2, iterations) {
  for (let it = 0; it < iterations; it++) {
    for (let color = 0; color < 2; color++) {  // 0=red, 1=black
      for (let i = 1; i < N - 1; i++) {
        for (let j = 1; j < N - 1; j++) {
          if ((i + j) % 2 !== color) continue;
          const idx = i * N + j;
          phi[idx] = (
            phi[(i+1) * N + j] + phi[(i-1) * N + j] +
            phi[i * N + j+1] + phi[i * N + j-1] +
            h2 * rho[idx]
          ) / 4;
        }
      }
    }
  }
}

5. Метод наддосягань SOR

SOR (Successive Over-Relaxation) покращує Гаусса-Зейделя введенням параметра релаксації ω ∈ (1, 2):

З ω = ω_opt кількість ітерацій зменшується з O(N²) до O(N) — принципово інший порядок! Для N=100 це зменшення з ~28 000 до ~300 ітерацій.

class PoissonSolverSOR {
  constructor(N, h = 1.0 / N) {
    this.N = N;
    this.h2 = h * h;
    this.omega = 2 / (1 + Math.sin(Math.PI / N)); // ω_opt
    this.phi = new Float64Array(N * N);  // потенціал
    this.rho = new Float64Array(N * N);  // ρ/ε₀
    this.fixed = new Uint8Array(N * N);  // маска ГУ Діріхле
  }

  setCharge(i, j, value) { this.rho[i * this.N + j] = -value; }
  setFixed(i, j, value) {
    const idx = i * this.N + j;
    this.phi[idx] = value;
    this.fixed[idx] = 1;
  }

  step(iterations = 1) {
    const { N, h2, omega, phi, rho, fixed } = this;
    for (let it = 0; it < iterations; it++) {
      for (let color = 0; color < 2; color++) {
        for (let i = 1; i < N - 1; i++) {
          for (let j = 1; j < N - 1; j++) {
            if ((i + j) % 2 !== color) continue;
            const idx = i * N + j;
            if (fixed[idx]) continue;
            const gs = (
              phi[(i+1)*N+j] + phi[(i-1)*N+j] +
              phi[i*N+j+1] + phi[i*N+j-1] - h2 * rho[idx]
            ) / 4;
            phi[idx] += omega * (gs - phi[idx]);
          }
        }
      }
    }
  }

  residual() {
    const { N, h2, phi, rho, fixed } = this;
    let max = 0;
    for (let i = 1; i < N - 1; i++) {
      for (let j = 1; j < N - 1; j++) {
        const idx = i * N + j;
        if (fixed[idx]) continue;
        const lap = (
          phi[(i+1)*N+j] + phi[(i-1)*N+j] +
          phi[i*N+j+1] + phi[i*N+j-1] - 4*phi[idx]
        ) / h2;
        max = Math.max(max, Math.abs(lap + rho[idx]));
      }
    }
    return max;
  }
}

6. Граничні умови Діріхле та Неймана

Ghost cells для умови Неймана

7. Ізолінії та силові лінії: JS реалізація

Ізолінії потенціалу (equipotential)

Ізолінія рівня c: множина точок де φ(x,y) = c. Алгоритм Marching Squares перебирає квадрати сітки та визначає конфігурацію перетину рівня:

function marchingSquaresIsolines(grid, N, levels) {
  const segments = [];
  for (const level of levels) {
    for (let i = 0; i < N - 1; i++) {
      for (let j = 0; j < N - 1; j++) {
        const v = [
          grid[i*N+j], grid[i*N+j+1],
          grid[(i+1)*N+j], grid[(i+1)*N+j+1]
        ];
        const code = (+(v[0]>level)) | (+(v[1]>level))<<1
                    | (+(v[2]>level))<<2 | (+(v[3]>level))<<3;
        if (code === 0 || code === 15) continue; // всередині або зовні
        const pts = msLookup(code, v, level, i, j); // таблиця переходів
        if (pts) segments.push(...pts);
      }
    }
  }
  return segments;
}

Силові лінії поля

Силова лінія — крива, дотична до вектора E = −∇φ у кожній точці. Інтеграція методом RK4 від початкової точки (або набору точок навколо зарядів):

function traceFieldLine(phi, N, startI, startJ, stepSize = 0.5, maxSteps = 2000) {
  const path = [[startI, startJ]];
  let [ci, cj] = [startI, startJ];

  const grad = (i, j) => {
    const gx = (phi[Math.round(i+1)*N + Math.round(j)]
               - phi[Math.round(i-1)*N + Math.round(j)]) / 2;
    const gy = (phi[Math.round(i)*N + Math.round(j+1)]
               - phi[Math.round(i)*N + Math.round(j-1)]) / 2;
    return [gx, gy];
  };

  for (let s = 0; s < maxSteps; s++) {
    const [gx, gy] = grad(ci, cj);
    const mag = Math.hypot(gx, gy);
    if (mag < 1e-6) break;
    ci -= stepSize * gx / mag; // E = -∇φ → пропускаємо знак (лінія від +)
    cj -= stepSize * gy / mag;
    if (ci < 0 || ci >= N || cj < 0 || cj >= N) break;
    path.push([ci, cj]);
  }
  return path;
}

8. Мультирозв'язок та застосування

Мультигрід (Multigrid)

Причина повільності SOR: видаляє лише "короткохвильові" помилки за ітерацію, "довгохвильові" (гладкі) зберігаються. Мультигрід виправляє це: переходить на грубіші сітки, де довгохвильові компоненти стають короткохвильовими і швидко пригнічуються (V-cycle або W-cycle):