스프링으로 연결되어 서로 평행한 직선 위에서 움직이는 일련의 질점

이번엔 멀쩡한 모델을 골랐기를 바랍니다. 

(얼핏 보면 되게 예쁘게 생김

 

구체적인 상황은 다음과 같습니다. 

직선은 4개가 있는 것으로 가정합니다. 따라서 질점도 4개입니다. 

중력장은 작용하지 않고, 각 직선이 서로 떨어진 거리는 d입니다. 

스프링의 원래 길이는 d이며, 스프링 상수는 k입니다. 인접한 질량 m의 질점을 서로 연결합니다. 

 

초기값으로 들어가는 것은 각 질점의 y방향 위치로 합니다. 

---

라그랑지안을 계산했을 때 매우 깔끔한 4가지 형태의 식을 얻습니다. 

(n개 직선일 상황으로 만들어보기도 쉬울 것 같습니다...만 4개만 하기로 합니다)

 

m=k=d=1을 가정합니다. 

초기값은 x, y, z, w로 하며, 그 속도 초기값은 모두 0으로 설정합니다. 

 

---

from matplotlib import animation, rc
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

m = 1
d = 1
k = 1

def dydt(u,t): 
    x = u[0]
    y = u[1]
    z = u[2]
    w = u[3]
    dx = u[4]
    dy = u[5]
    dz = u[6]
    dw = u[7]

    dydt0 = dx
    dydt1 = dy
    dydt2 = dz
    dydt3 = dw
    dydt4 = -k/m*( x-y - (x-y)/(1+(x-y)**2)**0.5 )
    dydt5 = -k/m*( y-x - (y-x)/(1+(x-y)**2)**0.5 ) -k/m*( y-z - (y-z)/(1+(z-y)**2)**0.5 )
    dydt6 = -k/m*( z-y - (z-y)/(1+(z-y)**2)**0.5 ) -k/m*( z-w - (z-w)/(1+(z-w)**2)**0.5 )
    dydt7 = -k/m*( w-z - (w-z)/(1+(w-z)**2)**0.5 )
    return [dydt0, dydt1, dydt2, dydt3, dydt4, dydt5, dydt6, dydt7]

frame = 200

y0 = np.array( [1, 0, -2, 0] )
y0 = np.hstack([y0, np.zeros(4)])
t = np.linspace(0, frame/10, frame)
sol = odeint(dydt, y0, t)
xx = sol[:,0]
yy = sol[:,1]
zz = sol[:,2]
ww = sol[:,3]

plt.plot(0 * np.ones(frame), xx, c='r')
plt.plot(1 * np.ones(frame), yy, c='g')
plt.plot(2 * np.ones(frame), zz, c='b')
plt.plot(3 * np.ones(frame), ww, c='k')
plt.show()

 

from matplotlib import animation, rc

fig, ax = plt.subplots(figsize=(5, 6))

ax.set_xlim(( -1, 4))
ax.set_ylim((-3, 3))


line1, = ax.plot([], [], color='gold', lw=6)
line2, = ax.plot([], [], color='gold', lw=6)
line3, = ax.plot([], [], color='gold', lw=6)
line4, = ax.plot([], [], color='gold', lw=6)
spring1, = ax.plot([], [], color='dimgray', lw=2)
spring2, = ax.plot([], [], color='dimgray', lw=2)
spring3, = ax.plot([], [], color='dimgray', lw=2)
mass1, = ax.plot([], [], mfc='w', marker='o', color='r', ms='15')
mass2, = ax.plot([], [], mfc='w', marker='o', color='g', ms='15')
mass3, = ax.plot([], [], mfc='w', marker='o', color='b', ms='15')
mass4, = ax.plot([], [], mfc='w', marker='o', color='k', ms='15')

def init():
    line1.set_data([], [])
    return (line1,)

dummy = np.linspace(0,49,50)
spring = 0.1*np.sin(dummy*2/7*np.pi)

def animate(i):
    x = xx[i]
    y = yy[i]
    z = zz[i]
    w = ww[i]
    mass1.set_data(0, x)
    mass2.set_data(1, y)
    mass3.set_data(2, z)
    mass4.set_data(3, w)
    
    dx = y-x
    lx = (1+dx**2)**0.5
    xsx = np.linspace(0,1,50)
    xsy = np.linspace(x,y,50)
    xsx -= spring*dx/lx
    xsy += spring/lx
    spring1.set_data(xsx, xsy)
    
    dy = z-y
    ly = (1+dy**2)**0.5
    ysx = np.linspace(1,2,50)
    ysy = np.linspace(y,z,50)
    ysx -= spring*dy/ly
    ysy += spring/ly
    spring2.set_data(ysx, ysy)
    
    dz = w-z
    lz = (1+dz**2)**0.5
    zsx = np.linspace(2,3,50)
    zsy = np.linspace(z,w,50)
    zsx -= spring*dz/lz
    zsy += spring/lz
    spring3.set_data(zsx, zsy)
    
    line1.set_data([0, 0], [-3, 3])
    line2.set_data([1, 1], [-3, 3])
    line3.set_data([2, 2], [-3, 3])
    line4.set_data([3, 3], [-3, 3])
    return (line1, line2, line3, line4, spring1, spring2, spring3, mass1, mass2, mass3, mass4)

anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=frame, interval=20, blit=True)
                               
rc('animation', html='html5')
anim

 

 

초기위치 1, 0, -2, 0

 

지그재그로 놓였을 때 다음과 같이 움직임을 알 수 있었습니다. 

 

 

 

초기위치 2, 0, 2, 0

위와 같이 조금 대칭되게 놓으면 모두 균일한 운동을 할 줄 알았더니, 이 상황에서 진동 양상은 고르지 않았습니다. 

복잡한 진동양상 : 다만 x와 w, y와 z는 대칭되는 운동

 

 

 

초기위치 3, 0, 0, 0

한쪽에만 변위가 주어졌을 때의 운동입니다. 

파동의 반사를 떠올리게 하네요. 

진동의 phase가 전달된다는 느낌의 구성이라 시간에 따른 변위 그래프도 상당히 재미있습니다.