9. 力学系 - linalg

9. 力学系

Pythonで学ぶ線形代数学

ニュートンの運動方程式

code: 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

https://gyazo.com/5ffcd45d35ced894d1765aba5714b88b

https://gyazo.com/ceadd79fd38cbd20a9e002b2044d8391

code: 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)

https://gyazo.com/1614acb51e1473449342700e2759a059

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

https://gyazo.com/6db77a64e75becabb6375e01f06ccdbc

線形の微分方程式

code: 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

https://gyazo.com/932fa7464dfd522f5e1eb747ddaa37ff

$ \lambda_1=1, \lambda_2=1

https://gyazo.com/0edcf09e79c2f8fdb5b10ff7242bdc4f

$ \lambda_1=1, \lambda_2=2

https://gyazo.com/4c6bfe06f8c96cc103863be291c8bda3

$ \lambda_1=1, \lambda_2=-1

https://gyazo.com/fd9c1a980106c58af9fb08e5b73515d7

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

$ \lambda=0

https://gyazo.com/dff0026c0d78ac1bdd14ac0dbed8854d

$ \lambda=1

https://gyazo.com/7268d5fed4fe11fd9a62169434450c60

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

$ a=0, b=1

https://gyazo.com/636de5b49de6c971b768992d20b0b3dc

$ a=1, b=1

https://gyazo.com/67282c47b29e41e94ec186d98ae349dc

マルコフ確率場

code: 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[Xi, j]

ij = self.matrixij

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 = [Xi0, j, Xi1, j, Xi, j0, Xi, j1]

nw, ne, sw, se = [Xi0, j0, Xi0, j1, Xi1, j0, Xi1, j1]

a = Xi, 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):

Xi, j = 1 - Xi, j

scrn.update(X)

ソースコード: gibbs.py

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

https://gyazo.com/3a7a8aa8233872ee6c1de15c42fbc319

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

https://gyazo.com/f3f35c98398f38a2c70477a358396802

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

https://gyazo.com/a1faceaff9c61b8ad74dff4c11ef9bc9

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

https://gyazo.com/6a00b68dda8a0ccd0a2128a74a54a7c6

1径数半群

code: 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:

Sx-1 += dt

else:

Sy.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(Sx), len(Sx)

print(s / tmax)

m = s / n

axsx.hist(Sx, bins=10)

t = arange(0, 3, 0.01)

axsx.plot(t, exp(-t / m) / m * s)

axsx.set_xlim(0, 3), axsx.set_ylim(0, 600)

plt.show()

ソースコード: semigroup.py

状態変化

https://gyazo.com/b17ad1342af39e08411cce359dc3b4fa

滞在時間の分布

https://gyazo.com/fdbfe6d18763da10af3061e116673adf