<jupyter_start><jupyter_text>OVモデル<jupyter_code>import numpy as np
import pathfollowing as pf
import numpy.random as rd
from scipy.integrate import ode, solve_ivp
from scipy.linalg import solve, eigvals
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set('poster', 'whitegrid', 'dark', rc={"lines.linewidth": 2, 'grid.linestyle': '-'})
Ncars = 10

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

def dV(x):
    return 1.0 / np.cosh(x-2)**2

def OV(t, x, a):
    y = np.zeros(2*(Ncars-1))
    for k in range(Ncars-2):
        y[2*k] = x[2*k+1]
        y[2*k+1] = V(x[2*(k+1)]) - V(x[2*k]) - x[2*k+1]
    y[-2] = x[-1]
    y[-1] = V(a - np.sum([x[::2]])) - V(x[-2]) - x[-1]
    return y

def f(x, a):
    y = np.zeros(2*(Ncars-1))
    for k in range(Ncars-2):
        y[2*k] = x[2*k+1]
        y[2*k+1] = V(x[2*(k+1)]) - V(x[2*k]) - x[2*k+1]
    y[-2] = x[-1]
    y[-1] = V(a[0] - np.sum([x[::2]])) - V(x[-2]) - x[-1]
    return y

def fx(x, a):
    A = np.zeros((2*(Ncars-1), 2*(Ncars-1)))
    for k in range(Ncars - 2):
        A[2*k, 2*k+1] = 1.0
        A[2*k+1, 2*k] = -dV(x[2*k])
        A[2*k+1, 2*k+1] = -1.0
        A[2*k+1, 2*(k+1)] = dV(x[2*(k+1)])
    A[-2, -1] = 1.0
    A[-1, ::2] = -dV(a[0] - np.sum([x[::2]]))
    A[-1, -2] += -dV(x[-2])
    A[-1, -1] = -1.0
    return A

def fa(x, a):
    y = np.zeros(2*(Ncars-1))
    y[-1] = dV(a[0] - np.sum([x[::2]]))
    return y
Npts = 64
def func(x, a):
    T = x[-1]
    h = 1.0 / Npts
    n = (len(x) - 1) // Npts
    ss = np.sin([2*np.pi*j*h for j in range(Npts)])
    y = np.zeros(len(x))
    # tau_0 = 0
    y[:n] = (x[n:2*n] - x[(Npts-1)*n:Npts*n])*0.5/h - T * f(x[:n], a)
    # tau_1 ... tau_{N-2}
    for j in range(1, Npts-1):
        y[j*n:(j+1)*n] = (x[(j+1)*n:(j+2)*n] - x[(j-1)*n:j*n])*0.5/h - T * f(x[j*n:(j+1)*n], a)
    # tau_{N-1} = (N-1) * h
    y[(Npts-1)*n:Npts*n] = (x[:n] - x[(Npts-2)*n:(Npts-1)*n])*0.5/h - T * f(x[(Npts-1)*n:Npts*n], a)
    # Phase condition
    y[-1] = x[:-1:n] @ ss
    return y

def dfdx(x, a):
    T = x[-1]
    n = (len(x)-1) // Npts
    h = 1.0 / Npts
    A = np.zeros((len(x), len(x)))
    I = np.identity(n)
    # j = 0
    A[:n, n:2*n] += I * 0.5 / h
    A[:n, (Npts-1)*n:Npts*n] -= I * 0.5 / h
    A[:n,:n] += -T * fx(x[:n], a)
    A[:n,-1] += -f(x[:n], a)
    # j = 1 ... N-2
    for j in range(1, Npts-1):
        A[j*n:(j+1)*n, (j+1)*n:(j+2)*n] += I * 0.5 / h
        A[j*n:(j+1)*n, (j-1)*n:j*n] -= I * 0.5 / h
        A[j*n:(j+1)*n, j*n:(j+1)*n] += -T * fx(x[j*n:(j+1)*n], a)
        A[j*n:(j+1)*n, -1] = -f(x[j*n:(j+1)*n], a)
    # j = N - 1
    A[(Npts-1)*n:Npts*n, :n] += I * 0.5 / h
    A[(Npts-1)*n:Npts*n, (Npts-2)*n:(Npts-1)*n] -= I * 0.5 / h
    A[(Npts-1)*n:Npts*n, (Npts-1)*n:Npts*n] += -T*fx(x[(Npts-1)*n:Npts*n], a)
    A[(Npts-1)*n:Npts*n, -1] += -f(x[(Npts-1)*n:Npts*n], a)
    # Phase condition
    A[-1, :-1:n] = np.sin([2*np.pi*j*h for j in range(Npts)])
    return A

def dfda(x, a):
    T = x[-1]
    n = (len(x) - 1) // Npts
    h = 1.0 / Npts
    b = np.zeros(len(x))
    for j in range(Npts):
        b[j*n:(j+1)*n] = -T*fa(x[j*n:(j+1)*n], a)
    return b<jupyter_output><empty_output><jupyter_text>一様交通流の追跡<jupyter_code>a=np.array([1.0])
x=np.array([a[0] / Ncars if m%2==0 else 0.0 for m in range(2*(Ncars-1))])
f(x,a)
a=np.array([1.0])
x=np.array([a[0] / Ncars if m%2==0 else 0.0 for m in range(2*(Ncars-1))])
bd,bp,lp,hb=pf.pathfollow(x, a, f, fx, fa, nmax=50, h=1.0, epsr=1.0e-10, epsb=1.0e-10, problem='equilibria', quiet=True)
print(hb)
print(bd[-1]['a'])
bd_r = np.array([bd[m]['a'][0] for m in range(len(bd))])
bd_x = np.array([bd[m]['x'][0] for m in range(len(bd))])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(bd_r, bd_r/Ncars)
ax.plot(bd_r, bd_x, '.')<jupyter_output><empty_output><jupyter_text>渋滞流への追跡の切り替え<jupyter_code>vr = bd[hb[0]]['vr']
vi = bd[hb[0]]['vi']
vr /= np.sqrt(vr@vr)
vi /= np.sqrt(vi@vi)
rr = np.hypot(vr, vi)
tt = np.array([np.angle(vr[i]+1.j*vi[i]) for i in range(len(vr))])
ev = eigvals(fx(bd[hb[0]+1]['x'], bd[hb[0]+1]['a']))
ev
ndim = 2*(Ncars-1)
x = np.zeros(ndim*Npts+1)
da = np.abs(bd[hb[0]]['a'][0] - bd[hb[0]+1]['a'][0])
for j in range(Npts):
    x[ndim*j:ndim*(j+1)] = bd[hb[0]+1]['x'] + 10 * np.sqrt(da) * rr * np.cos(2*np.pi*j/Npts + tt - tt[0])
x[-1] = 2*np.pi/np.imag(ev[0])
y = np.copy(x)
a = bd[hb[0]+1]['a']
for m in range(10):
    b = func(y, a)
    A = dfdx(y, a)
    y -= solve(A, b)
    print(np.linalg.norm(b))
# x = np.copy(y)
print(y)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x[1:-1:2*ndim],'o-')
ax.plot(x[3:-1:2*ndim],'o-')
ax.plot(x[5:-1:2*ndim],'o-')
x = np.copy(y)
bdc=pf.pathfollow00(x, a, func, dfdx, dfda, nmax=100, h=1.0, epsr=1.0e-10, epsb=1.0e-10)
bdc_r = np.array([bdc[m]['a'][0] for m in range(len(bdc))])
bdc_xmax = np.array([np.max(bdc[m]['x'][:-1:2*ndim]) for m in range(len(bdc))])
bdc_xmin = np.array([np.min(bdc[m]['x'][:-1:2*ndim]) for m in range(len(bdc))])
x = np.copy(y)
bdc2=pf.pathfollow00(x, a, func, dfdx, dfda, nmax=50, h=-1.0, epsr=1.0e-10, epsb=1.0e-10)
bdc2_r = np.array([pt['a'][0] for pt in bdc2])
bdc2_xmax = np.array([np.max(pt['x'][:-1:2*ndim]) for pt in bdc2])
bdc2_xmin = np.array([np.min(pt['x'][:-1:2*ndim]) for pt in bdc2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel("L")
ax.plot(bd_r, bd_x, '-')
ax.plot(bdc_r[:93], bdc_xmax[:93], '-k')
ax.plot(bdc_r[:93], bdc_xmin[:93], '-k')
ax.plot(bdc2_r[:31], bdc2_xmax[:31], '-k')
ax.plot(bdc2_r[:31], bdc2_xmin[:31], '-k')<jupyter_output><empty_output><jupyter_text>一様交通流から別の渋滞流に切り替える<jupyter_code>vr = bd[hb[1]]['vr']
vi = bd[hb[1]]['vi']
vr /= np.sqrt(vr@vr)
vi /= np.sqrt(vi@vi)
rr = np.hypot(vr, vi)
tt = np.array([np.angle(vr[i]+1.j*vi[i]) for i in range(len(vr))])
ev = eigvals(fx(bd[hb[1]+1]['x'], bd[hb[1]+1]['a']))
ev
ndim = 2*(Ncars-1)
x = np.zeros(ndim*Npts+1)
da = np.abs(bd[hb[1]]['a'][0] - bd[hb[1]+1]['a'][0])
for j in range(Npts):
    x[ndim*j:ndim*(j+1)] = bd[hb[1]+1]['x'] + 10 * np.sqrt(da) * rr * np.cos(2*np.pi*j/Npts + tt - tt[0])
x[-1] = 2*np.pi/np.imag(ev[2])
y = np.copy(x)
a = bd[hb[1]+1]['a']
for m in range(10):
    b = func(y, a)
    A = dfdx(y, a)
    y -= solve(A, b)
    print(np.linalg.norm(b))
# x = np.copy(y)
print(y)
x = np.copy(y)
bdc3=pf.pathfollow00(x, a, func, dfdx, dfda, nmax=100, h=1.0, epsr=1.0e-10, epsb=1.0e-10)
bdc3_r = np.array([bdc3[m]['a'][0] for m in range(len(bdc3))])
bdc3_xmax = np.array([np.max(bdc3[m]['x'][:-1:2*ndim]) for m in range(len(bdc3))])
bdc3_xmin = np.array([np.min(bdc3[m]['x'][:-1:2*ndim]) for m in range(len(bdc3))])
x = np.copy(y)
bdc4=pf.pathfollow00(x, a, func, dfdx, dfda, nmax=50, h=-1.0, epsr=1.0e-10, epsb=1.0e-10)
bdc4_r = np.array([pt['a'][0] for pt in bdc4])
bdc4_xmax = np.array([np.max(pt['x'][:-1:2*ndim]) for pt in bdc4])
bdc4_xmin = np.array([np.min(pt['x'][:-1:2*ndim]) for pt in bdc4])
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)
ax.set_xlim(0, 40)
ax.set_ylim(0,4)
ax.set_xlabel("L")
ax.plot(bd_r, bd_x, '-k')
ax.plot(bd_r[hb], bd_x[hb], 'xk', label="HB")
ax.plot(bdc_r[:93], bdc_xmax[:93], '--k')
ax.plot(bdc_r[:93], bdc_xmin[:93], '--k')
ax.plot(bdc2_r[:31], bdc2_xmax[:31], '--k')
ax.plot(bdc2_r[:31], bdc2_xmin[:31], '--k')
ax.plot(bdc3_r[:50], bdc3_xmax[:50], '--k')
ax.plot(bdc3_r[:50], bdc3_xmin[:50], '--k')
ax.plot(bdc4_r[:13], bdc4_xmax[:13], '--k')
ax.plot(bdc4_r[:13], bdc4_xmin[:13], '--k')
ax.legend()
# plt.savefig("bd_OV.pdf", bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.29315215928965427, 'fraction_numerical': 0.062962464684515, 'mean_word_length': 3.470234515935057, 'pattern_counts': {'":': 1, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '350348', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4318, 'n_tokens_neox': 4174, 'n_words': 688, 'path': 'notebooks/continuation/ov.ipynb', 'repo_name': 'tmiyaji/sgc164', 'size': '7433'}