Алгоритм Ланцоша (собственные значения 2D-лапласиана)
Нахождение нескольких наименьших собственных значений (eigenvalues) двумерного дискретного оператора Лапласа методом Ланцоша. Собственные значения -- это специальные числа, которые характеризуют матрицу: они показывают, как оператор растягивает или сжимает векторы в разных направлениях. Программа строит подпространство Крылова, находит собственные значения трёхдиагональной матрицы с помощью бисекции и последовательности Штурма и сравнивает с аналитическими формулами.
Задача
Дан дискретный двумерный оператор Лапласа $-\Delta_h$ на сетке $N \times N$ в квадрате $[0,1]^2$ с нулевыми граничными условиями (Дирихле). Это большая разреженная матрица размера $(N-2)^2 \times (N-2)^2$ (по числу внутренних узлов). Нужно найти несколько наименьших собственных значений $\lambda$, то есть таких чисел, что:
$A v = \lambda v, \quad v \neq 0$
где $A$ -- матрица дискретного лапласиана. Прямое вычисление всех собственных значений стоит $O(N^6)$ операций, что неприемлемо для больших сеток. Метод Ланцоша позволяет найти лишь несколько нужных собственных значений гораздо быстрее.
Аналитические собственные значения дискретного лапласиана на равномерной сетке:
$\lambda_{p,q} = \frac{4}{h^2}\left(\sin^2\frac{\pi p}{2(N-1)} + \sin^2\frac{\pi q}{2(N-1)}\right), \quad p,q = 1,2,\ldots,N-2$
Идея алгоритма
- Выбираем начальный вектор $v_0$ и нормируем его.
- Строим подпространство Крылова -- последовательность векторов $v_0, Av_0, A^2 v_0, \ldots$, ортогонализуя каждый новый вектор относительно всех предыдущих (полная реортогонализация).
- В процессе получаем трёхдиагональную матрицу $T$ размера $m \times m$ с коэффициентами $\alpha$ (диагональ) и $\beta$ (наддиагональ). Собственные значения $T$ приближают собственные значения исходной матрицы $A$.
- Находим собственные значения $T$ методом бисекции с использованием последовательности Штурма (рекуррентное соотношение, которое позволяет подсчитать, сколько собственных значений меньше заданного числа $\lambda$).
- Сравниваем результат с аналитическими формулами.
Разбор
Параметры и начальный вектор
Сетка $32 \times 32$, размер подпространства Крылова $m = 100$, ищем $5$ наименьших собственных значений. Начальный вектор -- комбинация синусов для лучшего «возбуждения» нескольких мод.
цел Nx, Ny
Nx := 32
Ny := 32
вещ h, h2, pi
pi := 3.141592653589793
h := 1.0 / (Nx - 1)
h2 := h * h
| Число внутренних узлов и длина вектора состояния: (Nx-2)*(Ny-2)
цел NI, NJ, N
NI := Nx - 2
NJ := Ny - 2
N := NI * NJ
| Параметры Ланцоша
цел m, want
m := 100 | размер подпространства Крылова (больше для лучшей точности)
want := 5 | cколько наименьших значений хотим
Умножение матрицы на вектор (matvec)
Оператор $A = -\Delta_h$ действует по пятиточечному шаблону на 2D-сетке: $y_{i,j} = (4 x_{i,j} - x_{i-1,j} - x_{i+1,j} - x_{i,j-1} - x_{i,j+1}) / h^2$. Граничные значения равны нулю (Дирихле).
алг matvec_A(цел Nx, цел Ny, цел NI, цел NJ, вещ h2, вещ таб x[1:NI,1:NJ], вещ таб y[1:NI,1:NJ])
нач
| y = A x, A = (4 u - sum_neigh)/h^2, Дирихле 0 на границе
цел i, j
вещ c
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
y[i,j] := 0.0
кц
кц
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
c := 4.0 * x[i,j]
если i > 1 то c := c - x[i-1,j] все
если i < NI то c := c - x[i+1,j] все
если j > 1 то c := c - x[i,j-1] все
если j < NJ то c := c - x[i,j+1] все
y[i,j] := c / h2
кц
кц
кон
Цикл Ланцоша
На каждом шаге $k$:
- Нормируем вектор $w$ предыдущего шага: $v_1 = w / \beta_{k-1}$.
- Умножаем на оператор: $w = A \cdot v_1$.
- Вычитаем проекцию на предыдущий вектор: $w = w - \beta_{k-1} \cdot Q_{k-1}$.
- Вычисляем $\alpha_k = (v_1, w)$ и $w = w - \alpha_k \cdot v_1$.
- Полная реортогонализация: вычитаем проекции на все предыдущие базисные векторы $Q_0, \ldots, Q_k$ для численной устойчивости.
- Вычисляем $\beta_k = |w|$.
нц для k от 1 до m-1 шаг 1
если beta[k-1] = 0 то
| ранг меньше k, досрочно завершаем заполнение
alpha[k] := 0
выход
все
inv := 1.0 / beta[k-1]
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
v1[i,j] := w[i,j] * inv
Q[k,i,j] := v1[i,j]
кц
кц
matvec_A(Nx, Ny, NI, NJ, h2, v1, w)
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - beta[k-1] * Q[k-1,i,j]
кц
кц
alpha[k] := dot2(NI, NJ, v1, w)
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - alpha[k] * v1[i,j]
кц
кц
| Однопроходная реортогонализация по всем предыдущим векторам Q[0..k]
нц для r от 0 до k шаг 1
s := 0.0
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
s := s + Q[r,i,j] * w[i,j]
кц
кц
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - s * Q[r,i,j]
кц
кц
кц
если k < m-1 то
beta[k] := norm2_2d(NI, NJ, w)
все
кц
Бисекция и последовательность Штурма
Для нахождения собственных значений трёхдиагональной матрицы $T$ используется бисекция. На каждом шаге бисекции с помощью последовательности Штурма подсчитывается, сколько собственных значений лежит ниже пробного значения $\lambda$. Рекуррентное соотношение:
$d_0 = \alpha_0 - \lambda, \quad d_i = \alpha_i - \lambda - \frac{\beta_{i-1}^2}{d_{i-1}}$
Число отрицательных $d_i$ равно числу собственных значений, меньших $\lambda$.
алг цел sturm_count(цел m, вещ таб alpha[0:m-1], вещ таб beta[0:m-2], вещ lam)
нач
| Счёт Штурма через LDL' рекурсию: d0 = a0 - lam; di = ai - lam - (bi-1^2)/d{i-1}
цел i, count
вещ d, prev
count := 0
d := alpha[0] - lam
если d < 0.0 то count := count + 1 все
нц для i от 1 до m-1 шаг 1
| избегаем деления на ноль: если |d| очень мало, слегка отодвинем
если abs(d) < 1e-14 то
если d >= 0.0 то d := 1e-14 иначе d := -1e-14 все
все
d := alpha[i] - lam - (beta[i-1] * beta[i-1]) / d
если d < 0.0 то count := count + 1 все
кц
знач := count
кон
Сама бисекция ищет $t$-е по возрастанию собственное значение, сужая интервал $[a, b]$ за 80 шагов:
нц для t от 1 до want шаг 1
a := lo
b := hi
| ищем t‑е по возрастанию собственное значение
нц для j от 0 до 80 шаг 1
mid := 0.5 * (a + b)
cnt := sturm_count(m, alpha, beta, mid)
если cnt >= t то
b := mid
иначе
a := mid
все
кц
evals[t-1] := 0.5 * (a + b)
кц
Сравнение с аналитикой
Точные собственные значения для первых пяти мод вычисляются по формуле:
$\lambda_{p,q} = \frac{4}{h^2}\left(\sin^2\frac{\pi p}{2(N-1)} + \sin^2\frac{\pi q}{2(N-1)}\right)$
вещ lam11, lam12, lam21, lam13, lam22
lam11 := (4.0 / h2) * (sin(pi * 1.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 1.0 / (2.0 * (Ny - 1)))**2)
lam12 := (4.0 / h2) * (sin(pi * 1.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 2.0 / (2.0 * (Ny - 1)))**2)
lam21 := (4.0 / h2) * (sin(pi * 2.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 1.0 / (2.0 * (Ny - 1)))**2)
lam13 := (4.0 / h2) * (sin(pi * 1.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 3.0 / (2.0 * (Ny - 1)))**2)
lam22 := (4.0 / h2) * (sin(pi * 2.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 2.0 / (2.0 * (Ny - 1)))**2)
вывод "Analytic few (unordered): ", lam11, ", ", lam12, ", ", lam21, ", ", lam13, ", ", lam22, "\n"
Полная программа
алг
нач
| Ланцош для симметричного оператора A = (−Δ_h) с Дирихле в 2D на [0,1]^2.
| Цель: получить несколько наименьших собственных значений (ритц‑значения) и сравнить с аналитикой.
цел Nx, Ny
Nx := 32
Ny := 32
вещ h, h2, pi
pi := 3.141592653589793
h := 1.0 / (Nx - 1)
h2 := h * h
| Число внутренних узлов и длина вектора состояния: (Nx-2)*(Ny-2)
цел NI, NJ, N
NI := Nx - 2
NJ := Ny - 2
N := NI * NJ
| Параметры Ланцоша
цел m, want
m := 100 | размер подпространства Крылова (больше для лучшей точности)
want := 5 | cколько наименьших значений хотим
| Вектора и рабочие буферы
вещ таб v0[1:NI, 1:NJ], v1[1:NI, 1:NJ], w[1:NI, 1:NJ]
вещ таб Q[0:m-1, 1:NI, 1:NJ] | база Ланцоша (для векторов на будущее)
вещ таб alpha[0:m-1], beta[0:m-2]
цел i, j, k, t, r
вещ s, nrm, inv, a, b, lo, hi, mid, cnt
| Инициализация начального вектора (случайный/гладкий) и нормализация
| Возьмём v0[i,j] = sin(π x) sin(π y) + добавки для лучшего возбуждения мод
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
v0[i,j] := sin(pi * (i * h)) * sin(pi * (j * h)) + 0.5 * sin(pi * (2.0 * i * h)) * sin(pi * (j * h)) + 0.5 * sin(pi * (i * h)) * sin(pi * (2.0 * j * h))
кц
кц
nrm := 0.0
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
nrm := nrm + v0[i,j] * v0[i,j]
кц
кц
nrm := sqrt(nrm)
inv := 1.0 / nrm
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
v0[i,j] := v0[i,j] * inv
Q[0,i,j] := v0[i,j]
кц
кц
| Ланцош цикл
| beta[-1] = 0, v_{-1} не определён
| w = A v0
matvec_A(Nx, Ny, NI, NJ, h2, v0, w)
alpha[0] := dot2(NI, NJ, v0, w)
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - alpha[0] * v0[i,j]
кц
кц
| Однопроходная реортогонализация для численной устойчивости
s := 0.0
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
s := s + Q[0,i,j] * w[i,j]
кц
кц
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - s * Q[0,i,j]
кц
кц
если m > 1 то
beta[0] := norm2_2d(NI, NJ, w)
все
нц для k от 1 до m-1 шаг 1
если beta[k-1] = 0 то
| ранг меньше k, досрочно завершаем заполнение
alpha[k] := 0
выход
все
inv := 1.0 / beta[k-1]
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
v1[i,j] := w[i,j] * inv
Q[k,i,j] := v1[i,j]
кц
кц
matvec_A(Nx, Ny, NI, NJ, h2, v1, w)
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - beta[k-1] * Q[k-1,i,j]
кц
кц
alpha[k] := dot2(NI, NJ, v1, w)
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - alpha[k] * v1[i,j]
кц
кц
| Однопроходная реортогонализация по всем предыдущим векторам Q[0..k]
нц для r от 0 до k шаг 1
s := 0.0
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
s := s + Q[r,i,j] * w[i,j]
кц
кц
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
w[i,j] := w[i,j] - s * Q[r,i,j]
кц
кц
кц
если k < m-1 то
beta[k] := norm2_2d(NI, NJ, w)
все
| v0 <- v1 для следующей итерации
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
v0[i,j] := v1[i,j]
кц
кц
кц
| Найдём наименьшие собственные значения тридиагонали T с помощью бисекции и счёта Штурма
| Оценки спектра (Гершгорин): [lo, hi]
a := alpha[0]
b := alpha[0]
s := 0.0
нц для k от 0 до m-1 шаг 1
если alpha[k] < a то a := alpha[k] все
если alpha[k] > b то b := alpha[k] все
если k < m-1 и beta[k] > s то s := beta[k] все
кц
lo := a - 2.0 * s
hi := b + 2.0 * s
вещ таб evals[0:want-1]
вещ eps
eps := 1e-10
нц для t от 1 до want шаг 1
a := lo
b := hi
| ищем t‑е по возрастанию собственное значение
нц для j от 0 до 80 шаг 1
mid := 0.5 * (a + b)
cnt := sturm_count(m, alpha, beta, mid)
если cnt >= t то
b := mid
иначе
a := mid
все
кц
evals[t-1] := 0.5 * (a + b)
кц
| Вывод результатов и сравнение с аналитикой для первых мод (1,1), (1,2), (2,1), (1,3), (2,2)
вывод "Lanczos2D: N=", Nx, "x", Ny, ", m=", m, ", want=", want, "\n"
нц для t от 0 до want-1 шаг 1
вывод " λ[", t+1, "] ≈ ", evals[t], "\n"
кц
| Аналитические значения для 2D Лапласиана: λ_{p,q} = (4/h^2)(sin^2(π p/(2(Nx-1))) + sin^2(π q/(2(Ny-1))))
вещ lam11, lam12, lam21, lam13, lam22
lam11 := (4.0 / h2) * (sin(pi * 1.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 1.0 / (2.0 * (Ny - 1)))**2)
lam12 := (4.0 / h2) * (sin(pi * 1.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 2.0 / (2.0 * (Ny - 1)))**2)
lam21 := (4.0 / h2) * (sin(pi * 2.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 1.0 / (2.0 * (Ny - 1)))**2)
lam13 := (4.0 / h2) * (sin(pi * 1.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 3.0 / (2.0 * (Ny - 1)))**2)
lam22 := (4.0 / h2) * (sin(pi * 2.0 / (2.0 * (Nx - 1)))**2 + sin(pi * 2.0 / (2.0 * (Ny - 1)))**2)
вывод "Analytic few (unordered): ", lam11, ", ", lam12, ", ", lam21, ", ", lam13, ", ", lam22, "\n"
кон
алг вещ dot(цел N, вещ таб x[0:N-1], вещ таб y[0:N-1])
нач
вещ s
цел i
s := 0.0
нц для i от 0 до N-1 шаг 1
s := s + x[i] * y[i]
кц
знач := s
кон
алг вещ norm2(цел N, вещ таб x[0:N-1])
нач
вещ s
цел i
s := 0.0
нц для i от 0 до N-1 шаг 1
s := s + x[i] * x[i]
кц
s := sqrt(s)
знач := s
кон
алг вещ dot2(цел NI, цел NJ, вещ таб x[1:NI,1:NJ], вещ таб y[1:NI,1:NJ])
нач
вещ s
цел i, j
s := 0.0
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
s := s + x[i,j] * y[i,j]
кц
кц
знач := s
кон
алг вещ norm2_2d(цел NI, цел NJ, вещ таб x[1:NI,1:NJ])
нач
вещ s
цел i, j
s := 0.0
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
s := s + x[i,j] * x[i,j]
кц
кц
s := sqrt(s)
знач := s
кон
алг matvec_A(цел Nx, цел Ny, цел NI, цел NJ, вещ h2, вещ таб x[1:NI,1:NJ], вещ таб y[1:NI,1:NJ])
нач
| y = A x, A = (4 u - sum_neigh)/h^2, Дирихле 0 на границе
цел i, j
вещ c
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
y[i,j] := 0.0
кц
кц
нц для i от 1 до NI шаг 1
нц для j от 1 до NJ шаг 1
c := 4.0 * x[i,j]
если i > 1 то c := c - x[i-1,j] все
если i < NI то c := c - x[i+1,j] все
если j > 1 то c := c - x[i,j-1] все
если j < NJ то c := c - x[i,j+1] все
y[i,j] := c / h2
кц
кц
кон
алг цел sturm_count(цел m, вещ таб alpha[0:m-1], вещ таб beta[0:m-2], вещ lam)
нач
| Счёт Штурма через LDL' рекурсию: d0 = a0 - lam; di = ai - lam - (bi-1^2)/d{i-1}
цел i, count
вещ d, prev
count := 0
d := alpha[0] - lam
если d < 0.0 то count := count + 1 все
нц для i от 1 до m-1 шаг 1
| избегаем деления на ноль: если |d| очень мало, слегка отодвинем
если abs(d) < 1e-14 то
если d >= 0.0 то d := 1e-14 иначе d := -1e-14 все
все
d := alpha[i] - lam - (beta[i-1] * beta[i-1]) / d
если d < 0.0 то count := count + 1 все
кц
знач := count
кон