FESADev/docs/hht-implicit-dynamics/numerical-review-probes.md

# HHT 독립 시간정확도 검산

[검토 보고서](numerical-review.md)의 상세 evidence다. 2026-09-22 Python 표준 라이브러리만 사용했다. FESA solver, Abaqus 또는 reference 비교 실행이 아니다. 아래 오차는 승인 tolerance를 새로 정하는 값이 아니라 정확해와의 수렴 관찰값이다.

## 모델과 결과

1. Smooth SDOF: `T=.5(1+q²)v²`, `c=qv²`, `q_exact=.4+sin(t)`, `f=(1+q_exact²)a_exact+q_exact v_exact²`, `0≤t≤1`. 초기 `q=.4,v=1,a=0`. 입력은 일관된 무차원 검산 단위다. `alpha=-.05,beta=.275625,gamma=.55`.
2. Moving-null: `x(q)=q/|q|`, `T=.5|Dx(q)v|²`, `U=.5(|q|-1)²`, `q∈R²\{0}`. `r=|q|`, `d=q/r`, `M=(I-ddᵀ)/r²`, `c=-2(q·v)Mv/r²`, `f_int=(r-1)d`. `ker M=span(q)`가 회전한다. 정확해 `q=(cos t,sin t)`, `v=(-sin t,cos t)`, `a=-q`; 초기 위치·속도·가속도 대수조건은 일관된다. Endpoint 반경 대수조건 `r=1`은 각도 parameterization으로 정확하게 강제하며 tangent 평형만 Newton으로 푼다. 이 모델은 shell/drilling 구현이나 일반 singular-DAE 정리의 대체물이 아니다.

| SDOF 방법 | `h=1/512` U 오차 | `h=1/1024` U 오차 | 관측 차수 |
| --- | ---: | ---: | ---: |
| Naive `M1 A1+(1+alpha)g1-alpha g0=0` | 1.0272820368e-5 | 5.2200068232e-6 | 0.976709 |
| Weighted mass `((1+alpha)M1-alpha M0)A1+(1+alpha)g1-alpha g0=0` | 3.2914014847e-7 | 8.2299941973e-8 | 1.999739 |
| Endpoint physical `M1 a1+g1=0`, `A1=(1+alpha)a1-alpha a0` | 2.7561411953e-7 | 6.8902116013e-8 | 2.000030 |

`g=c-f`. 마지막 방법의 physical acceleration 오차는 `2.3176116737e-7 → 5.7545108678e-8`로 줄지만 algorithmic `A` 오차는 `5.2612378749e-5 → 2.6344436092e-5`다. 따라서 `A`를 endpoint physical acceleration으로 보고하면 안 된다.

Moving-null 마지막 방법은 같은 두 h에서 U 오차 `4.0606413493e-7 → 1.0152120526e-7`, V 오차 `4.1505360817e-7 → 1.0376522277e-7`, physical a 오차 `4.1511310820e-7 → 1.0262747089e-7`를 보였다. `max|q·v|`는 `8.4282e-11 → 1.0594e-11`, `max|q·a+v²|`는 `7.9368e-7 → 1.9445e-7`다. Endpoint 위치 대수조건을 지켜도 미분된 속도·가속도 조건이 매 frame 정확히 0이 되는 것은 아니다. 이 결함의 시간수렴과 초기 일관성을 별도로 검증해야 한다.

## 독립 수학 대조

Smooth 물리식 `M(q)a+g(q,v,t)=0`에서 Newmark 내부 가속도의 정확해 대응은 `A=a+alpha*h*adot+O(h²)`다. Naive 식의 일차 결함은

`alpha*h*(M*adot+gdot)=-alpha*h*Mdot*a`.

따라서 `Mdot*a≠0`인 일반 문제에서 상수 M 증명을 그대로 사용할 수 없다. Weighted mass는 이 항을 상쇄하지만 `alpha∈(-1,0)`이면 서로 다른 두 시점의 양의 준정부호 질량의 양의 결합이다. 그 영공간은 두 영공간의 교집합이므로, 영공간이 회전하면 rank가 늘어날 수 있다. 단순 SDOF의 2차 결과가 이를 정당화하지 않는다.

Physical/algorithmic 가속도 분리는 [Brüls–Arnold 원문](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf)의 Eq. (5), (7), (8)에 `alpha_m=0,alpha_f=-alpha`를 대입한 것이다. 원문의 p.6 equivalent-model/index-3 가정은 비특이 mass 및 제약 정칙성을 사용한다. Singular drilling mass에 그 증명을 직접 적용하지 않는다. 이 문서의 제조해와 결함 계산은 독립 검산이다.

## 재현 코드와 명령

Workspace root에서 다음 명령으로 아래 Python block만 실행한다. 외부 package, 파일 생성 또는 solver 실행은 없다.

```powershell
python -c 'from pathlib import Path; s=Path("docs/hht-implicit-dynamics/numerical-review-probes.md").read_text(encoding="utf-8"); exec(s.split(chr(96)*3+"python\n",1)[1].split(chr(96)*3,1)[0])'
```

```python
import math

alpha = -.05
beta = (1-alpha)**2/4
gamma = .5-alpha

def exact(t):
    return .4+math.sin(t), math.cos(t), -math.sin(t)

def force(t):
    q, v, a = exact(t)
    return (1+q*q)*a+q*v*v

def scalar(h, variant):
    q, v, ap = exact(0)
    A = ap
    for n in range(round(1/h)):
        t = (n+1)*h
        m0 = 1+q*q
        g0 = q*v*v-force(n*h)
        oldap = ap
        pred = q+h*v+h*h*(.5-beta)*A
        vpred = v+h*(1-gamma)*A
        def calc(x):
            an = (x-pred)/(beta*h*h)
            vn = vpred+gamma*h*an
            m = 1+x*x
            if variant == 'physical':
                pn = (an+alpha*oldap)/(1+alpha)
                r = m*pn+x*vn*vn-force(t)
                dr = 2*x*pn+m/((1+alpha)*beta*h*h)+vn*vn+2*x*vn*gamma/(beta*h)
            else:
                me = m if variant == 'naive' else (1+alpha)*m-alpha*m0
                dm = 2*x if variant == 'naive' else (1+alpha)*2*x
                r = me*an+(1+alpha)*(x*vn*vn-force(t))-alpha*g0
                dr = dm*an+me/(beta*h*h)+(1+alpha)*(vn*vn+2*x*vn*gamma/(beta*h))
                pn = an
            return r, dr, vn, an, pn
        x = pred+beta*h*h*A
        for _ in range(20):
            r, dr, *_ = calc(x)
            dx = -r/dr
            x += dx
            if abs(dx) < 1e-14:
                break
        else:
            raise RuntimeError('scalar Newton failed')
        q = x
        _, _, v, A, ap = calc(q)
    qe, ve, ae = exact(1)
    return abs(q-qe), abs(v-ve), abs(ap-ae), abs(A-ae)

def moving_null(h):
    q, v, ap = [1.,0.], [0.,1.], [-1.,0.]
    A = ap[:]
    theta = 0.
    max_cv = max_ca = 0.
    dot = lambda x,y: sum(x[i]*y[i] for i in range(2))
    for _ in range(round(1/h)):
        pred = [q[i]+h*v[i]+h*h*(.5-beta)*A[i] for i in range(2)]
        vp = [v[i]+h*(1-gamma)*A[i] for i in range(2)]
        old = ap[:]
        theta += h
        def calc(angle):
            d = [math.cos(angle),math.sin(angle)]
            tangent = [-d[1],d[0]]
            an = [(d[i]-pred[i])/(beta*h*h) for i in range(2)]
            vn = [vp[i]+gamma*h*an[i] for i in range(2)]
            pn = [(an[i]+alpha*old[i])/(1+alpha) for i in range(2)]
            dv, tv = dot(d,vn), dot(tangent,vn)
            residual = dot(pn,tangent)-2*dv*tv
            derivative = 1/((1+alpha)*beta*h*h)-dot(pn,d)-2*(tv*tv+dv*(-dv+gamma/(beta*h)))
            return residual, derivative, d, vn, pn, an
        for _ in range(15):
            r, dr, *_ = calc(theta)
            dx = -r/dr
            theta += dx
            if abs(dx) < 1e-14:
                break
        else:
            raise RuntimeError('moving-null Newton failed')
        _, _, q, v, ap, A = calc(theta)
        max_cv = max(max_cv,abs(dot(q,v)))
        max_ca = max(max_ca,abs(dot(q,ap)+dot(v,v)))
    qe = [math.cos(1),math.sin(1)]
    ve = [-math.sin(1),math.cos(1)]
    ae = [-x for x in qe]
    errors = [math.hypot(*[x[i]-y[i] for i in range(2)]) for x,y in ((q,qe),(v,ve),(ap,ae))]
    return *errors, max_cv, max_ca

for variant in ('naive','weighted_mass','physical','moving_null'):
    print(variant)
    previous = None
    for k in range(5,11):
        h = 2**(-k)
        errors = moving_null(h) if variant == 'moving_null' else scalar(h,variant)
        order = math.log(previous/errors[0],2) if previous else float('nan')
        print(f'h={h:.9g} errors='+','.join(f'{x:.10e}' for x in errors)+f' pU={order:.6f}')
        previous = errors[0]
```

Newton correction cutoff `1e-14`는 이 독립 double-precision probe의 반복 종료값이다. Solver convergence 또는 reference acceptance tolerance로 채택하지 않는다. 미세 h의 acceleration 오차에는 displacement 차분의 roundoff가 포함되므로 표의 관측값을 정확한 수학 상수로 해석하지 않는다.

## 후속 full Cartesian moving-null 검산

추가 Numerical 작업에서 메인이 [full_coordinate_probe.py](../../.harness/build/hht-numerical/full_coordinate_probe.py)를 작성했고 수치 검토자가 residual/Jacobian 식을 읽고 아래 명령을 재실행했다. Exit `0`, 2026-09-22. 이 script는 ignored build evidence이며 FESA source가 아니다. 위 angle-reduced probe와 달리 q의 두 Cartesian 성분을 직접 Newton으로 풀며 radial projection·normalization을 하지 않는다. 같은 `x=q/|q|` 모델의 exact 2×2 Jacobian을 사용한다.

```powershell
python .harness/build/hht-numerical/full_coordinate_probe.py
```

오차는 각 h의 **모든 시점**에서 q/v/physical a의 2-norm 최대다. 앞 절의 종단 오차와 다른 측정량이다.

| h | max q error | max v error | max physical a error | 직전 h 대비 q/v/a 차수 |
| --- | ---: | ---: | ---: | --- |
| 1/128 | 6.4948801161e-6 | 6.6395874231e-6 | 1.2690547915e-5 | 1.99932 / 1.99962 / 2.01024 |
| 1/256 | 1.6240799623e-6 | 1.6601117645e-6 | 3.1667175260e-6 | 1.99968 / 1.99981 / 2.00269 |
| 1/512 | 4.0606299832e-7 | 4.1505196280e-7 | 7.9259171465e-7 | 1.99985 / 1.99992 / 1.99834 |

마지막 grid의 최대 `||q|-1|=2.22045e-16`, `|q·v|=8.39836e-11`, `|q·a+v²|=7.91118e-7`. 상쇄 전 `Ma`의 entry products·c·f_int를 사용한 최대 normwise residual은 `7.94824e-11`, Newton 최대 3 corrections, 마지막 correction 최대 `4.00220e-16`이었다. 이 probe의 종료는 correction `<=1e-13`이며 N2의 모든 AND gate를 구현한 solver가 아니다. 위 관찰을 N2–N7 전체 통과나 실제 shell 통과로 확대하지 않는다.

별도 trial `q=(.85,.52)`, predictor `(.83,.50)`, velocity predictor `(-.35,.71)`, old physical a `(-.7,-.3)`, h=.1, direction `(.6,-.8)`의 exact effective Jacobian을 중앙차분했다. N5의 j=8,10,12,14,16,18,20 순서에서 상대 오차는 `[3.13656e-5,1.96035e-6,1.22522e-7,7.65779e-9,4.77929e-10,3.25134e-11,1.04250e-11]`이다. 연속한 내부 sweep 점 j=12/14 등이 `1e-6`을 만족한다. 이는 이 2D trial/direction에 대한 미분 evidence이며 shell D3 oracle를 대체하지 않는다.

Full-coordinate 결과는 새로운 [C04b](verification-cases.md)의 moving-null 구조와 max-time 측정을 지지한다. C04b의 전체 prescribed grids, 초기화 실패 mutation, regularity·scale·반복오차 분리 등 모든 판정을 실행한 것은 아니다. 더 미세한 h에서 residual이 증가할 수 있으므로 h=1/1024까지 같은 residual 기준으로 통과했다는 주장을 하지 않는다.