反转奇异矩阵（拉普拉斯算子到逆拉普拉斯算子）

Question

我正在将用 Octave/Matlab 编写的二维 CFD 代码移植到 Fortran。域是周期性的，因此该方案基于FT。下面的矩阵，'laplacian'：

0 -1 -4 -9 -4 -1

-1 -2 -5 -10 -5 -2

-4 -5 -8 -13 -8 -5

-9 -10 -13 -18 -13 -10

-4 -5 -8 -13 -8 -5

-1 -2 -5 -10 -5 -2

表示 6 x 6 网格上 FT 的拉普拉斯算子。我想要逆矩阵，即使拉普拉斯算子是奇异的。在 matlab/octave、'inv(laplacian)' returns 所有 'Inf' 中，但是 '1./laplacian' returns 正确答案（尽管是 (1,1) 元素，返回作为 Inf，必须设置为零）。

问题是如何使用LAPACK翻译第二种形式。我通常使用的矩阵求逆序列 'DGETRF/DGETRI' 因信息=4 而失败，不足为奇。还有大约二十个其他 DxxTRF。有谁知道什么可能有机会做 Octave 所做的事情？

Answer 1

如果你正在做我认为你正在做的事情，那么你想将傅立叶变换的不同波数乘以从拉普拉斯算子的特征值导出的不同系数。

类似于

lambda(kx, ky, kz) = (kx**2 + ky**2 + kz**2)

注意 "matrix" 中的 1、4 和 9。它们是这些正方形 kx**2.

这不是矩阵求逆，实际上只是将 1.0 除以 table 形式的数字。 table 看起来像一个矩阵，因为你的代码是二维的，所以你只有 lambda(kx, ky).

拉普拉斯算子的整个实际矩阵将非常大（N 乘以 N，其中 2D 中的 N=nx*ny 和 3D 中的 N=nx*ny*nz）并且具有 lambdas 在对角线上，其他地方都为零。逆矩阵在对角线上有 1./labdas。所以你的操作在某种意义上是矩阵求逆，但与你想象的不同。

你要做的是

FT(kx,ky,kz) = FT(kx,ky,kz) / (kx**2 + ky**2 + kz**2)

也可以写成

FT = FT * inv_laplacian

哪里

inv_laplacian = 1. / laplacian

其中 laplacian 是系数（特征值）数组。

这不是矩阵求逆，它只是将 1.0 除以数字。

现在，特征值 0 有什么用？如果您不启用浮点异常捕获，我建议您不要启用它们，那么您只需这样做：

 FT = FT * inv_laplacian

因为 inv_laplacian(0,0,0) 为 0，FT(0,0,0) 除以 0 且未定义（NaN 或类似）。您可以将其设置为任何您想要的。 FT(0,0,0)的意思就是你这个字段的平均值，随意。所以只是

 FT(0,0,0) = 0 !or any number you want

就是这样。

顺便说一句，我在现实世界的科学代码中看到过极端的做法，例如：

for i =0, nkx-1
  for j =0, nkx-1
    for k =0, nkx-1
      FT(i, j, k) = FT(i, j, k) / (cos(i*ax) + cos(j*ay) + cos(k*az)
    end 
  end
end

计算了 AGES。这与你的情况非常相似，你的特征值不是余弦而是平方。

重点是系数在时间上是恒定的并且是可分离的。

一个人只能计算一次：

lambdax(i) = i**2 !or cos(ax*i)
lambday(i) = j**2 !or cos(ay*j)
lambdaz(i) = k**2 !or cos(az*k)

然后

FT(kx,ky,kz) = FT(kx,ky,kz) / (lambdax(kx) + lambday(ky) + lambdaz(kz))

您可以查看我的快速泊松求解器 PoisFFT 的源代码以查看示例 https://github.com/LadaF/PoisFFT/blob/master/src/poisfft-solvers-inc.f90 并找到合适的边界条件。您的情况可能 PoisFFT_Solver2D_FullPeriodic 第 152 行。

反转奇异矩阵（拉普拉斯算子到逆拉普拉斯算子）

Invert a singular matrix (Laplacian to inverse Laplacian)

matlab

fortran

matrix

octave

lapack