곡면 위의 운동 — 구속과 라그랑주 승수

구슬은 고리 위에서만 움직인다 — 구속조건이 자유도를 어떻게 줄이고, 라그랑주 승수가 어떻게 구속력의 정체를 알려주는가.

들어가며

1장의 뉴턴 역학은 입자가 3차원 공간을 자유롭게 떠다닌다고 가정했지만, 실제 문제 대부분은 그 가정과 어긋난다. 구슬은 철사 고리 위에만 머물고, 진자 추는 일정한 길이의 줄 끝에 매달려 있으며, 동전은 미끄러지지 않고 굴러간다. 산업 로봇의 손끝은 6 자유도 안에서만 움직이고, 자동차의 바퀴는 차체와 평행한 방향으로만 굴러간다. 이런 제한된 운동을 한꺼번에 부르는 이름이 구속(constraint) 이다.

본론 1 은 구속을 홀로노믹과 비홀로노믹 두 부류로 가른다. 본론 2 는 수직 원형 고리 위 구슬을 끝까지 풀어 자유도가 어떻게 줄어드는지를 본다. 본론 3 은 7 장에서 본격적으로 다룰 라그랑주 승수(Lagrange multiplier — 구속식을 운동방정식에 끌어들이기 위해 도입하는 미지수로, 그 값이 구속력의 진폭에 정비례한다) 가 구속력 그 자체를 측정한다는 사실을 한 발 먼저 본다.

이 책 전체에서 구속은 물리적으로 입자가 갈 수 있는 영역을 제한하는 식을 뜻한다. 줄의 인장력, 고리의 수직항력, 빙판의 수직항력 같은 구속력은 식 안에 나타나지만 외부에서 주어진 알려진 힘이 아니라 운동과 함께 결정되는 미지의 힘이다. 중력 같은 외부 힘이 알려진 양이라면, 구속력은 미지수에서 출발해 운동 안에서 풀려 나온다.

본론 1 — 두 종류의 구속

구속조건의 두 부류를 가르는 기준은 식이 좌표만의 식으로 적분 가능한가에 있다.

홀로노믹 구속(holonomic constraint — 좌표·시간만의 방정식으로 적히는 구속) 은

$$f(q, t) = 0$$

으로 적힌다. $q = (q^1, \dots, q^n)$ 은 입자(들)의 좌표 묶음이고, $f$ 는 그 좌표를 받아 실수 값을 돌려주는 함수다. ”$f = 0$ 을 만족하는 점만 허용된다”는 약속이 좌표 공간에서 허용 영역을 한 차원 좁힌다. 반지름 $R$ 인 원형 고리 위의 구슬은 평면에서

$$f(x, y) = x^2 + y^2 - R^2 = 0$$

을 만족하므로 2차원이 1차원 곡선으로 떨어진다. 길이 $\ell$ 인 실에 매달린 단진자는 $x^2 + y^2 = \ell^2$ 한 줄이 자유도를 2 에서 1 로 줄인다. 단위 반지름 구면 $S^2$ 위의 점은 $x^2 + y^2 + z^2 = 1$ 을 만족하며 3차원이 2차원 곡면으로 줄어든다. 빙판 위 스케이트 선수의 발끝이 빙판 위에만 머문다면 $z = 0$ 이라는 홀로노믹 구속이다.

홀로노믹 구속의 결정적인 미덕은 좌표 가운데 한 개를 해소해서 변수의 개수를 줄일 수 있다는 점이다. $x^2 + y^2 = R^2$ 에서 $x = R\cos\theta$, $y = R\sin\theta$ 라 두면 두 변수 $(x, y)$ 가 한 변수 $\theta$ 로 바뀌고 구속은 자동으로 만족된다. 1장에서 본 일반좌표의 도입이 바로 이 작업이다. 좌표가 곡면을 따라 흐르게끔 잡혀 있으니 구속력은 식 안에 등장조차 하지 않는다. 단진자에서 $(x, y)$ 대신 $\theta$ 하나만 쓰는 풀이가 이 미덕의 가장 작은 예다.

비홀로노믹 구속(non-holonomic constraint — 속도에 걸리지만 좌표만의 식으로 적분되지 않는 구속) 은

$$\sum_i a_i(q)\, \dot q^i = 0$$

으로 적히지만 어떤 매끄러운 함수 $f(q) = 0$ 형태로 다시 쓸 수가 없다. 미끄러지지 않는 동전의 구름 조건이 대표적이다. 동전이 평면 위에 서 있고 그 면이 평면과 수직이라 하자. 동전의 위치를 $(x, y)$ 로, 가리키는 방향각을 $\phi$ 로 적으면 미끄러지지 않는다는 조건은

$$\dot x \cos\phi + \dot y \sin\phi = 0$$

— 순간의 속도 방향을 동전이 굴러가는 쪽으로만 제한한다. 그러나 동전이 어디에 있을 수 있는가 자체는 제한하지 않는다. 시간을 충분히 들이면 평면 위 어디에든 동전을 가져다 놓을 수 있다 — 직진과 회전을 적당히 섞으면 된다. 자전거·스케이트 날·자동차 바퀴가 모두 이런 부류다. 좌표 공간 자체는 줄어들지 않는데 순간 속도의 방향만 제약된다.

만약 어떤 함수 $f(x, y, \phi)$ 가 있어서

$$\frac{\partial f}{\partial x} = \cos\phi, \quad \frac{\partial f}{\partial y} = \sin\phi, \quad \frac{\partial f}{\partial \phi} = 0$$

을 만족한다면, 동전의 구름 조건은 $df/dt = 0$ 이 되어 $f(x, y, \phi) = \text{const}$ — 홀로노믹 구속 한 줄로 다시 적힌다. 그러나 둘째 식의 $\partial f / \partial y = \sin\phi$ 가 $\phi$ 에 의존하는데 셋째 식은 $\partial f / \partial \phi = 0$ 이어야 한다 — 모순이다. 그러니 그런 $f$ 는 존재하지 않고, 구름 조건은 적분되지 않는다.

반례로 평면 위 입자가 $\dot x - \dot y = 0$ 을 따른다고 하자. 좌변을 시간에 대해 적분하면 $x - y = \text{const}$ — 좌표만의 식이다. 식의 모양은 속도에 걸린 조건이지만 본질은 홀로노믹이다. 한 단계 더 가서 $y\dot x - x\dot y = 0$ 같은 식은 $-x^2 \frac{d}{dt}(y/x) = 0$ 으로 다시 적을 수 있으므로 $y/x = \text{const}$ — 원점을 지나는 직선 위에 머문다는 홀로노믹 조건이다. 그러나 이런 적분 인자가 모든 속도 조건에 대해 존재하지는 않는다 — 그래서 동전의 구름 조건은 진짜 비홀로노믹이다. 판정 기준 — 프로베니우스 정리(Frobenius theorem — 속도에 걸린 조건들의 모임이 좌표 식으로 적분 가능한지 판정하는 미분기하 정리) — 는 7 장에서 다룬다.

비홀로노믹 구속이 있는 계는 자유도의 정의가 미묘하다. 좌표 공간 차원으로 본 자유도(위치 자유도)는 줄어들지 않지만 속도가 갈 수 있는 방향은 제한된다 — 자전거가 어디든 갈 수 있지만 옆으로 순간 이동은 못 하듯이. 동전의 위치 자유도는 3 ($x, y, \phi$), 순간 속도 자유도는 2.

구속은 시간 의존성에 따라서도 두 부류로 나뉜다. 경경적(scleronomic — 시간이 식에 명시적으로 등장하지 않는) 구속은 $f(q) = 0$, 유경적(rheonomic — 시간이 명시적으로 들어가는) 구속은 $f(q, t) = 0$. 회전하는 판 위의 한 홈에 놓인 구슬은 그 홈의 위치가 $t$ 에 따라 바뀌므로 유경적이다. 둘 다 홀로노믹이지만 유경적인 쪽은 에너지가 보존되지 않을 수 있다 — 구속이 시간에 따라 일을 해 줄 수 있기 때문이다.

구속에는 양방향과 단방향의 구분도 있다. 줄이 매단 진자에서 줄이 늘어나지 않는다는 조건은 양방향이다 — 추가 줄 안쪽으로도 바깥쪽으로도 갈 수 없다. 책상 위의 공은 책상에서 떨어질 수 는 있지만 책상을 뚫고 내려갈 수 는 없으므로 단방향이다. 단방향 구속은 부등식 $f(q) \geq 0$ 으로 적히고, 일반적으로 라그랑주 승수가 한쪽 부호만 허용되는 식으로 다루어진다. 이 책에서는 양방향 — 등식 $f(q) = 0$ — 인 경우만 다루지만, 단방향 사례는 충돌·접촉 역학에서 중요하다.

이 장은 홀로노믹·경경적·양방향 쪽만 깊이 다룬다. 본론 2 의 예제 — 수직 원형 고리 위의 구슬 — 도 정확히 이 부류 안에 있다.

본론 2 — 고리 위의 구슬, 처음부터 끝까지

질량 $m$ 인 구슬이 반지름 $R$ 인 수직 원형 고리 위를 미끄러진다. 중력가속도는 $g$. 데카르트 좌표로 적으면 위치는 $(x, y, z)$ 세 변수이지만 고리가 평면 위에 있다면 $z = 0$ 으로 자유도 하나가 떨어지고, 평면 안에서도 $x^2 + y^2 = R^2$ 가 또 하나를 떨어뜨린다. 3차원 공간에서 출발했지만 구속 두 개로 결국 1 자유도짜리 계가 된다.

자유도가 1 이므로 좌표 하나로 위치를 적는다. 가장 자연스러운 선택은 고리의 가장 아래쪽에서 잰 각도 $\theta$ 다. 위치는

$$x = R\sin\theta, \qquad y = R(1 - \cos\theta).$$

$\theta = 0$ 일 때 $(x, y) = (0, 0)$ — 가장 아래 — 이고 $\theta = \pi$ 일 때 $(x, y) = (0, 2R)$ — 가장 위. $\theta$ 는 $[0, 2\pi)$ 안에서 한 바퀴를 돈다 — 1 장에서 본 $S^1$ 의 단순한 사례다.

속도는

$$\dot x = R\cos\theta\,\dot\theta, \qquad \dot y = R\sin\theta\,\dot\theta.$$

두 성분의 제곱을 더하면

$$\dot x^2 + \dot y^2 = R^2 (\cos^2\theta + \sin^2\theta)\,\dot\theta^2 = R^2 \dot\theta^2,$$

피타고라스 항등식 덕분에 운동에너지가

$$T = \tfrac{1}{2} m (\dot x^2 + \dot y^2) = \tfrac{1}{2} m R^2 \dot\theta^2.$$

자유 입자의 $\tfrac{1}{2} m v^2$ 와 같은 모양 — 다만 질량 $m$ 이 관성 모멘트(moment of inertia — 회전 운동에서 질량의 역할을 하는 양으로, 질량에 회전축까지 거리의 제곱을 곱해 모은 것) $m R^2$ 로, 속도가 각속도 $\dot\theta$ 로 바뀐다.

위치에너지는 높이 $y = R(1 - \cos\theta)$ 에 $mg$ 를 곱해

$$V = mgR(1 - \cos\theta).$$

$\theta = 0$ 에서 $V = 0$, $\theta = \pi$ 에서 $V = 2mgR$ — 고리의 위·아래 차이만큼의 위치에너지 차다.

7 장의 라그랑지언 $L = T - V$ 를 결과만 빌려 쓴다. 상수항 $mgR$ 은 운동방정식에 영향을 주지 않으므로 떼어내면

$$L = \tfrac{1}{2} m R^2 \dot\theta^2 + m g R \cos\theta.$$

(부호가 바뀐 이유는 $-(1-\cos\theta) = \cos\theta - 1$ 에서 상수항 $-1$ 을 떼어 낸 결과다.) 오일러–라그랑주 방정식(Euler–Lagrange equation — 작용을 변분으로 멈추게 하는 운동방정식, $L$ 의 좌표·속도에 대한 편미분으로 적힌다)

$$\frac{d}{dt}\!\left( \frac{\partial L}{\partial \dot\theta} \right) - \frac{\partial L}{\partial \theta} = 0$$

에 넣는다. 부분 도함수는

$$\frac{\partial L}{\partial \dot\theta} = mR^2 \dot\theta, \qquad \frac{\partial L}{\partial \theta} = -mgR\sin\theta,$$

식은 $m R^2 \ddot\theta - (-mgR\sin\theta) = 0$, 즉 $m R^2 \ddot\theta + m g R \sin\theta = 0$. $mR$ 로 나누면

$$\ddot\theta + \frac{g}{R} \sin\theta = 0.$$

진자 방정식이다. 원래 $(x, y, z)$ 세 좌표에 구속 두 개가 걸린 3-DOF 문제였지만, 일반좌표 $\theta$ 하나로 단번에 1-DOF 문제가 됐다. 구속이 좌표 안으로 대수적으로 흡수되어 사라진 셈이다.

진폭이 작으면 $\sin\theta \approx \theta$ — 테일러 전개의 1차 항 — 로 근사할 수 있고 운동방정식은

$$\ddot\theta + \frac{g}{R}\theta = 0$$

— 각진동수 $\omega = \sqrt{g/R}$ 인 조화진동이다. 일반해는 $\theta(t) = A\cos(\omega t + \varphi)$ 이고 주기는 $T_0 = 2\pi/\omega = 2\pi\sqrt{R/g}$. $R = 0.2\,\mathrm{m}$ 이면 $\omega = \sqrt{9.81/0.2} \approx 7.00\,\mathrm{rad/s}$, 주기 $T_0 \approx 0.898\,\mathrm{s}$. $R = 1.0\,\mathrm{m}$ 이면 $\omega \approx 3.13\,\mathrm{rad/s}$, 주기 $T_0 \approx 2.01\,\mathrm{s}$ — 길이 1 미터짜리 진자가 한 번 흔들리는 데 약 2초.

진폭이 크면 1차 근사가 깨지면서 주기가 진폭에 의존한다. 진폭이 클수록 마루 근처에서 $\sin\theta < \theta$ 가 강하게 작동해 복원이 약해지고, 한 주기에 걸리는 시간이 길어진다. 주기를 완전 타원적분(complete elliptic integral — 진자의 정확한 주기 공식에 등장하는 적분으로, 진폭이 클수록 작은 진폭 주기에서 멀어진다) 으로 적을 수 있는데, 진폭 $\theta_0$ 에서의 주기는 대략 $T \approx T_0 (1 + \theta_0^2/16 + \cdots)$ 로 늘어난다. $\theta_0 = \pi/3$ 정도면 주기가 약 5% 늘어나고 $\theta_0 = \pi/2$ 이면 약 18% 늘어난다. 이 비조화성(anharmonicity — 진폭에 따라 주기와 모양이 달라지는, 조화진동에서 벗어난 성질) 효과를 본 절 마지막의 파이썬 코드에서 직접 본다.

같은 운동을 데카르트로 다시 적어 보면 왜 일반좌표 길이 그렇게 짧은지가 더 분명하다. 데카르트로는 변수가 $(x, y, z)$ 셋, 구속이 $z = 0$ 과 $x^2 + y^2 = R^2$ 둘이다. 운동방정식은 $m\ddot{\mathbf{r}} = m\mathbf{g} + \mathbf{N}$ — 중력에 더해 고리의 수직항력 $\mathbf{N}$ 이 들어간다. $\mathbf{N}$ 은 고리에 수직한 방향으로만 작용한다는 조건과 구속식을 두 번 미분해서 얻은 식까지 합쳐, 미지수 5 개($x, y, z, N_x, N_y$ — $N_z$ 는 $z = 0$ 에서 자동 풀린다)에 식 5 개를 푸는 셈이 된다. 자유도 1 짜리 계가 풀리는데 미지수는 5 개 — 1 장 본론 1 에서 본 어색함의 정량적 모습이다.

$\theta$ 한 변수로 적은 풀이는 식이 한 줄, 미지수가 하나, 구속력 $\mathbf{N}$ 은 식에 등장조차 하지 않는다 — 좌표가 구속을 자동으로 만족시키니 구속력은 대수적으로 사라진다.

그러면 $\mathbf{N}$ 은 어떻게 구하나? 고리의 수직항력은 공학적으로 중요한 양이다. 고리가 그 힘을 견딜 수 있는지, 어느 지점에서 가장 큰 힘이 걸리는지가 설계 문제로 직접 들어온다. 진자 줄에 걸리는 인장력이 줄의 항복 강도를 넘으면 줄이 끊어진다. 롤러코스터의 트랙 설계에서도 마루 꼭대기에서 차량이 트랙에서 떨어지지 않게 최소 수직항력을 확보해야 한다는 조건이 직접 들어간다. 일반좌표 길은 운동의 모양은 깔끔하게 주지만 구속력 그 자체에 대해서는 아무 말을 해 주지 않는다.

본론 3 — 라그랑주 승수가 가리키는 것

구속력을 알고 싶을 때는 좌표를 줄이지 않고 원래 좌표 $(x, y)$ 를 그대로 두면서 구속식

$$f(x, y) = x^2 + y^2 - R^2 = 0$$

을 별도로 부과한다. 운동방정식에 한 항을 더해

$$m \ddot{\mathbf{r}} = -\nabla V + \lambda\, \nabla f.$$

기호 $\nabla$ (나블라)는 그래디언트(gradient — 스칼라 함수에 작용해 그 함수가 가장 빨리 증가하는 방향을 가리키는 벡터를 돌려주는 연산자) 로, 평면 안에서 $\nabla V = (\partial V/\partial x,\ \partial V/\partial y)$. $\lambda$ (lambda) 는 라그랑주 승수 라는 새 미지수다.

운동방정식의 오른쪽에 외부 힘 $-\nabla V$ 와 구속력 $\lambda \nabla f$ 가 나란히 들어가 있다. 외부 힘은 알려진 양, 구속력은 미지수. $(x, y)$ 두 좌표에 $\lambda$ 한 미지수, 식은 운동방정식 두 줄에 구속식 한 줄 — 균형이 맞다.

평면 안에서 $f = x^2 + y^2 - R^2$ 의 그래디언트는

$$\nabla f = (2x,\ 2y) = 2(x, y).$$

구슬이 고리 위에 있다면 $(x, y)$ 는 원의 한 점 — 원 중심에서 그 점을 향하는 반지름 방향 벡터다. 크기 $2R$ 의 그 벡터가 곡면(여기서는 원)에 수직이다 — 원에 접하는 방향을 따라 움직이면 $f$ 가 변하지 않으니 $\nabla f$ 는 접하는 방향에 직교한다. 그래서 $\nabla f$ 의 방향은 곡면의 바깥 법선 방향(outward normal — $f > 0$ 쪽을 가리키는 곡면 수직 벡터) 이다.

추가된 항 $\lambda \nabla f$ 는 곡면에 수직한 방향의 힘이다. 그 방향이 바로 수직항력의 방향 — 고리가 구슬을 안쪽으로 밀거나 밖으로 잡아당기는 방향. 즉 $\lambda \nabla f$ 가 구속력의 정체이고 그 진폭(부호 포함)이 $\lambda$ 다.

$\lambda$ 의 값을 구하려면 구속식 $f = 0$ 을 시간에 대해 두 번 미분한다.

$$\frac{d}{dt}f = 2x\dot x + 2y\dot y = 0,$$ $$\frac{d^2}{dt^2}f = 2\dot x^2 + 2x\ddot x + 2\dot y^2 + 2y\ddot y = 0.$$

여기에 운동방정식 $m\ddot x = -\partial V/\partial x + 2\lambda x$, $m\ddot y = -\partial V/\partial y + 2\lambda y$ — 즉 $\ddot x = -(1/m)\partial V/\partial x + (2\lambda/m) x$, $\ddot y$ 도 마찬가지 — 를 넣고 정리하면 $\lambda$ 에 대한 한 줄의 식이 나온다. 이 식을 풀면 $\lambda$ 가 운동의 상태 $(x, y, \dot x, \dot y)$ 의 함수로 결정된다.

원형 고리에서 끝까지 밀어 보면 매끄럽게 떨어진다. $\theta$ 좌표로 옮겨 쓰면 $\dot x^2 + \dot y^2 = R^2 \dot\theta^2$ — 본론 2 그대로다. $x\ddot x + y\ddot y$ 는 위치에너지 항과 라그랑주 승수 항으로 분리되고, 정리하면

$$2R^2\dot\theta^2 + 2(x, y)\cdot \ddot{\mathbf{r}} = 0$$

이라는 한 줄에서 출발해, 운동방정식과 구속을 함께 풀면 $\lambda$ 가 $\dot\theta^2$ 와 고리의 자세에 따른 중력 성분의 결합으로 떨어진다. 중력이 식에서 빠지는 수평으로 놓인 매끈한 원형 고리에서는 구속력의 크기가 정확히 구심력(centripetal force — 원 중심을 향해 회전 운동을 유지시키는 힘으로 크기 $m R \dot\theta^2$)

$$|\lambda \nabla f| = 2|\lambda| R = m R \dot\theta^2$$

이 된다. 즉

$$\lambda = \frac{m \dot\theta^2}{2} \cdot (\pm 1).$$

부호는 곡면의 어느 쪽이 바깥인가에 따라 달라진다. 고리가 구슬을 안쪽으로 잡아당기는 상황(구슬이 고리 안쪽에서 미끄러지는 경우)이면 $\lambda < 0$, 바깥에서 미끄러지는 경우면 $\lambda > 0$ — 부호는 좌표 선택과 바깥의 정의에 따라 달라진다.

라그랑주 승수가 구심력의 진폭에 정비례한다. 구슬이 빠르게 돌수록($\dot\theta$ 가 클수록) 고리가 구슬을 잡아 두는 데 필요한 힘이 커지고 그에 비례해 $\lambda$ 도 커진다. 수직 원형 고리에서는 여기에 중력의 수직 성분 — $-mg\cos\theta$ 같은 항 — 이 더해져서 $\lambda$ 가 위치에 따라 변하는 더 일반적인 식이 된다. 어느 경우든 구속력의 크기가 $\lambda$ 한 숫자로 잡힌다.

흔한 혼동 셋. 첫째, 라그랑주 승수의 부호는 좌표 선택에 의존한다. 같은 구속을 $f = x^2 + y^2 - R^2$ 로 적든 $g = R^2 - x^2 - y^2$ 로 적든 물리적 구속은 똑같지만 $\nabla f$ 와 $\nabla g$ 의 방향이 정반대다. 한쪽에서 $\lambda > 0$ 인 상황이 다른 쪽에서는 $\lambda < 0$ 으로 나온다. 구속력 그 자체는 $\lambda \nabla f$ — 두 부호 변화가 상쇄되어 — 좌표 선택과 무관하지만, $\lambda$ 만 따로 떼어 보면 어느 방향을 바깥으로 잡았는가에 따라 부호가 바뀐다.

둘째, 벡터 $\nabla f$ 는 곡면의 바깥 법선 방향이다. 정의상 $f > 0$ 쪽이 바깥, $f < 0$ 쪽이 안쪽이라는 약속이다. $f = x^2 + y^2 - R^2$ 로 적으면 원의 밖이 $f > 0$, 안이 $f < 0$ 이므로 바깥 법선은 원 중심에서 멀어지는 방향이다. $g = R^2 - x^2 - y^2$ 로 적으면 같은 곡면이지만 바깥의 방향이 반대다.

셋째, 라그랑주 승수는 시간의 함수다. $\lambda = \lambda(t)$ 로, 구슬이 운동하는 동안 매 순간 다른 값을 가진다. 수평 고리에서는 $\lambda \propto \dot\theta^2$ 이고, 에너지 보존에 의해 $\dot\theta^2$ 이 위치에 따라 변하므로 $\lambda$ 도 같이 변한다. 수직 고리에서는 중력이 더해져서 $\lambda(t)$ 의 모양이 한층 복잡해진다. 가장 아래쪽($\theta = 0$)에서 구슬이 가장 빠르게 지나가니 $|\lambda|$ 가 가장 크고, 마루 근처에서는 $|\lambda|$ 가 가장 작다. 롤러코스터를 탔을 때 골짜기에서 더 많이 짓눌리고 마루에서 거의 둥둥 뜨는 그 감각이 이 식 한 줄에 압축되어 있다.

라그랑주 승수는 단순히 제약식을 만족시키기 위해 동원한 수학 도구가 아니라 구속력의 진폭 그 자체를 가리키는 물리량이다. 세 줄 요약: $\lambda$ 는 구속력의 크기, $\nabla f$ 는 곡면의 바깥 법선, $\lambda$ 는 시간의 함수. 정식 유도는 7 장 변분원리와 구속 에서 다룬다.

파이썬으로 확인

# 고리 위의 구슬: theta_ddot = -(g/R) sin(theta) 를 RK4로 적분.
import numpy as np
import matplotlib.pyplot as plt

g, R = 9.81, 0.2
dt, T_end = 1e-3, 4.0
N = int(T_end / dt)

def f(state):
    th, om = state
    return np.array([om, -(g / R) * np.sin(th)])

state = np.array([np.pi / 3, 0.0])   # theta0 = 60도, 정지에서 출발
ts = np.linspace(0, T_end, N + 1)
ths = np.empty(N + 1)
oms = np.empty(N + 1)
ths[0], oms[0] = state

for i in range(N):
    k1 = f(state)
    k2 = f(state + 0.5 * dt * k1)
    k3 = f(state + 0.5 * dt * k2)
    k4 = f(state + dt * k3)
    state = state + (dt / 6) * (k1 + 2*k2 + 2*k3 + k4)
    ths[i + 1], oms[i + 1] = state

fig, ax = plt.subplots(2, 1, sharex=True)
ax[0].plot(ts, ths); ax[0].set_ylabel(r'$\theta$ [rad]')
ax[1].plot(ts, oms); ax[1].set_ylabel(r'$\dot\theta$ [rad/s]')
ax[1].set_xlabel('t [s]')
plt.tight_layout(); plt.show()

상태 변수는 $(\theta, \dot\theta)$ 두 개 — 1 장의 어휘로 위상공간 차원 $2 \times 1 = 2$ 인 1-DOF 계의 상태 한 점. 함수 f(state) 는 그 상태에서의 시간 도함수 $(\dot\theta, \ddot\theta)$ 를 돌려준다 — 본론 2 의 $\ddot\theta = -(g/R)\sin\theta$ 를 그대로 옮긴 결과다. 적분기는 4차 룽게–쿠타(RK4 — Runge–Kutta 4th order — 한 시간 단계 안에서 네 점의 기울기 $k_1,\dots,k_4$ 를 평균해 갱신하는 표준 ODE 적분기). 명시적 오일러보다 한 단계 정확도가 높지만 심플렉틱은 아니다. 4 초 정도의 짧은 시간에서는 에너지 드리프트가 거의 보이지 않지만 적분 구간이 길어지면 누적된다 — 1 장 마지막 절의 명시적/심플렉틱 오일러 비교와 같은 사정이다.

초기조건은 $\theta_0 = \pi/3$ (60도), $\dot\theta_0 = 0$ — 정지 상태에서 60도 위에서 놓는다. $R = 0.2\,\mathrm{m}$, $g = 9.81\,\mathrm{m/s^2}$ 이면 본론 2 의 작은 진폭 주기 $T_0 = 2\pi\sqrt{R/g} \approx 0.898\,\mathrm{s}$. 4 초 동안 적분하면 대략 네다섯 번의 주기를 본다.

궤적은 한 주기를 닫고 다시 같은 진폭으로 돌아온다. 진폭이 작았다면 거의 사인 곡선이지만 $\theta_0 = \pi/3$ 처럼 큰 진폭에서는 마루 근처에서 약간 평평해지는 비조화성이 눈에 띈다. $\sin\theta$ 의 비선형성이 만들어내는 효과다. 정량적으로 큰 진폭의 주기가 작은 진폭의 주기보다 길어지므로 4 초 동안의 진동 횟수가 작은 진폭 가정으로 셈한 횟수보다 약간 적게 나온다. 진폭을 — $\theta_0 = 0.1,\ 0.5,\ 1.0,\ \pi/3,\ \pi/2$ — 바꿔 가며 돌려 보면 비조화성이 어느 진폭에서 눈에 들어오기 시작하는지를 손으로 확인할 수 있다.

이 코드에서 구속력은 식 자체에 없다. 좌표를 $\theta$ 하나로 잡았기 때문에 구속력은 운동방정식 안에 등장하지 않는다. 구속력까지 추적하려면 본론 3 의 라그랑주 승수 $\lambda$ 를 매 시간 단계에서 함께 계산하면 된다. 가장 단순한 형태로는 $\lambda$ 가 $\dot\theta^2$ 와 중력의 수직 성분의 결합으로 떨어지므로 state 와 f(state) 가 갱신될 때마다 한 줄을 더해 $\lambda(t)$ 를 함께 적는다. 정식 다룸은 7 장 이후로 미룬다.

총 에너지는 $E = \tfrac{1}{2} m R^2 \dot\theta^2 + m g R (1 - \cos\theta)$ 이고 시간 적분 동안 수학적으로 일정해야 한다. 매 시간 단계에서 $E(t)$ 를 계산해 그래프로 그리면, RK4 의 수치적 정확도가 4 초 동안 만드는 에너지 드리프트 — 보통 $10^{-8}$ 수준 — 가 눈에 들어온다. 1 장의 명시적 오일러와 비교하면 RK4 가 얼마나 잘 작동하는지 그대로 보이고, 동시에 심플렉틱이 아닌 적분기이므로 시간이 충분히 길어지면 드리프트가 누적된다는 한계도 분명해진다. 이 차이가 11 장에서 다시 심플렉틱 적분기가 등장하는 이유다.

다음 장으로

3장: 텐서와 공변미분에서는 이 장에서 곡면이라 부른 대상을 다양체라는 정식 이름으로 다시 만나고, 그 위에서 벡터를 어떻게 미분해야 좌표 선택에 의존하지 않는지를 본다. 라그랑지언이 “좌표를 바꿔도 같은 방정식을 준다”는 점이 우연이 아니라는 사실도 거기서 분명해진다.