有限要素法による非時間依存2次元シュレディンガー方程式

以前, 有限要素法による非時間依存1次元シュレディンガー方程式で1次元の場合を解いたので, 今回は2次元の場合を解きます.

2次元の時間に依存しないシュレディンガー方程式は

$$\left(-\frac{\hbar^2}{2m}\nabla^2+V(x,\,y)\right)\phi(x,\,y)=E\phi(x,\,y)$$

と表せます. ただし, $\displaystyle\nabla^2=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$です.

計算編

1.弱形式を求める

1次元のときと同様に弱形式を求めるところから始めます. 冒頭のシュレディンガー方程式の両辺に, 波動関数の複素共役 $\phi^*(x,\,y)$ を重み関数としてかけ, 領域 $\Omega$ 内で積分します.

$$-\frac{\hbar^2}{2m}\int_\Omega\phi^*\nabla^2\phi\;dS+\int_\Omega\phi^*V\phi\;dS=E\int_\Omega\phi^*\phi\;dS\;\;\;(1.1)$$

式(1.1)左辺第1項は

$$\nabla\cdot(f\nabla g)=(\nabla f)\cdot (\nabla g)+f\nabla^2 g$$

の関係式を用いて

$$-\frac{\hbar^2}{2m}\int_\Omega\phi^*\nabla^2\phi\;dS=-\frac{\hbar^2}{2m}\left(-\int_\Omega(\nabla\phi^*)\cdot(\nabla\phi)\;dS+\int_\Omega\nabla\cdot(\phi^*\nabla\phi)\;dS\right)$$

さらに, ガウスの発散定理

$$\int_\Omega\nabla\cdot \bm{A}\;dS=\int_{\partial\Omega}\bm{A}\cdot\bm{n}\;ds$$

を用いると（$\partial\Omega$ は領域 $\Omega$ の境界を表します）

$$-\frac{\hbar^2}{2m}\int_\Omega\phi^*\nabla^2\phi\;dS=-\frac{\hbar^2}{2m}\left(-\int_\Omega(\nabla\phi^*)\cdot(\nabla\phi)\;dS+\int_{\partial\Omega}\phi^*(\nabla\phi)\cdot\bm{n}\;ds\right)$$

と書けます. Dirichlet境界条件 $\forall(x,\;y)\in\partial\Omega,\;\phi(x,\;y)=0$ より, 境界上で $\phi^*=0$ なので

$$-\frac{\hbar^2}{2m}\int_\Omega\phi^*\nabla^2\phi\;dS=\frac{\hbar^2}{2m}\int_\Omega(\nabla\phi^*)\cdot(\nabla\phi)\;dS$$

となります. 1次元の場合と同じような形をしていますね. 従って, 以下の式(1.2)に示す弱形式が得られます.

$$\frac{\hbar^2}{2m}\int_\Omega(\nabla\phi^*)\cdot(\nabla\phi)\;dS+\int_\Omega\phi^*V\phi\;dS=E\int_\Omega\phi^*\phi\;dS\;\;\;(1.2)$$

2.形状関数

1次元のときと同じように形状関数を求めます. 今回は1つの要素を三角形1次要素（各頂点で近似する）とします. ちなみに三角形2次要素では頂点と各頂点の中点でも近似を行います.
本来であれば任意の三角形に対応できるようにするべきですが, 今回は簡単のために三角形の形状を固定します. $x$軸, $y$軸それぞれに $N+1$ 個の節点を設け, 領域を $N^2$ 個の要素に分割し, さらに各長方形について傾きが負の方の対角線で分割します. 従って, 領域は $2N^2$ 個の要素に分割されます. また, $x_0^i$, $x_1^i$, $y_0^j$, $y_1^j$ ($i,\,j=0,\,1,\,...\,,\,N-1$)によって作られる長方形のうち, $x_0^i$, $y_0^j$ を含む三角形要素を $(ij)_0$ 要素, もう一方を$(ij)_1$ 要素と呼ぶこととします.

領域の分割と形状関数

要素の形状が定まったので, ここからは近似関数を求めていきます. 1次元の場合は直線でしたが2次元の場合は平面になります. 平面の方程式は

$$ax+by+cz+d=0$$

で与えられます. まずは $(ij)_0$ 要素における形状関数を求めていきましょう. ここで得たい形状関数は

$$\phi=h^{(ij)_0}_{00}(x,\,y)\phi^{(ij)_0}_{00}+h^{(ij)_0}_{10}(x,\,y)\phi^{(ij)_0}_{10}+h^{(ij)_0}_{01}(x,\,y)\phi^{(ij)_0}_{01}\;\;\;(2.1)$$

を満たすような関数です. ただし, $\phi^{(ij)_0}_{kl}=\phi(x_k^i,\,y_l^j)$ です. 添え字が入り乱れてますが, 上付き添え字で要素の指定, 下付き添え字で要素の頂点を指定しています. 同様に, $h^{(ij)_0}_{kl}$ は形状関数の基底です. また, 要素が隣接していることから, $\phi$ に関して以下の式(2.2)に示すような関係があることに留意しましょう. あとで全体行列を構成するときに必要になります.

$$\phi^{(ij)_1}_{11}=\phi^{((i+1)j)_0}_{01}=\phi^{((i+1)j)_1}_{01}=\phi^{(i(j+1))_0}_{10}=\phi^{(i(j+1))_1}_{10}=\phi^{((i+1)(j+1))_0}_{00}\;\;\;(2.2)$$

式(2.1)を求める方法ですが, 今回は法線ベクトルを利用する方法を用います. 三点 $P_0$, $P_1$, $P_2$ を含む平面の法線ベクトル $\bm{n}$ は外積を用いて

$$\bm{n}=\overrightarrow{P_0P_1}\times\overrightarrow{P_0P_2}$$

で求められます. 法線ベクトルが平面に平行な任意のベクトルと直交している, つまり内積が0であることを利用すると, 平面の方程式

$$\bm{n}_x(x-x_0)+\bm{n}_y(y-y_0)+\bm{n}_z(z-z_0)=0$$

が得られます.

ここまでを実際に計算してみましょう. 法線ベクトルは

$$\begin{align*} \bm{n}&=\begin{pmatrix*} x^i_1-x^i_0\\ 0\\ \phi^{(ij)_0}_{10}-\phi^{(ij)_0}_{00} \end{pmatrix*}\times \begin{pmatrix*} 0\\ y^j_1-y^j_0\\ \phi^{(ij)_0}_{01}-\phi^{(ij)_0}_{00} \end{pmatrix*} \\ &=\begin{pmatrix*} -(y^j_1-y^j_0)(\phi^{(ij)_0}_{10}-\phi^{(ij)_0}_{00})\\ -(x^i_1-x^i_0)(\phi^{(ij)_0}_{01}-\phi^{(ij)_0}_{00})\\ (x^i_1-x^i_0)(y^j_1-y^j_0) \end{pmatrix*} \end{align*}$$

と計算できます. よって, 平面の方程式は

$$-(y^j_1-y^j_0)(\phi^{(ij)_0}_{10}-\phi^{(ij)_0}_{00})(x-x^i_0)-(x^i_1-x^i_0)(\phi^{(ij)_0}_{01}-\phi^{(ij)_0}_{00})(y-y^j_0)+(x^i_1-x^i_0)(y^j_1-y^j_0)(\phi^{(ij)_0}-\phi^{(ij)_0}_{00})=0$$

変形すると

$$\phi^{(ij)_0}=\left(1-\frac{x-x^i_0}{x^i_1-x^i_0}-\frac{y-y^j_0}{y^j_1-y^j_0}\right)\phi^{(ij)_0}_{00}+\frac{x-x^i_0}{x^i_1-x^i_0}\phi^{(ij)_0}_{10}+\frac{y-y^j_0}{y^j_1-y^j_0}\phi^{(ij)_0}_{01}$$

が得られます. さらに, この式を式(2.1)と見比べることで基底関数が得られます.

$$h^{(ij)_0}_{00}=\left(1-\frac{x-x^i_0}{x^i_1-x^i_0}-\frac{y-y^j_0}{y^j_1-y^j_0}\right)$$

$$h^{(ij)_0}_{10}=\frac{x-x^i_0}{x^i_1-x^i_0}$$

$$h^{(ij)_0}_{01}=\frac{y-y^j_0}{y^j_1-y^j_0}$$

$(ij)_1$ 要素の形状関数の計算は省略し, 以下に基底関数だけ示しておきます.

$$h^{(ij)_1}_{11}=\left(1+\frac{x-x^i_1}{x^i_1-x^i_0}+\frac{y-y^j_1}{y^j_1-y^j_0}\right)$$

$$h^{(ij)_1}_{10}=-\frac{y-y^j_1}{y^j_1-y^j_0}$$

$$h^{(ij)_1}_{01}=-\frac{x-x^i_1}{x^i_1-x^i_0}$$

3.離散化と要素行列

前節で求めた形状関数を使って式(1.2)の弱形式を離散化していきましょう. $(ij)_p$ 要素における $\phi(x,\,y)$ は以下の式(3.1)のように書けます.

$$\phi(x,\,y)=\sum_{kl\in m_p}\phi^{(ij)_p}_{kl}h^{(ij)_p}_{kl}(x,\,y)$$

ただし, $m_0=\{00,\,10,\,01\}$, $m_1=\{11,\,10,\,01\}$ です. ポテンシャル $V(x,\,y)$ に関しても同様の近似を行うこととします. また, $l^i_x=x^i_1-x^i_0$, $l^j_y=y^j_1-y^j_0$ とおくと, 式(1.2)左辺第1項の $\nabla\phi$ は

$$\nabla\phi^{(ij)_0}=\left( \frac{-\phi^{(ij)_0}_{00}+\phi^{(ij)_0}_{10}}{l_x^i}\;\; \frac{-\phi^{(ij)_0}_{00}+\phi^{(ij)_0}_{01}}{l_y^j} \right)^{\rm{T}}$$

$$\nabla\phi^{(ij)_1}=\left( \frac{\phi^{(ij)_0}_{11}-\phi^{(ij)_0}_{01}}{l_x^i}\;\; \frac{\phi^{(ij)_0}_{11}-\phi^{(ij)_0}_{10}}{l_y^j} \right)^{\rm{T}}$$

と表せます. よって, $(\nabla\phi^*)\cdot(\nabla\phi)$ は以下に示す行列 $L^{(ij)}$ を用いると

$$L^{(ij)}=\begin{pmatrix*} \displaystyle\frac{1}{(l_x^i)^2}+\frac{1}{(l_y^j)^2}&\displaystyle-\frac{1}{(l_x^i)^2}&\displaystyle-\frac{1}{(l_y^j)^2}\\ &&\\ \displaystyle-\frac{1}{(l_x^i)^2}&\displaystyle\frac{1}{(l_x^i)^2}&\displaystyle 0\\ &&\\ \displaystyle-\frac{1}{(l_y^j)^2}&\displaystyle 0&\displaystyle\frac{1}{(l_y^j)^2} \end{pmatrix*}$$

$$(\nabla\phi^{*(ij)_0})\cdot(\nabla\phi^{(ij)_0})=(\phi^{*(ij)_0}_{00}\;\;\phi^{*(ij)_0}_{10}\;\;\phi^{*(ij)_0}_{01})L^{(ij)}\begin{pmatrix*} \phi^{(ij)_0}_{00}\\ \phi^{(ij)_0}_{10}\\ \phi^{(ij)_0}_{01} \end{pmatrix*}$$

$$(\nabla\phi^{*(ij)_1})\cdot(\nabla\phi^{(ij)_1})=(\phi^{*(ij)_1}_{11}\;\;\phi^{*(ij)_1}_{01}\;\;\phi^{*(ij)_1}_{10})L^{(ij)}\begin{pmatrix*} \phi^{(ij)_1}_{11}\\ \phi^{(ij)_1}_{01}\\ \phi^{(ij)_1}_{10} \end{pmatrix*}$$

と表せます（$\phi$ の下付き添え字の順番に注意!）. これは $x$, $y$ に依存しないので, 式(1.2)左辺第1項は

$$\frac{\hbar^2}{2m}\sum_{(i,\,j)\in M}S^{(ij)}\left((\phi^{*(ij)_0}_{00}\;\;\phi^{*(ij)_0}_{10}\;\;\phi^{*(ij)_0}_{01})L^{(ij)}\begin{pmatrix*} \phi^{(ij)_0}_{00}\\ \phi^{(ij)_0}_{10}\\ \phi^{(ij)_0}_{01} \end{pmatrix*}+ (\phi^{*(ij)_1}_{11}\;\;\phi^{*(ij)_1}_{01}\;\;\phi^{*(ij)_1}_{10})L^{(ij)}\begin{pmatrix*} \phi^{(ij)_1}_{11}\\ \phi^{(ij)_1}_{01}\\ \phi^{(ij)_1}_{10} \end{pmatrix*} \right)$$

と表せます. ただし, $M=\{(i,\,j)\in\mathbb{Z}^2|0\le i,\,j\le N-1\}$, $S^{(ij)}=l_x^il_y^j/2$ （三角要素の面積）です.

次に式(1.2)左辺第2項を計算していきます. まずは変数を変換して簡単な形に変えましょう. $(ij)_0$ 要素の積分は

$$\begin{align*} \int_{\Omega^{(ij)_0}}f(x,\,y)\;dS&=\int_{x^i_0}^{x^i_1}\int_{y^j_0}^{-\frac{l^j_y}{l^i_x}(x-x^i_1)+y_0}f(x,\,y)\;dxdy \end{align*}$$

と表せますが, $s=\frac{x-x^i_0}{l^i_x}$, $t=\frac{y-y^j_0}{l^j_y}$ とおくと

$$\begin{align*} \int_{\Omega^{(ij)_0}}f(x,\,y)\;dS&=l^i_xl^j_y\int_0^1\int_0^{1-s}f(s,\,t)\;dsdt \end{align*}$$

となり, 計算が比較的容易になります. ポテンシャル $V$ も基底 $h^{(ij)_0}_{kl}(x,\,y)$ を用いて近似するので, 被積分関数には $(h_{00})^3=(1-s-t)^3$, $(h_{00})^2h_{10}=(1-s-t)^2s$, ... , $(h_{01})^3=t^3$ が現れます（3つの基底から重複を許して3つ選ぶ組合せ, つまり10種ある）. 頑張ってこれらをすべて積分すると（対称性も意識しながらやると楽になります） $(ij)_0$ 要素の積分は

$$(\phi^{*(ij)_0}_{00}\;\;\phi^{*(ij)_0}_{10}\;\;\phi^{*(ij)_0}_{01}) V^{(ij)_0} \begin{pmatrix*} \phi^{(ij)_0}_{00}\\ \phi^{(ij)_0}_{10}\\ \phi^{(ij)_0}_{01} \end{pmatrix*}$$

となります. ただし

$$V^{(ij)_0}=\frac{l_x^il_y^j}{120}\begin{pmatrix*} 6V_{00}+2V_{10}+2V_{01} & 2V_{00}+2V_{10}+V_{01} & 2V_{00}+V_{10}+2V_{01}\\ 2V_{00}+2V_{10}+V_{01} & 2V_{00}+6V_{10}+2V_{01} & V_{00}+2V_{10}+2V_{01}\\ 2V_{00}+V_{10}+2V_{01} &V_{00}+2V_{10}+2V_{01} & 2V_{00}+2V_{10}+6V_{01} \end{pmatrix*}$$

です（上付き添え字省略）. $(ij)_1$ 要素については計算せずとも, 対称性により, $00$を$11$に, $01$と$10$を入れ替えるだけでよいことが分かります.

$$V^{(ij)_1}=\frac{l_x^il_y^j}{120}\begin{pmatrix*} 6V_{11}+2V_{01}+2V_{10} & 2V_{11}+2V_{01}+V_{10} & 2V_{11}+V_{01}+2V_{10}\\ 2V_{11}+2V_{01}+V_{10} & 2V_{11}+6V_{01}+2V_{10} & V_{11}+2V_{01}+2V_{10}\\ 2V_{11}+V_{01}+2V_{10} &V_{11}+2V_{01}+2V_{10} & 2V_{11}+2V_{01}+6V_{10} \end{pmatrix*}$$

$$(\phi^{*(ij)_1}_{11}\;\;\phi^{*(ij)_1}_{01}\;\;\phi^{*(ij)_1}_{10}) V^{(ij)_1} \begin{pmatrix*} \phi^{(ij)_1}_{11}\\ \phi^{(ij)_1}_{01}\\ \phi^{(ij)_1}_{10} \end{pmatrix*}$$

以上より, 式(1.2)左辺第2項は

$$\sum_{(i,\,j)\in M}\left((\phi^{*(ij)_0}_{00}\;\;\phi^{*(ij)_0}_{10}\;\;\phi^{*(ij)_0}_{01}) V^{(ij)_0} \begin{pmatrix*} \phi^{(ij)_0}_{00}\\ \phi^{(ij)_0}_{10}\\ \phi^{(ij)_0}_{01} \end{pmatrix*}+ (\phi^{*(ij)_1}_{11}\;\;\phi^{*(ij)_1}_{01}\;\;\phi^{*(ij)_1}_{10}) V^{(ij)_1} \begin{pmatrix*} \phi^{(ij)_1}_{11}\\ \phi^{(ij)_1}_{01}\\ \phi^{(ij)_1}_{10} \end{pmatrix*}\right)$$

と書けます.

最後に式(1.2)右辺を計算しましょう. とはいえ, これは左辺第2項と計算の流れが同じなので省略し, 結果だけを以下に示しておきます.

$$B^{(ij)}=\frac{l_x^il_y^j}{24} \begin{pmatrix*} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{pmatrix*}$$

とおくと, $(ij)_0$ 要素は

$$(\phi^{*(ij)_0}_{00}\;\;\phi^{*(ij)_0}_{10}\;\;\phi^{*(ij)_0}_{01}) B^{(ij)} \begin{pmatrix*} \phi^{(ij)_0}_{00}\\ \phi^{(ij)_0}_{10}\\ \phi^{(ij)_0}_{01} \end{pmatrix*}$$

$(ij)_1$ 要素は

$$(\phi^{*(ij)_1}_{11}\;\;\phi^{*(ij)_1}_{01}\;\;\phi^{*(ij)_1}_{10}) B^{(ij)} \begin{pmatrix*} \phi^{(ij)_1}_{11}\\ \phi^{(ij)_1}_{01}\\ \phi^{(ij)_1}_{10} \end{pmatrix*}$$

となり, 式(1.2)右辺は

$$E\sum_{(i,\,j)\in M}\left( (\phi^{*(ij)_0}_{00}\;\;\phi^{*(ij)_0}_{10}\;\;\phi^{*(ij)_0}_{01}) B^{(ij)} \begin{pmatrix*} \phi^{(ij)_0}_{00}\\ \phi^{(ij)_0}_{10}\\ \phi^{(ij)_0}_{01} \end{pmatrix*}+ (\phi^{*(ij)_1}_{11}\;\;\phi^{*(ij)_1}_{01}\;\;\phi^{*(ij)_1}_{10}) B^{(ij)} \begin{pmatrix*} \phi^{(ij)_1}_{11}\\ \phi^{(ij)_1}_{01}\\ \phi^{(ij)_1}_{10} \end{pmatrix*} \right)$$

以上より

$$A^{(ij)_k}=\frac{\hbar^2}{2m}S^{(ij)}L^{(ij)}+V^{(ij)_k}$$

$$\bm{\phi}^{(ij)_0}=\left(\phi^{(ij)_0}_{00}\;\;\phi^{(ij)_0}_{10}\;\;\phi^{(ij)_0}_{01}\right)^{\rm{T}}$$

$$\bm{\phi}^{(ij)_1}=\left(\phi^{(ij)_1}_{11}\;\;\phi^{(ij)_1}_{01}\;\;\phi^{(ij)_1}_{10}\right)^{\rm{T}}$$

とおくと, 式(1.2)は

$$\sum_{(i,\,j)\in M}\sum_{k=0}^1(\bm{\phi}^{*(ij)_k})^{\rm{T}}A^{(ij)_k}\bm{\phi}^{(ij)_k}=E\sum_{(i,\,j)\in M}\sum_{k=0}^1(\bm{\phi}^{*(ij)_k})^{\rm{T}}B^{(ij)}\bm{\phi}^{(ij)_k}$$

と書けます.

4.全体行列を構成

全体行列を構成します. 最終的に得られる方程式は式(2.2)を考慮すると

$$\left(\phi^{*00}\;\;\phi^{*10}\;\;...\;\;\phi^{*NN}\right) A \begin{pmatrix*} \phi^{00}\\ \phi^{10}\\ \vdots\\ \phi^{NN} \end{pmatrix*} =E \left(\phi^{*00}\;\;\phi^{*10}\;\;...\;\;\phi^{*NN}\right) B \begin{pmatrix*} \phi^{00}\\ \phi^{10}\\ \vdots\\ \phi^{NN} \end{pmatrix*}$$

という形をしています. これに注意しながら要素行列の各要素を全体行列に対応させましょう.

さらに, 得られた方程式に変分原理を用いることにより

$$A\bm{\phi}=EB\bm{\phi}$$

を得ます.

5.境界条件処理

Dirichlet条件 $\forall(x,\;y)\in\partial\Omega,\;\phi(x,\;y)=0$ を用います. 1次元のときは0行0列とN行N列を削除していましたが, 今回は境界にあたる行, および列をすべて削除しましょう.

プログラミング編

6.Pythonで2次元有限要素法

import numpy as np
import matplotlib.pyplot as plt
from collections.abc import Callable
import scipy.sparse.linalg as spsl

class SchrodingerFEM2():

    DIRAC_CONSTANT: float = 1.054571817e-34
    DIRAC_CONSTANT_DOUBLE: float = DIRAC_CONSTANT**2

    EIG_VAL_FILE = 'eigenvalue.csv'
    EIG_VEC_FILE = 'eigenvector.csv'

    n: int  #要素数

    m: float    #質量

    matrix_genned: bool #行列は生成済みか否か

    x: np.ndarray   #節点のx座標
    y: np.ndarray   #節点のy座標
    lx: np.ndarray  #各要素のx軸に平行な辺の長さ
    ly: np.ndarray  #各要素のy軸に平行な辺の長さ
    A: np.ndarray   #A行列
    B: np.ndarray   #B行列

    eigen_val_vec: list #計算結果を保持する

    V: Callable[[float, float], float]  #ポテンシャル



    def __init__(self, n: int, m: float, x: np.ndarray, y: np.ndarray, V: Callable[[float, float], float]):
        self.n = n  
        self.m = m  
        self.x = x  
        self.y = y  
        self.V = V  

        self.matrix_genned = False  
        self.eigen_val_vec = []     

        if len(x) != n+1:
            raise ValueError(f'Array x size should be {n+1} but is {len(x)}.')
        if len(y) != n+1:
            raise ValueError(f'Array x size should be {n+1} but is {len(y)}.')
        
        self.lx = np.zeros(n)
        self.ly = np.zeros(n)
        for i in range(n):
            self.lx[i] = self.x[i+1]-self.x[i]
            self.ly[i] = self.y[i+1]-self.y[i]
        
        mat_size = (n+1)*(n+1)
        self.A = np.zeros((mat_size, mat_size), np.float64)
        self.B = np.zeros((mat_size, mat_size), np.float64)
    
    # ポテンシャル項の計算用の3×3行列が3つ
    V_base = np.array([[[6, 2, 2], [2, 2, 1], [2, 1, 2]], [[2, 2, 1], [2, 6, 2], [1, 2, 2]], [[2, 1, 2], [1, 2, 2], [2, 2, 6]]], np.float64)

    def gen_V_element(self, i: int, j: int)->tuple[np.ndarray, np.ndarray]:
        Vij0 = self.V(self.x[i], self.y[j])*self.V_base[0]+self.V(self.x[i+1], self.y[j])*self.V_base[1]+self.V(self.x[i], self.y[j+1])*self.V_base[2]
        Vij1 = self.V(self.x[i+1], self.y[j+1])*self.V_base[0]+self.V(self.x[i], self.y[j+1])*self.V_base[1]+self.V(self.x[i+1], self.y[j])*self.V_base[2]

        Vij0 *= (self.lx[i]*self.ly[j])/120
        Vij1 *= (self.lx[i]*self.ly[j])/120

        return Vij0, Vij1
    
    def gen_L_matrix(self, i: int, j: int)->np.ndarray:
        tmpx = 1/self.lx[i]**2
        tmpy = 1/self.ly[j]**2
        return np.array([[tmpx+tmpy, -tmpx, -tmpy], [-tmpx, tmpx, 0], [-tmpy, 0, tmpy]], np.float64)

    def gen_A_element(self, i: int, j: int)->np.ndarray:
        Aij0 = np.zeros((3, 3), np.float64)
        Aij1 = np.zeros((3, 3), np.float64)
        
        Vij0, Vij1 = self.gen_V_element(i, j)

        Sij = self.lx[i]*self.ly[j]/2
        tmpij = (self.DIRAC_CONSTANT_DOUBLE*Sij/(2*self.m))*self.gen_L_matrix(i, j)

        Aij0 += (tmpij+Vij0)
        Aij1 += (tmpij+Vij1)

        Aij = np.zeros((4, 4), np.float64)

        #2つの三角要素を事前に合体
        for k in range(3):
            for l in range(3):
                Aij[k][l] += Aij0[k][l]
                Aij[k+1][l+1] += Aij1[2-k][2-l]
        return Aij
    
    def gen_B_element(self, i: int, j: int)->np.ndarray:
        tmpBij = np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]])

        Bij = np.zeros((4, 4), np.float64)

        #2つの三角要素を事前に合体
        for k in range(3):
            for l in range(3):
                Bij[k][l] += tmpBij[k][l]
                Bij[k+1][l+1] += Bij[2-k][2-l]
        
        Bij *= (self.lx[i]*self.ly[j])/24

        return Bij
    
    #全体行列への対応のための配列
    idx = [(0, 0), (1, 0), (0, 1), (1, 1)]

    def gen_matrix(self)->None:
        for i in range(self.n):
            for j in range(self.n):
                Aij = self.gen_A_element(i, j)
                Bij = self.gen_B_element(i, j)
                for k in range(4):
                    for l in range(4):
                        self.A[i+self.idx[k][0]+(j+self.idx[k][1])*(self.n+1)][i+self.idx[l][0]+(j+self.idx[l][1])*(self.n+1)] += Aij[k][l]
                        self.B[i+self.idx[k][0]+(j+self.idx[k][1])*(self.n+1)][i+self.idx[l][0]+(j+self.idx[l][1])*(self.n+1)] += Bij[k][l]

        #境界条件処理
        self.A = self.A[self.n+1:-(self.n+1),self.n+1:-(self.n+1)]
        self.B = self.B[self.n+1:-(self.n+1),self.n+1:-(self.n+1)]
        for i in range(self.n-1):
            self.A = np.delete(self.A, i*(self.n-1), 0)
            self.A = np.delete(self.A, i*(self.n-1), 1)
            self.A = np.delete(self.A, (i+1)*(self.n-1), 0)
            self.A = np.delete(self.A, (i+1)*(self.n-1), 1)
            self.B = np.delete(self.B, i*(self.n-1), 0)
            self.B = np.delete(self.B, i*(self.n-1), 1)
            self.B = np.delete(self.B, (i+1)*(self.n-1), 0)
            self.B = np.delete(self.B, (i+1)*(self.n-1), 1)

        print('Matrix calculation has been done !\n\n', end='')
        self.matrix_genned = True
    
    def solve(self)->None:
        if not self.matrix_genned:
            raise ValueError('Matrix has not been generated.')
        
        val, vec = spsl.eigs(self.A, self.n**2, self.B, maxiter=self.n**2*500, which='SI', sigma=0.05)
        print('Eigen value has been calculated !\n\n', end='')

        for i in range(len(val)):
            self.eigen_val_vec.append((val[i], vec[:, i]))

        self.eigen_val_vec.sort(key=lambda x:x[0])
    
    #1次元の解の配列を2次元に変換
    def reshape_vec(self, i: int)->np.ndarray:
        vec = np.zeros((self.n+1, self.n+1))

        for k in range(self.n-1):
            for l in range(self.n-1):
                vec[l+1][k+1] = self.eigen_val_vec[i][1][l+k*(self.n-1)]
        
        return vec
    
    def get_matrix(self)->tuple[np.ndarray, np.ndarray]:
        return self.A, self.B
    
    def get_eigenvalue(self, i)->np.complex64:
        return self.eigen_val_vec[i][0]
    
    def get_eigenvector(self, i)->np.ndarray:
        return self.reshape_vec(i)
    
    def export_data(self)->None:
        val_file_handler = open(self.EIG_VAL_FILE, mode='w')
        vec_file_handler = open(self.EIG_VEC_FILE, mode='w')

        for val_vec in self.eigen_val_vec:
            val_file_handler.write(f'{val_vec[0]}\n')
            vec_file_handler.write(','.join(map(str, val_vec[1]))+'\n')
        
        val_file_handler.close()
        vec_file_handler.close()

# 無限井戸型ポテンシャル（とみなせる）
def V(x: float, y: float)->float:
    return 0

if __name__=='__main__':
    N = 25  #要素数
    m = SchrodingerFEM2.DIRAC_CONSTANT_DOUBLE/2

    L = 5
    x = np.linspace(-L, L, N+1)
    y = np.linspace(-L, L, N+1)

    fem = SchrodingerFEM2(N, m, x, y, V)
    fem.gen_matrix()
    fem.solve()

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    X, Y = np.meshgrid(x, y)
    
    print(fem.get_eigenvalue(1))
    phi = fem.get_eigenvector(1)
    ax.plot_surface(X, Y, phi, cmap='jet')

    plt.show()

7.例

無限井戸型ポテンシャルの2番目にエネルギーが小さい状態