Pythonによる遅延微分方程式の数値解析と可視化

# はじめに

遅延微分方程式を python で簡単に解く例を紹介します。

[遅延微分方程式](https://en.wikipedia.org/wiki/Delay_differential_equation)の解析解を見つけることは、一般には非常に困難です。

- 数値的解法 (方程式を離散化し、数値的に近似解を求める。この記事で紹介。)
- 定性的分析 (方程式の挙動、特に安定性や周期解などの性質を理解する)
- 近似的な解法 (特定の条件下で、方程式の近似解を求める)

などの方法で解かれ、数値的な手法や定性的な分析を通じて方程式の挙動を理解することが一般的です。

https://colab.research.google.com/drive/1edNB_T3625ulScT2VmdAqAeaM02TKWt5?usp=sharing

に、この記事の google colab をおいてます。

---

## 遅延微分方程式の定義

遅延微分方程式は通常、以下の形式で表されます:
``` math:
\frac{d x(t)}{d t}=f(x(t), x(t-\tau))
``` 

ここで、 $x(t)$ は時刻 $t$ でのシステムの状態、 $\tau$ は遅延時間、 $f$ は関数です。

---

## 遅延化学反応モデル

化学反応の例として、反応速度定数 $k_1$ と $k_2$ を持つ反応を考えます。遅延時間 $\tau$ を含む微分方程式は以下のようになります：

``` math:
\frac{d N(t)}{d t}=-k 1 \cdot N(t)+k 2 \cdot N(t-\tau)
```

ここで、 $N(t)$ は時刻 $t$ での化学物質の濃度です。

---

## python による数値解法

Python の `scipy.integrate.solve_ivp` 関数を用いて微分方程式を解くことはできますが、遅延項を直接扱うことはできません。そこで、遅延項をリストに保存し、それを参照することで遅延を考慮することにします。

---

### 離散化して時間発展を解く部分

``` python:
import numpy as np
import matplotlib.pyplot as plt

def solve_delayed_chemical_reaction(tfinal, k1, k2, tau):
    dt = 0.1  # 時間ステップ
    N = 1.0  # 初期濃度
    t = 0
    results = [(t, N)]  # 結果を保存

    while t < tfinal:
        # 遅延時間前の状態を取得
        N_tau = next((x[1] for x in reversed(results) if x[0] <= t - tau), N)

        # 微分方程式の計算
        dNdt = -k1 * N + k2 * N_tau
        N += dNdt * dt

        # 次の時間ステップ
        t += dt

        # 状態を記録
        results.append((t, N))

    return results
```

---

ここで、差分化して時間発展を計算しています。

``` python:
N_tau = next((x[1] for x in reversed(results) if x[0] <= t - tau), N)
```

の部分は、`(x[1] for x in reversed(results) if x[0] <= t - tau)` はジェネレータ式で、`results` リストを逆順に走査し、時刻 `t - tau` よりも前の時刻の要素（`x[0]`）の濃度（`x[1]`）を生成します。

`next()` はこのジェネレータから最初の要素を取り出します。python の組み込み関数の `next()` 関数は、オプショナルな第二引数としてデフォルト値を取ることができます。今の場合は、もし該当する要素がなければ、デフォルト値 `N` を返します。この問題の場合は、時刻遅れが小さいことを想定して、逆順に探索することで、遅延時間 `tau` 前の濃度 `N_tau` を効率的に求めます。

---

### パラメータを与えてプロットする関数

遅延時間 $\tau$ と反応速度定数 $k_1$、$k_2$ を変化させて、それらが化学反応に与える影響を観察します。異なるパラメータの組み合わせで濃度の変化をプロットすることで、遅延と反応速度が反応ダイナミクスにどのように影響を与えるかを視覚的に理解することができます。

``` python:
def explore_parameter_effects(taulist=[10,15], k1list=[0.3,0.4], k2list=[0.1,0.2]):
    for tau in taulist:
        for k1 in k1list:
            for k2 in k2list:
                results = solve_delayed_chemical_reaction(100, k1, k2, tau)
                t_values, N_values = zip(*results)
                plt.plot(t_values, N_values, label=f"tau={tau}, k1={k1}, k2={k2}")
    plt.legend(loc='center left', bbox_to_anchor=(1., .5))
    plt.yscale("linear")
    plt.xscale("linear")
    plt.show()
```

---

### 生成されるプロット

- デフォルトの引数 (taulist=[10,15], k1list=[0.3,0.4], k2list=[0.1,0.2])

``` python:
explore_parameter_effects()
```

![スクリーンショット 2024-01-23 18.25.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/8b8805dc-da4f-70e4-a560-0d4826be5d73.png)

- 引数を指定する例 (taulist=[15,30], k1list=[0.2], k2list=[0.1,0.15])

``` python:
explore_parameter_effects(taulist=[15,30], k1list=[0.2], k2list=[0.1,0.15])
```

![スクリーンショット 2024-01-23 18.25.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/fbeff837-bff2-7d93-2e55-8a040304abe5.png)


<details><summary>コード全文</summary>

``` python:
import numpy as np
import matplotlib.pyplot as plt

def solve_delayed_chemical_reaction(tfinal, k1, k2, tau):
    dt = 0.1  # 時間ステップ
    N = 1.0  # 初期濃度
    t = 0
    results = [(t, N)]  # 結果を保存

    while t < tfinal:
        # 遅延時間前の状態を取得
        N_tau = next((x[1] for x in reversed(results) if x[0] <= t - tau), N)

        # 微分方程式の計算
        dNdt = -k1 * N + k2 * N_tau
        N += dNdt * dt

        # 次の時間ステップ
        t += dt

        # 状態を記録
        results.append((t, N))

    return results

def explore_parameter_effects(taulist=[10,15], k1list=[0.3,0.4], k2list=[0.1,0.2]):
    for tau in taulist:
        for k1 in k1list:
            for k2 in k2list:
                results = solve_delayed_chemical_reaction(100, k1, k2, tau)
                t_values, N_values = zip(*results)
                plt.plot(t_values, N_values, label=f"tau={tau}, k1={k1}, k2={k2}")
    plt.legend(loc='center left', bbox_to_anchor=(1., .5))
    plt.yscale("linear")
    plt.xscale("linear")
    plt.show()

explore_parameter_effects()

explore_parameter_effects(taulist=[15,30], k1list=[0.2], k2list=[0.1,0.15])
```
</details>

---

# まとめ

遅延の存在は、システムの安定性、振動、周期的な挙動など、新しいダイナミクスの特性をもたらしますが、多くの遅延微分方程式は非線形であり、その解析は数学的に挑戦的ですが、python でプロットしてみると、少し親近感が湧く気がします。

---

## 補足１:電磁気学の遅延ポテンシャルとの違い

電磁気の遅延ポテンシャルは、あくまで一つの解です。

- 電磁気学で登場する遅延ポテンシャルは、マクスウェル方程式の特定の解の一つです。方程式は、普通の時間遅れの項がない微分方程式です。
- 遅延微分方程式は、遅延効果を直接方程式に組み込んだ形式で、過去の状態が現在の状態に影響を与えることを組み込んだ方程式です。

ちなみに、マクスウェル方程式に、無理やり時間遅れの項を入れて遅延微分方程式を作ることは数学的には可能なのか(は知りませんが..)、光速度が一定など物理的制約を破らずないことは必要です。

---