这至少是一个想法，它是否有效是一个不同的问题。

假设您对未知数进行排序，以便首先拥有域内部的未知数，然后是边界内的所有未知数。然后对应于您的问题的矩阵按以下方式分解： 其中，表示内部的形状函数，表示域的边界。只是拉普拉斯算子的狄利克雷边界条件矩阵，我们知道如何有效地反转它。然后可以考虑使用 Schur 补码方法，首先求解边界未知数，然后求解内部未知数，或者反过来。

$$A = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ C^{\partial\circ} & D^{\partial\partial} \end{pmatrix}$$ $\circ$$\partial$$A^{\circ\circ}$

如果您遵循 Stokes 系统的 Silvester-Wathen 预条件子的论点，您还可以基于 Schur 补集方法为整个矩阵构造一个好的预条件子；它可能看起来像这样： 其中是 Schur 补码。值得思考是否可以通过更简单的矩阵逼近，就像斯托克斯方程一样，可以通过质量矩阵逼近 Schur 补码。

$$P^{-1} = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ 0 & S^{\partial\partial} \end{pmatrix}^{-1}$$ $S^{\partial\partial}=D^{\partial\partial}-D^{\partial\circ}[A^{\circ\circ}]^{-1}B^{\circ\partial}$$S^{\partial\partial}$

正如其他人指出的那样，在这种情况下，（代数）多重网格实际上可以成为一个很好的预处理器。下面是scikit-fem和pyamg的概念验证实现。它表明，pyamg 的预处理器使 GMRES 迭代的次数或多或少独立于未知数的数量。

我只用几千个未知数的网格对此进行了测试。主要的计算成本是计算 的左侧零空间A，用于使右侧一致。（在保证一致性的实际应用中不需要的计算。）

import matplotlib.pyplot as plt
import meshzoo
import numpy as np
import pyamg
import scipy.linalg
import scipyx
import skfem as fem
from skfem.helpers import dot
from skfem.models.poisson import laplace

rng = np.random.default_rng(0)

tol = 1.0e-8

for n in range(5, 41, 5):
    # points, cells = meshzoo.rectangle_tri((0.0, 0.0), (1.0, 1.0), n)
    points, cells = meshzoo.disk(6, n)
    print(f"{n = }, {len(points) = }")

    @fem.BilinearForm
    def flux(u, v, w):
        return dot(w.n, u.grad) * v

    mesh = fem.MeshTri(points.T.copy(), cells.T.copy())
    basis = fem.InteriorBasis(mesh, fem.ElementTriP1())
    facet_basis = fem.FacetBasis(basis.mesh, basis.elem)

    lap = fem.asm(laplace, basis)
    boundary_terms = fem.asm(flux, facet_basis)

    A = lap - boundary_terms

    b = rng.random(A.shape[1])
    # make the system consistent by removing A's left nullspace components from the
    # right-hand side
    lns = scipy.linalg.null_space(A.T.toarray()).T
    for n in lns:
        b -= np.dot(b, n) / np.dot(n, n) * n

    ml = pyamg.smoothed_aggregation_solver(
        A,
        coarse_solver="jacobi",
        symmetry="nonsymmetric",
        max_coarse=100,
    )
    M = ml.aspreconditioner(cycle="V")

    _, info = scipyx.gmres(A, b, tol=tol, M=M, maxiter=20)
    # res = b - A @ info.xk

    num_unknowns = A.shape[1]
    plt.semilogy(
        np.arange(len(info.resnorms)), info.resnorms, label=f"N={num_unknowns}"
    )

plt.xlabel("step")
plt.ylabel("residual")
plt.legend()
plt.savefig("out.png", bbox_inches="tight")
plt.show()

如果你没有强制任何边界条件，你应该删除问题的零空间（为了确保有唯一的解决方案，在 PETSc 中可以使用 MatSetNullSpace 来实现），然后在问题上使用多重网格（例如，PCMG PETSc)。

令我惊讶的是，有 n+1 个零特征值。我希望只有一个。您是否正在求解表面上的拉普拉斯方程（如嵌入 3D 空间中的 2D 表面）并且边界是该表面的边界（如嵌入 3D 空间中的 1D 边界）？这是我能想象的唯一一种情况，即零空间人口众多。