Python에서 간단하게 진자 운동 시뮬레이션을 애니메이션으로 구현하기

요즘은 회사안에서 개인의 만족도에 대한 생각들을 많이 하는 편입니다. 물론 지금의 일이 아주아주 재미있고 멋진데 사실 저는 좀 더 많은 일을 할 수 있으면 좋겠다는 생각을 하고 있거든요. 뭔가 의견만 내면 의도대로 잘 되지 않을때가 있으니 차라리 내가 책임지더라도 한번 끝까지 밀어부쳐보고 싶다는 생각을 하게 되죠. 그런데 이게 문제가 되는 것은 조직내에서의 이런 돌출 행동에 대한 시선과 또 내 이름이 들어갈 이 로봇이 정말 멋지게 완성되었으면 좋겠다는 순수한 생각이 한 50%, 그러면서 내가 좀 많은 부분을 할 수 있는 능력이 (있다는 것이 아니라) 있으면 좋겠다는 약간은 불순한 생각이 또 한 50%가 되면서 분명 오해의 소지가 있는거죠. 그러니... 이런 저런 생각으로 머리가 복잡할때는 그러면서 답이 없는 듯한 고민일 때는... 약간 머리도 식힐겸.. 이런 기초 지식에 가까운 글을 포스팅하고 있는 거겠죠???ㅎㅎ

오늘의 주제는 희한하게 진자(펜들럼 pendulum)와 Python, simulation, Runge Kutta로 표현되는데요. 요근래 계속 비슷한 글을 올리고 있죠.^^. 일단 처음에 저는 미분방정식을 수치해석적으로 푸는 방법 중 유명한 Runge Kutta 방법을 소개

[The Robot/Prog.Lang.] - MATLAB에서 4차 Runge Kutta를 이용하여 1차 혹은 2차 미분방정식을 푸는 예제

했었구요. 그리고 그 Runge Kutta를 실습할 수 예제로 Pendulum의 자유 낙하 운동을 대상으로 했지요. 물론 MATLAB에서도 그 예제를 다루었지만,

[The Robot/Prog.Lang.] - MATLAB에서 직접 2차 미방을 풀어 진자 운동 구현하기

그건 그저 제가 가장 잘 사용하는 tool이라 검증용으로 MATLAB을 이용했던거구요. 실제로는 Processing이라는 꽤 Rapid하게 그래픽적인 결과물을 얻을 수 있는 언어로 한 번 더 보여드렸구요.

[The Robot/Prog.Lang.] - Processing에서 진자 운동을 애니메이션으로 시뮬레이션하기

이제 오늘은 이런 미분 방정식 Runge Kutta + pendulum 시리즈를 종결할 Python으로 구현하기가 되는거죠^^

이미 이전 글에서 다루었던

이렇게 천장에 매달려 자유운동을 하는 진자를 시뮬레이션 하는 것이 목적인거죠.^^.

물론 우리가 대상으로 하는 운동역학 방정식은 위와 같구요. 이제 구현된 코드를 볼까요^^.

import pygame
from pygame.locals import *
import numpy as np

h = 0.05
t = v = 0
x = 30*np.pi/180
pen_fm = 0.01
pen_m = 0.1
pen_l = 100*0.01
pen_J = 0.02
pen_g = 9.8
gndCenterX = 150
gndConterY = 20
penLength = pen_l*100*2

먼저 2D에서 애니메이션을 표현하기 위한 많은 방법이 있겠죠. matplotlib를 이용해서 그냥 쉽게 그려도 되지만 살짝 Processing 같은 느낌을 주고 싶어서 pygame이라는 모듈을 사용했습니다. 그리고 numpy는 삼각함수 때문에 사용한거구요. 나머지는 파라미터들이죠. 아 그런데 마지막의 penLength는 실제 진자의 길이(pen_l)에 100배해서 애니메이션에서 잘 보이게 한 것입니다.

def calcODEFunc(tVal, xVal, vVal):
	return -pen_fm/(pen_m*pen_l*pen_l+pen_J)*vVal-pen_m*pen_g*pen_l/(pen_m*pen_l*pen_l+pen_J)*xVal

def solveODEusingRK4(t, x, v):
	kx1 = v
	kv1 = calcODEFunc( t, x, v )
	 
	kx2 = v + h*kv1/2
	kv2 = calcODEFunc( t + h/2, x + h*kx1/2, v + h*kv1/2 )
	 
	kx3 = v + h*kv2/2
	kv3 = calcODEFunc( t + h/2, x + h*kx2/2, v + h*kv2/2 )
	 
	kx4 = v + h*kv3
	kv4 = calcODEFunc( t + h, x + h*kx3, v + h*kv3 )
	 
	dx = h*(kx1 + 2*kx2 + 2*kx3 + kx4)/6
	dv = h*(kv1 + 2*kv2 + 2*kv3 + kv4)/6

	return x+dx, v+dv

먼저 Runge Kutta를 다룬 이전 글들에서와 동일하게 풀고자 하는 미방을 표현한 함수가 calcODEusingRK4구요. Runge Kutta를 풀고 있는 것이 solveODEusingRK4입니다. 이 이야기는 이미 이전에 MATLAB을 통해 설명[바로가기]을 했었지요.

pygame.init()

srf = pygame.display.set_mode((300,300))

font = pygame.font.SysFont('Vernada.ttf', 25)
aurthorSrf = font.render('by PinkWink', True, (50,50,50))

loopFlag = True

while loopFlag:
	for event in pygame.event.get(QUIT):
		loopFlag = False

	srf.fill((255,255,255))

	t = t + h
	[x, v] = solveODEusingRK4(t,x,v)

	updatedX = gndCenterX + penLength*np.sin(x)
	updatedY = gndConterY + penLength*np.cos(x)

	pygame.draw.line(srf, (100,100,100), (gndCenterX, gndConterY), (updatedX, updatedY), 2)
	pygame.draw.circle(srf, (100,100,100), (int(updatedX), int(updatedY)), 10, 0)

	pygame.draw.line(srf, (0,0,0), (10,20), (290,20), 10)
	srf.blit(aurthorSrf, (180,270))

	pygame.time.delay(40)
	pygame.display.flip()

이제 실제 애니메이션을 그리는 부분은 위 부분인데요. 초기화하는 부분이 있고, while문을 무한으로 돌리는 부분이 있지요. 여기서 저도 잘은 모르겠는데, pygame.event를 이용해서 종료 버튼이 눌러졌음을 확인하고 종료를 해야합니다. for event~라는 저 구문을 사용안했더니 그냥 x 표를 눌러도 꺼지지 않는 기이한 현상이 생기더군요. ㅎㅎ 아무튼 이렇게 구현하고 나면,

이런 실행 화면을 볼 수 있습니다. 전체 코드는

import pygame
from pygame.locals import *
import numpy as np

h = 0.05
t = v = 0
x = 30*np.pi/180
pen_fm = 0.01
pen_m = 0.1
pen_l = 100*0.01
pen_J = 0.02
pen_g = 9.8
gndCenterX = 150
gndConterY = 20
penLength = pen_l*100*2

def calcODEFunc(tVal, xVal, vVal):
	return -pen_fm/(pen_m*pen_l*pen_l+pen_J)*vVal-pen_m*pen_g*pen_l/(pen_m*pen_l*pen_l+pen_J)*xVal

def solveODEusingRK4(t, x, v):
	kx1 = v
	kv1 = calcODEFunc( t, x, v )
	 
	kx2 = v + h*kv1/2
	kv2 = calcODEFunc( t + h/2, x + h*kx1/2, v + h*kv1/2 )
	 
	kx3 = v + h*kv2/2
	kv3 = calcODEFunc( t + h/2, x + h*kx2/2, v + h*kv2/2 )
	 
	kx4 = v + h*kv3
	kv4 = calcODEFunc( t + h, x + h*kx3, v + h*kv3 )
	 
	dx = h*(kx1 + 2*kx2 + 2*kx3 + kx4)/6
	dv = h*(kv1 + 2*kv2 + 2*kv3 + kv4)/6

	return x+dx, v+dv

pygame.init()

srf = pygame.display.set_mode((300,300))

font = pygame.font.SysFont('Vernada.ttf', 25)
aurthorSrf = font.render('by PinkWink', True, (50,50,50))

loopFlag = True

while loopFlag:
	for event in pygame.event.get(QUIT):
		loopFlag = False

	srf.fill((255,255,255))

	t = t + h
	[x, v] = solveODEusingRK4(t,x,v)

	updatedX = gndCenterX + penLength*np.sin(x)
	updatedY = gndConterY + penLength*np.cos(x)

	pygame.draw.line(srf, (100,100,100), (gndCenterX, gndConterY), (updatedX, updatedY), 2)
	pygame.draw.circle(srf, (100,100,100), (int(updatedX), int(updatedY)), 10, 0)

	pygame.draw.line(srf, (0,0,0), (10,20), (290,20), 10)
	srf.blit(aurthorSrf, (180,270))

	pygame.time.delay(40)
	pygame.display.flip()

입니다. 구현되는 모습은

입니다. 요즘 이렇게 저는 약간의 취미같은 느낌으로 Processing이랑 Python을 틈틈이 정말 틈틈이 일주일에 한 1~2시간 정도 이렇게 공부하는데요. 어라~ 정말 재미있네요. ㅎㅎ. 아무튼 뭐 별 볼일 없는 주제지만 살짝 마무리 합니다. 아마 예약 발행될 날짜가 10월 31일 일텐데요. 미리미리 새로운 11월을 활기차게 보내자구요. ㅎㅎ^^