generated at
9. 力学系

ニュートンの運動方程式

newton.py
from vpython import * Ball = sphere(color=color.red) Wall = box(pos=vec(-10, 0, 0), length=1, width=10, height=10) Spring = helix(pos=vec(-10, 0, 0), length=10) dt, x, v = 0.01, 2.0, 0.0 while True: rate(1 / dt) dx, dv = v * dt, -x * dt x, v = x + dx, v + dv Ball.pos.x, Spring.length = x, 10 + x
ソースコード: newton.py



newton2.py
from numpy import arange import matplotlib.pyplot as plt def taylor_1st(x, v, dt): dx = v * dt dv = -x * dt return x + dx, v + dv def taylor_2nd(x, v, dt): dx = v * dt - x / 2 * dt ** 2 dv = -x * dt - v / 2 * dt ** 2 return x + dx, v + dv update = taylor_1st # taylor_2nd dt = 0.1 # 0.01, 0.001 path = [(2.0, 0.0)] # (x, v) for t in arange(0, 500, dt): x, v = path[-1] path.append(update(x, v, dt)) plt.plot(*zip(*path)) plt.axis('scaled'), plt.xlim(-3.0, 3.0), plt.ylim(-3.0, 3.0) plt.show()
ソースコード: newton2.py

テーラー展開の1次近似(dt = 0.1, 0.01, 0.001)

テーラー展開の2次近似(dt = 0.1, 0.10, 0.001)


線形の微分方程式
phasesp.py
from numpy import array, arange, exp, sin, cos from numpy.random import uniform import matplotlib.pylab as plt def B1(lmd1, lmd2): return lambda t: array([[exp(lmd1 * t), 0], [0, exp(lmd2 * t)]]) def B2(lmd): return lambda t: exp(lmd * t) * array([[1, 0], [t, 1]]) def B3(a, b): return lambda t: exp(a * t) * array([ [ cos(b * t), sin(b * t)], [-sin(b * t), cos(b * t)]]) B = B1(1, 0) # B1(lmd1, lmd2), B2(lmd), B3(a, b) V = uniform(-1, 1, (100, 2)) T = arange(-10, 10, 0.01) [plt.plot(*zip(*[B(t).dot(v) for t in T])) for v in V] plt.axis('scaled'), plt.xlim(-1, 1), plt.ylim(-1, 1) plt.show()
ソースコード: phasesp.py

\begin{bmatrix}e^{t\lambda_1}&0\\0&e^{t\lambda_2}\end{bmatrix}の場合:

\lambda_1=1, \lambda_2=0

\lambda_1=1, \lambda_2=1

\lambda_1=1, \lambda_2=2

\lambda_1=1, \lambda_2=-1



\begin{bmatrix}e^{t\lambda}&0\\1&e^{t\lambda}\end{bmatrix}の場合:

\lambda=0

\lambda=1


\begin{bmatrix}a&b\\-b&a\end{bmatrix}の場合:

a=0, b=1

a=1, b=1


マルコフ確率場
gibbs.py
from numpy import random, exp, dot from tkinter import Tk, Canvas class Screen: def __init__(self, N, size=600): self.N = N self.unit = size // self.N tk = Tk() self.canvas = Canvas(tk, width=size, height=size) self.canvas.pack() self.pallet = ['white', 'black'] self.matrix = [[self.pixel(i, j) for j in range(N)] for i in range(N)] def pixel(self, i, j): rect = self.canvas.create_rectangle x0, x1 = i * self.unit, (i + 1) * self.unit y0, y1 = j * self.unit, (j + 1) * self.unit return rect(x0, y0, x1, y1) def update(self, X): config = self.canvas.itemconfigure for i in range(self.N): for j in range(self.N): c = self.pallet[X[i, j]] ij = self.matrix[i][j] config(ij, outline=c, fill=c) self.canvas.update() def reverse(X, i, j): i0, i1 = i - 1, i + 1 j0, j1 = j - 1, j + 1 n, s, w, e = [X[i0, j], X[i1, j], X[i, j0], X[i, j1]] nw, ne, sw, se = [X[i0, j0], X[i0, j1], X[i1, j0], X[i1, j1]] a = X[i, j] b = 1 - 2 * a intr1 = b intr20 = b * sum([n, s, w, e]) intr21 = b * sum([nw, ne, sw, se]) intr3 = b * sum([n * ne, ne * e, e * n, e * se, se * s, s * e, s * sw, sw * w, w * s, w * nw, nw * n, n * w]) intr4 = b * sum([n * ne * e, e * se * s, s * sw * w, w * nw * n]) return intr1, intr20, intr21, intr3, intr4 N = 100 beta = 1.0 I = [-4.0, 1.0, 1.0, 0.0, 0.0] scrn = Screen(N) X = random.randint(0, 2, (N, N)) while True: for i in range(-1, N - 1): for j in range(-1, N - 1): S = reverse(X, i, j) p = exp(beta * dot(I, S)) if random.uniform(0.0, 1.0) < p / (1 + p): X[i, j] = 1 - X[i, j] scrn.update(X)
ソースコード: gibbs.py

\beta = 1.0; I = [-4.0, 1.0, 1.0, 0.0, 0.0]

\beta = 2.0; I = [0.0, 1.0, -1.0, 0.0, 0.0]

\beta = 4.0; I = [-2.0, 2.0, 0.0, -1.0, 2.0]

\beta = 1.5; I = [-2.0, -1.0, 1.0, 1.0, -2.0]


1径数半群

semigroup.py
from numpy import arange, array, eye, exp from numpy.random import choice, seed from numpy.linalg import eig import matplotlib.pyplot as plt seed(2020) dt, tmax = 0.01, 1000 T = arange(0.0, tmax, dt) G = array([[-3, 4, 0], [ 1, -4, 5], [ 2, 0, -5]]) v = eig(G)[1][:, 0] print(v / sum(v)) dP = eye(3) + G * dt X = [0] S = [[dt], [], []] for t in T: x = X[-1] y = choice(3, p=dP[:, x]) if x == y: S[x][-1] += dt else: S[y].append(dt) X.append(y) plt.figure(figsize=(15, 5)) plt.plot(T[:2000], X[:2000]) plt.yticks(range(3)) fig, axs = plt.subplots(1, 3, figsize=(20, 5)) for x in range(3): s, n = sum(S[x]), len(S[x]) print(s / tmax) m = s / n axs[x].hist(S[x], bins=10) t = arange(0, 3, 0.01) axs[x].plot(t, exp(-t / m) / m * s) axs[x].set_xlim(0, 3), axs[x].set_ylim(0, 600) plt.show()
ソースコード: semigroup.py

状態変化

滞在時間の分布