비선형 이산 작용과 에너지 보정

정확도와 구조는 서로 다른 질문에 답한다

수치 궤적이 매끄럽고 에너지가 거의 일정해 보여도 위상 흐름을 왜곡할 수 있다. 반대로 symplectic 형식을 보존하는 방법에서 실제 물리적 에너지가 조금씩 진동할 수도 있다. 예쁜 궤도 하나나 유리한 진단값 하나로 방법의 적절성을 모두 판단하지 않는다.

이 단원은 선형 진동자의 상수 행렬을 다시 비교하는 대신 실제 비선형 사차 퍼텐셜을 적분한다. velocity Verlet, 고전 RK4, Verlet 뒤에 에너지를 보정하는 방법을 비교한다. 마지막 방법은 핵심 차이를 드러낸다. 에너지 값 하나를 고쳐도 사상의 미분은 달라지며, 실수인 다음 운동량을 만들지 못할 수도 있다.

그림의 수치 점은 실제 실행한 갱신이다. 중단한 궤적, 미정인 도함수, 계산하지 않은 조회 단계는 그대로 구분한다. 마지막 유효한 점을 복사해 계산이 끝난 것처럼 표시하지 않는다.

Hamiltonian과 물리 단위를 정한다

질량 $m>0$, 강성 $k\ge0$, 사차 계수 $a\ge0$에 대해

$$U(q)=\frac k2q^2+\frac a4q^4,\qquad H(q,p)=\frac{p^2}{2m}+U(q)$$

를 사용한다. 위치는 m, 운동량은 kg m/s, 강성은 N/m, 사차 계수는 N/m³다. 운동방정식은 $\dot q=p/m$, $\dot p=-kq-aq^3$이다. 위치와 운동량의 오차는 단위가 달라 원래 수치를 그대로 더한 상태 오차로 만들 수 없다.

퍼텐셜이 커지는 경우에는 비음수 에너지가 실제 운동을 제한한다. 두 퍼텐셜 계수가 모두 0이면 자유 입자로 돌아간다. 따라서 이 유계 역학 문제의 거친 보폭에서 수치 궤적이 커지는 현상을 연속 해 자체의 발산으로 해석하면 안 된다. 정확한 역학과 선택한 이산 사상을 먼저 구분한다.

Kick–drift–kick 갱신을 실제 실행한다

고정 보폭 $h$에서 velocity Verlet은

$$p_{1/2}=p-\frac h2U'(q),\qquad Q=q+\frac hm p_{1/2},\qquad P=p_{1/2}-\frac h2U'(Q)$$

로 계산한다. Kick은 위치를 고정한 채 운동량을 바꾸고 drift는 운동량을 고정한 채 위치를 옮긴다. 각각 정준 사상이며 그 합성도 이 성질을 유지한다. 같은 순서를 $-h$로 실행하면 정확 산술에서 원래 상태로 돌아온다.

시간 구간 $T$와 단계 수 $N$을 골라 $h=T/N$으로 둔다. 세분화할 때에는 보폭뿐 아니라 단계 수도 함께 바꾸어 최종 물리 시간을 유지한다. 서로 다른 종료 시각을 비교하면 이산화 오차에 다른 물리적 실험의 효과까지 섞인다.

이산 작용에서 Verlet을 유도한다

$h\ne0$에서 두 위치 사이의 작용을

$$L_d(q,Q;h)=\frac{m(Q-q)^2}{2h} -\frac h2\left[U(q)+U(Q)\right]$$

로 근사하자. 두 끝점의 이산 운동량 관계는

$$p=-D_1L_d=\frac mh(Q-q)+\frac h2U'(q),\qquad P=D_2L_d=\frac mh(Q-q)-\frac h2U'(Q)$$

다. 이 식을 풀면 앞의 kick–drift–kick이 나온다. 따라서 이 방법은 에너지 오차가 우연히 작게 측정된 계산만이 아니라 이산 변분 원리에서 얻은 사상이다.

실습은 같은 $L_d$로 두 이산 Legendre 잔차를 계산한다. Verlet에서는 반올림 오차 범위에서 사라지지만 RK4나 에너지 보정 뒤의 운동량은 일반적으로 이 관계를 만족하지 않는다. 이는 선언한 이산 작용이 어느 사상을 생성하는지 검사하는 것이며, 다른 가능한 모든 변분 정식화까지 배제하는 주장이 아니다.

수치 사상 전체를 미분한다

위상 순서를 $(q,p)$로 두고 한 단계 Jacobian을 $M=D\Phi_h$라 하자. 성분 단위는

$$[M]=\begin{pmatrix}1&\mathrm{s/kg}\\\mathrm{kg/s}&1\end{pmatrix}$$

다. 한 자유도에서는 $M^TJM=J$가 $\det M=1$로 줄어든다. 더 높은 차원에서는 행렬식 1만으로 충분하지 않다. 제품 코드는 각 대수 단계의 도함수를 전방 미분으로 함께 전달하며 에너지 보정의 모든 입력 의존성도 포함한다.

전체 궤적의 Jacobian은 단계별 행렬을 순서대로 곱한 값이다. 부동소수점으로 작은 행렬식 결함을 측정했다고 모든 초기 상태의 정리를 얻지는 않는다. 접선 성분이 커지면 행렬식의 두 곱을 빼는 계산 자체도 나쁜 조건을 가질 수 있다. 방법의 해석적 구조 성질과 수치로 계산한 행렬식 곡선을 따로 읽는다.

높은 차수와 symplectic 성질을 비교한다

고전 RK4는 힘·속도의 네 단계를 평가해 가중치 $1,2,2,1$로 결합한다. 충분히 매끄러운 해의 고정 구간에서는 4차 정확도를 갖는다. 같은 보폭이라면 2차 Verlet보다 짧은 시간의 상태 오차가 작은 경우가 많다.

하지만 이 비선형 진동자에서 RK4는 일반적으로 symplectic 방법이 아니다. 상태 오차와 Jacobian 결함을 나란히 비교하면 어느 방법이 모든 지표에서 동시에 이기는 것은 아님을 알 수 있다. 자유 입자는 예외로, 세 방법이 모두 동일한 정확한 전단으로 줄어든다.

방법의 차수는 정칙성·안정성 가정 아래 보폭이 작아질 때의 점근적 주장이다. 임의로 큰 보폭이나 모든 시간 구간, 반올림 오차가 지배하는 비교에서도 항상 명목상의 오차 비율이 나온다는 뜻이 아니다. 실험 결과와 이론의 적용 조건을 함께 기록한다.

임시 위치를 유지하면서 에너지를 보정한다

Verlet이 제안한 $(Q,P_v)$에서 위치를 그대로 두고

$$W=H(q,p)-U(Q),\qquad P=\operatorname{sign}(P_v)\sqrt{2mW}$$

로 운동량만 고르자. 원래 방정식의 물리적 충격량을 유도한 것이 아니라 특정 후처리 규칙을 정의한 것이다. 실수 가지가 존재하면 입력 에너지를 유지하지만 $W<0$이면 그 위치에서 목표 에너지를 가질 실수 운동량이 없다. 보폭을 몰래 줄이거나 음의 제곱근을 0으로 자르지 않고 중단한다.

$W=0$이면 유한한 $P=0$은 존재하지만 일반적인 제곱근 미분 공식은 사용할 수 없다. 상태를 반환하면서 도함수는 미정으로 둔다. 특별히 단순화되는 모든 영근 사상이 수학적으로 미분 불가능하다는 주장은 아니다. 자유 운동과 0 보폭은 알려진 전단·항등식으로 처리한다. $P_v=0$, $W>0$이면 이 부호 규칙이 가지를 결정하지 못한다.

보정에서 입력 에너지의 도함수도 포함한다

여기서 목표는 현재 입력의 $H(q,p)$다. 입력을 미소 변화시킬 때도 고정된 외부 상수 하나가 아니다. 사상을 미분할 때 이 차이를 빠뜨리기 쉽다.

$A=Q_q=1-h^2U''(q)/(2m)$, $B=Q_p=h/m$라 두면 정칙한 비영 운동량에서

$$P_q=\frac mP\left[U'(q)-U'(Q)A\right],\qquad P_p=\frac mP\left[p/m-U'(Q)B\right]$$

다. 따라서

$$\det D\Phi_h=\frac{Ap-hU'(q)}P$$

이며 일반적으로 1이 아니다. 미분하면서 목표 에너지를 고정해 버리면 하나의 고정 에너지 곡선에 대한 사영처럼 다른 사상을 검사하게 된다. 올바른 에너지 값과 올바른 위상 기하를 구별해야 하는 이유다. 같은 이름의 보정이라도 무엇을 입력에 의존시키는지부터 명시한다.

잘라 쓴 수정 Hamiltonian을 조심해서 해석한다

이 velocity Verlet 관습에서 선도 수정 Hamiltonian은

$$H_{\mathrm{trunc}}=H+h^2\left[ \frac{U''(q)p^2}{12m^2}-\frac{(U'(q))^2}{24m}\right]$$

이다. 기호 전개를 하면 한 단계 변화의 $h^4$항까지 상쇄되어 국소 결함이 $O(h^5)$부터 남는다. 작고 고정된 보폭으로 유계 영역을 따라갈 때 실제 에너지보다 이 값의 변화가 훨씬 작을 수 있는 이유를 설명한다.

하지만 잘라 쓴 점근식이지 정확한 보존량이나 무조건적인 안정성 증명이 아니다. 수정 방정식의 논의에는 보폭·정칙성·방문 영역에 대한 가정이 필요하다. 실습은 실제 잘린 식의 변화량을 표시하고 이를 0이라고 선언하지 않는다. 선형 진동자에서만 성립하던 정확한 이차 불변량을 사차 퍼텐셜에 그대로 가져오지도 않는다.

독립적인 물리 기준 궤적을 마련한다

브라우저의 기준 계산은 RK4 단계 이분과 Richardson 보정을 사용한다. 한 번의 전체 보폭과 두 번의 절반 보폭을 비교하며 위치·운동량을 각각의 척도로 검사한다. 절대 허용치는 $10^{-12}$ m와 $10^{-12}$ kg m/s, 상대 허용치는 $10^{-11}$이다. 이는 국소 오차 추정 조건이지 엄밀한 전역 오차 상한이 아니다.

독립 검증은 더 엄격한 허용치의 SciPy DOP853과 별도의 허용치 세분화 해를 비교한다. 마지막 점만이 아니라 기준 격자 전체를 확인한다. 단계 값과 Jacobian은 독립 기호 미분·고정밀 계산으로 검사하므로 기준 적분기의 결과를 symplectic 성질의 증명으로 사용하지 않는다.

무차원 상태 오차가 필요하면 위치와 운동량의 차이를 각각 보고한 물리적 척도로 나눈 뒤 결합한다. 원래 단위의 오차 곡선도 따로 남긴다. 기준 계산의 예산이 끝나는 상황 역시 명시하며, 그런 결과를 정확한 해라고 이름만 바꾸지 않는다.

같은 시간에서 세분화하고 관측 차수를 읽는다

같은 $T$에서 $N,2N,4N$단계를 비교한다. 정규화한 종료 오차가 $e_N,e_{2N}$이면 관측 차수는

$$r=\log_2(e_N/e_{2N})$$

다. 기본 비선형 예제에서는 Verlet이 약 2차, RK4가 약 4차를 보인다. 에너지 보정은 전환점이나 실수 가지 실패, 부호 규칙 변화 근처에서 균일한 관측 차수를 보장하지 않는다.

미완료 실행이나 선언한 비교 바닥값 이하의 오차에는 차수를 보고하지 않는다. 무차원 바닥값 $10^{-10}$은 기준 해·반올림 척도에 가까운 불안정한 비율을 피하기 위한 실험 규칙이다. 더 작은 오차가 수학적으로 0이라는 정리가 아니다. 실패를 제외한 결과만으로 방법 전체의 차수를 과장하지 않는다.

보존량뿐 아니라 역방향 실행도 검사한다

선택한 방법으로 $N$단계 진행한 뒤 $-h$로 $N$단계를 돌아온다. 처음과 마지막의 위치·운동량 차이를 따로 보고한다. Verlet은 대칭적인 합성이므로 반올림 오차 범위에서 돌아온다.

에너지 보정 규칙은 양방향에서 에너지를 유지해도 왕복 오차가 남을 수 있다. 순방향 계산을 마치고도 역방향 중간에서 실수 보정 운동량이 없어질 수도 있다. 그 경우 역방향 오차는 실패 이유와 함께 미정이다. 미정을 0으로 바꾸면 계산 실패를 성공으로 오인한다.

RK4 역시 일반적으로 왕복 결함이 있다. 가역성·symplectic 성질·물리 에너지 거동·수렴 차수는 서로 관련되지만 같은 검사가 아니다. 어느 결과를 실제 측정했는지, 어느 성질을 해석적으로 증명했는지 구분해 기록한다.

한 단계를 계산하고 제어의 의미를 읽는다

$m=k=a=1$, $q_0=1$ m, $p_0=0$, $T=2$ s, $N=16$을 선택한다. $h=1/8$ s이며 첫 반 kick은 $-1/8$, 위치는 $Q=63/64$ m, Verlet의 마지막 운동량은 약 $-0.246139288$ kg m/s다.

초기 에너지는 $3/4$ J, 한 단계 에너지는 약 $0.749526756$ J다. 에너지 보정은 운동량을 약 $-0.248054504$ kg m/s로 바꾸어 에너지를 복원하지만 한 단계 행렬식은 약 $1.007843$이다. 에너지에는 올바른 값이어도 주변 위상 상태의 기하에는 다른 결과를 만든다.

단계 수는 $N=4\,2^{\mathrm{level}}$이고 시간 코드는 2로, 초기 위치·운동량 코드도 각 SI 단위에서 2로 나눈다. 유한한 갱신이라도 위치 $10^6$ m, 운동량 $10^9$ kg m/s의 표시 한계를 넘으면 중단하며 이를 실제 특이성이라고 부르지 않는다. 조회는 정수 단계다. $N$을 넘거나 중단 이후의 단계는 미정이다. $T=0$에서는 모든 갱신이 항등식이지만 $h$로 나누는 $L_d$ 공식은 여전히 정의되지 않는다.

연구에서 어떤 성질을 따로 검사할지 정한다

학습한 보정이 에너지 손실을 줄이면서 위상 부피·가역성·상태 정확도는 바꿀 수 있다. 한 개의 스칼라 목적을 개선하고 나머지가 따라온다고 가정하지 말고 각 성질을 기록한다.

미분 가능한 시뮬레이터에서는 구현된 사상의 미분을 독립 기호식이나 고정밀 도함수와 비교한다. 보정 목표의 입력 의존성을 포함한다. 보존 주장을 할 때에는 정확한 모형, 방법, 보폭 규칙, 정밀도, 실제 완료 시간과 실패 정의역을 명시한다.

이 짧은 실험이 임의 신경망 적분기의 전역 오차 상한이나 장시간 안정성을 증명하지는 않는다. 비선형 모형에서도 주장을 이해하기 쉽게 정식화하고 반증할 수 있도록 검산한 예제다. 후속 논문 실험은 해당 문제의 추가 가정과 더 긴 검증을 별도로 제공해야 한다.

해설과 함께 여섯 가지를 확인한다

문제 1

기본 예제의 첫 Verlet 위치를 구하라.

해설: $U'(1)=2$ N이므로 반 kick은 $-1/8$ kg m/s다. $h/m=1/8$을 곱한 위치 이동은 $-1/64$ m이고 $Q=63/64$ m다. 이산 작용의 양 끝 운동량 관계에서도 같은 값을 얻는다.

문제 2

에너지를 복원한 보정이 symplectic이라고 바로 말할 수 없는 이유는?

해설: 사상의 행렬식은 임시 위치와 입력 에너지의 도함수까지 포함한다. 예제에서 약 $1.007843$으로 1이 아니다. 궤적이 에너지 수준에 남아도 주변 위상 상태들의 넓이는 왜곡될 수 있다.

문제 3

음의 $W$는 무엇을 뜻하는가?

해설: 고른 임시 위치의 퍼텐셜이 전체 가용 에너지보다 크다. 위치를 고정한 채 실수 운동량으로 고칠 수 없다. 연속 초기값 문제의 실패가 아니라 해당 갱신 규칙의 실패다. 모형이나 보폭을 몰래 바꾸지 말고 중단을 보고한다.

문제 4

세분화에서 보폭과 단계 수를 함께 바꾸는 이유는?

해설: $T=Nh$를 유지해야 이산화만 비교한다. 단계 수를 고정하고 보폭만 절반으로 줄이면 더 이른 시각에 끝나므로 종료 상태 차이에 다른 물리적 시간의 효과가 섞인다.

문제 5

작은 잘린 수정 에너지 변화가 정확한 에너지 보존을 증명하는가?

해설: 아니다. 선도 $h^2$ 보정까지만 포함한다. 기호 국소 결함은 원래 에너지보다 높은 차수에서 시작하지만 모든 차수에서 0은 아니다. 실제 변화량과 점근 논의의 가정을 계속 확인해야 한다.

문제 6

순방향을 마친 뒤 역방향 보정이 실패하면 어떤 오차를 보고하는가?

해설: 완료된 역방향 위치·운동량이 없으므로 둘 다 실패 이유와 함께 미정이다. 0이라고 하면 왕복 성공을 허위 인증한다. 역방향도 완료되는 다른 사례에서는 실제 비영 성분 오차를 보고할 수 있다.

읽을거리와 독립 검증

분할·이산 변분 원리는 Hairer의 Symplectic Integrators 1·7절, 수정 방정식과 그 가정은 Hairer·Lubich·Wanner의 Störmer–Verlet 논문을 참고한다. 독립 DOP853은 SciPy solve_ivp 문서를 따른다. 특정 에너지 보정과 그 도함수, 표시 경계, 수치 비교는 여기서 직접 유도하고 검산한다.