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

はじめに

有限要素法による非時間依存1次元シュレディンガー方程式で1次元, 有限要素法による非時間依存2次元シュレディンガー方程式で2次元のシュレディンガー方程式を有限要素法を用いて解いたので, 今回は3次元のシュレディンガー方程式を解く.

3次元のシュレディンガー方程式は

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

です（$dV$ は体積分を表す.）. ただし, $\displaystyle\nabla^2=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}$ です. また, 境界条件としてDirichlet境界条件 $\forall(x,\,y,\,z)\in \partial\Omega,\,\phi(x,\,y,\,z)=0$ を用います.

計算編

1.弱形式

2次元の場合と同じですが一応.
冒頭のシュレディンガー方程式の両辺に重み関数として $\phi$ の複素共役 $\phi^*$ をかけ, 両辺を領域 $\Omega$ で積分すると以下の式(1.1)を得ます.

$$-\frac{\hbar^2}{2m}\int_\Omega\phi^*\nabla^2\phi\;dV+\int_\Omega\phi^*V\phi\;dV=E\int_\Omega\phi^*\phi\;dV\;\;\;(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\;dV=-\frac{\hbar^2}{2m}\left(-\int_\Omega(\nabla\phi^*)\cdot(\nabla\phi)\;dV+\int_\Omega\nabla\cdot(\phi^*\nabla\phi)\;dV\right)$$

となります. さらに, ガウスの発散定理

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

を用いると

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

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

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

を得ます. 以上より弱形式は以下の式(1.2)のようになります.

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

2.形状関数

2次元のときはメッシュ形状を具体的に定めましたが, 一般的に書いた方が記述が見やすいと思ったので今回は一般的なメッシュを考えます. 用いる要素は四面体1次要素とし, 領域を 0から $N-1$ までの $N$ 個に分割するものとします.

$e$ 番目（$e=0,\,1,\,...\,,\,N-1$）の要素について考えましょう. この四面体要素の4頂点の座標を $(x^e_i,\,y^e_i,\,z^e_i)$, $\phi$ の値を $\phi^e_i$（$i=0,\,1,\,2,\,3$）としましょう. 形状関数として線形な関数

$$\phi=\alpha x+\beta y+\gamma z+\delta$$

を用いるとすると $\alpha$, $\beta$, $\gamma$, $\delta$ は以下の連立方程式を満たします.

$$\begin{pmatrix*} x_0^e&y_0^e&z_0^e&1\\ x_1^e&y_1^e&z_1^e&1\\ x_2^e&y_2^e&z_2^e&1\\ x_3^e&y_3^e&z_3^e&1 \end{pmatrix*} \begin{pmatrix*} \alpha^e\\ \beta^e\\ \gamma^e\\ \delta^e \end{pmatrix*} = \begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

ここで, 行列 $S^e$ を

$$S^e= \begin{pmatrix*} x_0^e&y_0^e&z_0^e&1\\ x_1^e&y_1^e&z_1^e&1\\ x_2^e&y_2^e&z_2^e&1\\ x_3^e&y_3^e&z_3^e&1 \end{pmatrix*}$$

とし, さらに行列 $M^e$ を

$$\begin{align*} M^e&= \begin{pmatrix*} x_0^e&y_0^e&z_0^e&1\\ x_1^e&y_1^e&z_1^e&1\\ x_2^e&y_2^e&z_2^e&1\\ x_3^e&y_3^e&z_3^e&1 \end{pmatrix*}^{-1}\\ &\\ &={S^e}^{-1} \end{align*}$$

とします. よって

$$\begin{pmatrix*} \alpha^e\\ \beta^e\\ \gamma^e\\ \delta^e \end{pmatrix*} =M^e \begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

となるので, 行列 $M$ の $i$ 行 $j$ 列要素を $M_{ij}$ （0-indexed）で表すと

$$\phi^e=\sum_{j=0}^3\left(M_{0j}^ex+M_{1j}^ey+M_{2j}^ez+M_{3j}^e\right)\phi_j^e$$

となります. これが1次四面体要素の形状関数です. また, $M_{0j}^ex+M_{1j}^ey+M_{2j}^ez+M_{3j}^e$ が形状関数の基底だと分かります. ここで, 基底を $L_j^e=M_{0j}^ex+M_{1j}^ey+M_{2j}^ez+M_{3j}^e$ としておきましょう.

3.離散化と要素行列

式(1.2)の弱形式は

$$\frac{\hbar^2}{2m}\sum_{e=0}^{N-1}\int_{\Omega^e}(\nabla{\phi^*}^e)\cdot(\nabla\phi^e)\;dV+\sum_{e=0}^{N-1}\int_{\Omega^e}{\phi^*}^eV^e\phi^e\;dV =E\sum_{e=0}^{N-1}\int_{\Omega^e}{\phi^*}^e\phi^e\;dV\;\;\;(3.1)$$

と離散化できます.

式(3.1)左辺第1項から計算していきましょう. $\nabla\phi$ は

$$\begin{align*} \nabla\phi&= \begin{pmatrix*} \alpha\\ \beta\\ \gamma \end{pmatrix*} \\ &= \begin{pmatrix*} M_{00}&M_{01}&M_{02}&M_{03}\\ M_{10}&M_{11}&M_{12}&M_{13}\\ M_{20}&M_{21}&M_{22}&M_{23}\\ \end{pmatrix*} \begin{pmatrix*} \phi_0\\ \phi_1\\ \phi_2\\ \phi_3 \end{pmatrix*} \end{align*}$$

となります （上付き添え字略）. よって

$$\begin{align*} (\nabla\phi^*)\cdot(\nabla\phi)&=(\nabla\phi^*)^{\rm{T}}(\nabla\phi)\\ &=\left(\phi_0^*\;\;\phi_1^*\;\;\phi_2^*\;\;\phi_3^*\right) \begin{pmatrix*} M_{00}&M_{10}&M_{20}\\ M_{01}&M_{11}&M_{21}\\ M_{02}&M_{12}&M_{22}\\ M_{03}&M_{13}&M_{23} \end{pmatrix*} \begin{pmatrix*} M_{00}&M_{01}&M_{02}&M_{03}\\ M_{10}&M_{11}&M_{12}&M_{13}\\ M_{20}&M_{21}&M_{22}&M_{23}\\ \end{pmatrix*} \begin{pmatrix*} \phi_0\\ \phi_1\\ \phi_2\\ \phi_3 \end{pmatrix*}\\ &\\ &= \left(\phi_0^*\;\;\phi_1^*\;\;\phi_2^*\;\;\phi_3^*\right) \begin{pmatrix*} M_{00}^2+M_{10}^2+M_{20}^2&M_{00}M_{01}+M_{10}M_{11}+M_{20}M_{21}&M_{00}M_{02}+M_{10}M_{12}+M_{20}M_{22}&M_{00}M_{03}+M_{10}M_{13}+M_{20}M_{23}\\ *&M_{01}^2+M_{11}^2+M_{21}^2&M_{01}M_{02}+M_{11}M_{12}+M_{21}M_{22}&M_{01}M_{03}+M_{11}M_{13}+M_{21}M_{23}\\ *&*&M_{02}^2+M_{12}^2+M_{22}^2&M_{02}M_{03}+M_{12}M_{13}+M_{22}M_{23}\\ *&*&*&M_{03}^2+M_{13}^2+M_{23}^2 \end{pmatrix*} \begin{pmatrix*} \phi_0\\ \phi_1\\ \phi_2\\ \phi_3 \end{pmatrix*} \end{align*}$$

となります（対称行列になるので下三角要素を省略）.

$$K =\frac{\hbar^2}{2m} \begin{pmatrix*} M_{00}^2+M_{10}^2+M_{20}^2&M_{00}M_{01}+M_{10}M_{11}+M_{20}M_{21}&M_{00}M_{02}+M_{10}M_{12}+M_{20}M_{22}&M_{00}M_{03}+M_{10}M_{13}+M_{20}M_{23}\\ *&M_{01}^2+M_{11}^2+M_{21}^2&M_{01}M_{02}+M_{11}M_{12}+M_{21}M_{22}&M_{01}M_{03}+M_{11}M_{13}+M_{21}M_{23}\\ *&*&M_{02}^2+M_{12}^2+M_{22}^2&M_{02}M_{03}+M_{12}M_{13}+M_{22}M_{23}\\ *&*&*&M_{03}^2+M_{13}^2+M_{23}^2 \end{pmatrix*}$$

とおくと, これは $x$, $y$, $z$ に依存しないので, 式(3.1)左辺第1項は

$$\frac{\hbar^2}{2m}\sum_{e=0}^{N-1}\int_{\Omega^e}(\nabla{\phi^*}^e)\cdot(\nabla\phi^e)\;dV = \sum_{e=0}^{N-1} \left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right) T^eK^e \begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

ただし, $T^e$ は $e$ 番目の四面体の体積であり

$$T^e = \frac{1}{6}|\det S^e|$$

で計算できます.

体積の計算について

4点 $P_0\;(x_0,\,y_0,\,z_0)$, $P_1\;(x_1,\,y_1,\,z_1)$, $P_2\;(x_2,\,y_2,\,z_2)$, $P_3\;(x_3,\,y_3,\,z_3)$ がなす四面体を考えます. 4点を平行移動しても移動後の4点がなす四面体の体積は変化しないので, $P_0$ が原点に来るように平行移動すると

$$\det S= \det \begin{pmatrix*} 0&0&0&1\\ x_1-x_0&y_1-y_0&z_1-z_0&1\\ x_2-x_0&y_2-y_0&z_2-z_0&1\\ x_3-x_0&y_3-y_0&z_3-z_0&1 \end{pmatrix*}$$

となります. これを1行目に関して余因子展開すると

$$\det S=-\det \begin{pmatrix*} x_1-x_0&y_1-y_0&z_1-z_0\\ x_2-x_0&y_2-y_0&z_2-z_0\\ x_3-x_0&y_3-y_0&z_3-z_0 \end{pmatrix*}$$

（絶対値を計算するので符号を無視しつつ）計算を進めると

$$\begin{align*} \det S&=(x_1-x_0)(y_2-y_0)(z_3-z_0)+(y_1-y_0)(z_2-z_0)(x_3-x_0)+(z_1-z_0)(x_2-x_0)(y_3-y_0)-(z_1-z_0)(y_2-y_0)(x_3-x_0)-(y_1-y_0)(x_2-x_0)(z_3-z_0)-(x_1-x_0)(z_2-z_0)(y_3-y_0)\\ &\\ &= \begin{pmatrix*} (x_2-x_0)(y_3-y_0)-(y_2-y_0)(x_3-x_0)\\ (y_1-y_0)(x_3-x_0)-(x_1-x_0)(y_3-y_0)\\ (x_1-x_0)(y_2-y_0)-(y_1-y_0)(x_2-x_0) \end{pmatrix*} \cdot \begin{pmatrix*} z_1-z_0\\ z_2-z_0\\ z_3-z_0 \end{pmatrix*}\\ &\\ &= \left( \begin{pmatrix*} x_1-x_0\\ x_2-x_0\\ x_3-x_0 \end{pmatrix*} \times \begin{pmatrix*} y_1-y_0\\ y_2-y_0\\ y_3-y_0 \end{pmatrix*} \right) \cdot \begin{pmatrix*} z_1-z_0\\ z_2-z_0\\ z_3-z_0 \end{pmatrix*} \end{align*}$$

となります. これはいわゆるスカラー三重積であり, 絶対値は3つのベクトル

$$\begin{pmatrix*} x_1-x_0\\ x_2-x_0\\ x_3-x_0 \end{pmatrix*}\;\;\; \begin{pmatrix*} y_1-y_0\\ y_2-y_0\\ y_3-y_0 \end{pmatrix*}\;\;\; \begin{pmatrix*} z_1-z_0\\ z_2-z_0\\ z_3-z_0 \end{pmatrix*}$$

が張る平行六面体の体積になっています （外積の部分で底面積を計算し, 外積のベクトルに対して内積をとる, つまり正射影することで高さを計算していると考えられます）. よって, これを6で割った値は四面体の体積になります.

次に式(3.1)左辺第2項を計算していきます. $\phi^{*e}$, および $\phi^e$ を基底 $L_j$ との線形結合で表すと

$$\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right) \int_{\Omega^e} \begin{pmatrix*} L_0^e\\L_1^e\\L_2^e\\L_3^e \end{pmatrix*} V^e \left(L_0^e\;\;L_1^e\;\;L_2^e\;\;L_3^e\right) \;dV \begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

となります. 例によってポテンシャルも基底 $L_j$ を用いて近似すると, 4つの基底のうちから3つを選んだ積を領域内で積分する必要があることが分かります. この積分は

$$\int_{\Omega^e}L_0^kL_1^lL_2^mL_3^n\;dV=6T\frac{k!l!m!n!}{(k+l+m+n+3)!}$$

で計算できます.

四面体1次要素の積分の導出

行列 $S$ の $ij$ 余因子を

$$\Delta_{ij}=(-1)^{i+j} \det \begin{pmatrix*} \ddots&\vdots&\vdots&\\ \cdots&S_{(i-1)(j-1)}&S_{(i-1)(j+1)}&\cdots\\ \cdots&S_{(i+1)(j-1)}&S_{(i+1)(j+1)}&\cdots\\ &\vdots&\vdots&\ddots \end{pmatrix*}$$

と表します（行列 $S$ から $i$ 行と $j$ 行を除いた行列の行列式に $(-1)^{i+j}$ をかけた値）. このとき

$$\begin{align*} M&=S^{-1}\\ &= \frac{1}{\det S} \begin{pmatrix*} \Delta_{00}&\Delta_{10}&\Delta_{20}&\Delta_{30}\\ \Delta_{01}&\Delta_{11}&\Delta_{21}&\Delta_{31}\\ \Delta_{02}&\Delta_{12}&\Delta_{22}&\Delta_{32}\\ \Delta_{03}&\Delta_{13}&\Delta_{23}&\Delta_{33}\\ \end{pmatrix*} \end{align*}$$

となります（行と列の添え字に注意）.

基底をここに改めて示しておきます.

$$L_j^e=M_{0j}^ex+M_{1j}^ey+M_{2j}^ez+M_{3j}^e$$

計算したい積分は

$$\int_{\Omega^e}L_0^kL_1^lL_2^mL_3^n\;dV$$

です. 頑張って計算すると $L_0+L_1+L_2+L_3=1$ になるので, $L_0=X$, $L_1=Y$, $L_2=Z$, $L_3=1-X-Y-Z$ とおきます. 変数を $x$, $y$, $z$ から $X$, $Y$, $Z$ に変換したいので頑張ってヤコビアンを計算すると

$$\begin{align*} \left|\det \begin{pmatrix*} \displaystyle\frac{\partial x}{\partial X}&\displaystyle\frac{\partial y}{\partial X}&\displaystyle\frac{\partial z}{\partial X}\\ &&\\ \displaystyle\frac{\partial x}{\partial Y}&\displaystyle\frac{\partial y}{\partial Y}&\displaystyle\frac{\partial z}{\partial Y}\\ &&\\ \displaystyle\frac{\partial x}{\partial Z}&\displaystyle\frac{\partial y}{\partial Z}&\displaystyle\frac{\partial z}{\partial Z} \end{pmatrix*} \right| &= \left|(\det S)^3\det \begin{pmatrix*} M_{11}M_{22}-M_{12}M_{21}&M_{21}M_{02}-M_{22}M_{01}&M_{01}M_{12}-M_{02}M_{11}\\ M_{12}M_{20}-M_{10}M_{22}&M_{22}M_{00}-M_{20}M_{02}&M_{02}M_{10}-M_{00}M_{12}\\ M_{10}M_{21}-M_{11}M_{20}&M_{20}M_{01}-M_{21}M_{00}&M_{00}M_{11}-M_{01}M_{10} \end{pmatrix*} \right|\\ &= \left|\det \begin{pmatrix*} \Delta_{11}\Delta_{22}-\Delta_{21}\Delta_{12}&\Delta_{12}\Delta_{20}-\Delta_{22}\Delta_{10}&\Delta_{10}\Delta_{21}-\Delta_{20}\Delta_{11}\\ \Delta_{21}\Delta_{02}-\Delta_{01}\Delta_{22}&\Delta_{22}\Delta_{00}-\Delta_{02}\Delta_{20}&\Delta_{20}\Delta_{01}-\Delta_{00}\Delta_{21}\\ \Delta_{01}\Delta_{12}-\Delta_{11}\Delta_{02}&\Delta_{02}\Delta_{10}-\Delta_{12}\Delta_{00}&\Delta_{00}\Delta_{11}-\Delta_{10}\Delta_{01} \end{pmatrix*} \right|\\ &= \left| \det S \right|\\ &=6T \end{align*}$$

となります（MATLABで計算しました）. よって

$$\int_{\Omega}L_0^kL_1^lL_2^mL_3^n\;dV=6T\int_{\Omega}X^kY^lZ^m(1-X-Y-Z)^n\;dXdYdZ$$

となります. 次に積分する範囲を考えましょう. $X=L_0$ は頂点 $(x_0,\,y_0,\,z_0)$ で1となり, 他の頂点では0であるような関数です. $Y=L_1$, $Z=L_2$ に関しても同様です. よって, 頂点 $(x_0,\,y_0,\,z_0)$ は $XYZ$ 空間の $X$ 軸上の $X=1$ の点に, $(x_1,\,y_1,\,z_1)$ は$Y$ 軸上の $Y=1$ の点に, $(x_2,\,y_2,\,z_2)$ は$Z$ 軸上の $Z=1$ の点に, そして頂点 $(x_3,\,y_3,\,z_3)$ は原点に写されたと考えられます. 従って

$$\int_{\Omega}L_0^kL_1^lL_2^mL_3^n\;dV=6T\int_0^1\int_0^{1-X}\int_0^{1-X-Y}X^kY^lZ^m(1-X-Y-Z)^n\;dXdYdZ$$

となります. さらに $\displaystyle\tilde{Z}=\frac{Z}{1-X-Y}$ と変数変換すると

$$\int_0^{1-X-Y}Z^m(1-X-Y-Z)^n\;dZ=(1-X-Y)^{m+n+1}\int_0^1\tilde{Z}^m(1-\tilde{Z})^n\;d\tilde{Z}$$

となります. $\tilde{Z}$ の積分は以下に示す第1種オイラー積分（その中でも特に $\alpha=0$, $\beta=1$ なのでベータ積分と呼ばれます. 調べたら証明は色々出てきます）です.

$$\int_\alpha^\beta(x-\alpha)^m(\beta-x)^n\;dx=\frac{m!n!}{(m+n+1)!}(\beta-\alpha)^{m+n+1}$$

これを用いると

$$\int_0^{1-X-Y}Z^m(1-X-Y-Z)^n\;dZ=(1-X-Y)^{m+n+1}\frac{m!n!}{(m+n+1)!}$$

と計算できます. 同じように $\displaystyle\tilde{Y}=\frac{Y}{1-X}$ と変数変換すると

$$\begin{align*} \int_0^{1-X}Y^l(1-X-Y)^{m+n}\;dY&=(1-X)^{l+m+n+2}\int_0^1\tilde{Y}^l(1-\tilde{Y})^{m+n+1}\;dY\\ &=(1-X)^{l+m+n+2}\frac{l!(m+n+1)!}{(l+m+n+2)!} \end{align*}$$

と計算できます. 従って

$$\begin{align*} \int_{\Omega}L_0^kL_1^lL_2^mL_3^n\;dV&=6T\frac{l!m!n!(m+n+1)!}{(m+n+1)!(l+m+n+2)!}\int_0^1X^k(1-X)^{l+m+n}\;dX\\ &=6T\frac{k!l!m!n!(m+n+1)!(l+m+n+2)!}{(m+n+1)!(l+m+n+2)!(k+l+m+n+3)!}\\ &=6T\frac{k!l!m!n!}{(k+l+m+n+3)!} \end{align*}$$

ポテンシャル $V$ を $L_j$ を用いて近似すると

$$V^e=\sum_{j=0}^{3}V^e_jL_j$$

です. 上に示した式を用いて式(3.1)左辺第1項の積分を計算すると

$$\int_{\Omega^e} \begin{pmatrix*} L_0^e\\L_1^e\\L_2^e\\L_3^e \end{pmatrix*} V^e \left(L_0^e\;\;L_1^e\;\;L_2^e\;\;L_3^e\right) \;dV=\frac{T^e}{120} \begin{pmatrix*} 6V_0^e+2V_1^e+2V_2^e+2V_3^e&2V_0^e+2V_1^e+V_2^e+V_3^e&2V_0^e+V_1^e+2V_2^e+V_3^e&2V_0^e+V_1^e+V_2^e+2V_3^e\\ *&2V_0^e+6V_1^e+2V_2^e+2V_3^e&V_0^e+2V_1^e+2V_2^e+V_3^e&V_0^e+2V_1^e+V_2^e+2V_3^e\\ *&*&2V_0^e+2V_1^e+6V_2^e+2V_3^e&V_0^e+V_1^e+2V_2^e+2V_3^e\\ *&*&*&2V_0^e+2V_1^e+2V_2^e+6V_3^e \end{pmatrix*}$$

となります（対称行列なので下三角要素は省略）. この行列を

$$V^e=\frac{T^e}{120} \begin{pmatrix*} 6V_0^e+2V_1^e+2V_2^e+2V_3^e&2V_0^e+2V_1^e+V_2^e+V_3^e&2V_0^e+V_1^e+2V_2^e+V_3^e&2V_0^e+V_1^e+V_2^e+2V_3^e\\ *&2V_0^e+6V_1^e+2V_2^e+2V_3^e&V_0^e+2V_1^e+2V_2^e+V_3^e&V_0^e+2V_1^e+V_2^e+2V_3^e\\ *&*&2V_0^e+2V_1^e+6V_2^e+2V_3^e&V_0^e+V_1^e+2V_2^e+2V_3^e\\ *&*&*&2V_0^e+2V_1^e+2V_2^e+6V_3^e \end{pmatrix*}$$

とおくと式(3.1)左辺第2項は

$$\sum_{e=0}^{N-1}\int_{\Omega^e}{\phi^*}^eV^e\phi^e\;dV=\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right) V^e \begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

となります.

最後に式(3.1)右辺を計算しますが, これは式(3.1)左辺第2項のときと同じように計算できるので途中を省略します. 式(3.1)右辺は行列

$$B^e=\frac{T^e}{20} \begin{pmatrix*} 2&1&1&1\\ 1&2&1&1\\ 1&1&2&1\\ 1&1&1&2 \end{pmatrix*}$$

を用いて

$$E\sum_{e=0}^{N-1}\int_{\Omega^e}\phi^{*e}\phi^e\;dV=E\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right) B^e \begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

以上より, 式(3.1)は

$$\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right)\left(T^eK^e+V^e\right)\begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*} =E\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right)B^e\begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

となり, $A^e=T^eK^e+V^e$ とおくことで

$$\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right)A^e\begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*} =E\sum_{e=0}^{N-1}\left(\phi_0^{*e}\;\;\phi_1^{*e}\;\;\phi_2^{*e}\;\;\phi_3^{*e}\right)B^e\begin{pmatrix*} \phi_0^e\\ \phi_1^e\\ \phi_2^e\\ \phi_3^e \end{pmatrix*}$$

と表されます.

4.全体行列を構成する

1次元, 2次元のときと同じです. $G$ 個の節点 $\phi^e_j$ に対して全節点での通し番号（グローバル節点番号）を割り当て, 最終的に

$$\left(\phi^{*0}\;\;\phi^{*1}\;\;\cdots\;\;\phi^{*G-1}\;\;\phi^{*G}\right)A \begin{pmatrix*} \phi^{0}\\\phi^{1}\\\vdots\\\phi^{G-1}\\\phi^{G} \end{pmatrix*} =E \left(\phi^{*0}\;\;\phi^{*1}\;\;\cdots\;\;\phi^{*G-1}\;\;\phi^{*G}\right)B \begin{pmatrix*} \phi^{0}\\\phi^{1}\\\vdots\\\phi^{G-1}\\\phi^{G} \end{pmatrix*}$$

となるように要素行列の各要素を全体行列に足し合わせましょう. 最後に変分原理を用いて

$$A \begin{pmatrix*} \phi^{0}\\\phi^{1}\\\vdots\\\phi^{G-1}\\\phi^{G} \end{pmatrix*} =EB \begin{pmatrix*} \phi^{0}\\\phi^{1}\\\vdots\\\phi^{G-1}\\\phi^{G} \end{pmatrix*}$$

と変形することで解くべき一般化固有方程式を得ます.

5.境界条件処理

これも1次元, 2次元のときと同じです. 対応する行, および列を削除しましょう.

プログラミング編

6.Pythonで有限要素法

PCスペックの都合でプログラムの正当性を確認できていないですがせっかくなので簡単な説明と併せて載せておきます.

要素クラスを定義します. 要素クラスは四面体の頂点の座標, M行列, 要素の体積を保持します.

class Element():
    # 要素クラス
    points: tuple  # 4頂点の座標
    global_numbers: tuple   # 4頂点のグローバル節点番号
    volume: float  # 要素の体積
    M: np.ndarray  # 要素のM行列

    def __init__(self, p0: tuple, p1: tuple, p2: tuple, p3: tuple, global_numbers: tuple) -> None:
        self.points = (p0, p1, p2, p3)
        self.global_numbers = global_numbers

        S = np.concatenate(
            [self.points, [[1], [1], [1], [1]]], axis=1, dtype=np.float64)
        self.volume = np.abs(npl.det(S))/6
        self.M = npl.inv(S)

        return None

メッシュクラスも定義します. 今回は以下に示すように直方体ごとに5つの四面体に分割します.

メッシュその1（3個）

メッシュその2（1個, 中心）

メッシュその3（1個）

class Mesh():
    # メッシュクラス
    N: int  # 1辺(2N-1)節点
    L: float  # 1辺2L
    X: np.ndarray  # x軸方向の節点座標
    Y: np.ndarray  # y軸方向の節点座標
    Z: np.ndarray  # z軸方向の節点座標
    elements: list  # 要素を保持
    indices: tuple = ((0, 1, 2, 4), (1, 4, 5, 7),
                      (1, 2, 3, 7), (2, 4, 6, 7), (1, 2, 4, 7))  # 要素生成用のインデックス

    def __init__(self, N: int, L: float) -> None:
        self.N = N
        self.L = L

        self.X = np.concatenate(
            [np.linspace(-L, 0, N-1, endpoint=False), np.linspace(0, L, N, endpoint=True)])
        self.Y = np.concatenate(
            [np.linspace(-L, 0, N-1, endpoint=False), np.linspace(0, L, N, endpoint=True)])
        self.Z = np.concatenate(
            [np.linspace(-L, 0, N-1, endpoint=False), np.linspace(0, L, N, endpoint=True)])
        self.elements = []
        return None

    def gen_mesh(self) -> None:
        # 直方体ごとに要素を生成する. 1つの直方体を5つの四面体に分割する.
        for i in range(2*self.N-2):
            for j in range(2*self.N-2):
                for k in range(2*self.N-2):
                    p = []
                    gn = []

                    # 直方体の8節点の座標とグローバル節点番号の割り当て
                    for l in range(2):
                        for m in range(2):
                            for n in range(2):
                                p.append(
                                    (self.X[i+l], self.Y[j+m], self.Z[k+n]))  # 節点を追加
                                if i+l == 0 or j+m == 0 or k+n == 0:
                                    # 境界のグローバル節点番号は-1
                                    gn.append(-1)
                                elif i+l == (2*self.N-2) or j+m == (2*self.N-2) or k+n == (2*self.N-2):
                                    # 境界のグローバル節点番号は-1
                                    gn.append(-1)
                                else:
                                    # グローバル節点番号を計算して割り当て
                                    gn.append((i+l-1)+(2*self.N-3) *
                                              (j+m-1)+((2*self.N-3)**2)*(k+n-1))
                    # 要素生成
                    for idx in self.indices:
                        points = (p[idx[0]], p[idx[1]], p[idx[2]], p[idx[3]])
                        gns = (gn[idx[0]], gn[idx[1]], gn[idx[2]], gn[idx[3]])
                        self.elements.append(Element(*points, gns))
        return None

残りは1次元, 2次元のときとほとんど同じです. 境界にあたる節点は初めから全体行列に含まないようにしています. また, 全体行列を疎行列クラス（SciPyのCSR形式）に変換しています.

class SchrodingerFEM3():

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

    N: int  # 要素数
    N_points: int  # 境界を除いた節点数

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

    A: np.ndarray   # A行列
    B: np.ndarray  # B行列

    elements: list  # 要素を保持
    eig: list  # 計算結果を保持

    # ポテンシャル項の計算用行列
    V_base: np.array = np.array([[[6, 2, 2, 2], [2, 2, 1, 1], [2, 1, 2, 1], [2, 1, 1, 2]], [[2, 2, 1, 1], [2, 6, 2, 2], [1, 2, 2, 1], [1, 2, 1, 2]],
                                 [[2, 1, 2, 1], [1, 2, 2, 1], [2, 2, 6, 2], [1, 1, 2, 2]], [[2, 1, 1, 2], [1, 2, 1, 2], [1, 1, 2, 2], [2, 2, 6, 6]]], np.float64)

    def __init__(self, N_points: int, m: float, V: Callable[[float, float, float], float], elements: list) -> None:
        self.N = len(elements)
        self.N_points = N_points
        self.m = m
        self.V = V
        self.elements = elements

        self.A = np.zeros((self.N_points, self.N_points))
        self.B = np.zeros((self.N_points, self.N_points))

        self.eig = []

        return None

    def gen_K_matrix(self, e: int) -> np.ndarray:
        Ke = np.zeros((4, 4), np.float64)
        for i in range(4):
            for j in range(4):
                for k in range(3):
                    Ke[i][j] += self.elements[e].M[k][i] * \
                        self.elements[e].M[k][j]
        Ke *= (self.DIRAC_CONSTANT_DOUBLE/(2*self.m))
        return Ke

    def gen_V_matrix(self, e: int) -> np.ndarray:
        # 要素の各頂点におけるポテンシャルの値
        V = [self.V(p[0], p[1], p[2]) for p in self.elements[e].points]
        Ve = np.zeros((4, 4), np.float64)
        for i in range(4):
            Ve += V[i]*self.V_base[i]
        Ve *= (self.elements[e].volume/120)
        return Ve

    def gen_A_element(self, e: int) -> np.ndarray:
        Ae = np.zeros((4, 4), np.float64)
        Ae += self.elements[e].volume*self.gen_K_matrix(e)+self.gen_V_matrix(e)
        return Ae

    def gen_B_element(self, e: int) -> np.ndarray:
        Be = np.array([[2, 1, 1, 1], [1, 2, 1, 1], [
                      1, 1, 2, 1], [1, 1, 1, 2]], np.float64)
        Be *= (self.elements[e].volume/20)
        return Be

    def gen_matrix(self) -> None:
        for e in range(self.N):
            Ae = self.gen_A_element(e)
            Be = self.gen_B_element(e)
            for i in range(4):
                i_gn = self.elements[e].global_numbers[i]
                if i_gn < 0:
                    continue
                for j in range(4):
                    j_gn = self.elements[e].global_numbers[j]
                    if j_gn < 0:
                        continue
                    self.A[i_gn][j_gn] += Ae[i][j]
                    self.B[i_gn][j_gn] += Be[i][j]
        self.A = sps.csr_matrix(self.A)
        self.B = sps.csr_matrix(self.B)
        return None

    def solve(self) -> None:
        val, vec = spsl.eigsh(self.A, self.N_points-2, self.B,
                              maxiter=self.N_points*600, which='SM', sigma=0.01)
        for i in range(val.shape[0]):
            self.eig.append((val[i], vec[:, i]))
        self.eig.sort(key=lambda x: x[0])
        return None

最後にまとめたものも載せておきます.

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


class Element():
    # 要素クラス
    points: tuple  # 4頂点の座標
    global_numbers: tuple   # 4頂点のグローバル節点番号
    volume: float  # 要素の体積
    M: np.ndarray  # 要素のM行列

    def __init__(self, p0: tuple, p1: tuple, p2: tuple, p3: tuple, global_numbers: tuple) -> None:
        self.points = (p0, p1, p2, p3)
        self.global_numbers = global_numbers

        S = np.concatenate(
            [self.points, [[1], [1], [1], [1]]], axis=1, dtype=np.float64)
        self.volume = np.abs(npl.det(S))/6
        self.M = npl.inv(S)

        return None


class Mesh():
    # メッシュクラス
    N: int  # 1辺(2N-1)節点
    L: float  # 1辺2L
    X: np.ndarray  # x軸方向の節点座標
    Y: np.ndarray  # y軸方向の節点座標
    Z: np.ndarray  # z軸方向の節点座標
    elements: list  # 要素を保持
    indices: tuple = ((0, 1, 2, 4), (1, 4, 5, 7),
                      (1, 2, 3, 7), (2, 4, 6, 7), (1, 2, 4, 7))  # 要素生成用のインデックス

    def __init__(self, N: int, L: float) -> None:
        self.N = N
        self.L = L

        self.X = np.concatenate(
            [np.linspace(-L, 0, N-1, endpoint=False), np.linspace(0, L, N, endpoint=True)])
        self.Y = np.concatenate(
            [np.linspace(-L, 0, N-1, endpoint=False), np.linspace(0, L, N, endpoint=True)])
        self.Z = np.concatenate(
            [np.linspace(-L, 0, N-1, endpoint=False), np.linspace(0, L, N, endpoint=True)])
        self.elements = []
        return None

    def gen_mesh(self) -> None:
        # 直方体ごとに要素を生成する. 1つの直方体を5つの四面体に分割する.
        for i in range(2*self.N-2):
            for j in range(2*self.N-2):
                for k in range(2*self.N-2):
                    p = []
                    gn = []

                    # 直方体の8節点の座標とグローバル節点番号の割り当て
                    for l in range(2):
                        for m in range(2):
                            for n in range(2):
                                p.append(
                                    (self.X[i+l], self.Y[j+m], self.Z[k+n]))  # 節点を追加
                                if i+l == 0 or j+m == 0 or k+n == 0:
                                    # 境界のグローバル節点番号は-1
                                    gn.append(-1)
                                elif i+l == (2*self.N-2) or j+m == (2*self.N-2) or k+n == (2*self.N-2):
                                    # 境界のグローバル節点番号は-1
                                    gn.append(-1)
                                else:
                                    # グローバル節点番号を計算して割り当て
                                    gn.append((i+l-1)+(2*self.N-3) *
                                              (j+m-1)+((2*self.N-3)**2)*(k+n-1))
                    # 要素生成
                    for idx in self.indices:
                        points = (p[idx[0]], p[idx[1]], p[idx[2]], p[idx[3]])
                        gns = (gn[idx[0]], gn[idx[1]], gn[idx[2]], gn[idx[3]])
                        self.elements.append(Element(*points, gns))
        return None


class SchrodingerFEM3():

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

    N: int  # 要素数
    N_points: int  # 境界を除いた節点数

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

    A: np.ndarray   # A行列
    B: np.ndarray  # B行列

    elements: list  # 要素を保持
    eig: list  # 計算結果を保持

    # ポテンシャル項の計算用行列
    V_base: np.array = np.array([[[6, 2, 2, 2], [2, 2, 1, 1], [2, 1, 2, 1], [2, 1, 1, 2]], [[2, 2, 1, 1], [2, 6, 2, 2], [1, 2, 2, 1], [1, 2, 1, 2]],
                                 [[2, 1, 2, 1], [1, 2, 2, 1], [2, 2, 6, 2], [1, 1, 2, 2]], [[2, 1, 1, 2], [1, 2, 1, 2], [1, 1, 2, 2], [2, 2, 6, 6]]], np.float64)

    def __init__(self, N_points: int, m: float, V: Callable[[float, float, float], float], elements: list) -> None:
        self.N = len(elements)
        self.N_points = N_points
        self.m = m
        self.V = V
        self.elements = elements

        self.A = np.zeros((self.N_points, self.N_points))
        self.B = np.zeros((self.N_points, self.N_points))

        self.eig = []

        return None

    def gen_K_matrix(self, e: int) -> np.ndarray:
        Ke = np.zeros((4, 4), np.float64)
        for i in range(4):
            for j in range(4):
                for k in range(3):
                    Ke[i][j] += self.elements[e].M[k][i] * \
                        self.elements[e].M[k][j]
        Ke *= (self.DIRAC_CONSTANT_DOUBLE/(2*self.m))
        return Ke

    def gen_V_matrix(self, e: int) -> np.ndarray:
        # 要素の各頂点におけるポテンシャルの値
        V = [self.V(p[0], p[1], p[2]) for p in self.elements[e].points]
        Ve = np.zeros((4, 4), np.float64)
        for i in range(4):
            Ve += V[i]*self.V_base[i]
        Ve *= (self.elements[e].volume/120)
        return Ve

    def gen_A_element(self, e: int) -> np.ndarray:
        Ae = np.zeros((4, 4), np.float64)
        Ae += self.elements[e].volume*self.gen_K_matrix(e)+self.gen_V_matrix(e)
        return Ae

    def gen_B_element(self, e: int) -> np.ndarray:
        Be = np.array([[2, 1, 1, 1], [1, 2, 1, 1], [
                      1, 1, 2, 1], [1, 1, 1, 2]], np.float64)
        Be *= (self.elements[e].volume/20)
        return Be

    def gen_matrix(self) -> None:
        for e in range(self.N):
            Ae = self.gen_A_element(e)
            Be = self.gen_B_element(e)
            for i in range(4):
                i_gn = self.elements[e].global_numbers[i]
                if i_gn < 0:
                    continue
                for j in range(4):
                    j_gn = self.elements[e].global_numbers[j]
                    if j_gn < 0:
                        continue
                    self.A[i_gn][j_gn] += Ae[i][j]
                    self.B[i_gn][j_gn] += Be[i][j]
        self.A = sps.csr_matrix(self.A)
        self.B = sps.csr_matrix(self.B)
        return None

    def solve(self) -> None:
        val, vec = spsl.eigsh(self.A, self.N_points-2, self.B,
                              maxiter=self.N_points*600, which='SM', sigma=0.01)
        for i in range(val.shape[0]):
            self.eig.append((val[i], vec[:, i]))
        self.eig.sort(key=lambda x: x[0])
        return None


if __name__ == '__main__':
    N = 8
    L = 3.
    m = SchrodingerFEM3.DIRAC_CONSTANT_DOUBLE

    mesh = Mesh(N, L)
    mesh.gen_mesh()

    fem = SchrodingerFEM3((2*N-3)**3, m, lambda x, y, z: 0., mesh.elements)
    fem.gen_matrix()
    fem.solve()

    for e in fem.eig[:30]:
        print(e[0])

参考

mutsumunemitsutan. Hatena Blog 数学とか語学とか楽しいよね. "【有限要素法】2次元有限要素法における面積座標（局所座標）の積分公式". 2018/10/08. https://mathlang.hatenablog.com/entry/2018/10/08/122539