# HHT 수치 제어·검증 설계 추가 근거

- Feature: `hht-implicit-dynamics`
- Owner: `research-agent`
- Date: `2026-09-22`
- Status: `Evidence and verification candidates; numerical thresholds not selected`
- 기준: [프로젝트 지침](../../AGENTS.md), [PRD](../PRD.md), [ADR](../ADR.md#adr-028), [기존 연구](research.md).

[관성 정식화](formulation-inertia.md)와 [독립 probe](numerical-review-probes.md)를 읽고
남은 질문만 조사했다. 기존 SDOF/moving-null 계산은 재실행하지 않았다.
**사실**은 열람한 일차 자료, **도출**은 명시한 수학적 계산, **후보**는 검증 제안이다.
수치 threshold·질량 채택·물리 범위·입출력 계약은 이 문서에서 결정하지 않는다.
현재 질량 선택은 사용자 결정으로 consistent mass가 확정되었다. Shell 속도 보간의
운동에너지를 그대로 적분하며, lumped는 이번 HHT 구현의 지원 요구에서 제외한다.

## 1. Physical acceleration, startup와 DAE 적용 경계

**사실 [N1].** Brüls–Arnold (2008) author postprint pp.3–4, Eq. (5)–(8)은
물리적 `q_ddot`와 보조 acceleration-like variable을 구분하고 보조 초기값을
물리적 초기 가속도와 같게 둔다. 이 문서의 `A`가 논문의 보조 `a`에 대응한다.
`alpha_m=0`, `alpha_f=-alpha`를 대입하면 현재 후보의
`A_(n+1)=(1+alpha)a_(n+1)-alpha a_n`이다. 같은 논문은 variable M일 때 이를
force-level 가중 평형과 동일시할 수 없다고 명시한다. 이 질문에는 기존 출처가 충분하다.

**사실 [N1]: 정리의 가정은 정확히 구분한다.** Printed p.6, “Equivalent model”은
`[[M,-L],[-h_qddot,I-h_y]]` block의 비특이성과 index-3 조건
`Phi_q g*_lambda`의 비특이성을 사용한다. 여기의 `h`는 논문 Eq. (4)의 algebraic
output function이며 시간간격이 아니다. Controller/algebraic coupling을 제거한
특수형에서는 첫 조건이 M의 비특이성으로 줄어든다. 논문 전체가 무조건 “M만
비특이이면 충분”하다고 읽어서는 안 된다. “Initial value condition”의 초기 오차 조건과
“Constrained case”의 정칙성 하에서 얻은 수렴 결과를 singular drilling에 자동 적용할 수 없다.

**도출/수치 제어에의 함의.** 현재 shell의 `psi=N(q)^T g=0`는 논문의 holonomic
`Phi=0`와 동일한 모델임이 아직 증명되지 않았다. 따라서 다음을 분리해서 관찰한다.

| 관찰 대상 | 필요한 검증 후보 | 남은 한계 |
| --- | --- | --- |
| 첫 시점 | `a0`는 우측 물리 평형·초기 대수조건으로 결정한 뒤 `A0=a0`; 첫 step과 이후 오차를 분리 | 작은 Newton residual로 불일관 초기값이 정당화되지 않음 |
| 물리 a와 A | 동일 시각 정확해에 대한 q/v/physical a 차수와 `A-a`를 별도 측정 | A를 물리 출력·각운동량 미분에 대입하면 안 됨 |
| Current null | 매 configuration에서 `MN`, `N^T I`, `psi`, 양의 부분공간 mass 및 `psi_q N` 정칙성 | 이전 시점 null, 고정 UR3, 전체 M의 역행렬로 대체 불가 |
| Hidden constraints | `psi_q v+psi_t`와 `psi_q a+psi_qq[v,v]+2psi_qt v+psi_tt`의 최대·종단 결함과 dt 수렴 | 위치 평형만 만족해도 두 결함이 0인 것은 아님 |
| 반례 coverage | 현재 움직이는 영공간 예제를 full redundant coordinates로 풀어 기존 reduced-angle probe와 대조 | 기존 probe는 위치조건을 정확히 강제했으므로 full-coordinate Newton 구현 검증이 아님 |
| 고정 간격 | 기존 후보에서 dt를 줄여 같은 물리 종단시각 비교 | adaptive dt·startup 재설정은 별도 scope이며 이번 근거로 추가하지 않음 |

**도출.** Chart 내부의 smooth q에서 projector와 null basis는 같은 물리 영공간을 나타내야
한다. Null vector의 부호·정규화 변경 자체와 rank 손실을 구분해야 한다.
`psi_q`에는 `DN`이 들어가므로 basis를 고정한 finite difference는 별도 식을 검증하게 된다.
투영 오차·hidden defect가 큰 결과를 발견하면 정확한 projection/condensation의 채택 여부는
formulation owner에게 돌린다. 여기서는 추가 projection을 자동 제안·채택하지 않는다.

## 2. Double precision, conditioning 및 미분 검증

**사실 [N2–N4].** IEEE double의 정상수 significand는 53 bits다. LAPACK DGECON은
선택한 1/∞ norm의 reciprocal condition number **추정값**을 반환한다. DGESVX는
equilibration 여부를 구별하고, factor pivot이 정확히 0인 경우와 reciprocal condition이
machine precision보다 작은 경우를 서로 다른 INFO로 보고한다. 후자는 해·오차 추정을
계속 계산할 수 있는 warning이다. 이 공개 library 동작은 FESA failure threshold를
자동 승인하지 않으며 vendor routine 채택 요구도 아니다.

**도출.** Binary64에서 1의 다음 표현수 간격 `eps=2^-52`, round-to-nearest 정상 범위의
unit roundoff `u=2^-53`를 구분한다. 문헌의 “machine epsilon” 정의와 library 반환 의미를
확인하지 않고 두 값을 혼용하면 안 된다. 작은 perturbation의 forward error는
대략 `condition × backward error`로 증폭되므로 residual 하나로 해의 정확도를 보장하지 못한다.

| 질문 | 근거에 따른 분리 및 검사 후보 |
| --- | --- |
| Physical mass-null인가? | `n=[0;J_l^-1 d]`, admissible free-space intersection 같은 분석적 운동학을 원본으로 삼는다. `MN`의 수치 결함과 양의 complement의 conditioning을 각각 검사한다. |
| 작은 양의 eigenvalue인가? | 단위·두께 때문에 작은 inertia를 rcond 하나로 drilling null로 바꾸지 않는다. 예상 추가 null/geometry 이상과 positive-space accuracy loss를 구분한다. |
| 어느 rcond인가? | Full singular M의 rcond는 positive inertia 품질을 설명하지 못한다. Positive-space mass, 초기화 algebraic block, effective Newton matrix를 각각 이름 붙이고 norm·scaling을 기록한다. |
| 단위를 바꾸면? | 병진·회전 단위가 섞인 원행렬 condition은 단위 선택에 의존한다. `q=D qhat`, `Rhat=D^T R/Escale`이면 `Khat=D^T K D/Escale`라는 virtual-work 일관 scaling을 대조 후보로 삼는다. 기준 scale의 수치 선택은 NR 소유다. |
| Library equilibration이면 충분한가? | Row/column equilibration은 solver conditioning 도구다. 물리적 nondimensionalization이나 예상 null proof와 동일하지 않다. 원·scaled 행렬 중 무엇을 추정했는지 보존한다. |

**사실 [N5].** GSL의 central differentiation은 truncation과 roundoff를 함께 추정하여
step을 고르며, 5-point와 3-point 차이를 오차 추정에 쓴다. 따라서 가장 작은 step을
항상 정답으로 간주하지 않는 독립 derivative 검사 설계가 공식 수치 library의 접근과 부합한다.

**도출/후보.** Dimensionless direction p에 대해 scaled 잔차의 중심차분과 `Khat p`를
비교하면 일반적인 smooth·비상쇄 구간에서 오차형은 `C_t delta²+C_r u/delta`다.
최적 step의 `u^(1/3)` 크기는 이 가정에서의 점근적 안내일 뿐 공통 고정값이 아니다.
Zero derivative, 큰 residual cancellation, branch 접근에서는 계수와 error floor가 달라진다.
여러 delta를 sweep하여 감소구간·최저 plateau·roundoff 증가구간을 기록하고, 모든 sample이
같은 허용 chart와 geometry 안에 있는지 확인한다. Unit-aware scale을 쓰고 임의
`max(1, norm)`로 작은 문제를 덮지 않는다.

Inertia partial derivative 검사에서는 q/v/physical a 중 하나만 바꾸고 나머지는 고정한다.
Effective tangent 검사에서는 같은 trial q에 연동되는 v, A, physical a를 모두 갱신한다.
두 검사는 서로 다른 미분이므로 섞지 않는다. 작은 dt의 acceleration 재구성에는
위치 차분 오차가 `1/dt²`로 증폭되므로 시간 차수 측정에서 너무 미세한 dt의 plateau를
수학적 차수 저하로 오인하지 않는다. 기준 oracle 정밀도와 nonlinear solve 오차도 분리한다.

## 3. 3차 회전 미분과 운동량·에너지 case 후보

**사실 [N6].** Higham–Relton의 Theorem 3.5/Algorithm 3.6은 고차 Fréchet derivative를
재귀 block matrix function의 우상단 block으로 구한다. Matrix exponential은 smoothness
조건을 만족한다. `P=[phi]_cross` 및 세 방향 행렬에 적용하면 24×24 exponential 하나의
3×3 우상단 block이 3차 혼합 미분을 준다.

**후보.** 이 방식에 `d0`를 곱한 결과를 Rodrigues/회전 Jacobian의 해석 미분과 비교하면
동일 scalar coefficient 도함수를 재사용하는 자기검증을 피할 수 있다. “Exact derivative”는
수학적 항등식을 뜻하며 부동소수점 exp 계산까지 오차가 없다는 뜻은 아니다.
Production 의존성 추가나 24×24 exp의 runtime 사용은 요구하지 않는다.

**자체 도출: 영회전의 유한 다항 oracle.** `H_p=[p]_cross` 등으로 쓰면

\[
D^3d(0)[p,s,r]=\frac16\sum_{\sigma\in S_3}
H_{\sigma(p)}H_{\sigma(s)}H_{\sigma(r)}d_0.
\]

이는 exponential Taylor series의 cubic 항에서 바로 얻는다. 방향 permutation 대칭,
각 방향의 선형성, 동일 방향 3회 `H_p³ d0`를 함께 검사할 수 있다. 비영 phi의
block-exponential oracle와 `D²d`의 방향 중심차분을 추가해 영점 특례만 맞는 오류를 검출한다.
Chart 경계 안쪽, 작은 phi, 서로 비평행한 방향을 포함하되 경계를 넘어 차분하지 않는다.

| ID | 실행 가능한 검증 설계 후보 | 판정에 필요한 한계 |
| --- | --- | --- |
| NR-E01 | 임의 유효 shell 상태에서 spatial rigid velocity를 부여하여 kinetic energy·총 P·총 L을 독립 volume integration과 대조 | Rigid velocity는 순간 운동학 검사다. 회전하는 탄성 shell이 자유운동에서도 강체를 유지한다는 가정은 하지 않음. Consistent/lumped 각각 자신의 measure를 사용 |
| NR-E02 | `I(q,v,a)`를 q/v/a별 미분하고 전체 effective tangent도 별도 directional check | 3차 director 항과 `D(H^T)` 누락을 분리 검출 |
| NR-E03 | 작은 smooth 회전 구간의 자유 shell에 일관된 initial velocity를 부여하거나 짧은 smooth loading 후 unload, P/L/총에너지·work와 dt 수렴 관찰 | 초기조건/시간하중이 CLI 범위 밖이면 kernel/procedure verification 후보로만 유지. 승인 input subset 확대 금지 |
| NR-E04 | 같은 spatial problem에 alpha=0 및 승인된 음수 alpha를 적용하고 저주파 해상 case와 높은 modal content case를 분리 | 일반 비선형 HHT의 총에너지 매 step 감소 또는 정확한 각운동량 보존을 무조건 요구하지 않음 |
| NR-E05 | 문헌 free-flying single-layer plate를 후속 benchmark 후보로 검토 | 아래 N7은 large rotation/energy 사례의 출처다. 현재 fixed chart를 넘어서는 tumble·장시간 운동은 그대로 채택할 수 없음 |

**사실 [N7].** Vu-Quoc–Tan (2003), §5.1.2/Figs.9–12는 free-flying single-layer plate의
큰 전체 운동과 P/L/에너지 및 dissipative/nondissipative integrator 결과를 비교한다.
해당 요소는 rotation DOF가 없는 solid-shell이므로 그 mass·응력 모델과 여기의 director shell은
다르다. Author-uploaded 본문의 관련 절을 확인했으나 그림의 하중 위치·수치 history를
독립 추출하지 않았다. 따라서 현 단계에서는 case family 근거이며 완성된 reference artifact가 아니다.

**추론.** 현재 shell 운동학에서는 `Pdot/Ldot`에 physical a를 사용하고 numerical drill
potential을 포함한 에너지·내력 일을 일관되게 대조해야 한다. N7의 보존형 integrator가
보이는 exact conservation을 HHT의 요구조건으로 옮길 수 없다. [기존 연구 S8](research.md#sources)의
비선형 Newmark 보존 한계 근거는 그대로 충분하다. 추가 shell 연구는 기하·하중·시간 구간이
승인된 chart와 입력 범위 안에 있는지 결정한 뒤에만 구체화할 필요가 있다.

## 출처·열람 한계와 인계

`accessed: 2026-09-22`. 문서 정적 검사만 수행하며 solver/build/test/reference 생성은 수행하지 않았다.

| ID | Stable citation / locator | 열람 범위 |
| --- | --- | --- |
| N1 | Brüls, O.; Arnold, M. (2008), *The generalized-alpha scheme as a linear multistep integrator: Toward a general mechatronic simulator*, J. Computational and Nonlinear Dynamics 3, 041007, [DOI:10.1115/1.2960475](https://doi.org/10.1115/1.2960475), [author PDF](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf) | Printed pp.3–6 Eq.5–8/Equivalent model, Initial value condition/Constrained case. Shell singular-mass 정리 아님 |
| N2 | Oracle, *Numerical Computation Guide*, [§2.2.3 Double Format](https://docs.oracle.com/cd/E37069_01/html/E39019/z4000ac019269.html) | 53-bit significand 설명. IEEE 표준 전체 적합성 시험을 수행한 것은 아님 |
| N3 | LAPACK, [DGECON](https://netlib.org/lapack/explore-html/d4/daf/group__gecon_ga4f9b830e19e12c7f082ddb497a57af18.html), Purpose/NORM/ANORM/RCOND | Normwise reciprocal condition estimate |
| N4 | LAPACK, [DGESVX](https://netlib.org/lapack/explore-html/d5/dbe/group__gesvx_ga82173a93234afc15d70b64233b3e5bc8.html), Description steps 1–6/INFO | Equilibration, singular pivot와 working-precision warning 구별 |
| N5 | GNU Scientific Library 2.8, [Numerical Differentiation](https://www.gnu.org/software/gsl/doc/html/diff.html), `gsl_deriv_central` | Truncation/roundoff 기반 step/error estimate. FESA tolerance 원본 아님 |
| N6 | Higham, N.J.; Relton, S.D. (2014), *Higher Order Fréchet Derivatives of Matrix Functions and the Level-2 Condition Number*, SIAM J. Matrix Anal. Appl. 35, 1019–1037, [DOI:10.1137/130945259](https://doi.org/10.1137/130945259), [author PDF](https://eprints.maths.manchester.ac.uk/2066/1/paper.pdf), Theorem 3.5/Algorithm 3.6 | 2013 author version printed pp.5–6; oracle identity와 계산법 |
| N7 | Vu-Quoc, L.; Tan, X.G. (2003), *Optimal solid shells for non-linear analyses of multilayer composites. II. Dynamics*, CMAME 192, 1017–1059, [DOI:10.1016/S0045-7825(02)00336-5](https://doi.org/10.1016/S0045-7825(02)00336-5), [author-uploaded text](https://www.researchgate.net/publication/222546721_Optimal_solid_shells_for_non-linear_analyses_of_multilayer_composites_II_Dynamics), §5.1.2/Figs.9–12 | Relevant text 확인; 하중 figure·곡선 digitization 미수행, exact benchmark reproduction 미정 |

Numerical Review에 인계할 결정은 scaling·conditioning/hidden-defect 판정·FD sweep과
시간 수렴 구간·case 최소 조합이다. Formulation에 남는 사항은 singular-mass regularity의
실제 shell 적용과 full-coordinate hidden constraint 검증이다. 사용자 결정에 따라 이번
HHT 구현은 consistent mass를 사용한다. 위 lumped 비교 근거는 연구 이력이며 현재 지원 요구가 아니다.