두개의 수평 용수철에 매달린 질점의 운동

뉴턴역학으로 하면 머리가 빠지므로 라그랑지안으로 쉽게쉽게 가도록 합니다. 

두 용수철이 초기 위치일 때의 길이를 l로 잡고 그 때 질점의 위치는 원점입니다. 

중력의 방향은 -y입니다. 

 

m=k=g=1, l=4를 기본적으로 가정합니다. 

초기 상태에서의 x, x', y, y'이 입력되는 초기값입니다.

 

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

k = 1
m = 1
g = 1
l = 4

def dydt(u,t): 
    x = u[0]
    y = u[1]
    dx = u[2]
    dy = u[3]
    dydt0 = dx
    dydt1 = dy
    dydt2 = -k/m*( 2*x - l*( (x+l)/((x+l)**2+y*y)**0.5 + (x-l)/((x-l)**2+y*y)**0.5 ) )
    dydt3 = -g -k/m*( 2*y - l*y*( 1/((x+l)**2+y*y)**0.5 + 1/((x-l)**2+y*y)**0.5 ) )
    return [dydt0, dydt1, dydt2, dydt3]

frame = 100

y0 = np.array( [1, 0, 0, 0] )
t = np.linspace(0, frame/10, frame)
y = odeint(dydt, y0, t)
xx = y[:,0]
yy = y[:,1]

plt.plot(xx, yy, c='r', mfc='w')
plt.vlines(-l, -10, 5, lw=5, color='gray')
plt.vlines(l, -10, 5, lw=5, color='gray')
plt.show()

 

위와 같이 미분방정식을 풀이하는 프로그램을 작성하면 아래와 같이 꽤나 준수한 모양의 그래프를 얻습니다. 

 

이제 이것을 용수철과 함께 묘화하여 애니메이션으로 넣어야 합니다. 

---

from matplotlib import animation, rc

fig, ax = plt.subplots(figsize=(4.1, 5))
ax.set_xlim(( -4.1, 4.1))
ax.set_ylim((-6.1, 4.1))

mass, = ax.plot([], [], mfc='w', marker='o', color='r', ms='15')
line_R, = ax.plot([], [], color='gray', lw=10)
line_L, = ax.plot([], [], color='gray', lw=10)
spring_L, = ax.plot([], [], color='k', lw=2)
spring_R, = ax.plot([], [], color='k', lw=2)

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

dummy = np.linspace(0,50,50)
spring = 0.3*np.sin(dummy*2/5*np.pi)

def animate(i):
    x = xx[i]
    y = yy[i]
    mass.set_data(x, y)
    
    dl = y/(x+4)
    ll = (1+dl**2)**0.5
    lsx = np.linspace(-4,x,50)
    lsy = np.linspace(0,y,50)
    lsx -= spring*dl/ll
    lsy += spring/ll
    spring_L.set_data(lsx, lsy)
    
    dr = y/(x-4)
    lr = (1+dr**2)**0.5
    rsx = np.linspace(x,4,50)
    rsy = np.linspace(y,0,50)
    rsx -= spring/lr*dr
    rsy += spring/lr
    spring_R.set_data(rsx, rsy)
    
    line_L.set_data([-4, -4], [-6, 4])
    line_R.set_data([4, 4], [-6, 4])
    return (line_L, line_R, spring_L, spring_R, mass)

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

(용수철의 묘사에는 조금의 수학적인 트릭이 들어가는데, 설명하지는 않겠습니다)

(x0,y0)=(1,0), (x'0,y'0)=(0,0)

초기 상태에서 오른쪽으로 조금 치우쳐있는 상황에서의 운동입니다. 

 

 

(x0,y0)=(0,1), (x'0,y'0)=(0,0)

초기상태에 조금 위쪽으로 치우친 경우의 운동입니다. 

수평방향의 운동은 없이 왕복 운동을 합니다. 

(주의 : 단진동이 아닙니다!)

 

 

(x0,y0)=(0,1), (x'0,y'0)=(4,4)

초기 상태에 오른쪽 위로 속도를 가지고 있는 상황입니다. 

 

이와 같이 다양한 모델들을 그려볼 수 있습니다. 

---

https://pinkwink.kr/1090

Matplotlib에서 그래프를 애니메이션으로 표현하기

애니메이션의 도입은 위 블로그에서 도움받았습니다. 

 

 

 

원점 시작, 초기 속도 (2,3)

TMI : 이 시스템의 궤적은 가끔 재미있는 패턴을 그립니다.

리사주 곡선이라고 착각할 수도 있겠지만 아니고, 보다 더 높은 복잡성을 띠는 곡선입니다. 

중력장이 없는, 원점 근방에서의 작은 운동만을 가정할 때 모든 해는 리사주 곡선으로 근사될 수 있습니다.