渋滞を表すOVモデル（最適速度モデル）をPython＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使って計算、GIFアニメにしてみた。又吉直樹のヘウレーカが渋滞特集だったので。

先日、

2/3の又吉直樹のヘウレーカ 「どうすれば“密”を避けられる?」（渋滞学で有名な西成活裕さんがゲスト）をTwitterでリアルタイムメモしてました。その記録＆昔、渋滞の微分方程式OVモデルを計算した結果。

というのを見てました。それで昔、OVモデル（最適速度模型、Optimal Velocity model）

dV_i/dt = a [U(X_i+1 - X_i) - V_i]

dX_i/dt = V_i

ここでU(x) = tanh(x-c) + tanh(c)

をExcelのVBAで計算したことを思い出した。今やるならPython＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)かな、ということでやってみた。

クリックでGIFアニメが始まります。

なんか逆回転するようにも見える、、、

ソースコードはこちら：

%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
import numba
import random
from scipy.stats import cauchy

INTERVAL = 0.05 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

Nt=10000
N=40
L=80
c=2.0
a=1.3
dx=L/N

tmin=0.0
tmax=500.0

t0=0
dt = (tmax-tmin)/Nt

t=np.linspace(0,tmax,Nt)
x0=np.empty(N*2)

def Velocity(x):
    return np.tanh(x-2)+np.tanh(2)

p=0
for i in range(N):
    x0[i]=p
    x0[N+i]=Velocity(dx)
    p+=dx
    p+=0.01*(random.uniform(0,1)-0.5)

sol= np.empty((Nt, N*2))
sol[0] = x0

def ov_model(t, x): #odeintのときとt,xの並びが逆
    f=np.zeros(N*2)

    for i in range(N):
        f[i]=x[N+i]
        if i< N-1:
            d=x[i+1]-x[i]
        else:
            d=x[0]-x[i]

        if d>L:
            d -= L

        if d<0:
            d +=L
            f[N+i]=a*(Velocity(d)-x[N+i])

    return f

solver=ode(ov_model)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆

k=1
while solver.successful() and solver.t < tmax:
    solver.integrate(t[k])

    sol[k] = solver.y
    for i in range(N):
        if sol[k,i]>L:
    sol[k,i] -=L
    k+= 1

theta=np.linspace(0,2*np.pi,100)

def update(i):
    plt.cla() #現在描写されているグラフを消去
    plt.xlim(-1.1,1.1)
    plt.ylim(-1.1,1.1) 
    plt.plot(np.cos(theta),np.sin(theta),c='Black',alpha=0.2)
    plt.scatter(np.cos(2*np.pi*sol[i,0:N]/L),np.sin(2*np.pi*sol[i,0:N]/L),c='Red')
    plt.xlabel("x")
    plt.ylabel("y")
plt.title("OV model by dop853")
plt.grid()

# Plot
fig = plt.figure(figsize=(8,8))
ani = FuncAnimation(fig, update, frames=range(0,Nt,15), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("ov_model.gif", writer = 'pillow',fps=FPS)