反轉奇異矩陣（拉普拉斯逆拉普拉斯算子）

Question

我將用Octave/Matlab編寫的二維CFD碼移植到fortran中。該域是週期性的，所以該方案基於FT。下面的矩陣， '拉普拉斯'：

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乘6網格上的FT的拉普拉斯算子。我想要逆矩陣，儘管拉普拉斯算子是單數的。在matlab/octave中，'inv（laplacian）'返回所有的'Inf'，然而'1./​​laplacian'返回正確答案（儘管返回Inf的（1,1）元素必須設置爲零）。

問題是如何使用LAPACK翻譯第二個表單。我常用的矩陣反轉序列'DGETRF/DGETRI'失敗，信息= 4，沒有什麼意外。還有另外20個DxxTRF。有沒有人知道什麼可能會做Octave做的事情？

Answer 1

如果你正在做我認爲你是的那麼你想要乘以不同的波數傅立葉變換乘以拉普拉斯算子的特徵值導出的不同係數。

喜歡的東西

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

通知的1,4和9中的 「矩陣」。他們是這些正方形kx**2。

這不是一個矩陣求逆，它實際上只是用表格形式寫出的數字除以1.0。該表看起來像一個矩陣，因爲你的代碼是2D的，所以你只有lambda(kx, ky)。

爲拉普拉斯算子，整個實際矩陣將是非常大的（N倍N，在2D其中N = nx*ny和N = nx*ny*nz在3D），希望對角和零其他地方的lambda秒。逆矩陣將在對角線上具有1./labda。因此，您的操作在某種意義上是矩陣求逆，但與您的想法不同。

你比什麼做是

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 

就是這樣。

---

BTW，我看到了一個真實的世界科學的代碼，甚至極端的做法，如：

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 

它經歷了千百年來計算。這與你的情況非常相似，你的特徵值不是餘弦，而是正方形。

問題是係數在時間上是恆定的，可分離的。

人們可以只計算一次：

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並找到適當的邊界條件。您的情況很可能是第152行的PoisFFT_Solver2D_FullPeriodic。