はじめまして、りょーつといいます。高専出身の大学院2年生です。研究の専門は力学や機構学で、Qiitaでは主に制御工学や数学に関する記事を書いています。
先日、恋人と量子力学の話題で盛り上がったことをきっかけに、ノリと勢いでシュレディンガー方程式を解いてみました。MATLABコードも載せておいたのでぜひ遊んでみてください。
目次
1.はじめに
2.シュレディンガー方程式の離散化
3.初期条件と正規化の扱い
4.シミュレーション
5.おわりに
1. はじめに
皆さんシュレディンガー方程式をご存じでしょうか？シュレディンガー方程式は、理学部や工学部では避けては通ることができない量子力学で頻繁に登場する、謎の偏微分方程式で、以下に示す見た目をしています。
i \hbar \frac{\partial }{\partial t} \psi(\boldsymbol{r},t)
=
\bigg \lbrack
\frac{-\hbar^2}{2m} \nabla^2 + V(\boldsymbol{r},t)
\bigg \rbrack
\psi(\boldsymbol{r},t)
\tag{1}
深い議論はさておき、物理でよく登場する運動方程式の量子力学バージョンという認識があれば十分だと思います(しらんけど)。ただし、$\psi(\boldsymbol{r},t)$が複素数であることに注意が必要です。
シュレディンガー方程式は、見た目通りとてつもなく複雑なので解析解を求めることが容易ではありません。大学の講義でも1次元かつ定常状態を仮定し、簡略化された以下のシュレディンガー方程式の解を求めることが多いかと思います。
\bigg \lbrack
\frac{-\hbar^2}{2m} \frac{\partial^2 }{{\partial x}^2} + V(x)
\bigg \rbrack
\psi(x)
=
0
\tag{2}
講義を受けたことがある方は共感いただけると思いますが、ちょっとでも変な条件を与えると、(2)式ですら解くのがとても大変になります(そのヤバさのおかげて解析学の級数解法や特殊多項式の分野が発展したらしい)。
上記のことからも、シュレディンガー方程式は人間の手に負える代物ではないと考えられます。
そこで本記事では、解析解を出すのを完全にあきらめて、MATLABに計算を丸投げする方法についてご紹介します。解析解をあきらめることで非定常状態のシミュレーションも比較的簡単に行うことができます。シュレディンガー方程式に関する深い理解は必要ありませんが、線形代数を多用するので数学の復習をしておくと良きです。
2. シュレディンガー方程式の離散化
では早速シュレディンガー方程式を離散化してみましょう。離散化ってどうやるの？と思われるかもしれませんが、内容はシンプルです。空間や時間をいくつかの格子で区切り、格子点どうしで偏微分の定義に沿って差分を取ってあげることで離散化ができます。
さすがに3次元非定常のシュレディンガー方程式は解いても可視化が難しいので、今回は1次元非定常のシュレディンガー方程式を離散化してみましょう。離散化する前の式は以下のように表されます。
i \hbar \frac{\partial }{\partial t} \psi(x,t)
=
\bigg \lbrack
\frac{-\hbar^2}{2m} \frac{\partial^2 }{{\partial x}^2} + V(x,t)
\bigg \rbrack
\psi(x,t)
\tag{3}
ここで、$i$は虚数単位、$x$は空間の位置、$t$は時刻、$\hbar$は正規化プランク定数、$m$は量子質量、$V(x,t)$は空間のポテンシャル、$\psi(x,t)$は波動関数を表します。まずは波動関数を離散化してみます。時刻$t$の刻み幅$\Delta t$で無限個に、空間$x$を刻み幅$\Delta x$で$K$個に離散化したとします。$\Delta$はラプラシアンではないのでご注意ください。このとき、離散化した後の各時刻と位置は以下のように計算できます。
t
=
n \Delta t \ \ \ \ \ (n = 0,1,2,\cdots)
\tag{4}
x
=
k \Delta x \ \ \ \ \ (k = 1,2,\cdots ,K)
\tag{5}
時刻と位置は$n$と$k$に対応するため、ある時刻・位置における波動関数の値を以下のように表記することにします。
\psi(k\Delta x, n \Delta t)
=
\psi^n_k
\tag{6}
さらに、数値計算を行う際は境界条件を考える必要があります。とりあえず簡単のために空間の端っこでは波動関数が$0$になると仮定しましょう。
\psi^n_0
=
\psi^n_{K+1}
=
0
\tag{7}
ここまでの準備を踏まえて偏微分を計算してみましょう。
まずは波動関数の時間偏微分を計算します。一番シンプルな後退差分法を使って離散化しましょう。となりあう格子の差を取ればいいので以下のようになるはずです。
\frac{\partial }{\partial t} \psi^n_k
=
\frac{\psi^n_k - \psi^{n-1}_k}{\Delta t}
\tag{8}
次に波動関数の$x$に関する2階偏微分を計算します。このとき、中心差分法を使用します。やり方は(8)式とほぼ同じですが、2回微分するので少しややこしいです。まずは1回微分してみましょう。
\frac{\partial }{\partial x} \psi^n_k
=
\frac{\psi^n_k - \psi^n_{k-1}}{\Delta x}
\tag{9}
1つズレた位置についても同様に計算しておきます。
\frac{\partial }{\partial x} \psi^n_{k+1}
=
\frac{\psi^n_{k+1} - \psi^n_k}{\Delta x}
\tag{10}
次にもう一度微分しましょう。ここで注意しないといけないのは2回目の微分を前進差分で行う点です。これは、計算結果の右辺について$\psi^n_k$が中心に来るようにしたいためです(計算結果を見るとなんとなくわかる)。
\frac{\partial^2 }{{\partial x}^2} \psi^n_k
=
\frac{\frac{\partial }{\partial x} \psi^n_{k+1} - \frac{\partial }{\partial x} \psi^n_k}{\Delta x}
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
=
\frac{\frac{\psi^n_{k+1} - \psi^n_k}{\Delta x} - \frac{\psi^n_k - \psi^n_{k-1}}{\Delta x}}{\Delta x}
\therefore
\frac{\partial^2 }{{\partial x}^2} \psi^n_k
=
\frac{\psi^n_{k+1} - 2 \psi^n_k + \psi^n_{k-1}}{{\Delta x}^2}
\tag{11}
ここで、
簡略化のために、ポテンシャル$V$やハミルトニアン$H$に対しても(6)式の添え字を適用しましょう。
V(x,t)
=
V^n_k
\tag{13}
H(x,t)
=
H^n_k
\tag{14}
さらに(11)式と(13)式を用いて$H^n_k \psi^n_k$を展開してみましょう。2階偏微分の結果にポテンシャルを足したものになるので以下のような結果が得られるはずです。
H^n_k \psi^n_k
=
-\frac{\hbar^2}{2m}
\frac{\psi^n_{k+1} - 2 \psi^n_k + \psi^n_{k-1}}{{\Delta x}^2} + V^n_k \psi^n_k
\tag{15}
少しイメージが付きづらいのでk=1から順に展開してみましょう。ただし、(7)式で与えた境界条件に注意します。
H^n_1 \psi^n_1
=
-\frac{\hbar^2}{2m}
\frac{\psi^n_{2} - 2 \psi^n_1}{{\Delta x}^2} + V^n_1 \psi^n_1
\tag{16}
H^n_2 \psi^n_2
=
-\frac{\hbar^2}{2m}
\frac{\psi^n_{3} - 2 \psi^n_2 + \psi^n_1}{{\Delta x}^2} + V^n_2 \psi^n_2
\tag{17}
H^n_3 \psi^n_3
=
-\frac{\hbar^2}{2m}
\frac{\psi^n_{4} - 2 \psi^n_3 + \psi^n_2}{{\Delta x}^2} + V^n_3 \psi^n_3
\tag{18}
\vdots
H^n_K \psi^n_K
=
-\frac{\hbar^2}{2m}
\frac{ - 2 \psi^n_K + \psi^n_{K-1}}{{\Delta x}^2} + V^n_K \psi^n_K
\tag{19}
ある程度規則性が見えてきたので(16)～(19)式を行列にまとめておきましょう。
\begin{bmatrix}
H^n_1 \psi^n_1 \\
H^n_2 \psi^n_2 \\
\vdots \\
H^n_K \psi^n_K \\
\\
\end{bmatrix}
=
-\frac{\hbar^2}{2m}
\begin{bmatrix}
-2 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
0 & 1 & -2 & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & 1 \\
0 & 0 & \cdots & 1 & -2
\end{bmatrix}
\begin{bmatrix}
\psi^n_1 \\
\psi^n_2 \\
\vdots \\
\psi^n_K \\
\\
\end{bmatrix}
+
\begin{bmatrix}
V^n_1 & 0 & \cdots & 0 \\
0 & V^n_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & V^n_K
\end{bmatrix}
\begin{bmatrix}
\psi^n_1 \\
\psi^n_2 \\
\vdots \\
\psi^n_K \\
\\
\end{bmatrix}
\tag{20}
(20)式をさらに簡略化したいと思います。そのためにハミルトニアン行列$H^n$と波動関数ベクトル$\boldsymbol{\psi}^n$を定義します。右上の添え字は時刻$n\Delta t$を意味しています。ハミルトニアン行列$H^n \in \mathbb{R}^K \times \mathbb{R}^K$を
H^n
=
-\frac{\hbar^2}{2m}
\begin{bmatrix}
-2 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
0 & 1 & -2 & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & 1 \\
0 & 0 & \cdots & 1 & -2
\end{bmatrix}
+
\begin{bmatrix}
V^n_1 & 0 & \cdots & 0 \\
0 & V^n_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & V^n_K
\end{bmatrix}
\tag{21}
波動関数ベクトル$\boldsymbol{\psi}^n \in \mathbb{R}^K$を
\boldsymbol{\psi}^n
=
\begin{bmatrix}
\psi^n_1 \\
\psi^n_2 \\
\vdots \\
\psi^n_K \\
\\
\end{bmatrix}
\tag{22}
とすれば、(20)式の計算、すなわちシュレディンガー方程式の右辺は以下の簡単な行列の積にまとめることができます。
\begin{bmatrix}
H^n_1 \psi^n_1 \\
H^n_2 \psi^n_2 \\
\vdots \\
H^n_K \psi^n_K \\
\\
\end{bmatrix}
=
H^n \boldsymbol{\psi}^n
\tag{23}
離散化したシュレディンガー方程式の左辺にあたる(8)式も(14)式を用いて離散化したいのですが、少し注意が必要です。(8)式の左辺には$n$番目の時刻の波動関数$\psi^n_k$と$n-1$番目の時刻の波動関数$\psi^n_{k-1}$が含まれています。そのため、右辺に$H^n_k \psi^n_k$を用いるか$H^n_{k-1} \psi^n_{k-1}$を用いるかが悩ましいのです。
結論を述べると、シンプルにこれらの平均を取った値を用いることにします。以上の議論を踏まえて(3)式を離散化すると以下のようになります。
i \hbar \frac{\psi^n_k - \psi^{n-1}_k}{\Delta t}
=
\frac{H^n_k \psi^n_k + H^{n-1}_k \psi^{n-1}_k}{2}
\therefore
\psi^n_k - \psi^{n-1}_k
=
-\frac{i\Delta t}{2\hbar}
(H^n_k \psi^n_k + H^{n-1}_k \psi^{n-1}_k)
\tag{24}
では(24)式も(16)～(19)式と同じように展開してみましょう。
\psi^n_1 - \psi^{n-1}_1
=
-\frac{i\Delta t}{2\hbar}
(H^n_1 \psi^n_1 + H^{n-1}_1 \psi^{n-1}_1)
\tag{25}
\psi^n_2 - \psi^{n-1}_2
=
-\frac{i\Delta t}{2\hbar}
(H^n_2 \psi^n_2 + H^{n-1}_2 \psi^{n-1}_2)
\tag{26}
\vdots
\psi^n_K - \psi^{n-1}_K
=
-\frac{i\Delta t}{2\hbar}
(H^n_K \psi^n_K + H^{n-1}_K \psi^{n-1}_K)
\tag{27}
なんだかとても綺麗な規則性が得られましたね。(25)～(27)式をベクトルにまとめてみましょう。
\begin{bmatrix}
\psi^n_1 \\
\psi^n_2 \\
\vdots \\
\psi^n_K \\
\\
\end{bmatrix}
-
\begin{bmatrix}
\psi^{n-1}_1 \\
\psi^{n-1}_2 \\
\vdots \\
\psi^{n-1}_K \\
\\
\end{bmatrix}
=
-\frac{i\Delta t}{2\hbar}
\begin{bmatrix}
H^n_1 \psi^n_1 \\
H^n_2 \psi^n_2 \\
\vdots \\
H^n_K \psi^n_K \\
\\
\end{bmatrix}
-\frac{i\Delta t}{2\hbar}
\begin{bmatrix}
H^{n-1}_1 \psi^{n-1}_1 \\
H^{n-1}_2 \psi^{n-1}_2 \\
\vdots \\
H^{n-1}_K \psi^{n-1}_K \\
\\
\end{bmatrix}
\tag{28}
さらに(22)式と(23)式を使って(28)式を簡略化すると以下の関係が得られます。
\boldsymbol{\psi}^n - \boldsymbol{\psi}^{n-1}
=
-\frac{i\Delta t}{2\hbar}
(H^n \boldsymbol{\psi}^n + H^{n-1} \boldsymbol{\psi}^{n-1})
\tag{29}
少しややこしい式変形や大きな行列が出てきて大変でしたが、(29)式は(3)式を離散化したものとなります。
最後に(29)式を変形して、シュレディンガー方程式の時間発展を求める計算式を作りましょう。移行して逆行列をかけるだけです。
\boldsymbol{\psi}^n
+
\frac{i\Delta t}{2\hbar}H^n \boldsymbol{\psi}^n
=
\boldsymbol{\psi}^{n-1}
-
H^{n-1} \boldsymbol{\psi}^{n-1}
\bigg( I + \frac{i\Delta t}{2\hbar}H^n \bigg)
\boldsymbol{\psi}^n
=
\bigg( I - \frac{i\Delta t}{2\hbar}H^{n-1} \bigg)
\boldsymbol{\psi}^{n-1}
\therefore
\boldsymbol{\psi}^n
=
{\bigg( I + \frac{i\Delta t}{2\hbar}H^n \bigg)}^{-1}
\bigg( I - \frac{i\Delta t}{2\hbar}H^{n-1} \bigg)
\boldsymbol{\psi}^{n-1}
\tag{30}
(30)式の表現行列にあたるものがシュレディンガー方程式の遷移行列になります。詳しくは第4章のシミュレーションで再度確認します。
3. 初期条件と正規化の扱い
第2章ではシュレディンガー方程式の離散化を行い、(30)式ではその遷移行列を導出しました。本章では(30)式の内容をMATLAB上で実装し、シュレディンガー方程式を実際に解いてみようと思います。
シュレディンガー方程式の時間発展を考える上で必要になるのは、波動関数の初期値$\boldsymbol{\psi}^{0}$です。この初期値の与え方が重要になります。
通常の数値計算であれば、初期値を適当に決めてもよいのですが、波動関数は確率密度関数と密接に関係しているため、適切に正規化しないといけません。
まずは正規化の話について少しまとめます。具体的には以下の制約があります。
\int_{-\infty}^{\infty}
\int_{-\infty}^{\infty}
\int_{-\infty}^{\infty}
{| \psi(\boldsymbol{r},t) |}^2
dx
dy
dz
=
1
\tag{31}
量子力学では、波動関数のノルムの2乗${|| \boldsymbol{\psi}(\boldsymbol{r},t) ||}^2$が確率密度に対応するため、全空間で積分した結果は1になるよ、ということを表しています。本記事では1次元のシュレディンガー方程式を扱うので、(31)式を簡略化し、以下の関係式を使用します。
\int_{-\infty}^{\infty}
{| \psi(x,t) |}^2
dx
=
1
\tag{32}
波動関数は複素数で表されるため、被積分関数であるノルムの計算には共役複素数を使用します。波動関数$\psi(x,t)$の共役複素数を$\psi^*(x,t)$としたとき、
\psi^* (x,t)
=
Re[\psi(x,t)] - i \ Im[\psi(x,t)]
\tag{33}
がなりたちます。こうすることで、ノルムの2乗${| \boldsymbol{\psi}(x,t) |}^2$が複素数の積で記述できるようになります。
\psi^*(x,t) \psi(x,t)
=
\big(Re[\psi(x,t)] - i \ Im[\psi(x,t)] \big)
\big(Re[\psi(x,t)] + i \ Im[\psi(x,t)] \big)
=
Re[\psi(x,t)]^2 + Im[\psi(x,t)]^2
\therefore
| \psi(x,t) |^2
=
\psi^*(x,t) \psi(x,t)
\tag{34}
本記事ではシュレディンガー方程式を離散化しているため、(32)式で表される正規化の操作も離散化して考える必要があります。ここでは積分を長方形近似で表してみましょう。ノルムが(34)式で与えられることを用いると、(32)式の離散化バージョンは以下のように行列の積で記述できると思います。
| \psi^n_1 |^2 \Delta x
+
| \psi^n_2 |^2 \Delta x
+
\cdots
+
| \psi^n_K |^2 \Delta x
=
1
\Delta x
\big(
{\psi^n_1}^* \psi^n_1
+
{\psi^n_2}^* \psi^n_2
+
\cdots
+
{\psi^n_K}^* \psi^n_K
\big)
=
1
\therefore
\Delta x
\begin{bmatrix}
{\psi^n_1}^* & {\psi^n_2}^* & \cdots & {\psi^n_K}^* \\
\end{bmatrix}
\begin{bmatrix}
\psi^n_1 \\
\psi^n_2 \\
\vdots \\
\psi^n_K \\
\\
\end{bmatrix}
=
1
\tag{35}
(35)式を見てお気づきの方もいるかと思いますが、エルミート転置を使うと非常に美しく表記することができます。行列のエルミート転置は、成分を全て共役複素数に変更し、転置を取ったものです。波動関数ベクトルのエルミート転置${\boldsymbol{\psi}^n}^*$は以下のように表されます。
{\boldsymbol{\psi}^n}^*
=
\begin{bmatrix}
{\psi^n_1}^* & {\psi^n_2}^* & \cdots & {\psi^n_K}^* \\
\end{bmatrix}
\tag{36}
以上のことから離散化した正規化条件は以下のようにまとめられます。
\Delta x \ {\boldsymbol{\psi}^n}^* \boldsymbol{\psi}^n
=
1
\tag{37}
初期条件の話に戻りましょう。(37)式は時間の情報$n$によらず常に成立する必要があるので、初期条件$\boldsymbol{\psi}^0$も
\Delta x \ {\boldsymbol{\psi}^0}^* \boldsymbol{\psi}^0
=
1
\tag{38}
を満たすように設定する必要があるのです。
最も簡単なのはインパルス応答を与える方法ですね。例えばある位置$j \in {1,2,\cdots ,K}$にインパルスを発生させる場合、以下のように設定するとよいです。
\psi^0_k
=
\frac{\delta_{jk}}{\sqrt{\Delta x}}
\tag{39}
(39)式中の$\delta_{jk}$はクロネッカーのデルタを意味しています。
また、数値計算特有の注意点として、常に(37)式がなりたつように波動関数$\boldsymbol{\psi}^n$を補正してあげる必要があります。具体的には以下に示す式で$\boldsymbol{\psi}^n$を補正します。
\boldsymbol{\psi}^n
\gets
\frac{\boldsymbol{\psi}^n}{\sqrt{\Delta x\ {\boldsymbol{\psi}^n}^* \boldsymbol{\psi}^n}}
\tag{40}
4. シミュレーション
第2章と第3章でいろいろ書きましたが、実際にプログラムで可視化したほうが分かりやすいと思い、MATLABで動くコードを書きました。変数名やコードの内容は、本記事の内容に対応させているので、数式と見比べてもらうと楽しいかなと思います。
ひとまず以下のコードをコピペして実行してみてください。MATLAB2024では問題なく動作しました。MATLABのすごいところは複素数を脳死で扱える点ですね(^▽^)
% 初期化動作---------------------------------------------------------------
clc
clear
% シュレディンガー方程式の設定-----------------------------------------------
hbar = 1; % 換算プランク定数\hbar
m = 80; % 量子質量m[kg]
% 数値計算の設定------------------------------------------------------------
K = 200; % 空間の最大値 (k = 1,2,・・・, K)
N = 10000; % 計算終了ループ回数 (n = 0,1,・・・,N)
Delta_x = 0.01; % 空間の格子幅Δx[m]
Delta_t = 0.001; % 時間の格子幅Δt[s]
% 変数の定義----------------------------------------------------------------
psi = zeros(K, 1); % 波動関数ベクトルψ^nの宣言
V = zeros(K,1); % ポテンシャル設定用ベクトルの宣言
I = eye(K); % 単位行列Iの定義
H_n = zeros(K); % ハミルトン行列H^nを定義
H_n_1 = zeros(K); % ハミルトン行列H^{n-1}を定義
% 初期状態設定--------------------------------------------------------------
psi(90) = 1/sqrt(Delta_x); % 波動関数ベクトルψ^0
for k = 1:K ,V(k) = Set_Potential(0,k); end % ポテンシャルV^0_k
H_n = -hbar*hbar/(2*m*Delta_x*Delta_x)*gallery('tridiag',K, 1, -2, 1) + diag(V); % ハミルトニアンH^0
% アニメーションウィンドウの設定----------------------------------------------
figure;
% 横軸の設定
x = (1:K)*1;
xlabel('k', 'FontSize', 10);
xlim([1 K]);
% 第1縦軸(ポテンシャル表示)の設定
yyaxis left
handle_v = plot(x,V, '-', 'Color', [1 0.7 0.7], 'LineWidth', 1.2);
ylabel('Potential', 'FontSize',10, 'Color', 'k');
ylim([0 500]);
set(gca, 'YColor', 'k');
% 第2縦軸(確立密度表示)の設定
yyaxis right
handle_d = plot(x,abs(psi).^2, 'k-', 'LineWidth', 1.5);
ylabel('Probability density', 'FontSize', 10,'Color', 'k');
ylim([0 10]);
set(gca, 'YColor', 'k');
% 時間展開------------------------------------------------------------------
for n = 1:N
H_n_1 = H_n; % ハミルトニアンのシフト操作H^{n-1} <- H^n
for k = 1:K ,V(k) = Set_Potential(n,k); end % ポテンシャル項V(x)を更新
H_n = -hbar*hbar/(2*m*Delta_x*Delta_x)*gallery('tridiag',K, 1, -2, 1) + diag(V); % ハミルトニアンH^nを計算
psi = (I + 1i*Delta_t/(2*hbar)*H_n)\(I - 1i*Delta_t/(2*hbar)*H_n_1)*psi; % psi^n = (I + (iΔt)/(2h)*H^n)^{-1} (I - (iΔt)/(2h)*H^{n-1})psi^{n-1}の計算
set(handle_v, 'YData', V); % ポテンシャルを赤線で描画
set(handle_d, 'YData', abs(psi).^2); % 確立密度を黒線で描画
drawnow; % アニメーションの更新
pause(0.005); % 一時停止
psi = psi / sqrt(Delta_x * (psi')*psi); % sqrt(Δx ψ*ψ)で波動関数を正規化(数値計算由来の誤差を修正)
end
% ポテンシャルの設定---------------------------------------------------------
function output = Set_Potential(n,k)
if k < 40+0.01*n, V(k) = 200 + 0.1*n; % 迫りくる井戸型ポテンシャルを設定
elseif k > 160-0.01*n, V(k) = 200 + 0.1*n; % 迫りくる井戸型ポテンシャルを設定
else , V(k) = 0;
end
output = V(k);
end
%--------------------------------------------------------------------------
実行すると、以下に示す図のようにゆっくり動く赤線と激しく動く黒線が描画されると思います。赤線がポテンシャル、黒線が各位置における確率密度を表します。
無事に実行できたら、プログラム内のいろいろな変数をいじってみてください。量子質量$m$を大きくすると拡散がゆっくりになったり、ポテンシャル$V^n_k$の値を大きくすると波が反射しやすくなったりして、シュレディンガー方程式のイメージが付きやすくなると思います！
5. おわりに
本稿ではシュレディンガー方程式を差分法で離散化し、MATLABで波動関数のふるまいをシミュレーションする方法について紹介しました。自身の専門分野ではないので、理解に苦しむ部分もありましたが、なんとか記事にまとまるくらいまで理解できました。数値計算の手法として有限要素法を使ったものもよく用いられているようなので、やる気が出た時にそっちにも挑戦してみようと思います。
今週も最後まで読んでいただき、ありがとうございました！