Lotka–Volterra方程式を逐次線形化＋離散化


---

## ✅ 1. 非線形連続時間モデル（原型）

被食者 \$x(t)\$、捕食者 \$y(t)\$：

$$
\begin{cases}
\displaystyle \frac{dx}{dt} = \alpha x - \beta x y \\[6pt]
\displaystyle \frac{dy}{dt} = \delta x y - \gamma y
\end{cases}
$$

* \$\alpha\$：被食者の自然増殖率
* \$\beta\$：捕食による被食者減少係数
* \$\delta\$：捕食により捕食者が増える係数
* \$\gamma\$：捕食者の自然死率

---

## ✅ 2. ヤコビ行列による逐次線形化（任意点まわり）

状態ベクトル \$\mathbf{z}(t) = \begin{bmatrix} x(t) \ y(t) \end{bmatrix}\$ に対し：

$$
\mathbf{f}(\mathbf{z}) =
\begin{bmatrix}
\alpha x - \beta x y \\
\delta x y - \gamma y
\end{bmatrix}
$$

ヤコビ行列（Jacobian）：

$$
\mathbf{J}(x, y) =
\begin{bmatrix}
\alpha - \beta y & -\beta x \\
\delta y & \delta x - \gamma
\end{bmatrix}
$$

よって、逐次線形近似：

$$
\dot{\mathbf{z}} \approx \mathbf{f}(x_0, y_0) + \mathbf{J}(x_0, y_0) \cdot 
\begin{bmatrix}
x - x_0 \\
y - y_0
\end{bmatrix}
$$

---

## ✅ 3. Euler法による数値離散化（非線形のまま）

時間刻み \$\Delta t\$ による陽的Euler法：

$$
\begin{cases}
x_{k+1} = x_k + \Delta t \cdot (\alpha x_k - \beta x_k y_k) \\[6pt]
y_{k+1} = y_k + \Delta t \cdot (\delta x_k y_k - \gamma y_k)
\end{cases}
$$

---

## ✅ 4. 線形化＋Euler離散：逐次線形状態更新（拡張可能構造）

各時刻 \$k\$ におけるヤコビアン：

$$
\mathbf{A}_k =
\begin{bmatrix}
\alpha - \beta y_k & -\beta x_k \\
\delta y_k & \delta x_k - \gamma
\end{bmatrix}
$$

状態差分更新：

$$
\mathbf{z}_{k+1} = \mathbf{z}_k + \Delta t \cdot \left[
\mathbf{f}(\mathbf{z}_k) + \mathbf{A}_k (\mathbf{z}_k - \mathbf{z}_k)
\right] = \mathbf{z}_k + \Delta t \cdot \mathbf{f}(\mathbf{z}_k)
$$

※非線形項による更新は維持しつつ、ヤコビ行列を随時更新することで**線形制御への展開可能な構造**となる。

---

## ✅ 5. 状態空間形式への整理（平衡点まわり）

平衡点 \$(x^*, y^*)\$ まわりでの線形化：

$$
\frac{d}{dt} \Delta \mathbf{z}(t) = \mathbf{A} \cdot \Delta \mathbf{z}(t)
$$

* \$\Delta \mathbf{z} = \mathbf{z} - \mathbf{z}^\*\$
* \$\mathbf{A} = \mathbf{J}(x^*, y^*)\$

---

## ✅ 6. ラプラス変換と伝達関数表現

$$
s \Delta \mathbf{Z}(s) = \mathbf{A} \cdot \Delta \mathbf{Z}(s)
\quad \Rightarrow \quad
\Delta \mathbf{Z}(s) = (s\mathbf{I} - \mathbf{A})^{-1} \cdot \Delta \mathbf{Z}(0)
$$

伝達関数形式による安定性・応答解析の出発点となる。

---

## ✅ 7. まとめ表

| 項目       | 数式・内容                                                                                        |
| -------- | -------------------------------------------------------------------------------------------- |
| 非線形モデル   | \$\dot{x} = \alpha x - \beta x y\$, \$\dot{y} = \delta x y - \gamma y\$                      |
| ヤコビ行列    | \$\begin{bmatrix} \alpha - \beta y & -\beta x \ \delta y & \delta x - \gamma \end{bmatrix}\$ |
| 線形近似     | \$\dot{z} \approx f(x\_0, y\_0) + J \cdot (z - z\_0)\$                                       |
| Euler法更新 | \$z\_{k+1} = z\_k + \Delta t \cdot f(z\_k)\$                                                 |
| 状態方程式    | \$\frac{d}{dt} \Delta z = A \Delta z\$                                                       |
| 伝達関数表現   | \$\Delta Z(s) = (sI - A)^{-1} \Delta Z(0)\$                                                  |

---

```python
# Program Name: lotka_volterra_linear_vs_nonlinear.py
# Creation Date: 20250617
# Overview: Compare nonlinear Lotka–Volterra solution with sequential linearization (discretized)
# Usage: Run script to simulate and plot nonlinear vs linearized dynamics

# --- ライブラリのインストール / Install required libraries ---
!pip install numpy matplotlib

# --- ライブラリのインポート / Import libraries ---
import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ定義 / Define parameters ---
alpha = 1.1   # 被食者の自然増殖率 / Natural growth rate of prey
beta = 0.4    # 捕食係数 / Predation rate
delta = 0.1   # 捕食による成長率 / Growth rate of predator from eating prey
gamma = 0.4   # 捕食者の自然死率 / Natural death rate of predator
T = 30        # シミュレーション時間 / Simulation duration
dt = 0.01     # 時間刻み幅 / Time step size
steps = int(T / dt)
x0, y0 = 10, 5  # 初期状態 / Initial state

# --- 非線形モデル定義 / Nonlinear Lotka–Volterra model ---
def f(x, y):
    dx = alpha * x - beta * x * y
    dy = delta * x * y - gamma * y
    return np.array([dx, dy])

# --- 線形化関数 / Jacobian-based linearization at given point (x, y) ---
def jacobian(x, y):
    return np.array([
        [alpha - beta * y, -beta * x],
        [delta * y, delta * x - gamma]
    ])

# --- 非線形オイラー解法 / Euler method for nonlinear system ---
x_nl, y_nl = [x0], [y0]
x, y = x0, y0
for _ in range(steps):
    dx, dy = f(x, y)
    x += dt * dx
    y += dt * dy
    x_nl.append(x)
    y_nl.append(y)

# --- 線形化による逐次オイラー解法 / Euler method with local linearization ---
z_lin = [np.array([x0, y0])]
z = np.array([x0, y0])
for _ in range(steps):
    J = jacobian(z[0], z[1])  # 現在の点でヤコビ行列を計算 / Compute Jacobian at current point
    dz = J @ z  # 線形近似ベクトル場で更新 / Linear vector field
    z = z + dt * dz
    z_lin.append(z)

# --- 時間軸 / Time vector ---
t = np.linspace(0, T, steps + 1)

# --- グラフ描画 / Plot nonlinear vs linearized results ---
plt.figure(figsize=(12, 5))

# --- 個体数の時間変化 / Population over time ---
plt.subplot(1, 2, 1)
plt.plot(t, x_nl, label="Prey (nonlinear)")
plt.plot(t, y_nl, label="Predator (nonlinear)")
plt.plot(t, [z[0] for z in z_lin], '--', label="Prey (linearized)")
plt.plot(t, [z[1] for z in z_lin], '--', label="Predator (linearized)")
plt.xlabel("Time")
plt.ylabel("Population")
plt.title("Lotka–Volterra Dynamics Over Time")
plt.legend()

# --- 相平面図 / Phase portrait ---
plt.subplot(1, 2, 2)
plt.plot(x_nl, y_nl, label="Nonlinear trajectory")
plt.plot([z[0] for z in z_lin], [z[1] for z in z_lin], '--', label="Linearized trajectory")
plt.xlabel("Prey")
plt.ylabel("Predator")
plt.title("Phase Portrait: Nonlinear vs Linearized")
plt.legend()

plt.tight_layout()
plt.show()

```
結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/e7fe2924-db17-49f9-9e66-3754754079a4.png)