02 連続モデルの数値解と解析解

連続時間の指数増殖モデル $\frac{dN}{dt} = rN$ とロジスティック成長モデル $\frac{dN}{dt} = rN(1 - N/K)$ をオイラー法 (Euler method)で解き，解析解 (analytical solution)と比較する．

解析解のプロット¶

指数増殖¶

指数増殖モデル $\frac{dx}{dt} = a x$ の解析解は $x(t) = x_0 e^{at}$ になる．増殖率 $a$ はマルサス係数 (Malthusian coefficient)と呼ばれる．

math.exp を用いて時系列を計算する．

# 02-01. 指数増殖（解析解）
import math

import matplotlib.pyplot as plt

初期値 $x_0 = 10$ ，増殖率 $a = 0.2$ ，刻み幅 $\Delta t = 0.1$ で $t = 0$ から $t = 100$ まで計算する．

a = 0.2
x0 = 10
dt = 0.1

t_list = []
x_list = []
for i in range(1001):
    t = dt * i
    x = x0 * math.exp(a * t)
    t_list.append(t)
    x_list.append(x)

plt.plot(t_list, x_list)

Solution to Exercise 1

r = 0.2
K = 1000
N0 = 10
dt = 0.1

t_list = []
N_list = []
for i in range(1001):
    t = dt * i
    N = K / (1 + ((K - N0) / N0) * math.exp(-r * t))
    t_list.append(t)
    N_list.append(N)

plt.plot(t_list, N_list)
plt.title("Logistic growth (analytical)")
plt.xlabel("Time (t)")
plt.ylabel("Pop. size (N)")

数値解（オイラー法）¶

指数増殖¶

オイラー法では $x_{t+\Delta t} = x_t + a x_t \, \Delta t$ で次の時刻の値を逐次計算する．離散モデルの実装と形が同じであることに注目する．

# 02-02. 指数増殖（数値解）
a = 0.2
x0 = 10

dt = 0.1

t = 0
x = x0
t_list = [t]
x_list = [x]
for i in range(1000):
    t = dt * (i + 1)
    x = x + a * x * dt

    t_list.append(t)
    x_list.append(x)

plt.plot(t_list, x_list)

解析解と数値解の比較¶

同じパラメータで解析解と数値解を計算し，1つの図に重ねて誤差の出方を観察する．

# 02-03. 指数増殖（解析解と数値解の比較）
a = 0.2
x0 = 10
t_end = 100

dt_a = 0.1
i_end_a = int(t_end / dt_a)

dt_n = 0.1
i_end_n = int(t_end / dt_n)

# 解析解
x_a_list = []
t_a_list = []
for i in range(i_end_a):
    t = dt_a * i
    x = x0 * math.exp(a * t)
    t_a_list.append(t)
    x_a_list.append(x)

# 数値解
t = 0
x = x0
t_n_list = [t]
x_n_list = [x]
for i in range(i_end_n):
    t = dt_n * (i + 1)
    x = x + a * x * dt_n
    t_n_list.append(t)
    x_n_list.append(x)

# 可視化
plt.plot(t_a_list, x_a_list)
plt.plot(t_n_list, x_n_list)

Solution to Exercise 2

r = 0.2
K = 1000
N0 = 10
dt = 0.1
t_end = 100
i_end = int(t_end / dt)

# 解析解
t_a_list = []
N_a_list = []
for i in range(i_end):
    t = dt * i
    N = K / (1 + ((K - N0) / N0) * math.exp(-r * t))
    t_a_list.append(t)
    N_a_list.append(N)

# 数値解（オイラー法）
t = 0
N = N0
t_n_list = [t]
N_n_list = [N]
for i in range(i_end):
    t = dt * (i + 1)
    N = N + r * N * (1 - N / K) * dt
    t_n_list.append(t)
    N_n_list.append(N)

plt.plot(t_a_list, N_a_list, label="Analytical")
plt.plot(t_n_list, N_n_list, label="Numerical (Euler)")
plt.title("Logistic growth")
plt.xlabel("Time (t)")
plt.ylabel("Pop. size (N)")
plt.legend()