モデルの概要
イギリスの数学者、気象学者であったL.F.Richardsonが提案したモデル。２か国間の戦力の増強の過程を予測するモデル。
$$ {\frac{dx}{dt} = ky-\alpha x+g} $$$${\frac{dy}{dt} = lx-\beta y+h}$$
x,y : X国、Y国の戦力
k,l : 防衛係数
α,β : 疲労係数
モデルの意味
ky、lxは自国の軍備の増加率が相手国の軍備の大きさに比例することを表す。
自国の戦力に比例して増加率が下がることを表している。
軍事増強への肯定的、否定的な民意などによる影響を表している。
両国の戦力の増強率は、以下の点$(\bar{x},\bar{y})$で、0となる。すなわち平衡点である。
$$\bar{x}=\frac{hk+\beta g}{\alpha \beta-kl}$$ $$\bar{y}=\frac{gl+\alpha h}{\alpha \beta-kl}$$
解の特性
特性根を$\lambda_1,\lambda_2$とすると、以下のような結論が得られる。
(Ⅰ)$\alpha \beta > kl$ のとき、$\lambda_1 <0,\lambda_2 <0$
(Ⅱ)$\alpha \beta > kl$ のとき、$\lambda_1\lambda_2 <0$
(Ⅰ)の場合、解は互いに負となり平衡点に収束する。
一方（Ⅱ）の場合、特性根の符号が異なるため、双曲線上の発散となる。
シミュレーションコード
係数を$k=l=1,\alpha = \beta = 2,g=2,h=3$として設定して、計算する。
scipyのsolve_ivpメソッドを用いて計算する。
まず、各微係数を計算する関数を定義する。
計算する期間:Tendを設定し、solve_ivpに関数、計算時間、初期値、定数値を入力し、しsolに代入する。
常微分方程式の計算
from scipy.integrate import solve_ivp
def odeRichar(t,z,l,alpha,beta,g,h):
dx=k*z[1]-alpha*z[0]+g
dx=l*z[0]-beta*z[1]+h
return [dx,dy]
k,l = 1,1
alpha, beta=2,2
x0 , y0 =10.0, -4.0 #初期の戦力
Tend = 7.0 #期間
#ODEを解く
"""
常微分方程式の初期値問題を解く関数。
"""
sol = solve_ivp(odeRichard, t_span=[0, Tend], y0 = [x0, y0], args=(k,l,alpha, beta, g, h), dense_output=True)
tc1 = np.linspace(0, Tend, 100) #時間配列
yc1 = sol.sol(tc1) #解の配列
結果の可視化
yc1とtc1を戦力均衡までの間の、戦力増強過程を可視化する。
fig, ax = plt.subplots(ncols=2, figsize=(14,4))
ax[0].plot(tc1, yc1[0].T, label='x')
ax[0].plot(tc1, yc1[1].T, label='y')
].T, label='y')
ax[0].grid()
ax[0].legend()
ax[0].set_xlabel('time')
ax[0].set_ylabel('amount of armaments')
ax[1].plot(tc1, yc1[0].T, label='x-y')
ax[1].grid()
ax[1].legend()
ax[1].set_xlabel('x')
ax[1].set_ylabel('y')
plt.show()
平衡点の検証
#理論の平衡点
barx = (h*k + beta*g)/(alpha*beta -k*l)
bary = (g*l + alpha*h)/(alpha*beta -k*l)
print('barx=', 'bary=',bary)
#数値計算結果
len = yc1[0].size
print('barx=', yc1[0][len-1], 'bary=', yc1[1][len-1])
微係数のベクトル場の可視化
$(dx,dy)$をベクトル場として、x-y平面に描画することで、収束発散を可視化する。
収束する場合
係数を$k=l=1,\alpha = \beta = 2,g=2,h=3$として設定して、計算する。
$$\alpha\times\beta = 4 > k \times l = 1$$
k, l = 1, 1 #防衛係数
alpha ,beta = 2, 2 #疲労係数
g, h = 2, 3 #不平因子
Xmin, Xmax = -5,12
Ymin, Ymax = -5,12
gridwidth = 2
X,Y=np.meshgrid(np.arange(Xmin,Xmax,gridwidth),np.arange(Ymin,Ymax,gridwidth))
fig = plt.subplots(figsize=(6,6))
x0 = (h*k + beta*g)/(alpha*beta -k*l)
y0 = (g*l + alpha*h)/(alpha*beta -k*l)
plt.plot(x0,y0,marker='o',color='red',markersize=12)
dx = k*Y - alpha*X + g
dy = l*X - beta*Y + h
U = dx
V = dy
plt.quiver(X,Y,U,V,color='black')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()
どのような戦力の組み合わせであっても、平衡点に向かって変化することがわかる。
発散する場合
係数を$k=l=2,\alpha = \beta = 1,g=2,h=3$として設定して、計算する。
$$\alpha\times\beta = 1 < k \times 4 = 1$$
同様のコードを用いて計算する。
xが正yが負の領域と、yが正xが負の領域では平衡点に向かっているが、平衡点付近では、x,yがともに正か負の方向に向かって変化するため、際限のない軍拡競争になることがわかる。
まとめ
Pythonと科学技術計算用ライブラリのScipy、グラフ描画ライブラリ matplotlibを用いて、リチャードソンの軍拡モデルを可視化した。