
"""
OVモデルの本体
"""
import numpy as np
from scipy.integrate import odeint

class ov_model:
    def __init__(self, carNum, L, a, dt, initialize=True):
        self.carNum = carNum
        self.dt =dt
        self.L = L
        self.set_params(L,a)
        if initialize:
            self.initialize_x_v()
        
    def set_params(self,L,a):
        self.a = a
        self.L = L
       
    def initialize_x_v(self): # must call this after set_params
        self.x = np.arange(0, self.L, self.L/self.carNum)        # initial position
        self.v = np.random.random(self.carNum)     # initial velocity
        self.x_v = np.concatenate((self.x,self.v))
        self.x_v_trace = np.empty((self.carNum*2))
        self.t_trace = np.zeros(1) # array of time for plot
        self.present_time = 0.0
        self.x_v_trace[:] = self.x_v       # data for plot
        
    def ov_eq(self, x_v, t, carNum, a, L):
        '''
        r.h.s of the ov equation
        '''
        x = x_v[:carNum]
        v = x_v[carNum:]
        x_fw = np.roll(x,-1)
        delx = (x_fw- x + L) % L
        return np.concatenate((v, a*(np.tanh(delx - 2.0) +np.tanh(2.0) -v)),axis=0)

    def update_unit(self, n_steps, store_trace=False):
        '''
        solving the ov equation by "n_steps" iterations
        '''
        t = np.arange(0, n_steps*self.dt, self.dt)
        new_x_v = odeint(self.ov_eq,self.x_v,t, args=(self.carNum,self.a,self.L))
        self.x_v = new_x_v[-1] % self.L
        self.present_time = self.present_time+self.dt*n_steps
        if store_trace:
            self.x_v_trace = np.vstack((self.x_v_trace, self.x_v))
            self.t_trace = np.vstack((self.t_trace,np.array([self.present_time])))
    
    def shift_x_v_t(self, n):
        if n > self.t_trace.size:
            n=-1
        self.t_trace = self.t_trace[n:]
        self.x_v_trace = self.x_v_trace[n:]


%matplotlib nbagg
from ipywidgets import interact_manual, Checkbox, widgets
from IPython.display import clear_output, display, HTML
model = ov_model(50, 100, 1.5, 0.01, True)

def show_trace_ov(carNum, rho, a, initialize, max_t):
    global model
    unit_steps =10
    dt=0.01
    if initialize:
        model = ov_model(carNum, carNum/rho, a, dt, initialize)
    else:
        model.shift_x_v_t(int(max_t/dt))
    for i in range(int(max_t/dt/unit_steps)):
        model.update_unit(unit_steps, True)
    # plot the traces
    import matplotlib.pyplot as plt
    fig = plt.figure(figsize=(10,5))
    axes = fig.add_subplot(1,1,1)
    axes.set_xlabel("time $t$")
    axes.set_ylabel("position $x$")
    axes.set_title("$\\rho=" + str(rho) + "$, $a=" + str(a) + "$")
    for i in range(carNum):
        axes.plot(model.t_trace, model.x_v_trace[:,i], '.', markersize=3)
    plt.show()
    
interact_manual(show_trace_ov, carNum=50, rho=widgets.FloatSlider(min=0.1,max=1.5, step=0.05, value=0.6),
                a=widgets.FloatSlider(min=0.5, max=3.0, step=0.1, value=1.5),
                initialize=Checkbox(value=True),max_t=(100,300))

<function __main__.show_trace_ov(carNum, rho, a, initialize, max_t)>


"""
OVモデルのアニメーション表示　（円形の道路上を車が移動する）
新しいウィンドウが開かないときは、FileメニューからNew Notebookを開き、
上の「OVモデルの本体」と書かれたセルと、このセルをコピーして実行してみてください。
違うノートブック間のコピー＆ペーストは、ブラウザのコピー、ペースト機能を使ってください。
（jupyterのメニューにあるコピー、ペーストは同じノートブックないだけで機能します。）
"""
import matplotlib.pyplot as plt
from matplotlib import animation
%matplotlib nbagg

carNum=50
rho=0.5
a=1.2
L=carNum/rho
model = ov_model(carNum, L, a, 0.01, initialize=True)

# create the animation frame and execute
radius = L/np.pi/2.0
fig = plt.figure(figsize=(6,6))
axes = fig.add_subplot(1,1,1)
pt, = axes.plot([], [], "o")
frame_size = radius*1.1
axes.set_xlim((-frame_size, frame_size))
axes.set_ylim((-frame_size, frame_size))
axes.set_title('a='+ str(a)+ ', $\\rho=$'+ str(rho))
    
def init_anim():
    pt, = axes.plot([], [], "o")
    model.initialize_x_v()
    return pt,
    
def anim_update(i):
    n_steps = 20
    model.set_params(L,a)
    model.update_unit(n_steps, False)
    arg_fac = 2.0*np.pi/L
    x = radius*np.cos(model.x_v[:carNum]*arg_fac)
    y = radius*np.sin(model.x_v[:carNum]*arg_fac)
    pt.set_xdata(x)
    pt.set_ydata(y)
    return pt,
    
anim=animation.FuncAnimation(fig, anim_update, init_func=init_anim,
                             frames=3000, interval=20, blit=False, repeat=False) # framesは終了までの繰り返し数
# Set blit=False when you use the following "save" function
# anim.save('ov_anim1.mp4', fps=15, extra_args=['-vcodec', 'libx264'])
plt.show()


import numpy as np
import matplotlib.pyplot as plt
%matplotlib nbagg
x = np.arange(0, 4.0, 0.02)
plt.plot(x, 2./(np.cosh(x-2.0)**2))
plt.xlabel("$\Delta x = \\rho^{-1}=L/N$")
plt.ylabel("$a$")
plt.text(0.3, 1.5, "stable", fontsize=12)
plt.text(1.5, 0.5, "unstable", fontsize=12)
plt.show()


'''
最適速度関数の描画
'''
import numpy as np
import matplotlib.pyplot as plt
%matplotlib nbagg
x = np.arange(0, 8.0, 0.02)
plt.plot(x, np.tanh(x-2.0)+np.tanh(2.0))
plt.xlabel("$\Delta x$")
plt.ylabel("$V(\Delta x)$")
plt.title("$V(\Delta x) = \\tanh (\Delta x - 2.0) + \\tanh (2.0)$")
plt.show()


class ov_model_ex(ov_model):
    def __init__(self, carNum, L, a, dt, initialize=True):
        super().__init__(carNum, L, a, dt, initialize)
        
    def get_headway(self):
        x = self.x_v[:carNum]
        x_fw = np.roll(x,-1)
        delx = (x_fw- x + self.L) % self.L
        return delx
    
import matplotlib.pyplot as plt
from matplotlib import animation
%matplotlib nbagg
carNum=200
rho=0.5
a=2.1
L=carNum/rho
model = ov_model_ex(carNum, L, a, 0.01, initialize=True)

# create the animation frame and execute
fig = plt.figure(figsize=(6,6))
axes = fig.add_subplot(1,1,1)
pt, = axes.plot([], [], "o")

axes.set_xlim((0, 4.0))
axes.set_ylim((0.0, 1.0+np.tanh(2.0)))
axes.set_title('a='+ str(a)+ ', $\\rho=$'+ str(rho))
axes.set_xlabel("$\Delta x$")
axes.set_ylabel("$v$")
# plot optimal velocity
x = np.arange(0,4.0,0.01)
y = np.tanh(x-2.0) + np.tanh(2.0)
axes.plot(x,y,"r-")

def init_anim():
    pt, = axes.plot([], [], "o")
    model.initialize_x_v()
    return pt,
    
def anim_update(i):
    n_steps = 20
    model.set_params(L,a)
    model.update_unit(n_steps, False)
    arg_fac = 2.0*np.pi/L
    x = model.get_headway()
    y = model.x_v[carNum:]
    pt.set_xdata(x)
    pt.set_ydata(y)
    return pt,
    
anim=animation.FuncAnimation(fig, anim_update, init_func=init_anim,
                               frames=3000, interval=10, blit=True, repeat=False)
# Set blit=False when you use the following "save" function
#anim.save('ov_anim.mp4', fps=15, extra_args=['-vcodec', 'libx264'])
plt.show()


%matplotlib inline
import numpy as np
from ipywidgets import interact, widgets
from IPython.display import clear_output, display, HTML

def show_4potential(a):
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(1,1,1)
    ax.set_xlim(( -4, 4))
    ax.set_ylim((-10, 10))
    ax.set_xlabel("$\\phi$")
    ax.set_ylabel("$f(\\phi)$")
    ax.set_title("$\\phi^4$ポテンシャル $y=a \\phi^2 +\\phi^4$, a=" + str(a))
    ax.grid()
    x = np.linspace(-4.0,4.0,100)
    y = a* x*x + x**4
    ax.plot(x,y)
    plt.show()

interact(show_4potential, a=widgets.FloatSlider(min=-6.0, max=2.0))

<function __main__.show_4potential(a)>

Table of Contents

概要説明¶

前処理¶

pythonによるプログラムとjupyter notebook¶

jupyter notebookの使い方¶

サンプルプログラムのダウンロードとjupyterの立ち上げ¶

jupyterの基本的なメニュー　編集、実行など¶

ノートブックの読み込み、新規作成¶

メニュー¶

編集上のTips¶

レポートについて¶

作成及び提出方法¶

実習¶

Tips¶

レポートへの画像の貼りこみ¶

Markdownの記法¶

最適速度モデル¶

OVモデルのシミュレーション¶

OVモデルの数値計算本体のクラス定義¶

軌跡の描画¶

車の動きのアニメーション¶

参考¶

各自のパソコンでの実習¶

説明図の作成¶

安定・不安定領域の図作成プログラム¶

最適速度関数の描画¶

車間距離-速度のアニメーション¶

不安定性の直感的理解¶

Table of Contents

概要説明¶

前処理¶

pythonによるプログラムとjupyter notebook¶

jupyter notebookの使い方¶

サンプルプログラムのダウンロードとjupyterの立ち上げ¶

jupyterの基本的なメニュー 編集、実行など¶

ノートブックの読み込み、新規作成¶

メニュー¶

編集上のTips¶

レポートについて¶

作成及び提出方法¶

実習¶

Tips¶

レポートへの画像の貼りこみ¶

Markdownの記法¶

最適速度モデル¶

OVモデルのシミュレーション¶

OVモデルの数値計算本体のクラス定義¶

軌跡の描画¶

車の動きのアニメーション¶

参考¶

各自のパソコンでの実習¶

説明図の作成¶

安定・不安定領域の図作成プログラム¶

最適速度関数の描画¶

車間距離-速度のアニメーション¶

不安定性の直感的理解¶

jupyterの基本的なメニュー　編集、実行など¶