중심력장 하에서 입자의 운동 (1) : 수치 미분 방정식 해법을 이용하여 문제 풀기 + 파이썬 시뮬레이션 코드

 

지난 포스팅 https://studyingrabbit.tistory.com/79

[고전역학 문제 풀이] 원형 빗면(콘, 깔때기)에서 구르는 입자의 운동

을 통하여 간략하게 중심력장 하에서 입자의 운동에 대해서 알아보았습니다. 지난 포스팅에서는 실제 시간에 따른 입자의 위치 $(r(t), \theta(t))$를 구하기 보다는 $r$에 대한 운동 방정식의 $V_{eff}(r)$의 개형을 통해 정성적인 입자의 움직임을 살펴보았습니다. 이번 포스팅에서는 지난 포스팅의 결과들을 활용하여 실제로 수치 미분 방정식의 해법을 이용하여 임의의 중심력장 하에서 입자의 시간에 따른 위치 $(r(t), \theta(t))$를 구해 보도록 하겠습니다. 

 

우리가 사용할 수치 미분 방정식 풀이 방법론은 지난 포스팅에서 여러번 언급한 바 있는 룽게-쿠타 방법입니다. 룽게-쿠타 방법에 대해서 잘 모르시는 분이라면, 이 포스팅의 뒷 부분을 읽기 보다는 룽게-쿠타 방법론에 대해서 다룬 이전 포스팅 https://studyingrabbit.tistory.com/55

룽게-쿠타 방법(Runge-Kutta method) 을 활용하여 보다 빠르고 정확하게 미분 방정식 문제를 수치 해법

을 먼저 보고 오는 것을 권장합니다.

 

지난 포스팅에서는 공간 자유도가 1차원인 시스템을 다루었는데, 주로 단조화 진동자 문제를 풀었습니다, 이번 문제는 공간 자유도가 2차원인 시스템입니다. 따라서 지난 포스팅에서 만들었던 파이썬 코드에서 공간에 대한 자유도를 하나 더 추가 시켜주어야 합니다. 그 방법은 단순 반복인데, 

k1x = velocity(t, x0, v0)
k1v = acceleration(t, x0, v0)

k2x = velocity(t+h/2, x0+h*k1x/2, v0+h*k1v/2)
k2v = acceleration(t+h/2, x0+h*k1x/2, v0+h*k1v/2)

지난 포스팅에서 위와 같은 코드를 이번 포스팅에서는 

k1x = velocity_x(t, x0, vx0, y0, vy0)
k1vx = acceleration_x(t, x0, vx0, y0, vy0)
k1y = velocity_y(t, x0, vx0, y0, vy0)
k1vy = acceleration_y(t, x0, vx0, y0, vy0)

k2x = velocity_x(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)
k2vx = acceleration_x(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)
k2y = velocity_y(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)
k2vy = acceleration_y(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)

위와 같이 $x, v_x$ 뿐 아니라, $y, v_y$에 대한 식을 추가해 주어야 합니다. 프로그램의 전체 코드 및 코드의 부분에 대한 설명은 이 포스팅의 맨 아래 부분에 있습니다.

 

위 코드에서는 $(r, \theta)$가 아니라 $(x, y)$를 이용하여 입자의 위치를 기술하였습니다. 중심력장 하에서 입자의 운동을 해석적으로 기술하기 위해서는 당연히 극좌표를 도입하는 것이 편리하지만, 단순 편한 코딩을 위해서는 직교좌표를 이용하는 것이 더 편리합니다. 물론, 보다 효율적이고 정교한 프로그래밍을 하기 위해서는 극좌표를 도입하는 것이 더 좋습니다. 그리고, 위 코드에서는 문제 풀이의 모든 부분은 수치 해법으로 처리하였습니다. 즉, $F_x = ma_x, F_y = ma_y$를 시작 부터 곧이 곧대로 수치적인 방법으로 푼 것 입니다. 역시나 보다 효율적이고 정교한 프로그래밍을 위해서는 해석적으로 간단히 할 수 있는 만큼 최대한 간단히 한 다음에, 해석적으로 풀 수 없는 부분만 수치적인 방법으로 푸는 코딩을 하는 것이 더 좋습니다. 특히 계산이 복잡하여 계산 시간이 오래 걸리는 문제일수록 효율적인 프로그램을 만들어야 합니다. 그러나 이 문제는 매우 간단한 문제이며, "머리를 쓰는(해석적인 방법으로 문제를 간단히 정리하는)" 과정을 최소화 하기 위해서 "처음 부터 곧이 곧대로" 수치적인 방법으로 문제를 푸는 코딩을 하였습니다. 

 

계산 조건

 

편의를 위해서 질량은 $1$로 두었습니다. 질량 $m$을 써야 하는 자리에 그냥 $1$을 쓴 것 입니다. 

결과를 표현한 에니메이션의 시간 스케일과 실제 문제의 시간 스케일은 일치하지 않습니다.

 

1. $\vec{F}(r) = -\frac{1}{r^2} \hat{r}$

 

가장 간단하면서도 중요한 문제인 케플러 문제를 풀어 보도록 하겠습니다. 케플러 문제는 (태양의) 중력으로 인한 행성의 움직임을 구하는 문제입니다. 매우 오래된 학문의 분과인 물리학 중에서도 가장 오래된 문제 중에 하나 입니다. 

 

1. 1 원운동

가장 간단한 원운동의 해 입니다. $x(0) = 1.0, y(0) = 0.0, v_x(0) = 0, v_y(0) = 1.0$ 입니다. 초기 위치의 원점으로 부터의 거리가 1이고, 원의 접선 방향 ($\hat\theta$) 방향의 속도가 1이기 때문에, 이 조건은 중력 = 구심력이 됩니다. 위 에니메이션에서 보는 바와 같이, 입자는 등속 원운동을 합니다. 

 

1.2 타원 운동

입자이 각운동량은 변화시키지 않고, 속도만 증가시켰습니다. 즉, 이번에는 같은 시작 위치에서, 초기 속도를 $v_y = 1.0, v_x = 0.2$와 같이 원점에서 멀어지는 방향($\hat{r}$)으로의 속도를 부여하였습니다. 원의 접선 방향 속도는 변화가 없기에, 각 운동량은 앞에서 본 원운동의 경우와 동일 합니다. 지난 포스팅에서 $V_{eff}(r)$를 이용한 정성적인 궤도의 분석의 결과를 활용하면, 이 조건에서의 운동은 $r_{min} \lt r \lt r_{max}$ 가 되고, 여기서 $r_{min} \ne r_{max}$이므로 타원 운동이 됩니다. 타원에는 두개의 초점이 있는데, 그 중 한 초점이 원점(힘의 중심)이 됩니다. 

 

케플러 문제에서 입자가 중력에 의해 구속된 상태에서의 궤적은 원 혹은 타원이라는 것을 아마 알고 계실 것 입니다(이걸 모르신다면, 이 포스팅을 읽기 보다는 케플러 문제에 대해서 다룬 고전역학 교과서를 보는 편이 더 좋습니다) 따라서 초기 속도를 바꾼다거나, 초기 위치를 바꾼다거나 했을 때, 입자의 역학적 에너지가 음수라서 중력에 의해 속박된 상태라면 항상 원 혹은 타원형의 궤적을 얻게 됩니다. 

 

원 혹은 타원의 특징은 닫힌 곡선이라는 것 입니다. 입자가 원점을 한 바퀴 돌게 되면, 즉 입자의 각도가 $0$에서  증가하여 $2\pi$로 되돌아 오면 입자의 위치 및 속도는 값은 항상 같은 값을 갖게 되는 것 입니다. 이는 케플러 문제의 아주 중요한 성질로, 이에 대한 해석적인 증명은 중심력장 운동의 해석적인 성질에 대해서 다룰 다음 포스팅에서 하도록 하겠습니다. 

 

입자의 초기 속도를 $v_y = 1.4$ 로 증가시켰습니다. 원운동에 비해서 입자의 운동에너지가 약 2배($1.4^2 = 1.96$) 정도 커진 조건이기에 원점에서 매우 멀리 떨어진 지점까지 도달할 수 있습니다. 그러나, 여전이 입자의 역학적에너지는 음수이기에 ($E = K + V = \frac{1}{2}1.4^2 + (-1) = -0.02$) 입자는 중심력장을 벗어 날 수 없고, 닫힌 곡선 궤도(타원) 운동을 하게 됩니다. 역학적 에너지가 커질 수록 타원이 이심률이 커지게 됩니다. 

 

1.3 중심력장을 벗어난 경우 : 포물선 혹은 쌍곡선 운동

입자의 초기 속도를 $v_y = 1.7$ 로 증가하여 역학적 에너지가 양수가 되도록 하였습니다. 이 경우 입자는 중심력장을 벗어나 무한히 먼 곳에 도달할 수 있습니다. 케플러 문제의 경우, 역학적 에너지가 정확히 0이 될 때의 궤적은 포물선, 역학적 에너지가 양수가 될 때의 궤적은 쌍곡선 입니다. 쌍곡선의 경우 원점에서 멀어지면 쌍곡선은 직선에 수렴하게 되는데, 중심력장의 영향은 거의 받지 않아서 입자는 등속 직선 운동에 가깝게 움직이게 됩니다. 

 

케플러 문제에서 입자의 궤적은 원추곡선으로 주어진다것을 해석적으로 증명할 수 있기에, 위에서 살펴본 경우는 해석적인 해를 수치적인 방법으로 재현한 것에 불과할 수 도 있습니다. 그러나, 단순히 해석적인 식을 종이에 써 보는 것과 실제로 시간에 따른 입자의 운동이 어떻게 되는지를 에니메이션을 통해서 보는 것은 분명한 차이가 있습니다. 이미 해석적인 답을 알고 있는 문제를 수치적인 방법을 통해서 확인하는 과정을 통해서, 우리의 코딩이 틀리지 않았다는 것을 간접적으로 보일 수 있습니다. 

 

2. $\vec{F}(r) = -\frac{1}{r^{3-\lambda^2}} \hat{r}$, $\lambda = \frac{p}{q} \lt 1$

 

케플러 문제를 일반화 하여, 거리의 일반적인 거듭 제곱에 비례하는 인력을 받는 경우의 문제를 풀어 보도록 하겠습니다. 위 식에서 $\lambda = 1$이 케플러 문제 입니다. 편의상 $\lambda \lt 1$로 제한하였는데, 이 힘은 중력(거리의 제곱에 반비례)에 비해서 원점에서 멀어질 수록 힘의 크기가 빠르게 0으로 수렴하는 힘이 됩니다. 

 

2.1 $\lambda = \frac{1}{2} \rightarrow \vec{F}(r) = -\frac{1}{r^{2.75}}\hat{r}$

 

$\lambda = \frac{1}{2}$ 경우 입니다. 초기 위치는 이번에도 역시 $x = 1, y = 0$으로 하겠습니다. 

 

2.1.1 원운동

 

$v_x = 0, v_y = 1$ 인경우, 원점 방향으로의 힘의 크기는 1이고, 원운동의 구심력 역시 1이므로 등속 원운동을 할 수 있습니다. 일반적으로 중심력장의 경우, 주어진 각 운동량에서 거리(위치)를 잘 잡으면 등속 원운동을 할 수 있습니다. 원운동의 시간에 따른 입자의 움직임은 케플러 문제의 경우 1.1와 완전히 동일하기 때문에 에니메이션은 생략합니다. 

 

2.1.2 원운동에서 약간의 운동 에너지가 더 해진 경우 : 원점을 두 바퀴 돈 후 제자리로 돌아오는 운동

반지름이 1, 속도가 1인 원운동과 비교했을 때, 각운동량의 크기는 같고, 에너지는 조금 높은 경우 입니다. 위 에니메이션과 같이 궤적은 원도 아니고 타원도 아닙니다. 입자는 원 주변을 2바퀴 돌고 (극좌표로 표현하면 $\theta$가 0 부터 $4\pi$만큼 증가하고 나서) 거의 제자리로 돌아와서 같은 운동을 반복합니다. 위 에니메이션이 궤적의 라인이 굵어지는 것을 봐서는 정확히 같은 궤적을 따라서 반복 운동하는 것은 아니고, 첫 번째 2바퀴 돌 때와 비교해서 거의 비슷한, 그러나 조금은 다른 궤적을 따라서 운동하게 됩니다. 중심력장을 기술하는데 쓴 $\lambda$가 $\frac{1}{2}$인데, 분모가 2이기 때문에 (거의) 원점을 2바퀴 돌고나서 제 자리로 돌아오는 궤적이 형성 됩니다. 

 

2.2 $\lambda = \frac{1}{3} \rightarrow \vec{F}(r) = -\frac{1}{r^{2.888}}\hat{r}$

$\lambda = \frac{1}{3}$ 으로 변경하였습니다. 이번에는 원점을 3바퀴 돌고나서 거의 같은 지점으로 돌아와서 반복 운동을 합니다. 역시 정확히 같은 궤적을 따라 움직이는 것은 아니고, 처음 3바퀴를 돌때와는 약간 다른 궤적입니다. 그러나 정성적으로 본다면 거의 같은 궤적을 따라 움직인다고 할 수 있습니다. 

 

일반적으로 $\lambda = \frac{p}{q}$ 와 같이 $\lambda$가 유리수 값을 갖고, 원운동에 비해서 아주 약간 어긋난 초기 조건 하에서 움직인다면, 입자는 원점을 $q$ 바퀴 돌고나서 초기 위치와 아주 유사한 곳으로 돌아오게 됩니다. 이 결과는 해석적인 방법으로 증명할 수 있으며, 이에 대해서는 다음 포스팅에서 하도록 하겠습니다. 2차원 조화진동자의 리사주 그림과 유사한 느낌이기도 합니다. 

 

3. $\vec{F}(r) = -\frac{1}{r^2} \hat{r} - \alpha \frac{1}{r^4} \hat{r}$

 

이번에는 $\frac{1}{r^2}$의 인력과 $\frac{1}{r^4}$에 비례하는 척력이 합해진 힘에 대해서 생각해 보겠습니다. 굳이 왜 이런 형태의 힘에 대해서 생각하는지 그 이유에 대해서 궁금해 하실 텐데, $\frac{1}{r^2}$는 중력에 의한 인력이고, $\frac{1}{r^4}$이 비례하는 인력 일반상대성이론에 의한 뉴턴의 중력의 보정항 입니다(가장 간단한 형태의 일반상대성이론의 해를 가정한 것 입니다). 즉, 중력이 매우 강한 (예를들어 태양-수성 정도의 짧은 거리) 영역에서는 상대성이론을 사용하여 정확한 중력을 계산해야 합니다. 뉴턴의 중력은 상대성이론에 따라 구해지는 정확한 중력의 첫 항이며, $\frac{1}{r^4}$이 비례하는 인력 상대론적 효과 입니다. 위에서 $\frac{1}{r^4}$의 계수를 그냥 $\alpha$로 놓았는데, 실제로는 $\alpha = \frac{3GML^2}{mc^2}$ 입니다. 여기서 $L$은 입자의 각 운동량 입니다.

 

3.1 $\alpha = 0.1$

 

편의상 $\alpha = 0.1$로 두고 계산해 보겠습니다. 이 경우에는 $r \gt 1$에서 첫 번째 항이 두 번째 항 보다 훨씬 더 커지게 되고, 두 번째 항을 일종의 섭동항으로 놓을 수 있습니다. 

$F = \frac{1}{r^n}$ 와 같이 간단한 경우에 비해서 궤적이 복잡해 졌습니다. 타원 비스무리 한 궤적을 따르는데, 하나의 닫힌 타원이 아니라 타원의 축이 이동하는 것 처럼 보입니다. 따라서 매 회전 마다 원점과 가장 가까워지는 점의 위치가 이동하게 되는데, 이를 세차 운동 이라고 합니다. 세차 운동의 주기와 원점을 한 바퀴 도는 운동의 주기의 비율이 무리수가 된다면, 입자의 궤적은 $r_{min} \lt r \lt r_{max}$ 영역의 공간을 빼곡하게 채우게 됩니다. 위 에니메이션에서는 에니메이션 파일의 용량 제한 상, 입자가 원점을 8바퀴 밖에 돌지 않습니다. 만일 오랜 시간 시뮬레이션을 한다면, $r_{min} \lt r \lt r_{max}$ 공간이 파란색으로 모두 채워질 것 입니다. 

 

실제로 태양과 매우 가까이에 있는 행성인 수성이 태향으로 부터 받는 중력을 계산할 때는 뉴턴의 중력 뿐 아니라, 일반 상대성이론에 대한 보정을 반드시 고려해야 합니다. 그 경우 위 에니메이션에서와 같이 세차 운동을 하게 됩니다. 실제로 수성은 세차 운동을 하는데, 세차 운동의 주기(크기)는 100년에 43초입니다. 여기서 초는 시간의 단위 초가 아니라 각도의 단위 초 인데, 각도의 단위 초는 $\frac{1}{3600}$ 도 입니다. 약 0.012도 입니다. 

 

4. $\vec{F}(r) = -r\hat{r}$

 

이번에는 인력이 거리 $r$에 비례하는 경우를 생각해 보겠습니다. 이는 2차원 단조화진동자 입니다. 2차원 단조화진동자 문제는 주로 직교 좌표계를 도입하여 $F_x = -kx, F_y = -ky$와 같이 높고, 각 차원 별로 자유도에 대해서 $x(t), y(t)$를 따로 구하는데, 이번에는 이를 $(r, \theta)$와 같이 극좌표계로 생각하여 문제를 푸는 것이라고 생각하면 됩니다. 

1, 2, 3에서 살펴본 거리가 멀어질수록 인력이 작아지는 힘의 경우, 충분한 운동 에너지가 주어지면 입자는 인력의 속박을 벗어날 수 있습니다. 하지만, 단조화진동자와 같이 원점에서 거리가 멀어질 수 록 힘이 커지는 경우에는 입자는 인력의 속박을 벗어 날 수 없습니다. 원점으로 부터 무한이 먼 곳으로 벗어나기 위해서는 무한이 큰 운동에너지가 필요하기 때문입니다.

 

위 에니메이션은 $(1,0)$의 초기 위치, $(0.3, 1)$의 초기 속도에 의한 운동입니다. 궤적은 타원입니다. 중력에 의한 운동에서와 같이 단조화 진동자의 운동의 궤적은 원 혹은 타원이 됩니다. 이는 간단하게 증명할 수 있으니, 설명은 생략합니다. $(x, y)$ 두 방향(임의의 수직한 두 방향)에 대해서 힘의 크기가 다르다면 리사쥬 도형과 같이 다양한 궤적이 가능하지만, 중심력장을 만족하는 단조화진동자의 경우, 항상 원 혹은 타원 궤도만 가능합니다. 

 

5. $\vec{F}(r) = -r^2\hat{r}$

 

이번에는 거리의 제곱에 비례하는 힘을 생각해 보겠습니다. 이번에도 역시 원 운동에서 약간 역학적 에너지가 증가한 경우의 궤적에 대해 계산해 보겠습니다.

앞에서와 같이 $(1,0)$의 초기 위치, $(0.3, 1)$의 초기 속도에 의한 운동의 궤적입니다. 앞에서 세차 운동의 경우 처럼 타원 비스무리한 궤적을 따라 움직이지만, 원점을 한 바퀴 돌 때 마다 타원의 축이 회전합니다. 또한 궤적은 $r_{min} \lt r \lt r_{max}$ 영역의 공간을 가득 채우려고 합니다. 

 

6. 정리

 

다양한 형태의 중심력장에 대한 입자의 궤적을 수치적 방법을 통해서 시뮬레이션 해 봤습니다. 원운동은 너무 단순하고, 매우 일반적인 운동은 또한 너무 일반적이고 복잡하기 때문에, 원운동의 조건에서 (각운동량은 고정 한채) 역학적 에너지를 조금 더 증가 시킨 초기 조건에 대해서 궤적을 알아 보았습니다. 중심력장 힘의 형태에 따라서 궤적은 타원이 되기도 하고, 세차 운동(혹은 좀 더 일반적인 운동)을 하기도 했습니다. 좀 더 다양한 초기 조건에 대해서 계산을 한다면 더욱 더 신기한 궤적이 얻어질 수 도 있습니다. 

 

위 방법은 단순 수치적인 계산이기 때문에, 매우 구체적인 조건에서의 궤적을 구하는 것은 가능하지만, 보다 일반적인 정량적인 혹은 정성적인 정리를 증명하는 것은 불가능 합니다. 예를들어 중력에 의한 속박된 운동의 궤도은 원과 타원만이 가능하다는 것을 해석적인 증명을 통해서 우리는 이미 알고 있습니다. 하지만, 이를 시뮬레이션만을 이용해서 증명하는 것은 불가능합니다. 앞서 언급한대로, 시뮬레이션은 매우 특정한 조건에 의한 운동의 궤도를 계산할 수 있고, 귀납적인 방법으로는 일반적인 정리를 증명할 수 없기 때문입니다. 

 

따라서, 중심력장의 운동에서 의미 있는 결과를 생각하고, 이를 증명하기 위해서는 반드시 해석적인 문제의 풀이가 필요합니다. 다음 포스팅에서는 해석적인 방법으로 이번 포스팅 시뮬레이션에서 관찰한 몇 가지 일반적이고 재미난 결과를 증명해 보도록 하겠습니다. 

 

7. 시뮬레이션에 사용한 파이썬 스크립트

 

룽게-쿠타 방법을 이용한 운동 방정식의 풀이 파이썬 코드 입니다. rungekutta_method 모듈에서 궤적 $(x(t), y(t))$를 계산하고, matplotlib 라이브러리를 이용하여 $(x(t), y(t))$를 그렸습니다. 잘 짜여진 코드는 아니지만, 이 코드를 활용하여 다양한 조건에서의 시뮬레이션을 해 볼 수 있습니다. 

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation

def velocity_x(t, x, vx, y, vy):
    return vx

def velocity_y(t, x, vx, y, vy):
    return vy

def acceleration_x(t, x, vx, y, vy):
    r = np.sqrt(x**2.0 + y**2.0)
    return -(x/r) *  1/ r**2
    
def acceleration_y(t, x, vx, y, vy):
    r = np.sqrt(x ** 2.0 + y ** 2.0)
    return -(y/r) * 1 / r**2
    
def rungekutta_method(x0, y0, vx0, vy0, t0, t1):
    ### time interval ###
    h = 0.01
    N = int((t1 - t0) / h)

    ### initial ###
    xs, vxs, ys, vys, ts = [x0], [vx0], [y0], [vy0], np.linspace(t0, t1, N)

    for i in range(1, N):
        ### equation for x ### : dx/dt = v

        t = t0 + h*(i-1)

        k1x = velocity_x(t, x0, vx0, y0, vy0)
        k1vx = acceleration_x(t, x0, vx0, y0, vy0)
        k1y = velocity_y(t, x0, vx0, y0, vy0)
        k1vy = acceleration_y(t, x0, vx0, y0, vy0)

        k2x = velocity_x(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)
        k2vx = acceleration_x(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)
        k2y = velocity_y(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)
        k2vy = acceleration_y(t+h/2, x0+h*k1x/2, vx0+h*k1vx/2, y0+h*k1y/2, vy0+h*k1vy/2)

        k3x = velocity_x(t+h/2, x0+h*k2x/2, vx0+h*k2vx/2, y0+h*k2y/2, vy0+h*k2vy/2)
        k3vx = acceleration_x(t+h/2, x0+h*k2x/2, vx0+h*k2vx/2, y0+h*k2y/2, vy0+h*k2vy/2)
        k3y = velocity_y(t+h/2, x0+h*k2x/2, vx0+h*k2vx/2, y0+h*k2y/2, vy0+h*k2vy/2)
        k3vy = acceleration_y(t+h/2, x0+h*k2x/2, vx0+h*k2vx/2, y0+h*k2y/2, vy0+h*k2vy/2)

        k4x = velocity_x(t+h, x0+h*k3x, vx0+h*k3vx, y0+h*k3y, vy0+h*k3vy)
        k4vx = acceleration_x(t+h, x0+h*k3x, vx0+h*k3vx, y0+h*k3y, vy0+h*k3vy)
        k4y = velocity_y(t+h, x0+h*k3x, vx0+h*k3vx, y0+h*k3y, vy0+h*k3vy)
        k4vy = acceleration_y(t+h, x0+h*k3x, vx0+h*k3vx, y0+h*k3y, vy0+h*k3vy)

        x1 = x0 + h*(k1x + 2*k2x + 2*k3x + k4x)/6
        vx1 = vx0 + h*(k1vx + 2*k2vx + 2*k3vx + k4vx)/6
        y1 = y0 + h*(k1y + 2*k2y + 2*k3y + k4y)/6
        vy1 = vy0 + h*(k1vy + 2*k2vy + 2*k3vy + k4vy)/6

        xs.append(x1)
        ys.append(y1)
        vxs.append(vx1)
        vys.append(vy1)

        x0 = x1
        y0 = y1
        vx0 = vx1
        vy0 = vy1

    return np.array(ts), np.array(xs), np.array(ys), np.array(vxs), np.array(vys)

def update(num, Z, point, line):
    point.set_offsets(Z[num])
    line.set_data(Z[:num, 0], Z[:num, 1])


if __name__ == '__main__':
    fig = plt.figure(figsize=(6, 6))
    ax = plt.axes()
    ax.grid()
    
    T, X, Y, VX, VY = rungekutta_method(1, 0, 0.3, 1.0, 0.0, 50)

    X, Y = X[::10], Y[::10]
    Z = np.hstack([X.reshape(-1,1), Y.reshape(-1,1)])
    N = len(X)

    xmax = max([np.abs(X).max(), np.abs(Y).max()])
    ax.set_xlim(-xmax*1.1, xmax*1.1)
    ax.set_ylim(-xmax*1.1, xmax*1.1)

    point = ax.scatter([], [], s=50)
    line, = ax.plot([], [])
    ax.scatter([0], [0], s= 100)

    ani = animation.FuncAnimation(fig, update, N, fargs=(Z, point, line), interval=10, blit=False)
    ani.save('a.gif', writer=animation.PillowWriter(fps=30))
    plt.show()