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 の固有値問題

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

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

ターミナルで LaplacianEigenvalues.edp を入手して実行 ($ はプロンプトなので入力しない)
$ curl -O https://m-katsurada.sakura.ne.jp/program/fem/LaplacianEigenvalues.edp
$ FreeFem++ LaplacianEigenvalues.edp
(領域の三角形分割の後、20個の固有関数の等高線を描くので、 \fbox{enter} を20回打つ。最後はもちろん \fbox{esc} で終了する。)

\includegraphics[width=5.5cm]{freefem++/eigenfunc-fig/fig0.eps} \includegraphics[width=5.5cm]{freefem++/eigenfunc-fig/fig5.eps} \includegraphics[width=5.5cm]{freefem++/eigenfunc-fig/fig19.eps}

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

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


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

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

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

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

すなわち

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

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

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

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

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

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

ゆえに

$\displaystyle \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加筆) 説明する必要が生じて、昨日久しぶりにここを見て、 抜本的な書き換えはしないけれど (時間がない)、少し補足することに。

桂田 祐史