import numpy as np
import math
g = 9.81
N = 1000
tau = 20.0
dt = tau/float(N-1)
SHO
state[0] : 늘어난거리, 위치
state[1] : 속도
dv_dt 유도식
F = -mg -kx (탄성력)
F = ma (뉴턴 제2법칙 가속도법칙)
ma = -mg -kx
a = -g -(kx/m)
x(늘어난 거리)를 이용해서 가속도를 구해주는 것
euler 함수에 들어갈 [속도, 가속도]를 구해준다
def SHO(state, time):
dx_dt = state[1]
dv_dt = -k/m * state[0] - g
return np.array([dx_dt, dv_dt])
euler
1차 일반 미분 방정식의 해를 근사시켜준다…
결국 Yn+1을 Yn을 이용해서 구해준다는 이야기
[Xn + v dt, Vn + a dt] == [Xn+1, Vn+1]
def euler(y, t, dt, f):
return y + f(y, t) * dt
k : 탄성계수
m : 질량
x0 : 초기위치
v0 : 초기속도
k = 3.5
m = 0.2
x0 = 0
v0 = 0
t = np.linspace(0.0, tau, N)
y = np.zeros((N, 2))
y[0, 0] = x0
y[0, 1] = v0
for i in range(N-1):
y[i+1] = euler(y[i], t[i], dt, SHO)
xdata = [y[i, 0] for i in range(N)]
vdata = [y[i, 1] for i in range(N)]
from pylab import *
add plot position
subplot(2,1,1)
plot(t, xdata)
xlabel("time")
ylabel("position")
<matplotlib.text.Text at 0x22a3a2c0780>
add plot velocity
subplot(2,1,2)
plot(t, vdata)
xlabel("time")
ylabel("velocity")
show()