見出し画像

「線形代数の半歩先」をPythonで写経 ～ ４章 時間発展方程式、微分方程式、偏微分方程式

ネイピア DS

第４部「ならべた数と移りゆく世界」

書籍の著者　大久保 潤 先生

この記事は、書籍「線形代数の半歩先」の 第４部「ならべた数と移りゆく世界」に掲載の「行列と時間発展系」に関する Python写経活動 のドキュメンタリーです。

第４部は「時間発展方程式」という分野になります。
この記事では、第19話、第21話、第23話の さまざまな時間発展方程式 に取り組みます！

第４部の世界は私にとって未知ワールド。
書籍のページをパラパラめくって直ぐに 目が点に なりました。

よくよく読んで気がついた「時間発展方程式」の特徴は次の２点です。

  • 時間経過とともに何かの法則に従って変化するデータを扱う領域

  • 何かの法則は「微分」のかたち

さらに目が点。。。

数学素人なので、どうぞお手柔らかにお願いいたします。
では書籍を開いて線形代数の旅に出発です🚀

いろいろな目の表情のイラスト：「いらすとや」さんより

はじめに

このブログシリーズは、書籍「線形代数の半歩先 データサイエンス・機械学習に挑む前の30話」（講談社サイエンティフィク、「テキスト」と呼びます）の Python 写経の実践を通じて得た個人的な知見を書きます。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

第４部 ならべた数と移りゆく世界

コトバンクさんの「発展方程式」の定義をお借りします。

時間とともに状態が変化していく物理現象があって，その変化のしかたが過去の履歴には無関係に，その瞬間の状態によって定まるならば，その変化のしかたは一つの微分方程式で記述される。このように状態の時間的発展を記述する方程式を，一般に発展方程式という。

https://kotobank.jp/word/%E7%99%BA%E5%B1%95%E6%96%B9%E7%A8%8B%E5%BC%8F

テキストも「時間 $${t}$$ の関数」「理学系、苦学系、経済学系など、さまざまな分野で使われます」と記載しています。
私はこの分野にほとんど触れたことがないので、超初歩的な部分で躓きながら、ほふく前進で写経を進めてまいりました。

ほふく前進のイラスト（自衛隊）：「いらすとや」さんより

この記事で取り扱うテーマは以下のとおりです。

  • 濃度の時間発展（図 19.3）

  • バネの時間発展

    • 固有値を使って解く方法（図 21.1、図 21.2）

    • 数値的に解く方法：オイラー法・ホイン法（図 21.3）

  • オルンシュタイン・ウーレンベック過程

    • エルミート関数を基底関数にして解く方法（図 23.2、図 23.4 [1]）

    • 空間を細かく分割して解く方法（図 23.4 [2]）

では、この記事のコードで利用するライブラリのインポートからスタートしましょう。

### インポート

# 数値計算
import numpy as np

# # 数学演算用の個別関数
from scipy.special import hermite       # エルミート多項式
from scipy.special import factorial     # 階乗

# # 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

１．濃度の時間発展

① 概要
テキスト第19話「移り変わりを数式で表現する」の最初に出現する微分方程式が「濃度の時間発展」です。
こちらのWebサイト（大阪大学大学院）がテキストに近い内容（に見える）ので、参考になるかもです。

テキストの数式をお借りします。
数式を起点にして図を写経いたします。

$$
\begin{align*}
\cfrac{d}{dt}\ x(t) &= -\alpha x(t) \tag{19.1} \\
x(t) &= x(0) e^{-\alpha t} \tag{19.2} \\
\end{align*}
$$

テキスト p.170, 171 の式（19.1）と（19.2）を引用

式（19.1）が微分方程式、式（19.2）は（19.1）の解です。
$${\alpha}$$ は減少の速さに関係するパラメータであり、式（19.2）では時間 $${t=0}$$ の初期条件 $${x(0)}$$ の肩に乗って指数関数を形作っています。
「肩に乗る」が大切なキーワードです。

② 実装
$${\alpha=3}$$、初期条件 $${x(0)= [0.1, 1, 3]}$$ を用いて、図 19.3「時間発展の様子：見た目が違っても、視点を変えれば本質は同じ」を実装します。

### p.171 図19.3  時間発展の様子：見た目が違っても、視点を変えれば本質は同じ

## 設定と準備
# 減少パラメータα
alpha = 3
# 初期状態 x0のリスト
x0_list = [3, 0.1, 1]
# 時間ｔ
times = np.linspace(0, 3, 1001)

## 描画処理
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.5), tight_layout=True)

# 濃度xの初期値0ごとにxの計算と可視化を繰り返し処理
for x0 in x0_list:
    # xの計算結果を格納する配列の初期化
    x = np.zeros(len(times), dtype=float)
    # 時間ｔごとに時間発展方程式の計算を繰り返し処理
    for i, t in enumerate(times):
        # 時間発展方程式の計算 ※式（19.2）
        x[i] = x0 * np.exp(-alpha * t)
    # [1] 時間発展方程式の計算結果の曲線の描画
    ax1.plot(times, x, label=f'$\\alpha={{{alpha}}}$, $x_0 = {{{x0}}}$')
    # [2] 時間発展方程式の計算結果を対数軸で描画
    ax2.plot(times, x, label=f'$\\alpha={{{alpha}}}$, $x_0 = {{{x0}}}$')

# 共通修飾
for ax in [ax1, ax2]:
    ax.set_xlabel('$t$', fontsize=12)
    ax.set_ylabel('$x(t)$', fontsize=12)
    ax.set(xticks=(0, 1, 2, 3))
# 個別修飾
ax1.set(title='[1] 初期値に応じて違う振る舞い？', yticks=(0, 1, 2, 3))
ax2.set_yscale('log')
ax2.set(title='[2] 縦軸を対数にするとどれも直線')
ax2.legend(bbox_to_anchor=(1.5, 1))

plt.show()

【実行結果】
指数的な減少を対数変換すると直線的になるのって面白いですね。

次のコードで時間発展方程式の解（19.2）を表現しています。

x[i] = x0 * np.exp(-alpha * t)

対数変換はチャート描画ライブラリである matplotlib に任せています。
次のように設定をするだけで、作図できます。

ax2.set_yscale('log')

２．バネの時間発展

① 概要
第19話にはもう一つの微分方程式（連立微分方程式）である「バネの時間発展」のイントロダクションが含まれています。
シリンダ内でバネの付いたピストンが動く仕組みを考えます。
バネの自然長を出発点にして、ピストンを引っ張って手を離すと、バネに引かれてピストンが左右に移動します。
テキストの図 19.4 を見ると数式と物理的な内容の繋がりがわかると思います。
こちらは数式で。

$$
\begin{align*}
m \cfrac{d^2}{dt^2}(t) &= -kx(t) - \gamma v(t) \tag{19.4} \\
\cfrac{d}{dt}x(t) &= v(t) \tag{19.5}  \\
\end{align*}
$$

テキスト p.172 式（19.4）（19.5）を引用

$${x(t)}$$ は時間 $${t}$$ におけるピストンの自然長からのずれの座標です。
式（19.4）は $${x(t)}$$ の加速度（２階微分）、式（19.5）の $${v(t)}$$ はピストンの速度（１階微分）です。
$${m}$$ はピストンの質量、$${k}$$ はバネ定数、$${\gamma}$$ は摩擦係数です。

えーっと、この数式から何をどうすればいいの！？

ということで、第21話「いくつかの時間発展を一度に解く」にジャンプして２つの解法を見に行きましょう！

◆ ◆ ◆

② 固有値を使って解く方法
３章で触れた「固有値」「固有ベクトル」を活用する解き方です。
テキストにはたくさんの数式が書かれています。
プログラミングに使えそうな数式をピックアップしました。

$$
\begin{align*}
L &= \begin{bmatrix} 0 & 1 \\ -(k/m) & -(\gamma/m) \end{bmatrix}  \tag{19.8} \\
\\
\cfrac{d}{dt}\bm x(t) &= \begin{bmatrix} \frac{d}{dt} x_1(t) \\ \frac{d}{dt} x_2(t)\end{bmatrix} \\
&= L \bm x(t) \tag{19.9} \\
\\
\bm x(t) &= c_1 e^{\lambda_1 t} \bm v_1 + \cdots, c_D e^{\lambda_D t} \bm v_D \\
&= \sum_{d=1}^D c_d e^{\lambda_d t} \bm v_d \tag{21.6} \\
\\
\bm x(0) &= c_1 \bm v_1 + \cdots + c_D \bm v_D \\
&= \begin{bmatrix} \bm v_1 & \cdots & \bm v_D \end{bmatrix} \begin{bmatrix} c_1 \\ \vdots \\ c_D \end{bmatrix} \\
&= P \begin{bmatrix} c_1 \\ \vdots \\ c_D\end{bmatrix} \tag{21.7}
\end{align*}
$$

テキスト p.173 式（19.8）（19.9）、p.189 式（21.6）（21.7）を引用

複雑さに直面して、変な汗が吹き出てきました。。。

計算の流れに整理し直すと、おそらく次のようになるでしょう。

  • 式（19.9）は時間発展を連立微分方程式の形で表しています。

  • 連立微分方程式の $${L}$$ は式（19.8）の要素を持つ係数行列です。

  • $${L}$$ の固有ベクトルが式（21.6）（21.7）の $${\bm v_d}$$ となります。

  • 式（21.7）の 初期状態 $${\bm x(0)}$$ を係数 $${c_d}$$ について解き、式（21.6）で時間発展方程式の解 $${\bm x(t)}$$ を得ます。

ということで、時間発展方程式の求解と、図 21.1「減衰と振動は、固有値から理解可能」の描画を実行します！

初期値はテキスト掲載の次の値を使います。

  • ピストンの質量 $${m=1}$$

  • バネ定数 $${k=3}$$

  • 摩擦係数 $${\gamma=1/2}$$

  • 時間 $${t=0}$$ のときの座標 $${x(0)=[1, 0]}$$

### p.191 図21.1 減衰と振動は、固有値から理解可能（図19.5と同じ図）

## 設定と準備
# パラメータの設定 m: ピストンの質量、k: バネ定数、gamma: 摩擦係数 p.174, p.190
m, k, gamma = 1, 3, 1/2
# ピストンの自然長からのずれの座標xの初期値 x(0)
x0 = np.array([1, 0])
# 時間発展を与える　D x D 行列 L の設定　※式（19.8）・（21.13）関連
L = np.array([[0, 1],            
              [-k/m, -gamma/m]])
print('L:')
print(L.round(3))

## 行列Lの固有値λ, 固有ベクトルvの算出
# 固有値・固有ベクトルの算出
λ, v = np.linalg.eig(L)
# 固有値の表示 ※式（21.14）関連 ※jが虚部
print('λ:')
print(λ.round(3))
print('v:')
# 固有ベクトルの表示 ※式（21.15）関連
print(v.round(3))

## 係数cの算出：numpyを用いて式（21.7）の連立方程式をcについて解く
c = np.linalg.solve(v, x0)
print('c:')
print(c)

## 時間発展方程式の解x(t)の算出
# 時間tの値の設定
times = np.linspace(0, 20, 1001)
# 時間発展方程式の解 x(t) の算出 ※式（21.6）
xt = np.real([sum([c[d] * np.exp(λ[d] * t) * v[:, d] for d in range(len(c))])
              for t in times]
)

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 時間発展方程式の解x1(t)の描画
ax.plot(times, xt[:, 0], label='$x_1(t)=x(t)$')
# 時間発展方程式の解x2(t)の描画
ax.plot(times, xt[:, 1], label='$x_2(t)=v(t)$')
# 修飾
ax.set(title='図 21.1 減衰と振動は、固有値から理解可能', xticks=range(0, 21, 5))
ax.set_xlabel('$t$', fontsize=12)
ax.set_ylabel('$x(t)$', fontsize=12)
ax.legend()
plt.show()

【実行結果】

行列 $${L}$$ の固有値 $${\lambda}$$（式 21.14 に相当）と 固有ベクトルを列にもつ行列 $${v}$$（式 21.15 に相当）、および係数 $${c}$$ は「j」のつく「虚数」を含む複素数になっています。

テキストによるとチャートの振動と減衰は次のように説明できるようです。

  • 固有値 $${\lambda}$$ が複素数であり、オイラーの公式から $${\sin, \cos}$$ が出てきて、振動の振る舞いに寄与しています。

  • 固有値 $${\lambda}$$ の実部 $${1/4=0.25 < 0}$$ が減衰に寄与しています。

■ コードの補足説明

  • $${L}$$ の固有ベクトルは「行列Lの固有値λ, 固有ベクトルvの算出」の箇所で変数 v に設定しています。

  • 式（21.7）の 係数 $${c_d}$$ の求解は「係数cの算出：numpyを用いて式（21.7）の連立方程式をcについて解く」の箇所で、numpy.linalg の solve （連立方程式の求解）を使いました。
    第１引数が固有ベクトル v です。式（21.7）の$${P}$$ に当たります。
    第２引数が初期条件 $${x(0)=[1, 0]}$$ です。
    固有ベクトルと変数 $${c}$$ の線形結合が 初期条件 $${x(0)}$$ となるように、変数 $${c}$$ の解を求めます。

  • 式（21.6）の時間発展方程式の解は「時間発展方程式の解x(t)の算出」の箇所で式（21.6）のとおり「xt = ～」を作っています。

◆ ◆ ◆

続いて減衰しない例として、テキストの図 21.2「減衰しない場合は固有値の実部がゼロ」を実装します。
係数行列 $${L'=\begin{bmatrix} 0 & 1 \\ -3 & 0 \end{bmatrix}}$$ を適用します。

### p.193 図21.2 減衰しない場合は固有値の実部がゼロ

## 設定と準備 図21.1と同じパラメータは省略
# 時間発展を与える　D x D 行列 L の設定　※式（21.20）関連
L_dash = np.array([[0, 1],            
                   [-3, 0]])
print("L':")
print(L_dash.round(3))

## 行列Lの固有値λ, 固有ベクトルvの算出
# 固有値・固有ベクトルの算出
λ, v = np.linalg.eig(L_dash)
# 固有値の表示 ※式（21.14）関連 ★実部がゼロ！
print('λ:')
print(λ.round(3))
print('v:')
# 固有ベクトルの表示 ※式（21.15）関連
print(v.round(3))

## 係数cの算出：numpyを用いて式（21.7）の連立方程式をcについて解く
c = np.linalg.solve(v, x0)
print('c:')
print(c)

## 時間発展方程式の解x(t)の算出
# 時間tの値の設定
times = np.linspace(0, 20, 1001)
# 時間発展方程式の解 x(t) の算出 ※式（21.6）
xt = np.real([sum([c[d] * np.exp(λ[d] * t) * v[:, d] for d in range(len(c))])
              for t in times]
)

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 時間発展方程式の解x1(t)の描画
ax.plot(times, xt[:, 0], label='$x_1(t)=x(t)$')
# 時間発展方程式の解x2(t)の描画
ax.plot(times, xt[:, 1], label='$x_2(t)=v(t)$')
# 修飾
ax.set(title='図 21.2 減衰しない場合は固有値の実部がゼロ', xticks=range(0, 21, 5))
ax.set_xlabel('$t$', fontsize=12)
ax.set_ylabel('$x(t)$', fontsize=12)
ax.legend(bbox_to_anchor=(1, 1))
plt.show()

【実行結果】
確かに固有値 $${\lambda}$$ の実部はゼロになっていて、チャートの曲線は減衰していません。

◆ ◆ ◆

③ 数値的に解く方法
こちらの解法は、数値計算手法である「オイラー法」と「ホイン法」を使います。
こちらも数式を手がかりにしてコードに落とし込みます。
上の式がオイラー法、下の式がホイン法です。

$$
\begin{align*}
\bm x(t + \Delta t) &= \bm x(t) + L \bm x(t) \Delta t \tag{21.27 }\\
\bm x(t + \Delta t) &= \bm x(t) + \cfrac{\Delta t}{2} \left( L \bm x(t) + L (\bm x(t) + L \bm x(t) \Delta t) \right) \\
&= \bm x(t) + L \bm x(t) \Delta t + L^2 \bm x(t) \cfrac{(\Delta t)^2}{2} \tag{21.28}  \\
\end{align*}
$$

テキスト p.195 式（21.27）（21.28）を引用

では図 21.3「手法を変えると大幅に精度が変わる」の実装に移ります！
最初に 左側の [1] のチャートです。
オイラー法、ホイン法、固有値展開（固有値を使って解く方法）の結果を比べます！
時間間隔 $${\Delta t}$$ は $${0.1}$$ としています。

### p.196 図21.3 手法を変えると大幅に精度が変わる

## 設定と準備
# パラメータの設定 m: ピストンの質量、k: バネ定数、gamma: 摩擦係数 p.174, p.190
m, k, gamma = 1, 3, 1/2
# ピストンの自然長からのずれの座標xの初期値 x(0)
x0 = np.array([1, 0])
# 時間発展を与える　D x D 行列 L の設定　※式（19.8）・（21.13）関連
L = np.array([[0, 1],            
              [-k/m, -gamma/m]])

## 固有値展開（時間発展方程式）のパラメータの算出
# 行列Lの固有値λ, 固有ベクトルvの算出
λ, v = np.linalg.eig(L)
# 係数cの算出：numpyを用いて式（21.7）の連立方程式をcについて解く
c = np.linalg.solve(v, x0)

## 時間発展方程式の解・近似解x(t)の算出
# 時間tの値の設定
delta_t = 0.1
times = np.arange(0, 20.1, delta_t)
# 固有値展開（時間発展方程式） ※式（21.6）
xt_eig = np.real([sum([c[d] * np.exp(λ[d] * t) * v[:, d] for d in range(len(c))])
                 for t in times])
# オイラー法 ※式（21.27）
xt_euler = [x0]
for t in range(1, len(times)):
    xt_euler.append(xt_euler[t-1] + (L @ xt_euler[t-1]) * delta_t)
xt_euler = np.array(xt_euler)
# ホイン法 ※式（21.28）
xt_huen = [x0]
for t in range(1, len(times)):
    xt_huen.append(
        xt_huen[t-1] + L @ xt_huen[t-1] * delta_t
        + L @ L @ xt_huen[t-1] * delta_t**2 / 2)
xt_huen = np.array(xt_huen)

## [1] の描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# ホイン法による解の近似解 x1(t) の描画
ax.plot(times, xt_huen[:, 0], color='tab:blue', label='ホイン法')
# オイラー法による解の近似解 x1(t) の描画
ax.plot(times, xt_euler[:, 0], color='tab:green', label='オイラー法')
# 固有値展開（時間発展方程式）の解 x1(t) の描画
ax.plot(times, xt_eig[:, 0], ls='--', color='tab:orange', label='固有値展開')
# 修飾
ax.set(title='図 21.3 [1] 手法を変えると大幅に精度が変わる', xticks=range(0, 21, 5),
       yticks=np.arange(-0.5, 1.1, 0.5))
ax.set_xlabel('$t$', fontsize=12)
ax.set_ylabel('$x_1(t)$', fontsize=12)
ax.legend()
plt.show()

【実行結果】

テキストのとおり、固有値展開とホイン法は重なっています。
ホイン法の方が精度が高いのでしょう。

■ コードの補足説明
オイラー法とホイン法の実装は「オイラー法 ※式（21.27）」、「ホイン法 ※式（21.28）」の箇所で、式（21.27）と（21.28）の通りに書いています。

◆ ◆ ◆

続いて 右側の [2] のチャートです。
時間間隔 $${\Delta t}$$ を $${[0.1, 0.01, 0.001, 0.0001]}$$ にして、固有値展開の方法とオイラー法・ホイン法との誤差を計算し、チャートで可視化します。
誤差計算部分は関数「compare_methods」で特定の時間間隔 delta_t の誤差を計算するようにしています。

## [2] の描画処理

# 固有値展開との誤差を算出する関数の定義
def compare_methods(L, λ, c, x0, delta_t, end_t=0.1):
    # 時間tの値の設定
    times = np.arange(0, end_t + delta_t, delta_t)
    # 固有値展開（時間発展方程式） ※式（21.6）
    xt_eig = np.real([sum([c[d] * np.exp(λ[d] * t) * v[:, d]
    for d in range(len(c))])
                     for t in times])
    # オイラー法 ※式（21.27）
    xt_euler = [x0]
    for t in range(1, len(times)):
        xt_euler.append(xt_euler[t-1] + (L @ xt_euler[t-1]) * delta_t)
    # ホイン法 ※式（21.28）
    xt_huen = [x0]
    for t in range(1, len(times)):
        xt_huen.append(
            xt_huen[t-1] + L @ xt_huen[t-1] * delta_t
            + L @ L @ xt_huen[t-1] * delta_t**2 / 2)
    # 戻り値：オイラー法の誤差、ホイン法の誤差
    return abs(xt_eig[-1] - xt_euler[-1]), abs(xt_eig[-1] - xt_huen[-1])

# オイラー法とホイン法の誤差の算出
# 算出するΔtのリスト
delta_ts = [10**-i for i in range(1, 5)]
# オイラー法とホイン法の誤差を格納するリストの初期化
euler_diffs, huen_diffs = [], []
# Δtごとにオイラー法とホイン法の誤差を算出
for delta_t in delta_ts:
    # オイラー法とホイン法の誤差を算出
    euler_diff, huen_diff = compare_methods(L, λ, c, x0, delta_t, end_t=0.1)
    # オイラー法とホイン法の誤差を格納
    euler_diffs.append(euler_diff)
    huen_diffs.append(huen_diff)

# 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# オイラー法の誤差の描画
ax.plot(delta_ts, np.array(euler_diffs)[:, 0], '-s', color='tab:green',
        label='オイラー法')
ax.text(x=10**-3.7, y=10**-5, s='$\propto \Delta t$', fontsize=12,
        color='tab:green')
# ホイン法の誤差の描画
ax.plot(delta_ts, np.array(huen_diffs)[:, 0], '-s', color='tab:blue',
        label='ホイン法')
ax.text(x=10**-1.5, y=10**-5, s='$\propto (\Delta t)^2$', fontsize=12,
        color='tab:blue')
# 修飾
ax.set(xscale='log', yscale='log', yticks=[10**-i for i in range(2, 11, 2)],
       title='図 21.3 [2] 手法を変えると大幅に精度が変わる')
ax.set_xlabel('$\\Delta\ t$', fontsize=12)
ax.set_ylabel('固有値展開からのずれ（誤差）', fontsize=12)
ax.legend()
plt.show()

【実行結果】

対数表示にすると直線的に見えます！
ホイン法の方が固有値展開との誤差が小さいようです。
テキストは「オイラー法は１次のオーダー、ホイン法は２次のオーダー」と記載しています。

途中何度も冷や汗をかきましたが、なんとか図の描画にたどり着けました！

３．オルンシュタイン・ウーレンベック過程

① 概要
水面にとても小さな粒子をそっと落とし、粒子が時間とともに動く様子を表現する「関数の時間発展」に移ります。
第22話「関数の時間変化を考える」でフォッカー・プランク方程式が紹介され、第23話「偏微分方程式を解く」でフォッカー・プランク方程式の１変数系を考える「オルンシュタイン・ウーレンベック過程」が登場します。
テキストおよび Wikipedia を参照してイメージを掴んでください！

水の中で眠る会社員たちのイラスト：「いらすとや」さんより

ということで、オルンシュタイン・ウーレンベック過程の数式（23.1）をお借りします。

$$
\cfrac{\partial}{\partial t} p(x,t) = \gamma \cfrac{\partial}{\partial x}\ (x p(x,t)) + \cfrac{D}{2} \cfrac{\partial^2}{\partial x^2}\ p(x,t) \tag{23.1}
$$

テキスト p.210 式（23.1）を引用

時刻 $${t}$$ で粒子が $${x}$$ の位置にいる確率と関係する関数として関数 $${p(x, t)}$$ が導入されています。
そして、$${\partial}$$ の記号は偏微分です！
偏微分方程式を解くにはどうしたらいいの？ということで、テキストの「基底関数を使う方法」と「空間を細かく分割する方法」にトライします。

◆ ◆ ◆

② 基底関数を使う方法
基底関数に「エルミート関数」を用いる方法で解きます！
エルミート関数の数式（23.4）をお借りします。

$$
\psi_n(x) = (2^n n! \sqrt{\pi})^{-1/2} e^{-x^2/2} H_n(x)  \tag{23.4}
$$

テキスト p.210 式（23.4）を引用

第２章で登場したエルミート多項式 $${H_n(x)}$$が使われます！！

さっそく図 23.2「エルミート多項式ではなく、エルミート関数」を描画しましょう。

### p.210 図23.2 エルミート多項式ではなく、エルミート関数

## 設定と準備
# x軸の値
xs = np.linspace(-5, 5, 1001)
# 次数のリスト
ns = list(range(3)) # [0, 1, 2] 
# エルミート関数の定義 ※式（23.4）
def psi(n, x):
    return (2**n * factorial(n, exact=True) * np.sqrt((np.pi)))**(-1/2) \
           * np.exp(-x**2 / 2) * hermite(n)(x)

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 次数ごとにエルミート多項式の算出とチャート描画を繰り返し処理
for n in ns:
    # エルミート関数の算出と曲線の描画
    ax.plot(xs, psi(n, xs), label=f'$\psi_{n}(x)$')
# 修飾
ax.set(xlim=(-5, 5), ylim=(-0.7, 0.8), title='図 23.2 エルミート関数')
ax.set_xlabel('$x$', fontsize=12)
ax.legend();

【実行結果】

テキストによるとエルミート関数は正規直交性があるとのこと。

◆ ◆ ◆

続けて図 23.4「計算方法は違っても、固有関数は大体一致する」の左側 [1] のエルミート関数を使って求めたチャートを描画しましょう。
こちらは関数の解ではなく、「固有関数」を描画するものです。
テキストの式（23.4）,（23.8）,（23.8）,（23.11）,（23.12）,（23.17）を使って、「何とか形になるように」たどり着きました。

### p.218 図23.4 計算方法は違っても、固有関数はだいたい一致する
# [1] エルミート関数を使って求めたもの

## 設定と準備
# 時間発展方程式（23.1）のパラメータ γ, D
γ, D = 1, 1
# エルミート関数の最大次数
N = 99
# 格子点の刻み幅 Δx
delta_x = 0.05
# 格子点の算出：x軸の値
xs = np.arange(-5, 5+delta_x, delta_x)

## 関数定義
# エルミート関数の定義 ※式（23.4）
def psi(n, x):
    return (2**n * factorial(n, exact=True) * np.sqrt((np.pi)))**(-1/2) \
           * np.exp(-x**2 / 2) * hermite(n)(x)
# 記号αの定義 ※式（23.8）
α = lambda n: 1/2 * np.sqrt(n * (n - 1))

## 行列Lの算出 ※式（23.12）,（23.11）
L = np.zeros((N + 1, N + 1))
for n in range(N + 1):
    L[n, n] = γ / 2 + D/2 * (-n - 1/2)
    if n >= 2:
        L[n, n - 2] = -γ * α(n) + D/2 * α(n)
    if n <= N - 2:
        L[n, n + 2] = γ * α(n + 2) + D/2 * α(n + 2)

## 行列Lの固有値λ, 固有ベクトルvの算出
λ, v = np.linalg.eig(L)
sort_idx = λ.argsort()[::-1]  # 固有値の昇順にソートするためのインデックス
λ = np.real(λ[sort_idx])      # 実数部のみ取り出し
v = np.real(v[:, sort_idx])   # 実数部のみ取り出し

## 固有関数φの算出 ※式（23.17）
phi = np.zeros((N + 1, len(xs)))
for m in range(N + 1):
    # エルミート関数 ψm(x) の算出
    psi_m = psi(m, xs)
    # 式（23.17）右辺の sum(v_im ･ ψm(x)) の計算
    for i in range(N + 1):
        phi[i] += v[m, i] * psi_m

## 固有関数φの可視化
# 描画対象の固有関数のインデックスの設定
i = 1
# 描画領域の設定
fig, ax = plt.subplots()
# 固有関数φの描画
ax.plot(xs, phi[i])
# 修飾
ax.set(yticks=np.arange(-0.8, 0.9, 0.4),
       title=f'[1] エルミート関数を使って求めたもの\n'
             f'$\phi_{i}(x),\ \lambda_{i}={λ[i]:.4f}$')
ax.set_xlabel('$x$', fontsize=12);

【実行結果】

■ コアとなる数式とコード実装の補足説明
線形作用素 $${\mathcal{L}}$$ に対して、スカラ $${\lambda_i}$$ と対応する関数 $${\varphi_i(x)}$$ が存在して、以下の式を満たすとき、$${\varphi_i(x)}$$ を固有関数、$${\lambda_i}$$ を固有値と呼びます。コードは「行列 L の固有値 λ, 固有ベクトル v の算出」で実装しています。

$$
\mathcal{L} \varphi_i(x) = \lambda_i \varphi(x) \tag{23.16}
$$

テキスト p.213 式（23.16）を引用

無限個のエルミート関数 $${\psi_m(x)}$$ を考えて形式的に行列 $${L}$$ を作り、$${L}$$ の固有ベクトルを $${\bm v_i}$$ とすると、固有関数 $${\varphi_i(x)}$$ は以下の式となります。コードは「固有関数 φ の算出 ※式（23.17）」で実装しています。

$$
\varphi_i(x) = \sum_{m=0}^M v_{im} \ket{\psi_m} = \sum_{m=0}^M v_{im} \psi_m(x) \tag{23.17}
$$

テキスト p.213 式（23.17）を引用

そして何よりも、一番難儀したのが、係数行列 $${L}$$ です。
式（23.11）の係数を拾い出して、コードの「行列Lの算出 ※式（23.12）,（23.11）」の箇所に実装しました。

◆ ◆ ◆

③ 空間を細かく分割する方法
この記事の最難関に到達しました。
ざっくりこの方法は、位置 $${x}$$ を $${\Delta x}$$ 刻みで格子状にして、離散的な点々を追いかける感じです。
大筋を含め、ぜひともテキストで考え方や数式をご確認ください！

実装にあたってテキストよりピックアップした数式を引用いたします。
■ 格子点上での関数の値 $${{p(x_m, t)}}$$ を使ったベクトル

$$
\bm p(t) = \begin{bmatrix}p(x_0, t) \\ \vdots \\ p(x_M, t)\end{bmatrix} \tag{23.28}
$$

テキスト p.216 式（23.28）を引用

■ 一階微分の中心差分

$$
\cfrac{d}{d x} p(x_m, t) \simeq \cfrac{p(x_m + \Delta x, t) - p(x_m - \Delta x, t)}{2 \Delta x} \tag{23.31}
$$

テキスト p.216 式（23.31）を引用

■ 二階微分の中心差分
こちらのWebサイトの記事を参考にさせていただきました。
ありがとうございます！

$$
\cfrac{d^2}{d x^2} p(x_m, t) \simeq \cfrac{p(x_m + \Delta x, t) -2p(x_m, t) + p(x_m - \Delta x, t)}{(\Delta x)^2}
$$

https://note.com/enjoy_simulation/n/n3b09adff43b3

■ 式（23.28）のベクトル $${\bm p(t)}$$ に対する時間発展方程式

$$
\cfrac{d}{dt}\ \bm p(t) = \widetilde{L} \bm p(t) \tag{23.37}
$$

テキスト p.217 式（23.37）を引用

数式の迷宮を彷徨い、いろいろモヤモヤしながらの実装となりました。

道に迷った人のイラスト（男性）：「いらすとや」さんより

そして何と言っても、テキストの図の固有値 $${\lambda_1}$$ と異なる値 になってしまいました。。。
ひとまず、実行結果を見てやってください。

### p.218 図23.4 計算方法は違っても、固有関数はだいたい一致する
### [2] 固有関数の可視化：数値微分・刻み幅0.05の格子を使って求めたもの
# ★★★ ＮＧ：書籍と固有値λの値が相違している（惜しいけど） ★★★

## 設定と準備
# 時間発展方程式（23.1）のパラメータ γ, D
γ, D = 1, 1
# 格子点の刻み幅 Δx
delta_x = 0.05
# 格子点の算出：x軸の値
xs = np.arange(-5, 5+delta_x, delta_x)
# 格子点の数
M = len(xs)

## 行列L_tildeの算出 ※式（23.31）
# 二階微分の数値微分の参考サイト: https://note.com/enjoy_simulation/n/n3b09adff43b3
L_tilde = np.zeros((M, M))
for m in range(M):  # 拡散項のx_m部分
    L_tilde[m, m] = D/2 * -2 / delta_x**2
    if m >= 1:      # 第1項：ドリフト項のx_m-1部分、第2項：拡散項のx_m-1部分
        L_tilde[m, m - 1] = - γ * xs[m - 1] / (2*delta_x) + D/2 / delta_x**2
    if m <= M - 2:  # 第1項：ドリフト項のx_m+1部分、第2項：拡散項のx_m+1部分
        L_tilde[m, m + 1] = + γ * xs[m + 1] / (2*delta_x) + D/2 / delta_x**2

## 行列L_tildeの固有値λ, 固有ベクトルv=固有関数の算出
λ, v = np.linalg.eig(L_tilde)
sort_idx = λ.argsort()[::-1]  # 固有値の昇順にソートするためのインデックス
λ = λ[sort_idx]               # 固有値λ
v = v[:, sort_idx]            # 固有ベクトルv

## 固有ベクトルv=固有関数の可視化
# 描画対象の固有関数のインデックスの設定
i = 1
# 描画領域の設定
fig, ax = plt.subplots()
# 固有ベクトルvの描画
ax.plot(xs, v[:, i], color='tab:orange')
# 修飾
ax.set(yticks=np.arange(-0.15, 0.151, 0.05),
       title=f'[2] 刻み幅{delta_x}の格子を使って求めたもの\n'
             f'$v_{i}(x),\ \lambda_{i}={λ[i]:.4f}$')
ax.set_xlabel('$x$', fontsize=12);

【実行結果】

テキストは固有値 $${\lambda_1=-1.0050}$$ ですが、こちらの実装では $${-1.0000}$$ となっています。
ご注意ください。
一応、[1] エルミート関数の固有値とこちらの実装の固有値が同じになっているので、何か御縁があるかもしれません。

■ コードの補足説明
何と言っても 係数行列 $${\widetilde{L}}$$（L_tilde）の計算が超絶難しかったです。
記述箇所は「行列L_tildeの算出 ※式（23.31）」です。
テキストに具体的な計算方法が見当たらなかったので、めっちゃ試行錯誤を繰り返しました。
現在のコードでは、式（23.1）のオルンシュタイン・ウーレンベック過程の方程式に、１階微分の中心差分と２階微分の中心差分を当てはめた結果となっています。

冷や汗は止みません。偏微分で偏頭痛になりそうです。。。

今回の写経は以上です。

冷や汗のイラスト：「いらすとや」さんより

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

いいなと思ったら応援しよう！

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します！

この記事が参加している募集

#

この春やりたいこと

2.0万件