Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# coding: utf-8
# In[1]:
import numpy as np
import scipy as sp
# matplotlib:
from matplotlib.pylab import plt
get_ipython().magic('matplotlib inline')
# In[2]:
# 単純な二階微分方程式の状態空間:
# In[3]:
# 微分方程式: y'' = ay' + by
# 一階微分化: x2' = a*x2 + b*x1; x1=y, x2=y', x2=x1'
# 状態: [x1', x2'] = [[0,1], [a,b]] * [x1, x2]
# 出力: y = [1,0]*[x1, x2]
a,b = 0.2, 1.3
x1, x2 = 1.3, 0.6
vec = np.array([x1, x2])
mat = np.array([[0,1], [a,b]]) # x2が次期にx1に行く。今のx2にはax1+bx2。
y = []
for i in range(10):
vec = mat.dot(vec)
y.append(vec[0]) #x1をyとして出力。
plt.plot(y)
# In[4]:
# クラス化:
class SpaceState(object):
def __init__(self):
self.a = 0.2
self.b = 1.3
self.x1 = 1.3
self.x2 = 0.6
def calc(self):
vec = np.array([self.x1, self.x2])
mat = np.array([[0,1], [self.a, self.b]]) # x2が次期にx1に行く。今のx2にはax1+bx2。
y = []
for i in range(10):
vec = mat.dot(vec)
y.append(vec[0]) #x1をyとして出力。
return y
# In[5]:
aho1 = SpaceState()
plt.plot(aho1.calc())
aho2 = SpaceState()
aho2.x2 = 2.3
plt.plot(aho2.calc())
# In[6]:
# バネの振動:
# In[7]:
# 微分方程式: my’’(t) = -cy’(t) -ky(t) + u(t)
# 一階微分化: x2' = -c/m*x2 + -k/m*x1; x2=x1'
# 状態: [x1’,x2’] = [[0,1], [-k/m,-c/m]]*[x1,x2] + [0,1/m]*u
# 出力: y = [1,0]*[x1,x2]
k = 0.2
c = 1.3
m = 3
x1, x2 = 1.3, 0.6
vec = np.array([x1, x2])
u = 2.3
mat = np.array([[0,1], [-k/m, -c/m]])
u_mat = np.array([0, 1/m])
y = []
for i in range(10):
vec = mat.dot(vec) + u_mat.dot(u)
y.append(vec[0])
plt.plot(y)
# In[8]:
# クラス化:
class Spring(object):
def __init__(self):
self.time = 10
self.k = 0.2
self.c = 1.3
self.m = 3
self.x1 = 1.3
self.x2 = 0.6
self.u = 2.3
def calc(self):
vec = np.array([self.x1, self.x2])
mat = np.array([[0,1], [-self.k/self.m, -self.c/self.m]])
u_mat = np.array([0, 1/self.m])
y = []
for i in range(self.time):
vec = mat.dot(vec) + u_mat.dot(self.u)
y.append(vec[0])
return y
# In[9]:
aho1 = Spring()
aho1.m = 1.3
plt.plot(aho1.calc())
aho2 = Spring()
aho2.x2 = 2.1
plt.plot(aho2.calc())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.