본문 바로가기
[硏 Ⅰ] 연구하다 - 교양/수치해석

벽에 용수철로 연결된 그릇 속, 두개의 용수철로 이어진 물체

by 천매 2021. 7. 9.

 

이번 상황은 매우 간단합니다. 

용수철의 원래 길이는 l로, 그릇의 두께는 없이, 폭은 d로 합니다. 

 


중간에 지운 흔적 : 처음 풀고 코드를 돌렸는데 해가 이상하게 나와서 식을 다시 수정...

실제로 해의 형태가 매우 좋게 나옵니다. 

조금의 시간만 들이면 수치해석적 방법 없이도 해석적인 해를 충분히 구할 수 있습니다.

(4계미분방정식 하나로 환원될 수 있기 때문)

 

하지만 수치해석적으로 구해봅시다. 

 


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

m = 1
k = 1
l = 2
D = 4

def dydt(u,t): 
    x = u[0]
    y = u[1]
    dx = u[2]
    dy = u[3]

    dydt0 = dx
    dydt1 = dy
    dydt2 = -k/m*(x -2*y -l +D)
    dydt3 = -k/m*( -x +4*y +l -2*D)
    return [dydt0, dydt1, dydt2, dydt3]

frame = 500

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

yy = xx + yy0

plt.plot(xx, 0 * np.ones(frame), c='r')
plt.plot(yy, 1 * np.ones(frame), c='g')

plt.show()
from matplotlib import animation, rc

fig, ax = plt.subplots(figsize=(1+l+D, 3.5))

ax.set_xlim(( -0.5, l+D+0.5))
ax.set_ylim((-0.5, 3))


line1, = ax.plot([], [], color='silver', lw=20)
line2, = ax.plot([], [], color='silver', lw=20)
line3, = ax.plot([], [], color='silver', lw=20)
line4, = ax.plot([], [], color='silver', lw=20)
line5, = ax.plot([], [], color='gray', lw=40)
spring1, = ax.plot([], [], color='dimgray', lw=2)
spring2, = ax.plot([], [], color='dimgray', lw=2)
spring3, = ax.plot([], [], color='dimgray', lw=2)
mass, = ax.plot([], [], mfc='khaki', marker='s', color='gold', ms='34')

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

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

def animate(i):
    x = xx[i]
    y = yy[i]

    mass.set_data(y, 1)
    
    spring1.set_data(np.linspace(0+h,x-h,50), np.ones(50)+spring)
    spring2.set_data(np.linspace(x+h,y-h,50), np.ones(50)+spring)
    spring3.set_data(np.linspace(y+h,x+D-h,50), np.ones(50)+spring)
    
    line1.set_data([0, 0], [0, 2.5])
    line2.set_data([x, x], [0.5, 1.5])
    line3.set_data([x, x+D], [0.5, 0.5])
    line4.set_data([x+D, x+D], [0.5, 1.5])
    line5.set_data([-0.5,l+D+0.5], [0, 0])
    return (line1, line2, line3, line4, spring1, spring2, spring3, mass)

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

이 모델의 재미있는 점은 여러가지가 있는데, 그 중 하나는 용수철의 길이 l과 접시의 폭 D를 잘 조절할 수 있다는 점입니다. 

참고로 D와 l을 어떻게 조절하든 접시 가운데가 안정평형이 나옵니다. 

다만 초기값이 모두 정지된 안정점이 아닌 이상, 평형에 도달하는 일은 없습니다... 

 

 

바로 아래 영상(위의 코드와 동일)에서는 용수철의 길이를 2로, 접시의 폭을 4로 했습니다. 

초기 상태에 접시와 벽을 잇는 용수철의 길이를 1, 접시 왼쪽 벽에서 물체에 달린 용수철의 길이는 3입니다. 

 


위 상황에서는 용수철의 원래 길이가 1입니다. 접시 내부 용수철이 거의 항상 팽팽한 상태에 있기는 합니다만, 그렇다고 해서 접시 중간에 물체가 고정되거나 하는 일은 없습니다...

 

(일부 값들을 테스트해보고 있는데 대부분 벽을 뚫습니다. 제대로 된 모델을 만들기 위해서는 벽에 닿았을 때 되돌아 탄성충돌을 하도록 식을 조절할 필요가 있습니다만... )


 

 

위 모델은 잔잔한 운동입니다.

용수철의 길이를 3으로 하였는데 평형 상태에서 얼마 떨어지지 않은 값을 초기치로 했습니다. 

 


 

 

위의 모델도 잔잔한 운동입니다. 위와 용수철의 길이는 같지만 조금 더 큰 진동을 합니다. 

 


마지막으로 가면 심심하니까 k=4로 두고(더 민감한 진동) 모델 하나의 푸리에 분석을 해봅니다. 

 

 

매우 정신없는 운동입니다. 

 

접시 왼쪽 끝과 물체의 변위만을 나타내면 다음과 같습니다.

 

이것을 fast Fourier transform 분석합니다. 

(x의 운동을 분석합니다)

xfft = np.fft.fft(xx)
f = np.fft.fftfreq(frame, 1/10)

xfft = xfft[f>0]
f = f[f>0]

plt.plot(f,np.abs(xfft)**2, c='plum', marker='o', \
         mfc='darkviolet', ms=4)
plt.loglog()
plt.xlabel('frequency')
plt.ylabel('amplitude')
plt.show()

원래 라그랑지안을 풀이한 연립미분방정식의 해의 형태를 볼 때 웬만해서는 봉우리가 단 2개만 나와야 합니다. 

다른 말로 하면 x와 y의 해는 서로 다른 진동수를 가진 삼각함수 2개의 합성이어야 했습니다. 

다른 말로 하면 위의 fft 분석 그래프에서 값 2개만 볼록 나와있고 나머지 값은 0이어야 했지만... 

 

봉우리의 진동수에서 주기값을 구해보면 대충 5초 주기와 1.5초 주기가 나옵니다. 

(4.444를 버리는 이유는 그것이 사실상 5.00에 흡수될 특이값이기 때문...)

이게 무엇을 의미하는고 하면, 위에 제시된 x와 y를 시간에 따라 그린 그림의 큰 형태를 결정하는 것은 주기가 5인 사인파이고, 그 다음으로 크게 변동을 결정하는 성분의 주기는 1.5이고... 라는 뜻입니다. 


접시 안의 물체 y의 운동도 마찬가지로 분석합니다. 

 

우선 접시의 위치를 기준으로 한 물체의 위치를 푸리에 분석합니다.

 

주요 봉우리는 1.5 부근임을 알 수 있습니다. 

 

여기에서 생각할 수 있는 점은, 접시 내부 용수철이 형성하는 진동 시스템의 주기는 1.5초 정도라는 것입니다. 

따라서 앞에서 분석한 접시의 운동에서 1.5초라는 것은, 바로 접시 내부 스프링에 의한 영향이 큽니다. 


 

그러면 y의 벽면에 대한 절대적인 위치를 푸리에 분석해봅니다. 

접시의 거시적인 움직임에 크게 지배받고 있다는 것을 알 수 있습니다. 

왜냐하면 처음에 분석한 그것의 형태와 푸리에 분석 결과가 동일하기 때문입니다.. 

 

아무튼 오늘은 여기서 마칩니다. 

 

 

참고로 올바른 node는 다음과 같이 나와야 합니다...

댓글