三つの積分器のエネルギーを比べる

한국어 원문으로 표시합니다.

목표

같은 물리 법칙을 세 가지 방법으로 이산화해 굴려 보고, 그 차이가 그림으로 드러나게 만듭니다. 이 실습이 끝나면 물리 엔진이 왜 특정 적분기를 고르는지 근거를 가지고 말할 수 있게 됩니다.

왜 중요한가

물리 엔진의 겉모습은 충돌과 마찰이지만, 그 아래에서 매 걸음 돌고 있는 것은 적분기입니다. 그리고 적분 방법은 오차만 내는 것이 아니라 계의 성질을 바꿉니다. 명시적 오일러로 용수철을 굴리면 에너지가 매 걸음 정확히 (1 + dt²) 배씩 늘어나 몇 초 만에 발산합니다. 힘 계산이 틀린 것이 아니라 순서가 틀린 것입니다.

두 줄의 순서만 바꾼 반암시적 오일러는 같은 비용으로 안정적입니다. 게임 물리 엔진이 이것을 쓰는 이유이고, 이 사실은 숫자를 직접 뽑아 그래프로 보면 잊히지 않습니다. 그래서 이 실습은 결과를 CSV 로 남기고 한 장의 PNG 에 겹쳐 그립니다.

단계

/root/integrator 에 도구 상자를 놓습니다.
/root/integrator/out/euler.csv — 명시적 오일러.
/root/integrator/out/semi.csv — 반암시적 오일러.
/root/integrator/out/verlet.csv — 속도 베를레.
/root/integrator/out/energy.png — 세 에너지 곡선을 한 장에.
/root/integrator/out/06-dt.txt — dt 를 바꾼 증가율.
/root/integrator/out/damped.csv 와 damped.png — 감쇠를 힘으로.

참고

그림은 nohup python3 -m http.server 8080 -d /root/integrator/out & 로 띄우고 웹 미리보기 에서 http://localhost:8080/energy.png 를 엽니다.
CSV 는 머리글 한 줄 + 401줄입니다. 0걸음째(초기 상태)를 빠뜨리지 마십시오.
흔한 실수 하나: 명시적 오일러를 두 줄로 나눠 쓰다가 새 x 를 속도 계산에 써 버리는 것입니다. 그 순간 그것은 반암시적 오일러가 되고, 두 CSV 가 똑같아집니다.
흔한 실수 둘: 베를레에서 반복문 전에 초기 가속도를 구해 두지 않는 것입니다. 첫 걸음이 틀립니다.

그리기 도구 상자를 놓는다

/root/integrator/gfxlib.py 를 예시 그대로 저장하고, /root/integrator/check.py 로 시험 무늬를 그려 /root/integrator/out/00-check.png 를 만드십시오. 무늬는 64x64 검은 배경에 (0,0)에서 (63,63)까지 흰색(255,255,255) 대각선을 긋고, 그 뒤에 (0,32)에서 (63,32)까지 빨강(255,0,0) 가로선을 겹쳐 그은 것입니다.

이 실습부터는 PNG 인코더를 다시 만들지 않습니다. 첫 실습에서 손으로 만든 것과 같은 코드를 도구로 내려 줍니다 — 여기서 배울 것은 파일 형식이 아니기 때문입니다.

실습 파드에는 볼륨이 없어서 앞 실습에서 만든 파일이 남아 있지 않습니다. 그래서 실습마다 도구 상자를 다시 놓는 것으로 시작합니다.

Canvas(w, h, bg) 를 만들고 line(x0, y0, x1, y1, rgb) 로 두 선을 그은 뒤 write_png(path) 로 저장하십시오. 가로선을 나중에 그어야 교차점 (32,32) 가 빨강이 됩니다.

이 실습에서는 에너지 곡선을 그리는 데 씁니다. 숫자만 보아서는 발산이 잘 안 보입니다.

명시적 오일러

/root/integrator/euler.py 로 용수철(질량 1, 강성 1, 가속도 a = -x)을 x=1.0, v=0.0 에서 시작해 dt=0.05 로 400걸음 굴리고, /root/integrator/out/euler.csv 에 step,x,v,energy 머리글과 401줄을 적으십시오. 에너지는 0.5*v*v + 0.5*x*x 입니다.

명시적 오일러는 한 걸음의 시작 시점 값만 씁니다. 위치와 속도를 동시에 갱신한다고 생각하면 됩니다.

x, v = x + dt * v, v - dt * x    # 오른쪽은 둘 다 옛 값

한 줄로 쓰면 실수하지 않습니다. 두 줄로 나눠 쓰면 새 x 를 v 계산에 써 버리기 쉬운데, 그러면 그것은 반암시적 오일러가 됩니다.

0걸음째(초기 상태)도 한 줄로 적어야 401줄이 됩니다. 마지막 에너지가 처음의 2.7배쯤 되어야 정상입니다 — 그게 이 단계에서 보려는 것입니다.

반암시적 오일러

/root/integrator/semi.py 로 같은 조건을 속도를 먼저 갱신하고 그 새 속도로 위치를 갱신하는 방식으로 굴려 /root/integrator/out/semi.csv 를 만드십시오. 형식은 앞 단계와 같습니다.

순서만 바뀝니다.

v = v - dt * x        # 속도 먼저
x = x + dt * v        # 방금 구한 새 속도를 쓴다

첫 걸음의 값으로 두 방법을 구분할 수 있습니다. 명시적 오일러는 x 가 그대로 1.0 이지만, 반암시적은 1 - dt*dt = 0.9975 가 됩니다. 이 한 자리가 두 방법의 차이 전부입니다.

에너지가 작은 폭으로 오르내리기만 하고 늘어나지 않아야 정상입니다. 계산 비용은 명시적 오일러와 똑같은데 결과는 완전히 다릅니다.

속도 베를레

/root/integrator/verlet.py 로 같은 조건을 속도 베를레로 굴려 /root/integrator/out/verlet.csv 를 만드십시오. 한 걸음은 x += v*dt + 0.5*a*dt*dt, a_new = -x, v += 0.5*(a + a_new)*dt 이고 a 를 a_new 로 바꿔 다음 걸음으로 갑니다.

가속도를 걸음마다 두 번 씁니다. 걸음 시작 시점의 a 와 위치를 옮긴 뒤의 a_new 를 평균 내서 속도를 갱신하는 것이 핵심입니다.

반복문에 들어가기 전에 a = -x 로 초기 가속도를 구해 두어야 합니다. 이걸 빠뜨리면 첫 걸음이 틀립니다.

첫 걸음의 x 는 1 - 0.5*0.0025 = 0.99875 가 나옵니다. 세 방법이 첫 걸음부터 서로 다른 값을 내므로 여기서 확인할 수 있습니다.

에너지가 거의 그대로 유지되어야 정상입니다 — 400걸음 뒤에도 처음의 1000분의 1 안쪽으로 붙어 있습니다.

세 곡선을 한 장에

/root/integrator/plot.py 로 세 CSV 의 에너지를 /root/integrator/out/energy.png(256x256, 검은 배경)에 그리십시오. 걸음 i 의 점은 px = int(i*255/400), py = 255 - int(E/2.0*255) 이고(0~255 밖은 자릅니다), 색은 오일러 (255,80,80), 반암시적 (80,255,80), 베를레 (80,120,255) 입니다. 그리는 순서도 이 순서입니다.

y 축은 에너지 0 을 화면 맨 아래(255), 에너지 2.0 을 맨 위(0)에 놓는 것입니다. 처음 에너지 0.5 는 255 - int(0.5/2*255) = 192 근처에 찍힙니다.

선을 이을 필요 없이 점만 찍어도 됩니다. 401개의 점이 256개 열에 들어가므로 자연스럽게 이어져 보입니다.

CSV 는 open(path).read().splitlines() 로 읽고 첫 줄(머리글)을 건너뛴 뒤 쉼표로 자르면 됩니다.

빨간 곡선만 위로 치솟고 초록과 파랑은 아래쪽에 붙어 있어야 합니다. 그 그림이 이 실습의 결론입니다.

dt 를 바꾸면 얼마나 나빠지나

명시적 오일러로 같은 시뮬레이션 시간 T=20 을 dt=0.01(2000걸음), dt=0.05(400걸음), dt=0.2(100걸음)로 굴려 마지막 에너지가 처음의 몇 배인지 /root/integrator/out/06-dt.txt 에 dt001=, dt005=, dt02= 세 줄로 적으십시오(소수점 아래 여섯 자리).

걸음 수는 T / dt 입니다. 세 경우 모두 같은 20초를 흉내 내지만 결과가 전혀 다릅니다.

이 계에 대해서는 닫힌 식이 있습니다. 한 걸음마다 에너지가 정확히 (1 + dt²) 배가 되므로 n걸음 뒤에는 (1 + dt²)^n 배입니다. 직접 굴린 값이 이 식과 맞는지 확인해 보십시오 — 맞으면 구현이 옳다는 강한 증거입니다.

dt 를 4배로 키웠을 뿐인데 증가율이 수십 배가 됩니다. 그래서 물리 엔진은 dt 를 키우는 대신 걸음을 여러 번 도는 쪽을 고릅니다.

감쇠는 힘으로 넣는다

반암시적 오일러에 감쇠를 힘으로 더해(a = -x - 0.3*v, 이때 v 는 갱신 전 값) dt=0.05 로 400걸음 굴리고 /root/integrator/out/damped.csv(형식은 앞과 같음)와 /root/integrator/out/damped.png(에너지 곡선, 5단계와 같은 좌표 변환, 색은 (255,200,80))를 만드십시오.

한 걸음은 이렇습니다.

a = -x - 0.3 * v      # v 는 아직 갱신하기 전 값
v = v + dt * a
x = x + dt * v

에너지가 한 번도 늘지 않고 줄어들어야 합니다. 400걸음 뒤에는 처음의 1퍼센트 아래로 내려갑니다.

속도에 매 걸음 0.99 를 곱하는 흔한 방식과 비교해 보십시오. 그 방식은 간단하지만 감쇠량이 dt 에 따라 달라져서, 프레임 속도가 다른 기계에서 물체가 다르게 움직입니다. 힘으로 넣으면 dt 가 달라져도 같은 물리를 흉내 냅니다.