「線形代数の半歩先」をPythonで写経～５章ファン・デル・ポル方程式、スナップショット・ペア、コルモゴロフの後退方程式

ネイピア DS

2025年4月17日 11:51

第５部「ならべた数のさらなる発展」

書籍の著者　大久保潤先生

この記事は、書籍「線形代数の半歩先」の第５部「ならべた数のさらなる発展」に掲載の「非線形系における線形性」に関する Python写経活動のドキュメンタリーです。

第５部に突入します。
第４部に引き続き「時間発展方程式」を取り扱います。
この記事では、第25話の ファン・デル・ポル方程式 に取り組みます！
数学素人なので、どうぞお手柔らかにお願いいたします。

◆◆ ◆

実は第５部を始めて直ぐに、混迷の沼に足を取られて身動きできない、みたいな感じになり、太刀打ちできない状況に陥りました…

「難しすぎて頭に入ってこないんです」

でもせっかく取り組み始めた書籍をリタイアしたくないです。
そこで私が採った方策は…

書籍の読み込みを一時停止
ChatGPTをフル活用する学習方法に切り替え
ChatGPT のサポートで書籍の図表を Python 実装
Python コードを読み解くことで書籍の数式にリーチ

ChatGPT活用型学習に切り替えたことによって、前章までの書きっぷりから大幅に変わっていますが、引き続きご愛顧のほどよろしくお願い申し上げます。

では書籍とChatGPTを開いて線形代数の旅に出発です🚀

はじめに

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

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

第５部ならべた数のさらなる発展

前説

第25話「時間発展データのために」では、時間が飛び飛びの時系列データを取り扱います。
飛び飛びとは、時間間隔 $${\Delta t_{\text{obs}}}$$ でのみデータを取得できることです。
正の整数 $${k}$$ で飛び飛びの離散時刻 $${t_k=t_0 + k \Delta t_{\text{obs}}}$$ を扱うのです。

また、スナップショット・ペアと呼ばれるデータを扱います。
スナップショット・ペアとはイメージ的には、状態 $${\bm x_k}$$ と、対応する $${\Delta t_{\text{obs}}}$$ 時間後の状態 $${\bm x_{k+1}}$$ をペアにして、$${(\bm x(t_0), \bm x(t_1)), (\bm x(t_1), \bm x(t_2)), \ldots}$$ のようなデータセットになります。

そして強敵「ファン・デル・ポル方程式」に取り組みます。
この記事の最終目標は、図 25.1「非線形な例題の軌跡をたどる」の描画です。
ただし、ファン・デル・ポル方程式の意味が分からないし、ファン・デル・ポル方程式から飛び飛びの時系列データを生成する方法も分かりません。
テキストから読み取れないのです。。。

そこで、第５部からはChatGPTに全面依拠して進めます！
ひとまずテキストから目をそむけます。。。
というわけで、ChatGPTの回答が正しいか、間違っているのかは判断できません！

ではでは生成AI劇場のはじまりはじまり～

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

### インポート

# 数値計算
import numpy as np

# 数学演算用の個別関数
from scipy.integrate import solve_ivp   # 常微分方程式の求解

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

１．ChatGPTにファン・デル・ポル方程式を尋ねる

① はじめの一歩
ChatGPTに数式の画像を投げるために、Markdown で２つの数式を書きます。

■ スナップショット・ペア（時間発展前後のペア）

$$
{(\bm x_n, \bm y_n) | n=1, 2, \ldots, N} \tag{25.1}
$$

テキスト p.229 式（25.1）を引用

■ データを生み出した方程式

$$
\begin{align*}
\cfrac{d}{dt} x_1(t) &= x_2(t) \tag{25.2} \\
\cfrac{d}{dt} x_2(t) &= \left( \epsilon x_2(t)(1 - x_1(t)^2) -x_1(t) \right) \tag{25.3} \\
\end{align*}
$$

テキスト p.230 式（25.2）（25.3）を引用

続けてChatGPTに投げます。
モデルは ChatGPT 4o です。

◆ ◆ ◆

【ご案内】ChatGPT画面について
ChatGPT画面で行ったChatGPTとのやりとりは、ChatGPT画面を画像貼り付けしています。
右寄せがユーザー（つまり私）の指示であり、グレー吹き出しの指示文章とグレー吹き出しの上にある画像等のアップロード表示が該当します。
左寄せがChatGPTの回答です。
やりとりの一部を切り取って掲示していますので、ご留意ください。

◆ ◆ ◆

方程式の画面ショットをアップロードして「これ何ですか？」から始めました！

② Pythonコードの制作
しめしめ、ChatGPTが反応しました！
Pythonコードを書きましょうか？と尋ねてきたので、テキスト p.230 に掲載の「$${\epsilon=1.0}$$ としてシミュレーションしました」を投げます。

んー、テキストの図 25.1 とは程遠い出来です。。。。
図 25.1 の見た目で、時刻 $${t=0}$$ の初期条件が $${x_1(0)=0.2, x_2(0)=0.8}$$ くらいなので、伝えます。

まだまだ図 25.1 には程遠いです。。。
テキストの図は時間の範囲が $${[0, 3]}$$ なので、これを伝えます。

ついに図 25.1 によく似たチャートがでてきました！
やったね！
Pythonコードを所望します。

③ 知識の拡張
コードをざっと見したところ、知らない関数 solve_ivp が含まれていたので質問します。

なんと！常微分方程式を解く関数ではないですか！
ファン・デル・ポル方程式の常微分方程式を解けるのです！
Python も scipy も ChatGPT も全員すごいです！！！

solve_ivp の説明一式は、記事の後半で紹介します。

ChatGPTのコードで図 25.1 が書けそうなことが分かって、小躍り気分を満喫するのもつかの間。。。
そうそう、当初の大問題「ファン・デル・ポル方程式」の意味がまだ理解できていないことを思い出して、おもむろに聞きます。

ふむふむ、なんとなく分かった気になります。
ファン・デル・ポル方程式の説明一式も、記事の後半で紹介します。

ChatGPTのおかげで全くの素人がちょっぴりファン・デル・ポル方程式に近づけました。
ありがとうGPT。

（注意書き）
テキストは「ファン・デル・ポル方程式」と呼んでいますが、ChatGPTは「ファン・デル・ポール方程式」と呼んでいます。
伸ばす棒が入っています。
この記事では、両方の書き方が併存しますので、ご了承ください。

２．図 25.1 のファン・デル・ポル方程式の描画

ChatGPT が書いた図 25.1 的なPythonコードを手がかりにして、テキストの形式に手作業で整えたものがこちらです！

### p.230 図25.1 非線形な例題の軌跡をたどる

## パラメータ設定
# 非線形性の強さ 𝜖
epsilon = 1.0
# 時間ステップ Δt_obs
dt_obs = 0.1
# 時間範囲
t_span = (0, 3)
# 結果を格納する時刻の配列
t_eval = np.arange(*t_span, dt_obs)

## ファン・デル・ポール方程式の関数定義
def van_der_pol(t, y):
    x1, x2 = y                               # y = [x1, x2]
    dx1dt = x2                               # 式（25.2）
    dx2dt = epsilon * x2 * (1 - x1**2) - x1  # 式（25.3）
    return [dx1dt, dx2dt]                    # 戻り値：x1, x2のtでの微分

## 初期条件の設定
y0 = [0.2, 0.9]

## 常微分方程式（ODE）の初期値問題を数値的に解く（ルンゲ・クッタ法） ※scipy利用
#  x1とx2の解がsol.yに格納される
sol = solve_ivp(van_der_pol,    # 微分方程式の関数
                t_span,         # 時間範囲（積分区間）
                y0,             # 初期値
                method='RK45',  # 数値解法：5次(4)の明示的ルンゲ・クッタ法
                t_eval=t_eval,  # 解を保存する時間
            )

## 可視化
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
# 色のリスト
cs = ['tab:red', 'lightgreen', 'violet', 'tab:blue','tab:orange']

# [1] 時間経過をプロット
# x1の折れ線グラフの描画 x1 = sol.y[0]
ax1.plot(sol.t, sol.y[0], '-o', color=cs[3], label='$x_1$')
# x2の折れ線グラフの描画 x2 = sol.y[1]
ax1.plot(sol.t, sol.y[1], '-o', color=cs[4], label='$x_2$')
# インデックス0, 25, 26のデータの色を変えて点と破線を描画
for i, t in enumerate([0, 25, 26]):
    ax1.plot(sol.t[t], sol.y[0][t], 'o', ms=8, color=cs[i])
    ax1.plot(sol.t[t], sol.y[1][t], 'o', ms=8, color=cs[i])
    ax1.plot([sol.t[t], sol.t[t]], sol.y[:, t], ls='--', color=cs[i])
# テキストの表示
ax1.text(x=0.8, y=-1.6, s='連続的な時間発展から\n離散的にデータを取得', fontsize=12)
# 修飾
ax1.set(xticks=(0, 1, 2, 3), title='[1] 時間経過をプロット')
ax1.set_xlabel('時間 $t$', fontsize=14) 
ax1.legend(loc='lower left', fontsize=14)

# [2] 平面座標にプロット
# x1, x2の散布図の描画
ax2.plot(sol.y[0], sol.y[1], 'o', color='navy')
# インデックス0, 25, 26のデータの色を変えて点を描画
for i, t in enumerate([0, 25, 26]):
    ax2.plot(sol.y[0][t], sol.y[1][t], 'o', ms=8, color=cs[i])
# 楕円の描画 matplotlib.paches.Ellipse
ellipse = Ellipse(xy=(sol.y[:, 26] + sol.y[:, 25]) / 2, 
                  height=0.2, width=0.45, angle=70, ls='--',
                  fc='none', ec='black')
ax2.add_patch(ellipse)
# テキストの表示
ax2.text(x=sol.y[0, 0], y=0.4, ha='center', s='最初の\nデータ', color=cs[0],
         fontsize=12)
ax2.text(x=0.68, y=-1.8, s='離散時間発展の前後の組が\nスナップショット・ペア',
         fontsize=12)

# 修飾
ax2.set(title='[2] 平面座標にプロット')
ax2.set_xlabel(r'$x_1$', fontsize=14)
ax2.set_ylabel(r'$x_2$', fontsize=14)

plt.show()

【実行結果】
図 25.1 に近いものを描くことができました！

■ コードの補足説明

① ファン・デル・ポル方程式の関数化
solve_ivp で微分方程式を解くにはまず、微分方程式を関数形式で記述する必要があります。
コードの「ファン・デル・ポール方程式の関数定義」の箇所に記述しています。
こちらはChatGPTのコードをそのまま使い、コメントを手作業で追加しました。
式（25.2）,（25.3）が実装されています。

② ファン・デル・ポル方程式の求解
solve_ivp でファン・デル・ポル方程式を求解する箇所は「常微分方程式（ODE）の初期値問題を数値的に解く（ルンゲ・クッタ法） ※scipy利用」です。
scipy.integrate（積分）の solve_ivp 関数です。
微分方程式を積分して原始関数を導出するようです。
こちらもChatGPTのコードをそのまま使い、コメントを手作業で追加しました。

③ チャート描画
描画はほぼ手作業で書き換えました。
がんばって楕円も追加しました！

ChatGPTによる解説コーナー

１．ChatGPTに習う「ファン・デル・ポル方程式」

ではChatGPTが教えてくれたファン・デル・ポル方程式のイロハに進みましょう。
生成AIのリアリティを体感する目的で、ChatGPTの回答をほぼ転記しています。

🌀 ファン・デル・ポール方程式（Van der Pol Equation）とは？

ファン・デル・ポール方程式は、非線形な振動現象を表す微分方程式で、物理・電子工学・生体システムなどさまざまな分野で使われる基本的なモデルです。

✅ 標準形（2階常微分方程式）

$$
\frac{d^2x}{dt^2} - \mu (1 - x^2) \frac{dx}{dt} + x = 0
$$

$${x(t)}$$：時間 $${t}$$ における変位
$${\mu}$$：非線形性と減衰の強さを表す定数（$${\mu > 0}$$）

🔁 一階の連立方程式への変形

新たに $${x_1 = x}$$, $${x_2 = \frac{dx}{dt}}$$ と置くことで、一階の連立常微分方程式に変形できます：

$$
\begin{align*}
\frac{dx_1}{dt} &= x_2 \\
\frac{dx_2}{dt} &= \mu \cdot x_2 (1 - x_1^2) - x_1 \\
\end{align*}
$$

この形は、数値的に解きやすく、シミュレーションにもよく使われます。

🎯 特徴

$$
\begin{array}{l:l}
特徴 & 内容 \\
\hline
\\
非線形性 & 減衰項が x に依存し、非線形のふるまいをする \\
\\
自励振動 & 外力がなくても、内部の非線形性により \\
& 振動が続く（リミットサイクル） \\
\\
非対称な振動 & 小さい振幅では増幅し、大きい振幅では \\
& 減衰されるため非対称になる \\
\\
応用範囲が広い & 電子回路、生体リズム、神経モデル \\
& などで使われる \\
\end{array}
$$

📈 振動のイメージ

$${\mu = 0}$$ のとき：単なる線形の調和振動（ばね・マス系と同じ）
$${\mu > 0}$$ のとき：非線形で安定な周期運動（リミットサイクル）を形成
$${\mu \gg 1}$$ のとき：鋭いスパイクのような波形（緩やかな上昇＋急激な減衰）

🧪 応用分野

❤️ 心臓の拍動モデル（Van der Pol 発振器）
🔌 三極真空管などの電子回路モデル
🧠 神経細胞のスパイク信号の近似
📊 非線形力学系やカオス理論の研究

💡 なぜ重要か？

最も基本的な「非線形・自己振動モデル」の1つとして知られている
リミットサイクルという現象を持つ典型例であり、教材にも頻繁に登場
シンプルながらも複雑で興味深い振る舞いを示す

🧰 参考：Pythonでの数値解法例（scipy）

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

mu = 1.0
def vdp(t, y):
    x1, x2 = y
    dx1 = x2
    dx2 = mu * x2 * (1 - x1**2) - x1
    return [dx1, dx2]

sol = solve_ivp(vdp, t_span=(0, 20), y0=[0.2, 0.8], t_eval=np.linspace(0, 20, 1000))

plt.plot(sol.y[0], sol.y[1])
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.title("Van der Pol 方程式の位相空間軌道")
plt.grid(True)
plt.show()

このように、ファン・デル・ポール方程式は数理モデルとして非常に重要で、非線形・周期・生物・電子などのキーワードに関心がある方にはとてもおすすめの題材です。

ゆるーい解説が沁みます。
なお、式中の $${\mu}$$ はテキストの $${\epsilon}$$ のことです。

２．ChatGPTに習う「solve_ivp」

続いて常微分方程式を解く関数 solve_ivp の解説です。

🧠 solve_ivp とは？

solve_ivp は、SciPy の 常微分方程式（ODE）の初期値問題を数値的に解く関数です。

✅ 基本の考え方

solve_ivp は "solve Initial Value Problem" の略で、連立常微分方程式を、初期条件から指定した時間範囲で 数値的に解いてくれる関数です。

📘 使用例（基本構文）

from scipy.integrate import solve_ivp

def f(t, y):
    return -y + t

t_span = (0, 5)
y0 = [1.0]

sol = solve_ivp(f, t_span, y0)

🔢 解いている方程式の意味

上記のコードは、次の微分方程式の初期値問題を解いています：

$$
\frac{dy}{dt} = -y + t, \quad y(0) = 1.0
$$

このように、`solve_ivp` は任意の1階微分方程式や連立微分方程式を初期条件付きで解くための関数です。

📦 主な引数の説明

$$
\begin{array}{l:l}
引数 & 意味 \\
\hline
\text{fun} & 微分方程式の関数 f(t, y) \\
\text{t\_span} & 時間範囲（例：(0, 5)） \\
\text{y0} & 初期値（y(t_0)） \\
\text{method} & 数値解法（例：\text{RK45}, \text{RK23}, \text{BDF}） \\
\text{t\_eval} & 結果を評価する時刻のリスト（任意） \\
\end{array}
$$

✅ 戻り値（`sol` オブジェクト）

solve_ivp の戻り値は OdeResult 型のオブジェクトです。主な属性は以下の通りです：

$$
\begin{array}{l:l}
属性名 & 内容 \\
\hline
\text{sol.t} & 出力された時刻の配列 \\
\text{sol.y} & 各時刻における状態（変数）の配列 \\
\text{sol.status} & 解のステータス（0 なら成功） \\
\text{sol.success} & 解が成功したかどうか（\text{True} / \text{False}） \\
\text{sol.message} & 解法のメッセージ（成功/失敗の説明） \\
\end{array}
$$

🧠 なぜ便利なのか？

自動で時間刻み（ステップサイズ）を調整してくれる（adaptive step size）
非線形や stiff な方程式にも対応可能
明示的・暗黙的な数値解法を選べる
複雑な連立方程式でも直感的に書ける
イベント検出や精度制御などの高機能オプションもあり

📌 ファン・デル・ポール方程式への応用

ファン・デル・ポール方程式のような連立微分方程式：

$$
\begin{align*}
\frac{dx_1}{dt} &= x_2 \\
\frac{dx_2}{dt} &= \varepsilon x_2 (1 - x_1^2) - x_1
\end{align*}
$$

これも solve_ivp で簡単に解くことができます：

def van_der_pol(t, y):
    x1, x2 = y
    dx1 = x2
    dx2 = epsilon * x2 * (1 - x1**2) - x1
    return [dx1, dx2]

sol = solve_ivp(van_der_pol, t_span=(0, 3), y0=[0.2, 0.8])

✅ まとめ

solve_ivp は、SciPy 標準の常微分方程式ソルバとして非常に汎用的でパワフルです：

初期値問題（ODE）の標準解法
シンプルな1変数ODEから複雑な非線形連立ODEまで対応
数値解法・出力時刻・精度などを細かく設定可能

🧰 主な method オプション

$$
\begin{array}{l:l}
\text{method} & 説明 \\
\hline
\text{RK45} & 明示的 \text{Runge-Kutta} 法（4/5次、デフォルト）\\
\text{RK23} & より粗い精度だが高速な \text{Runge-Kutta} 法 \\
\text{BDF} & 暗黙的多段法（\text{stiff} 系に強い） \\
\text{LSODA} & \text{stiff} / 非\text{stiff} 自動切替（\text{FORTRAN}ベース） \\
\end{array}
$$

📎 その他便利な引数

t_eval を指定すると、任意のタイミングの値を出力できます
rtol, atol で数値精度の制御が可能です
events パラメータで条件付き停止やイベント検出も可能です

とても詳しく説明してくれました。

コルモゴロフの後退方程式

まだまだ足りない・・・

じつは、テキスト p.235 のコルモゴロフの後退方程式にも興味が湧いたので、欲を出して、ChatGPTに聞いてみました。
テキストを読む中で、なんとなくコルモゴロフの後退方程式が将来予測に関係していそうだと感じたからです。
そうです、私は 25.1 節のタイトル「データから、少し先の未来を予測する」に興味を持ち、時間発展データの将来予測をしたいのです。
しかし、テキストを読んでも将来予測実践の糸口が見つからず、藁をも掴む気持ちで、ChatGPTに頼ろうと思ったのです。

ChatGPTとのやり取りの図は、ChatGPT画面の一部の切り取りです。
後半で整理した情報を書きますので、まずはChatGPTとのやりとりの雰囲気をお楽しみくださいな。

なんとなくコルモゴロフの後退方程式で未来予測ができそうな気がして、続けます。

そしてChatGPTがPython実装を提案してきました！
やったね！

コードを見て目眩がしました。。。
「このコードが何やってるのか、全く想像ができない。。。」
思わずつぶやきました。

そしてChatGPTは続けます。
私に優しく語りかけてくれるのです。

素人に難しいことはNGです。
簡単な事例で糸口を掴みに行きましょう。
どうやらブラウン運動の簡単な例を出してくれるようです。

ここでChatGPTはブラウン運動に対するコルモゴロフの後退方程式の数値解のPythonコードを書いてくれました。
早速、自環境で実行すると、不思議なチャートが現れました。

「おいおい、君が書いたコードにおかしな点があるんかい！」と思わず突っ込みたくなりますｗ

どうやら今回のコードは良さげなようです。
ChatGPTが書いたコードを鵜呑みにせず、「チェックして」と自己チェックをお願いすることで、より適切なコードが得られることを実感しました。

コードは後ほど紹介いたします。

ただコードの内容が理解できなかったので、コードの意味・内容をChatGPTに聞き返しました。

この質問は効果的でした。
コードの目的や動きを「とてもわかり易く」まとめてくれました。
コードの意味がわからないときは、どんどんChatGPTに質問しましょう！

ブラウン運動に対するコルモゴロフの後退方程式の数値解

ChatGPTのPythonコードを紹介します。
既にインポート済のライブラリを手作業で削除しましたが、これ以外の部分はChatGPTが書いた内容そのままです。

# 追加インポート
from scipy.sparse import diags

# パラメータ設定（修正版）
sigma = 1.0
T = 1.0         # 予測時間
dt = 0.0005     # 時間刻みを小さくして安定性を確保
n_steps = int(T / dt)

x_min, x_max = -5, 5
N = 201
x = np.linspace(x_min, x_max, N)
dx = x[1] - x[0]

# 初期条件: [-1, 1] にいるときのみ 1、それ以外は 0
u = np.where(np.abs(x) <= 1.0, 1.0, 0.0)

# 拡散項: (sigma^2 / 2) * ∂²u/∂x²
main_diag = -2 * np.ones(N)
off_diag = np.ones(N - 1)
L = (sigma**2 / (2 * dx**2)) * diags(
    [off_diag, main_diag, off_diag], offsets=[-1, 0, 1], format='csr'
)

# オイラー法で後退方程式を時間発展（u_t = L u）
for _ in range(n_steps):
    u += dt * (L @ u)
    u[0] = u[1]     # Neumann境界条件（端の値は内側と同じ）
    u[-1] = u[-2]

# 可視化
plt.figure(figsize=(8, 5))
plt.plot(x, u, label='u(x, T)')
plt.xlabel('$x$（現在の位置）')
plt.ylabel('確率（将来に $[-1, 1]$ にいる確率）')
plt.title('コルモゴロフ後退方程式による予測（ブラウン運動）')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

【実行結果】

私が欲しかったのは「時系列データの将来時点の予測値」ですが、コルモゴロフ後退方程式によるブラウン運動の予測は「現在の位置を踏まえて、将来、位置 $${[-1, 1]}$$ にいる確率」の形式で提供されました、がっくし。

ま、将来予測の観点で、コルモゴロフの後退方程式をこれ以上深追いしなくてもいいことが分かっただけでも良しとしましょう。

将来予測の続きは、次回以降の記事で紹介いたします！

コードの解説

続いて先程のコードに関するChatGPTの解説を掲載いたします。

✅ コードの目的

今、位置 $${x}$$ にいるブラウン粒子が、将来 $${t = 1}$$ に区間 $${[-1, 1]}$$ にいる確率を、コルモゴロフの後退方程式を使って数値的に求めるものです。

💻コード全体の構成

1. 必要なライブラリをインポート

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags

1️⃣ パラメータ設定

sigma = 1.0         # 拡散係数（ノイズの強さ）
T = 1.0             # 予測する時間（t=1）
dt = 0.0005         # 時間刻み（小さいほど安定）
n_steps = int(T / dt)  # 時間ステップ数

$${\frac{\partial u}{\partial t} = \frac{\sigma^2}{2} \frac{\partial^2 u}{\partial x^2}}$$ という後退方程式を、オイラー法で時間 $${t}$$ を増やしながら解きます。

2️⃣ 空間グリッドの作成

x_min, x_max = -5, 5
N = 201
x = np.linspace(x_min, x_max, N)
dx = x[1] - x[0]

$${x}$$ の範囲：$${[-5, 5]}$$
$${N}$$ 点に分割 → 空間刻み幅 $${dx}$$

3️⃣ 初期条件の定義（観測量 f(x)）

# 将来 X_t が [-1, 1] に入っているかを表す関数
u = np.where(np.abs(x) <= 1.0, 1.0, 0.0)

$${u(x, 0) = f(x)}$$ を設定
これは「将来 $${X_t \in [-1, 1]}$$ なら 1、それ以外なら 0」という観測関数

4️⃣ 差分法による L の行列化

行列 $${L}$$ に $${\frac{\partial^2 u}{\partial x^2}}$$ を設定

main_diag = -2 * np.ones(N)
off_diag = np.ones(N - 1)
L = (sigma**2 / (2 * dx**2)) * diags(
    [off_diag, main_diag, off_diag], offsets=[-1, 0, 1], format='csr'
)

2階微分を3点中心差分で離散化
これが後退方程式に出てくる演算子 $${\mathcal{L} u}$$

5️⃣ オイラー法で時間発展（後退方向）

for _ in range(n_steps):
    u += dt * (L @ u)
    u[0] = u[1]     # Neumann境界条件の簡易対応
    u[-1] = u[-2]

差分方程式：$${u^{n+1} = u^n + \Delta t \cdot L u^n}$$
毎ステップ、Neumann 境界条件（端の変化なし）として、
- 端の値を内側の値と等しくする：$${u[0] = u[1], u[-1] = u[-2]}$$

6️⃣ 可視化（結果表示）

plt.figure(figsize=(8, 5))
plt.plot(x, u, label='u(x, T)', color='orange')
plt.xlabel('$x$（現在の位置）')
plt.ylabel('確率（将来に $[-1, 1]$ にいる確率）')
plt.title('コルモゴロフ後退方程式による予測（ブラウン運動）【修正版】')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

横軸：今の位置 $${x}$$
縦軸：1秒後に $${[-1, 1]}$$ にいる確率

✅ 結果の意味まとめ

$${x = 0}$$ にいるとき：確率 ≈ 0.7 → 一番中心だから
$${x = \pm 2}$$ にいるとき：確率 ≈ 0.2 → 遠くて入りにくい
$${x = \pm 4}$$ にいるとき：確率 ≈ 0 → ほぼ到達不能

📘 補足：この計算のすごさ

通常なら何千回ものシミュレーションが必要な期待値を、偏微分方程式を1回解くだけで得られる！

📘 解説

🧠 コードの構造と意味

１．初期条件：観測関数 $${f(x)}$$

「将来 $${X_t}$$ が $${[-1,1]}$$ に入っていれば $${1}$$、それ以外は $${0}$$」

$$
f(x) =
\begin{cases}
1 & \text{if}\ x \in [-1, 1] \\
0 & \text{otherwise} \\
\end{cases}
$$

２．PDE（後退方程式）の形

これはブラウン運動の後退方程式

$$
\frac{\partial u}{\partial t} = \frac{\sigma^2}{2} \frac{\partial^2 u}{\partial x^2}
$$

３．数値的に解く方法

空間を差分（有限差分法）で離散化
時間は陽的オイラー法で進める
Neumann 境界条件として $${u[0] = u[1], u[-1] = u[-2]}$$ を使って端の暴走を防止

４．結果の解釈

$${u(x, T)}$$ は、今 $${x}$$ にいるときに、$${T}$$ 秒後に $${[-1,1]}$$ にいる確率
中心ほど高く（1に近く）、遠くほど低く（0に近い）

✅ この方法の意義

通常なら何千回もシミュレーションして平均をとる必要があるが、後退方程式を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の教科書です。
よかったらぜひ、お試しくださいませ。

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

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

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

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

この春やりたいこと

2.0万件