blog/content/posts/2021-05-21-pyphy-mulrienjin...

3.9 KiB

+++ date = 2021-05-21T13:18:27Z description = "" draft = false slug = "pyphy-mulrienjin-1pyeon-gibon-aidieo" title = "파이썬으로 구현한 뉴턴의 방정식과 적분방정식" math = true +++

우리의 세계에서 물체의 역학적 운동을 지배하는, 지배방정식이 있죠.뉴턴의 법칙, 방정식입니다.뉴턴의 법칙은 총 세개의 방정식이 있죠. 이 세가지 법칙만을 이용해서 물체의 운동을 계산할 수 있습니다.이것을 바탕으로 만들어진 역설이 라플라스의 역설이죠. 만약 모든 분자의 가속도, 속도, 위치 세가지 정보를 알고 있다면, 뉴턴의 법칙에 근거하여 그 다음 상황을 예측할 수 있고, 따라서 이 세상, 아니면 그것을 넘어서서 모든 분자의 정보를 알 수 있습니다.우리가 여기서 분석할 것은 가장 쉬운 법칙, 뉴턴의 제 2법칙, 가속도의 법칙입니다.뉴턴의 가속도 법칙은 F = \frac{dp}{dt} 로 표현할 수 있습니다. 이 법칙에 의하면, 질량의 변화가 없는 물체에서, a F 에 비례합니다. 따라서 각 순간 물체에 작용하는 힘의 크기를 안다면, 그 물체의 가속도를 알 수 있겠죠?

이때, 가속도와 속도, 변위의 정의를 이용하면 a = \frac{dv}{dt} 이고, v = \frac{dx}{dt}입니다.가속도를 안다면 적분을 통해 속도를 알 수 있고, 속도를 안다면 변위를 알 수 있죠.정상상태에서는 물체에 작용하는 힘의 방정식을 구할 수 있습니다. 이것을 질량으로 나눠 가속도를 얻을 수 있고, 속도와 변위의 방정식을 얻을 수 있습니다.

하지만, 컴퓨터는 대수적으로 적분할 수 없습니다. (물론 할 수 있는 알고리즘도 있지만 불완전하죠)그렇다면 어떻게 각 시점에서 가속도의 값을 알고, 속도와 변위의 값을 알 수 있을까요?이렇게 미분방정식을 통해서 각 시점의 이산적인 변수값을 알 수 있도록 하는 것이 오일러 근사법입니다.오일러 근사법을 통틀어 이렇게 수치적으로 방정식을 해석하는 학문을 통틀어 수치해석학이라고 합니다.이렇게 물리를 시뮬레이션 하거나 컴퓨터로 해를 구할 때는 수치해석학을 사용합니다.

오일러 방정식이란,

 \frac{dy}{dx} = f(x) 
 y_{n+1} = y_n + \frac{dy}{dx}|_{x_n} \Delta x 

이 식을 이용하면 특정 원점으로부터 \Delta x만큼의 일정한 간격을 띄어가며 y 를 근사할 수 있습니다.

이 방법으로 물체에 작용하는 힘을 알 때 변위와 속도를 구할 수 있습니다. 특정 \Delta t를 기준으로 \Deltan배만큼의 시점에서 구할 수 있습니다.

 v_{n+1} = v_n + \frac{dv}{dt}|_{t_n} \Delta t 
 x_{n+1} = x_n + \frac{dx}{dt}|_{t_n} \Delta t 

으로 말이죠..

이제 파이썬으로 오일러 근사법을 구현해 봅시다.

이렇게 데이터의 배열끼리의 연산을 취급할 때에는 파이썬의 수치분석 모듈인 Numpy를 사용합니다. 하지만, 간단한 시뮬레이션에서는 파이썬 기본으로도 충분합니다. 그래프를 표현하기 위해 matplotlib과 수식을 위해 math를 이용합니다.

import matplotlib.pyplot as plt 
import math 
x0 = 0.0
y0 = 1.0
dx = 0.01
step = 1000
dydx = lambda x, y: 1/y 
x = [x0 + dx*i for i in range(step)] 
y = [y0] 

for n in range(step-1): 
	y.append(y[n] + dydx(x[n],y[n]) * dx) 

plt.plot(x, y) 
plt.grid() 
plt.xlabel("x") 
plt.ylabel("y") 
plt.legend() 
plt.show()

{{< figure src="https://blog.kakaocdn.net/dn/KMPAl/btq9fmV0ntZ/HPdiTkaknSISR4yGvxKUSK/img.png" caption="\frac{dy}{dx} = \frac{1}{y}" >}}