« 四川料理 洛楽(京都駅 近鉄みやこみち)で酸辣湯麺+ミニ中華丼を食す。 | トップページ | 「標高と水の沸点の関係」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUPしてますが、すごいコメントがついていた! »

2021年2月 6日 (土)

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

先日、

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

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

dVi/dt = a [U(Xi+1 - Xi) - Vi]

dXi/dt = Vi

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

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

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

Ov_model

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

ソースコードはこちら:

%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)

 

« 四川料理 洛楽(京都駅 近鉄みやこみち)で酸辣湯麺+ミニ中華丼を食す。 | トップページ | 「標高と水の沸点の関係」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUPしてますが、すごいコメントがついていた! »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« 四川料理 洛楽(京都駅 近鉄みやこみち)で酸辣湯麺+ミニ中華丼を食す。 | トップページ | 「標高と水の沸点の関係」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUPしてますが、すごいコメントがついていた! »

最近の記事

2021年3月
  1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30 31      
フォト
無料ブログはココログ