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

HHT 독립 시간정확도 검산

검토 보고서의 상세 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 원문의 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 실행은 없다.

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])'

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를 작성했고 수치 검토자가 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을 사용한다.

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