Чисельне рівняння Шредінгера у 1D: розв'язувач Крank-Ніколсон

1. Дискретизація простору

Поділіть область [x_min, x_max] на N вузлів сітки з кроком dx = (x_max − x_min)/(N−1). Комплексна хвильова функція ψⱼ = ψ(xⱼ, t) зберігається як два паралельних буфери Float64Array, дійсна та уявна частини — ніколи як єдиний комплексний тип, оскільки в JavaScript такого немає:

const N  = 400;
const xMin = -20, xMax = 20;
const dx = (xMax - xMin) / (N - 1);
const x  = Float64Array.from({ length: N }, (_, j) => xMin + j * dx);
const psiRe = new Float64Array(N);
const psiIm = new Float64Array(N);

// Початковий гаусів хвильовий пакет: psi(x,0) = exp(-(x-x0)^2/(4*sigma^2)) * exp(i*k0*x)
function initGaussian(x0, sigma, k0) {
  let norm = 0;
  for (let j = 0; j < N; j++) {
    const env = Math.exp(-(x[j] - x0) ** 2 / (4 * sigma * sigma));
    psiRe[j] = env * Math.cos(k0 * x[j]);
    psiIm[j] = env * Math.sin(k0 * x[j]);
    norm += (psiRe[j] ** 2 + psiIm[j] ** 2) * dx;
  }
  const invSqrtNorm = 1 / Math.sqrt(norm);
  for (let j = 0; j < N; j++) { psiRe[j] *= invSqrtNorm; psiIm[j] *= invSqrtNorm; }
}

2. Гамільтоніан скінченних різниць

Друга просторова похідна в Ĥ = −ℏ²/(2m)∂²/∂x² + V(x) наближується стандартною трьохточковою центральною різницею:

Записаний як матриця, Ĥ — тридіагональний: −α на двох поза-діагоналях, (2α + Vⱼ) на діагоналі. Ця розрідженість — уся причина, чому вся схема лишається швидкою: прямий розв'язок з щільною матрицею коштував би O(N³) на крок часу; тридіагональна структура знижує це до O(N).

3. Чому явний метод Ейлера не працює

Очевидна перша спроба — ψ(t+dt) = ψ(t) − i(dt/ℏ)Ĥψ(t) — це явний (прямий) метод Ейлера, застосований до iℏ∂ψ/∂t = Ĥψ. Він безумовно нестабільний: оскільки Ĥ ермітовий з дійсними власними значеннями E, точний фактор еволюції в часі e^−iEdt/ℏ має модуль рівно 1, але його наближення Тейлора першого порядку 1 − iEdt/ℏ має модуль √(1+(Edt/ℏ)²) > 1 для будь-якого ненульового dt. Кожна мода зростає експоненційно, яким би малим не був крок часу.

4. Схема Крank-Ніколсон

Крank-Ніколсон усереднює дію гамільтоніана на старому та новому кроках часу, даючи схему, яка другого порядку точності за dt і точно унітарна (зберігає норму) для будь-якого кроку часу:

Права частина — простий тридіагональний добуток матриця-вектор (явний, дешевий). Ліва частина вимагає розв'язання тридіагональної лінійної системи на кожному кроці — все ще O(N), але справжнє розв'язання, а не пряме обчислення.

5. Розв'язання тридіагональної системи: алгоритм Томаса

Тридіагональна система Ax = d (піддіагональ a, діагональ b, наддіагональ c) розв'язується за один прохід вперед-назад, алгоритмом Томаса — методом Гаусса, спеціалізованим для ширини смуги 1:

// Комплексний тридіагональний розв'язок алгоритмом Томаса.
// a, b, c, d — масиви комплексних чисел {re, im}; b/d перезаписуються.
function thomasSolve(a, b, c, d, N) {
  // Прямий прохід: усунення піддіагоналі
  for (let j = 1; j < N; j++) {
    const w = cdiv(a[j], b[j - 1]);         // w = a[j] / b[j-1]
    b[j] = csub(b[j], cmul(w, c[j - 1]));  // b[j] -= w * c[j-1]
    d[j] = csub(d[j], cmul(w, d[j - 1]));  // d[j] -= w * d[j-1]
  }
  // Зворотна підстановка
  const x = new Array(N);
  x[N - 1] = cdiv(d[N - 1], b[N - 1]);
  for (let j = N - 2; j >= 0; j--) {
    x[j] = cdiv(csub(d[j], cmul(c[j], x[j + 1])), b[j]);
  }
  return x;
}
// cmul/cdiv/csub: тривіальні хелпери комплексної арифметики на об'єктах {re, im}
function cmul(p, q) { return { re: p.re*q.re - p.im*q.im, im: p.re*q.im + p.im*q.re }; }
function csub(p, q) { return { re: p.re - q.re, im: p.im - q.im }; }
function cdiv(p, q) { const d2 = q.re*q.re + q.im*q.im; return { re: (p.re*q.re + p.im*q.im)/d2, im: (p.im*q.re - p.re*q.im)/d2 }; }

6. Граничні умови та потенціали

Два вибори меж поширені. Жорсткі стінки фіксують ψ(x_min) = ψ(x_max) = 0 для всіх t (еквівалентно нескінченній прямокутній ямі, що містить усю сітку) — найпростіше реалізувати, але хвильові пакети нефізично відбиваються від країв області. Поглинаючі межі додають плавний уявний потенціал −iΓ(x) поблизу країв («комплексний поглинаючий потенціал», CAP), що гасить ψ до того, як вона досягне межі, імітуючи відкриту, нескінченну область:

// Потенціал: скінченна прямокутна яма/бар'єр + поглинаючі краї
function buildPotential(x, N, wellDepth, wellWidth, gamma0 = 2.0, absorbWidth = 3) {
  const V = new Float64Array(N);
  const Vi = new Float64Array(N);  // уявна частина: поглинаючий шар
  for (let j = 0; j < N; j++) {
    V[j] = Math.abs(x[j]) < wellWidth / 2 ? -wellDepth : 0;
    const distFromEdge = Math.min(x[j] - x[0], x[N-1] - x[j]);
    if (distFromEdge < absorbWidth) {
      Vi[j] = -gamma0 * (1 - distFromEdge / absorbWidth) ** 2;
    }
  }
  return { V, Vi };
}

7. Повний робочий код

Зібравши все разом: побудуйте тридіагональний Ĥ, зберіть матриці (I ± iĤdt/2ℏ), і рухайтесь на крок вперед кожен кадр.

function createSolver({ N, dx, hbar, mass, V, dt }) {
  const alpha = hbar * hbar / (2 * mass * dx * dx);
  const r = { re: 0, im: dt / (2 * hbar) };   // i*dt/(2*hbar) розкладено на re/im

  // Діагональ H (тридіагональна): 2*alpha + V[j]; поза-діагоналі: -alpha
  const diagH = Array.from({ length: N }, (_, j) => 2 * alpha + V[j]);

  function step(psiRe, psiIm) {
    // Права частина: (I - i*H*dt/2hbar) psi, тридіагональний matvec
    const dRe = new Float64Array(N), dIm = new Float64Array(N);
    for (let j = 0; j < N; j++) {
      let hRe = diagH[j] * psiRe[j], hIm = diagH[j] * psiIm[j];
      if (j > 0)     { hRe -= alpha * psiRe[j-1]; hIm -= alpha * psiIm[j-1]; }
      if (j < N - 1) { hRe -= alpha * psiRe[j+1]; hIm -= alpha * psiIm[j+1]; }
      // psi - i*(dt/2hbar)*H*psi = psi_re + H_im*(dt/2hbar), psi_im - H_re*(dt/2hbar)
      dRe[j] = psiRe[j] + r.im * hIm;
      dIm[j] = psiIm[j] - r.im * hRe;
    }
    // Ліва матриця (I + i*H*dt/2hbar): тридіагональна, розв'язується алгоритмом Томаса
    const a = Array.from({ length: N }, () => ({ re: 0, im: -r.im * -alpha }));
    const b = Array.from({ length: N }, (_, j) => ({ re: 1, im: r.im * diagH[j] }));
    const c = a.slice();
    const d = Array.from({ length: N }, (_, j) => ({ re: dRe[j], im: dIm[j] }));
    const x = thomasSolve(a, b, c, d, N);
    for (let j = 0; j < N; j++) { psiRe[j] = x[j].re; psiIm[j] = x[j].im; }
  }
  return { step };
}

8. Перевірка точності: нормування та енергія

Дві дешеві діагностики виловлюють майже кожну помилку дискретизації до того, як вона стане візуальним артефактом:

• Нормування: Σ|ψⱼ|²dx має лишатися надзвичайно близьким до 1 кожен кадр (дрейф понад ~10⁻⁶ за крок без поглинаючої межі означає помилку знаку, несиметричну матрицю, або dt завеликий відносно найтоншого енергетичного масштабу в Ĥ).
• Енергія: ⟨Ĥ⟩ = Σψⱼ*(Ĥψ)ⱼdx має бути сталою в часі для замкненої системи (без поглинаючого шару, без явної залежності V від часу) — Крank-Ніколсон точно зберігає це в точній арифметиці, тож дрейф тут сигналізує про помилку коду, а не обмеження схеми.