Verlet 적산법(integration, 계속해서 상태 변화가 있고 그 값을 누적해 더한다는 말)은 Wikipedia 글의 첫 줄에 밝힌 것처럼 뉴턴(Newton)의 운동 방정식을 적산하기 위한 수치 해석법이다. 이미 18세기부터 사용되었지만 1960년대에 프랑스 사람인 Verlet가 분자 운동학에서 많이 사용했다 해서 Verlet 적산법이라고 많이 부른다. (Verlet의 발음은 프랑스어라 익숙치 않지만 [vɛʁˈlɛ], 즉 베흘레와 유사한 것 같다. 외래어 표기법상으로는 베를레가 맞을 것 같다.)
이미 HTML5가 나오면서 canvas에서 애니메이션 등 다양한 그림 기법이 시도됐는데 최근에는 천이 찢어지는 애니메이션, 거미줄의 흔들림 등 물리적 운동 효과를 canvas에서 보이고 있다. 세상에나! JavaScript로 물리 운동을 시뮬레이션하다니!
기본 개념
그래서! 한 번 공부해보기로 했다. 우선 중학교만 나오면 다 아는 뉴턴의 운동 방정식을 단순하게 쓰자면 다음과 같다.
$$ x = vt $$
즉, 어떤 입자의 위치 $x$ 는 속도와 시간의 곱으로 구할 수 있다. 그런데 이걸 컴퓨터에 맞게 적산하려면, 즉 시간을 잘개 쪼개 매순간 $\Delta t$씩 흘러가는 동안 $i$번째 시간의 위치를 기준으로 $i + 1$번째 시점의 위치는 수치 해석적으로는 다음과 같이 쓴다.
$$ x_{i+1} = x_i + v_i \Delta t $$
반복해서 이 다음 $i + 2$번째 시점에는 다음과 같이 될 것이다.
$$ x_{i+2} = x_{i+1} + v_{i+1} \Delta t $$
Verlet 적분법은 여기서 속도를 구하기가 힘들거나 번거롭다는 데 주목했다. 대부분 중력 가속도와 같은 상수 값에 의해 속도가 결정되는데 매순간 속도를 구해야 한다는 건 번거로운 일이 아닐 수 없다. 그래서 $v_{i+1} = v_i + a_i\Delta t$, 즉 속도는 가속도에 시간을 곱해 구하므로 다음과 같이 전개할 수 있다.
$$ x_{i+2} = x_{i+1} + v_i\Delta t + a_i\Delta t^2 $$
여기서 $v_i$를 제거하기 위해 $x_{i+1} = x_i + v_i\Delta t$를 양변에서 빼주면 최종적으로 다음 공식이 도출된다.
$$ x_{i+2} = 2x_{i+1} - x_i + a_i\Delta t^2 $$
즉 i+1번째를 현재 시점으로 본 다면 그 다음 시점인 i+2번째의 위치는 현재 시점과 지난 시점의 위치, 가속도만 가지고 구할 수 있다.
사실 이것은 수치해석법이라 정확하게 미적분을 적용한 게 아니므로 위 식이 수학적으로 정확한 것은 아니다. 그러나 위 식에 대한 증명을 찾아보면 그 오차는 4차수(4th order) 오차이므로 큰 영향이 없고 컴퓨터에서 구현할 때 그 정도는 무시할 수 있게 된다.
구속 앨거리듬
구속은 입자간의 거리에 대한 조건으로 말하자면 입자간에 끈이 연결돼서 입자 운동을 제한하는 것이라 보면 된다. 가장 단순한 경우는 끈이고 좀더 복잡하게 한다면 스프링으로 연결한 것이 된다. 단순 연결이라고 본다면 이 Verlet 방법에서 구속에 따른 입자의 위치 계산은 단순히 위 (1)에서 두 점의 위치를 계산했을 때 그 두 점 사이의 거리가 원래의 거리와 차이가 난 것을 두 점에 "반땅"해주는(노나주는) 것이다.
좀더 수학적으로 말하자면 p1, p2라는 두 입자가 있을 때 위 (1)에 따른 거리가 d, 원래 거리가 l이라면 d - l 값을 반으로 갈라 p1, p2의 위치에 더하거나 빼서 조정해주는 것이다. p2의 좌표 값이 p1보다 크고 d > l이라면 p2의 좌표를 예로 들 경우 다음과 같이 조정하게 된다.
$$ x_{i+2} = x_{i+2} + \frac{d - l}{2} $$
지금 말한 건 1차원인 경우고 2차원인 경우는 $d - l$ 값을 그대로 노나줄 수는 없고 가로 좌표값과 세로 좌표값을 분리해서 비율적으로 노나줘야 한다.
분자 구조에 대한 그림에서 입자간의 연결을 흔히 스프링으로 해놓은 것을 기억할지 모르겠다. 만약 단순한 끈이 아니라 스프링으로 연결된 것으로 하자면 위에서 노나주는 값을 스프링의 탄성 비율에 따라 조정해야 할 것이다.
예시
자, 이제 실제로 뭔가 만들어봐야겠는데 밤 늦도록 공부만 실컷했고 하니 그 전에 좀 감을 잡기 위해 위 식에 따라 값들이 어떻게 변하는지 한번 계산해보자. 구속은 일단 다음으로 미루고 위치 변화만 살펴본다.
대부분의 물리 운동 시뮬레이션에서 가속도는 중력가속도면 충분하겠지만 Verlet 적산법이 얼마나 정확한지 보기 위해 가속도가 시간에 따라 변하는 함수인 경우를 적용해보자. 예를 들어 가속도를 $a_t = 10\cos(t)$라고 정의한다면 입자의 위치의 정확한 해는 두 번 적분하면 나온다.
- 가속도 $a_t = 10\cos(t)$
- 속도 $v_t = 10\sin(t)$
- 변위 $x_t = 10 - 10\cos(t)$
위의 공식을 0.1초씩 변경해가면서 구한 정확한 값과 Verlet 적산법으로 구한 값을 비교해보면 다음과 같다. 마지막 값을 보면 큰 차이가 없음을 볼 수 있다.
| 시간 $t$ | 가속도 $a_t$ | 속도 $v_t$ | 위치 $x_t$ | Verlet 위치 $x_{i+2}$ |
|---|---|---|---|---|
| 0 | 10 | 0 | 0 | 0 |
| 0.1 | 9.950041653 | 0.998334166 | 0.049958347 | 0.049958347 |
| 0.2 | 9.800665778 | 1.986693308 | 0.199334222 | 0.199417111 |
| 0.3 | 9.553364891 | 2.955202067 | 0.446635109 | 0.446882532 |
| 0.4 | 9.21060994 | 3.894183423 | 0.78939006 | 0.789881603 |
| 0.5 | 8.775825619 | 4.794255386 | 1.224174381 | 1.224986773 |
| 0.6 | 8.253356149 | 5.646424734 | 1.746643851 | 1.747850199 |
| 0.7 | 7.648421873 | 6.442176872 | 2.351578127 | 2.353247186 |
| 0.8 | 6.967067093 | 7.173560909 | 3.032932907 | 3.035128393 |
| 0.9 | 6.216099683 | 7.833269096 | 3.783900317 | 3.78668027 |
| 1 | 5.403023059 | 8.414709848 | 4.596976941 | 4.600393144 |
| 1.1 | 4.535961214 | 8.912073601 | 5.464038786 | 5.468136248 |
| 1.2 | 3.623577545 | 9.32039086 | 6.376422455 | 6.381238965 |
| 1.3 | 2.674988286 | 9.635581854 | 7.325011714 | 7.330577457 |
| 1.4 | 1.699671429 | 9.8544973 | 8.300328571 | 8.306665832 |
| 1.5 | 0.707372017 | 9.974949866 | 9.292627983 | 9.299750922 |
| 1.6 | -0.291995223 | 9.99573603 | 10.29199522 | 10.29990973 |
| 1.7 | -1.288444943 | 9.916648105 | 11.28844494 | 11.29714859 |
| 1.8 | -2.272020947 | 9.738476309 | 12.27202095 | 12.281503 |
| 1.9 | -3.232895669 | 9.463000877 | 13.23289567 | 13.24313719 |
| 2 | -4.161468365 | 9.092974268 | 14.16146837 | 14.17244244 |
$(14.17244244 - 14.16146837) / 14.16146837 = 0.000775$, 즉 약 $0.08%$ 오차 정도가 난다. 다른 가속도로도 시험해봤는데 가속도가 선형이거나 상수인 경우는 아예 오차가 없다(앞서 말한 것처럼 오차는 4차수 오차이므로 선형이거나 상수인 경우 4차 미분 값은 0이기 때문이겠지). 일반적으로 사용하는 경우 중력가속도처럼 상수 가속도를 사용하므로 Verlet 적산법이면 물리 운동 시뮬레이션으로는 아주 충분하다는 걸 알 수 있다.
초기값 계산
그런데 위에서 사실은 Verlet 계산 결과라고 한 열의 첫번째, 두번째 값은 상수를 넣은 것이다. 위 (1)번 공식에서 현재 위치와 그 전 위치를 알아야 다음 위치를 구할 수 있는데 첫번째, 두번째 위치는 그 전 위치가 없으므로 구할 수 없는 것이다.
그래서 일반적으로는 두번째 위치값을 첫번째 위치값과 같은 값으로 넣어준다(즉, 두 값 다 0으로 시작). 그러면 위 예시에서의 오차가 7%가 나오긴 하는데 다시 0.1초 후에는 그 위치를 따라 잡게 되고 더욱이 실제 애니메이션에서는 시간 간격이 0.1초보다 훨씬 작거나(30~60분의 1초 정도) 전체 시간이 상대적으로 긴 경우에는 이 오차는 더 느낄 수 없을 정도가 되기 때문에 이렇게 근사치로 적용한다. 즉, 위 예시는 2초 후의 오차인데 0.02초 간격일 때의 오차나 10초, 100초 후의 오차는 더 작아진다는 말이다.
그럼 여기까지 하고 이제 다음에는 꼭 실제 JavaScript 예제를 만들어보자. 근데 이거 시간이 많이 걸리네…