python
1import numpy as np 2 3_list = [-1, 2, 0, -3, 5]4 5def zerobound(x):6 return max(x, 0)7 8 9_list = list(map(zerobound, _list))10 11 12# 各定数13V1, V2, C1, C2, Lc = 6.75, 7.91, 0.13, 1.57, 5.014λ = 0.515 16 17# -> dydx18def fvdm(v, s_diff):19 # 各定数20 V1, V2, C1, C2, Lc = 6.75, 7.91, 0.13, 1.57, 5.021 λ = 0.822 """ 23 微分方程式の定義(odeint) 24 f = [x, v] 25 """26 #インデックススライス27 # L = [0,1,2,3,4,5]28 # print(L[2])29 # print(L[:2])２まで（終了ポイント）30 # print(L[1:])１以降（スタートポイント）31 # print(L[::2])2個ごと（等間隔）32 33 def V(s):34 # V1, V2, C1, C2, Lc = 6.75, 7.91, 0.13, 1.57, 5.035 return V1 + V2 * np.tanh(C1 * (s - Lc) - C2)36 37 diff_v = np.append(v, v[0])38 # print(v[-1])39 diff_v = np.diff(diff_v)40 # print('diff_v', diff_v)41 42 """ 43 三項演算子:ifを1行で使える。 44 boolean（ブール）の頭文字は大文字(python) 45 46 result = True if 1 > 0 else Folse 47 1>0ならTrue、他ならFolse 48 49 L = [i for i in range(6)] 50 0~5までのリスト 51 52 """53 54 # s_diff = np.array([i if i > 0 else L + i for i in s_diff])55 56 57 zero_bound_v = V(s_diff) - v 58 # print(zero_bound_v)59 # print('boud v', zero_bound_v)60 # print(V(s_diff))61 # print(v)62 # zero_bound_v = np.array(list(map(lambda x: max(x, 0), zero_bound_v)))63 64 65 v_dot = K * zero_bound_v + λ*diff_v 66 return v_dot 67 68def fvdm_t(v, s_diff, t):69 dvdt = fvdm(v, s_diff)70 71 72 73 74def diff(s):75 s_diff = np.append(s, s[0])76 s_diff = np.diff(s_diff)77 s_diff %= L 78 s_diff = np.array([i if i > 0 else L + i for i in s_diff])79 return s_diff 80 81# t s v82def runge_kutta(v, s, h):83 # if h == k and h is None:84 # x, y = x + h/2, y + k /285 86 # k187 # dydx = func(x, y)88 s_diff = diff(s)89 dydx = fvdm(v, s_diff)90 # y2 = dydx91 k1 = dydx * h 92 dydx = fvdm(v + (k1/2), s_diff + (k1 / 2))93 k2 = dydx * h 94 k3 = fvdm(v + k2/2, s_diff + k2 /2) * h 95 k4 = fvdm(v + k3/2, s_diff + k3 /2) * h 96 97 y2 = (k1 + 2* k2 + 2* k3 + k4)98 y2 = y2 / 699 y2 = v + y2 100 # for j, i in enumerate(y2):101 # if i > 15:102 # print(j, s_diff[j], i)103 # print((k1 + 2* k2 + 2* k3 + k4) / 6)104 105 106 return y2 107 108 109def fvdm_solver(x, y, s, h=0.001):110 y2 = np.array([y])111 s2 = np.array([s])112 x2 = np.array([0])113 114 for i, _ in enumerate(np.arange(0, x, h)): # 0.1115 dy = runge_kutta(y2[i], s2[i], h)116 # print(dy)117 y = y + dy 118 119 120 s_diff = diff(s)121 # print(s_diff)122 # print(y)123 for j, items in enumerate(zip(y, s_diff)):124 if items[0] > items[1]:125 # print(items[0], items[1])126 if items[0] - items[1] < 0.0:127 print(items)128 raise Exception 129 130 131 # if items[i] < items[i] + s_diff132 if items[1] >= 0:133 y[j] = items[1]134 else:135 y[j] = 0136 137 138 if items[0] > 30:139 y[j] = 30140 if (items[0] < 0.0):141 y[j] = 0142 143 # if sum(y) > L:144 # y = (y / sum(y)) * L145 # print(y)146 # print(sum(y))147 148 # print(y)149 150 # print(y)151 for x, j in enumerate(y):152 if j < 0:153 print(x,j)154 raise Exception 155 156 s = s + y * h 157 s = np.array([i if i < L else i - L for i in s])158 159 if all(map(lambda x: x >= 0, diff(s))) == False:160 raise Exception 161 162 x += h 163 y2 = np.vstack([y2, y])164 s2 = np.vstack([s2, s])165 x2 = np.append(x2, x)166 return (x2, y2, s2)167 168#台数、、、距離169N, K, c = 20, 0.85, 2170L = N * 1171 172# 初期位置173y = np.arange(N) * L / N 174 175v = [0] + [10] * (N-1)176# v = np.array([np.random.randint(0, 10) for i in range(N)])177 178t = 100179 180h = 0.005181 182x2, y2, s2= fvdm_solver(t, v, y, h)183 184# for i, y in enumerate(s2):185# print(i, y[0], y[-1])186# for i, y in enumerate(y2):187# # print(max(y))188# print(i, y[0])189 190# for i in range(100):191# print(V1 + V2 * np.tanh(C1 * (i - Lc) - C2))192 193import matplotlib.animation as animation 194import matplotlib.pyplot as plt 195import seaborn as sbn 196 197r = L / (2 * np.pi)198fig = plt.figure(figsize=(6, 6))199sbn.set_style('ticks')200ims = []201 202for index, x in enumerate(s2):203 theta = x / r 204 im = plt.plot(r * np.sin(theta), r * np.cos(theta), '.k', markersize=24)205 ims.append(im)206 207Writer = animation.writers['ffmpeg']208# writer = Writer(fps=15, bitrate=1800)209writer = Writer(fps=60)210ani = animation.ArtistAnimation(fig, ims, interval=10)211 212plt.legend()213ani.save("output.mp4", writer=writer)214plt.show()
0 コメント