B. 2次元 Laplacian の固有値問題を FreeFem++ で解く

(ここは書きなおす必要があることが判明したけれど、 しばらくは実行する時間が取れない…書いてあることを信用しないように。)

F. Hecht, O. Pironneau, FreeFem++ a software to solve PDE に分かりやすい解説がある (リンク切れ…と思ったら、マニュアルに入った…と安心していたらなくなった …と思ったら復活した)。 次に掲げるプログラムはそこに載っているもの (ただし変数名を少し書き換えた)である。

LaplacianEigenvalues.edp

// LaplacianEigenvlaues.edp
// original: https://doc.freefem.org/models/eigen-value-problems.html
// renamed some variables by mk (2024/8/29)

verbosity=0;
real shift = 20; //value of the shift (near 2 pi^2)
int nev = 20; //number of computed eigenvalues close to shift

// Mesh
mesh Th = square(20, 20, [pi*x, pi*y]);

// Fespace
fespace Vh(Th, P2);
Vh u, v;

// Problem
// AmsB = A - shift B ; // the shifted matrix (A minus shift B)
varf op (u, v)
    = int2d(Th)(
          dx(u)*dx(v)
        + dy(u)*dy(v)
        - shift* u*v
    )
    + on(1, 2, 3, 4, u=0)
    ;

// from the right hand side of -Laplacian u=lambda u
varf b (u, v) = int2d(Th)(u*v); //no boundary condition

matrix AmsB = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
matrix B = b(Vh, Vh, solver=CG, eps=1e-20);

// important remark:
// the boundary condition is make with exact penalization:
// we put 1e30=tgv on the diagonal term of the lock degree of freedom.
// So take Dirichlet boundary condition just on $a$ variational form
// and not on $b$ variational form.
// because we solve $ w=AmsB^-1*B*v $

// Solve
real[int] ev(nev); //to store the nev eigenvalue
Vh[int] eV(nev); //to store the nev eigenvector

int k = EigenValue(AmsB, B, sym=true, sigma=shift, value=ev, vector=eV,
    tol=1e-10, maxit=0, ncv=0);

// Display & Plot
for (int i = 0; i < k; i++){
    u = eV[i];
    real gg = int2d(Th)(dx(u)*dx(u) + dy(u)*dy(u));
    real mm = int2d(Th)(u*u) ;
    cout << "lambda[" << i << "] = " << ev[i] << ", err= " << int2d(Th)(dx(u)*dx(u) + dy(u)*dy(u) - (ev[i])*u*u) << endl;
    plot(eV[i], cmm="Eigenvector "+i+" value ="+ev[i], wait=true, value=true);
}

正方形領域 $\Omega=(0,1)\times(0,1)$ における、 Dirichlet 境界条件下の Laplacian の固有値問題

$$-\Laplacian u=\lambda u \Omega ,\quad u=0 \rd\Omega$$ (5)

について、 小さい方から20個の固有値を計算し、 対応する固有関数の等高線を描くプログラムである。

ターミナルで LaplacianEigenvalues.edp を入手して実行 ($ はプロンプトなので入力しない)

$ curl -O https://m-katsurada.sakura.ne.jp/program/fem/LaplacianEigenvalues.edp
$ FreeFem++ LaplacianEigenvalues.edp
(領域の三角形分割の後、20個の固有関数の等高線を描くので、 を20回打つ。最後はもちろん で終了する。)

なお、8行目の +on(1,2,3,4,u=0) を削除すると、 Dirichlet 境界条件の代わりに、 自然境界条件である Neumann 境界条件となる。

$$-\Laplacian u=\lambda u \Omega ,\quad \frac{\rd u}{\rd n}=0 \rd\Omega .$$ (6)

さて、 Dirichlet 境界条件の場合 (5)に、弱形式を導こう。

$$-\Laplacian u=\lambda u$$   $$\mbox{in $\Omega$}$$$$,\quad u=0$$   $$\mbox{(on $\rd\Omega$)}$$

の両辺に試験関数 $v$ をかけて積分し、部分積分すると (Greenの積分公式を用いて変形すると)

$$\dint_\Omega\nabla u\cdot\nabla v\;\DxDy=\lambda\int_\Omega u v\;\DxDy.$$

すなわち

$$\dint_\Omega\left(u_x v_x+u_y v_y\right)\DxDy=\lambda\dint_\Omega u v\;\DxDy.$$

これを有限要素法で離散化すると、(行列の)一般化固有値問題と呼ばれる

$$A \Vector{u}=\lambda B\Vector{u}$$

という形をした方程式が得られる。 ここで $B$ は正値対称行列、$A$ は実対称行列である。

「シフト法で解く」と簡単に書いてあるだけで、詳しいことは書いてない。 次のようなことかと推測する。 一般化固有値に近いと考えられる実数 $\sigma$ を与えたとき、 $A\Vector{u}=\lambda B\Vector{u}$ の両辺から $\sigma B\Vector{u}$ を引くと、

$$(A-\sigma B)\Vector{u}=(\lambda-\sigma)B\Vector{u}.$$

ゆえに

$$\left(A-\sigma B\right)^{-1} B\Vector{u}=\frac{1}{\lambda-\sigma}\Vector{u}.$$

これは $\Vector{u}$ が、行列 $\left(A-\sigma B\right)^{-1}$ の、 固有値 $\dfrac{1}{\lambda-\sigma}$ に属する固有ベクトルであることを示す。 ゆえに行列 $\left(A-\sigma B\right)^{-1}B$ について冪乗法を実行すれば、 絶対値が最大であるような固有値 $\mu=\dfrac{1}{\lambda-\sigma}$ が得られる。

細かいことだが、重要な注意がある。 ここでは square() を用いてメッシュ分割をしているが、 それを用いずに、自分で境界曲線を定義してメッシュ分割をすると、 真の固有値が重複(縮重)する場合であっても、 数値計算で得られる近似固有値は重複(縮重)しない。 これは自動メッシュ分割で得られる三角形分割は、 正方形の二面体群 (合同変換全体のなす群) $D_4$ の対称性を持たないため、 近似固有値が重複しないためである。

(2021/1/9加筆) 説明する必要が生じて、昨日久しぶりにここを見て、 抜本的な書き換えはしないけれど (時間がない)、少し補足することに。

• 行列の固有値問題について、 “一般化固有値問題” というキーワードを知っている必要がある。 また有限要素法により $A\bm{u}=\lambda B\bm{u}$, $A=((\nabla\phi_i, \nabla\phi_j))$, $B=\left((\phi_i, \phi_j)\right)$ という行列の一般化固有値問題に帰着される、 ということも分かっている必要がある。
• 行列の固有値問題についてシフト法 (FreeFem++ の Documentation の表現では、 the shifted-inverse method) について知っている必要がある。
• FreeFem++ の EigenValue() という関数については、まず、 現在の Documentation の §§4.7.46 EigenValue を見ると良い。 そこにParameters (引数と訳すべきか) の説明がある。ただし、 maxit, ncv については説明がない。 “Other parameters are available for more advanced use: see the ARPACK documentation” とあるので、 固有値計算ライブラリィ ARPACK を調べに行く必要がある (ちょっと面倒…サボります)。