newton.pyfrom 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
newton2.pyfrom 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()
phasesp.pyfrom 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()
gibbs.pyfrom 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)
semigroup.pyfrom 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()