diff --git a/docs/shell-transient-heat-transfer/formulation.md b/docs/shell-transient-heat-transfer/formulation.md index c7b8082..0fb4312 100644 --- a/docs/shell-transient-heat-transfer/formulation.md +++ b/docs/shell-transient-heat-transfer/formulation.md @@ -310,41 +310,57 @@ raw 초기 상태의 projection으로 바꾸지 않는다. 기본값 100이며 허용 증분 수의 상한이다. `d≥τ`이면 한 증분, `δ_0=τ`, 저장 시각 `τ`다. 그 외에는 **반복 덧셈 없이** 정수 `j`의 곱으로 정상 시각을 만든다. -정수비 근방은 다음 연산 bound로 한 번만 판정한다. 먼저 비율이 `INC+1` 이내인지 -확인하고 `τ`를 `[1/2,1)`에 놓는 같은 exact power-of-two로 `τ,d`를 재척도한다. -이 범위에서는 재척도된 `d`도 normal이며 아래 `p`의 overflow를 피할 수 있다. -아래 predicate는 재척도된 변수로 평가하며 결과와 시각·δ는 원래 단위를 유지한다. -`fl`은 binary64 nearest-even이고, -`rad(x)=max(x-nextafter(x,-∞), nextafter(x,+∞)-x)/2`는 해당 값의 보수적인 -rounding-cell 반폭이다. 반폭이 subnormal로 표현되지 않으면 powers-of-two scale에서 -계산한다. 여기서는 `q=fl(τ/d)`, `m=nearest_integer(q)`, `p=fl(m d)`를 쓴다. -`q`가 overflow하거나 `q>INC+1`이면 증분 수 초과다. 나머지 경우 `m≤2³²` 부근의 -정수는 double에 정확히 표현된다. Tie의 정수 선택도 nearest-even으로 고정한다. +정수비 근방의 동등성은 **입력 binary64로 반올림될 수 있는 실제 시간 구간**으로 판정한다. +입력 중심값을 오차 없는 물리량으로 간주하지 않는다. `fl`은 nearest-even이며 먼저 +`q=fl(τ/d)`를 구한다. `INC+1`은 uint64 범위에서 계산한다. +`q`가 overflow하거나 `q>INC+1`이면 증분 수 초과다. +나머지는 `m=nearest_integer(q)`를 후보로 쓴다. Tie는 nearest-even으로 고정하며 +`m≤2³²`는 정확히 표현된다. `INC`에 맞추려고 다른 정수 후보로 바꾸지 않는다. + +원래 입력 `x>0`의 predecessor/successor를 `x−,x+`라 하면 rounding cell은 \[ -B=d\,\operatorname{rad}(q)+\operatorname{rad}(p),\qquad -\widehat B=\operatorname{nextafter}(\operatorname{fl}(B),+\infty),\qquad -D=|\operatorname{fma}(-m,d,\tau)|,\qquad -\text{endpoint-close}\iff m\ge1\ \land\ D\le\widehat B. +I(x)=\left\langle (x^-+x)/2,\ (x+x^+)/2\right\rangle, +\qquad \text{endpoint-close}\iff m\ge1\ \land\ I(\tau)\cap mI(d)\ne\varnothing. \] -이는 임의 시간 epsilon이 아니다. 나눗셈 반올림의 dimensional uncertainty는 -`d rad(q)`, 곱셈의 uncertainty는 `rad(p)`이므로 합이 두 연산에 대한 보수적 bound다. -특히 quotient가 정수 `m`으로 반올림되는 경우 -`|τ-md|≤d rad(q)`이고, 정상 시각 곱 자체의 표시 불확실성은 `rad(p)`다. -Fused remainder로 곱의 선행 반올림을 제거하므로 `τ-fl(md)`의 이중 오차를 비교하지 않는다. -이 bound 안의 endpoint는 정상 마지막 증분과 같은 계수를 사용하는 명시적 정책이다. -실제 remainder가 이 범위 밖이면 작더라도 보존한다. 외부 입력 숫자의 임의 오차나 물리 -시간 tolerance를 추가하지 않는다. 재척도된 `rad(q)`는 power-of-two이고 `d rad(q)`는 -normal 범위의 exact binary scaling이다. 따라서 `fl(B)`는 덧셈 한 번의 반올림뿐이며 -`nextafter` 한 번이 보수적 상계를 준다. FMA 한 번의 잔여 반올림까지 포함하면 -accepted endpoint의 정확한 거리 bound는 `|τ-md|≤Bhat+rad(D)`다. 이렇게 **표현 가능한 -threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를 고정한다. FMA를 두 개의 -일반 연산으로 분해하거나 compiler 재결합에 맡기지 않는다. +꺾쇠의 두 끝점은 `x`의 binary64 significand 최하위 bit가 0(even)이면 포함하고, +1(odd)이면 제외한다. Subnormal도 같은 규칙이며 최소 양수의 predecessor는 0이다. +최대 finite의 successor는 부동소수점 `+∞` 대신 **수학적인 다음 lattice 값 `2^1024`**로 +둔다. 최대 finite significand는 odd이므로 overflow midpoint를 cell에 포함하지 않는다. +Midpoint를 double로 반올림하지 않고 아래 정수·지수 표현으로 정확히 구성한다. +입력 변환은 correctly-rounded nearest-even을 전제로 하며 임의 입력 오차는 더하지 않는다. -1. Endpoint-close이면 `N=m`. 아니면 `N=floor(q)+1=ceil(τ/d)`다. Quotient 반올림으로 - 정수 경계를 넘을 수 있는 경우는 위 endpoint-close bound 안에 이미 포함되므로 - 이 분기에서 floor는 안전하다. 먼저 nearest-integer branch를 판정하는 순서를 유지한다. +`I(τ)=(a_L,a_U)`, `mI(d)=(b_L,b_U)`의 끝점 포함 flag를 보존하여 +`L=max(a_L,b_L)`, `U=min(a_U,b_U)`를 구한다. `LU`면 nonclose다. +`L=U`면 그 접점이 **두 cell 모두에 포함될 때만** close다. 나눗셈·곱셈 오차와 입력 cell +폭을 중복 가산하는 epsilon은 없다. `q`는 정수 후보만 만들고, 겹침 판정에는 반올림된 +`fl(md)`를 사용하지 않는다. 같은 endpoint를 허용하지 않는 cell은 잔여시간이 작아도 구별한다. + +정확한 predicate는 임의정밀도 dependency 없이 다음 고정 폭 연산으로 구현할 수 있다. + +- `e=frexp(τ).exponent`, `s=2^(−e)`로 두어 `sτ∈[1/2,1)`로 만든다. **원래 입력의 + cell을 먼저 구성한 뒤** `s`를 곱한다. Normal화된 값의 prev/next로 cell을 다시 만들면 + 원래 subnormal 입력의 표현 범위를 줄이므로 금지한다. `s` 자체를 double로 만들 필요는 없다. +- `d<τ`, 유한 `q≤INC+1≤2³²`이면 정확한 `τ/d`는 `2³²+2^(−21)` 이하이므로 + `sd>2^(−34)`이다. Scaled 중심값과 원래 cell의 midpoint는 모두 lattice `2^(−88)`의 + 정수배다. Normal 입력의 최소 반간격은 이 lattice보다 작지 않으며, subnormal은 더 + 성긴 원래 간격을 유지한다. Binade 경계의 비대칭 간격도 따로 보존한다. +- `X=τs2^88`, `Y=ds2^88`와 네 cell endpoint를 원래 binary64 significand/exponent에서 + shift·덧셈으로 정확히 구성한다. Scaled cell 끝점은 2보다 작으므로 각 정수는 89 bit 이하다. + `m` 및 floor 교정용 `floor(q)+1`은 최대 `2³²+1`이므로 곱은 **122 bit 이하**다. +- 두 `uint64` 상당의 high/low word로 저장·비교한다. High word를 먼저 비교하고, 같으면 + low word를 비교한다. 덧셈/뺄셈의 carry/borrow와 shift의 word 경계를 명시적으로 처리한다. + 곱은 네 base-`2³²` digit와 정수 승수의 최대 두 digit로 계산할 수 있다. 각 + `a_i*b_j+result_digit+carry≤2^64−1`이므로 중간값도 uint64에 들어가며 다음 digit으로 + carry를 넘긴다. 최상위 128 bit 밖 carry가 없음을 검사한다. `__int128`, extended float, + compiler 연산 재결합 또는 arbitrary-precision library에 의존하지 않는다. + +1. Endpoint-close이면 `N=m`. 아니면 중심값의 정확한 `N=ceil(τ/d)`를 쓴다. + `f=floor(q)`에서 정수 비교 `fY≤X<(f+1)Y`를 확인하고 필요하면 `f`를 한 번 ±1 + 교정한다. `q≤2³²`의 division 오차는 1보다 작으므로 한 번이면 충분하다. + `X=fY`이면 `N=f`, 아니면 `N=f+1`이다. Division 반올림을 증분 수에 그대로 전달하지 않는다. + `N`은 `INC`와 비교하기 전 uint32로 축소하지 않는다. 2. `N>INC`이면 solve 전에 실패한다. 정상 저장 시각은 `t_j=fl(j d)` (`j=1,…,N−1`), 마지막은 **입력 그대로** `t_N=τ`다. 모든 시각이 유한하고 엄격히 증가하며 `t_{N−1}<τ`인지 preflight에서 확인한다. `t_0=0`이다. @@ -352,7 +368,9 @@ threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를 마지막 증분만 `δ=fl(τ-(N−1)d)`로 만든다. 이 remainder는 `fma(-(N−1),d,τ)`처럼 곱-차를 한 번에 반올림한다. `δ>0`·유한성을 검사한다. 저장된 `fl(jd)`끼리 다시 빼서 정상 `δ`를 만들지 않는다. 개념적인 격자 `j d`와 - 저장 표현의 차이 및 endpoint-close 조정은 위 rounding bound로 제한된다. + 저장 표현의 차이는 각 곱의 rounding cell로 제한된다. Endpoint-close 조정은 입력 cell + overlap으로 제한된다. 원래 단위의 `|τ−md|`는 해당 방향의 `τ` 반간격과 `m`배 `d` + 반간격 합 이내이고, 계산된 FMA remainder에는 한 번의 최종 rounding만 남는다. 4. 해당 증분의 `A`, RHS, full residual, reaction, `ΔE` 수지, BE 소산은 모두 **같은 실제 `δ`**를 사용한다. Endpoint-close는 동일 계수를 재사용하며 실제 shortened 마지막 `δ`가 `d`와 다를 때 refactorize한다. @@ -363,6 +381,107 @@ threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를 `nextafter`로 endpoint-close 경계를 양쪽에서 이동한 사례, 아주 작고 큰 시간 척도, `INC=N`/`N−1`, uint32 overflow 방지를 독립 검증한다. 이 격자는 adaptive time stepping이 아니다. +### 7.2 입력 cell 수정의 실제 검산 기록 + +2026-09-29 Python 3.11 `Fraction`을 독립 exact-dyadic oracle로 사용한 아래 scratch를 +실행했다(exit 0). 이는 문서 수식 검산이며 solver build/test나 수치 correctness 승인이 아니다. +기존 `Bhat=nextafter(fl(d rad(q)+rad(fl(md))),+∞)`는 입력 cell을 누락하여 첫 사례를 +nonclose로 오판했다. Scratch의 임의정밀도는 검산용이며 production dependency가 아니다. + +```text +old-failure: m=31, exactD=8.326672684688674e-16 > Bhat=7.507374033534679e-16 +corrected: close=True, N=31; INC=31 accept, INC=30 reject +tau=nextafter(.9,+inf), d=.3: close=False, N=4, delta=1.6653345369377348e-16 +tau=.9, d=.3: close=True, N=3 +tau=20*minsub, d=6*minsub: even/even midpoint included, close=True, N=3 +tau=17*minsub, d=5*minsub: odd/odd midpoint excluded, close=False, N=4 +power-of-two scales [-900,-400,0,400,900]: 5/5 +maxfinite, original-subnormal, binade cases: 5/5 +seed41 input-cell versus fixed-lattice predicate: 10000/10000 +uint32 INC boundary: N=4294967295 accept; N=4294967296 reject +``` + +
+실행한 재현 scratch (Python 표준 라이브러리) + +```python +import math, random, struct +from fractions import Fraction as F +def cell(x): + a=F(x); p=F(math.nextafter(x,0.0)) + nx=math.nextafter(x,math.inf) + n=F(nx) if math.isfinite(nx) else a+(a-p) + even=(struct.unpack('>Q',struct.pack('>d',x))[0]&1)==0 + return (a+p)/2,(a+n)/2,even +def overlap(a,b): + lo=max(a[0],b[0]); hi=min(a[1],b[1]) + return lo=t: return 1 + if predicate(t,d): return round(t/d) + z=F(t)/F(d); return -(-z.numerator//z.denominator) +def count_fixed(t,d): + if d>=t: return 1 + if fixed(t,d): return round(t/d) + e=math.frexp(t)[1]; s=F(2)**(-e+88) + x,y=int(F(t)*s),int(F(d)*s); f=math.floor(t/d) + if f*y>x: f-=1 + elif (f+1)*y<=x: f+=1 + assert f*y<=x<(f+1)*y + return f if x==f*y else f+1 +def rad(x): + return max(F(x)-F(math.nextafter(x,-math.inf)), + F(math.nextafter(x,math.inf))-F(x))/2 +t,d=5.351455173733616,.17262758624947144 +m=round(t/d); q=t/d; p=m*d +old=math.nextafter(float(F(d)*rad(q)+rad(p)),math.inf) +distance=abs(F(t)-m*F(d)) +assert distance>F(old) and predicate(t,d) and fixed(t,d) +print('old-failure:',m,float(distance),old,distance<=F(old)) +print('corrected:',predicate(t,d),'N=',count(t,d),'INC31=',count(t,d)<=31,'INC30=',count(t,d)<=30) +t2=math.nextafter(.9,math.inf) +print('real-remainder:',predicate(t2,.3),'N=',count(t2,.3),'delta=',float(F(t2)-3*F(.3))) +print('decimal:',predicate(.9,.3),count(.9,.3)) +tiny=math.ulp(0.0) +for a,b in [(20,6),(17,5)]: + print('midpoint:',a,b,cell(a*tiny)[2],cell(b*tiny)[2],predicate(a*tiny,b*tiny),count(a*tiny,b*tiny)) +for e in [-900,-400,0,400,900]: + tt,dd=math.ldexp(t,e),math.ldexp(d,e) + assert predicate(tt,dd) and fixed(tt,dd) and count(tt,dd)==31 +print('power-of-two-scales: 5/5') +maxf=float.fromhex('0x1.fffffffffffffp+1023') +for tt,dd in [(maxf,maxf/31),(20*tiny,6*tiny),(17*tiny,5*tiny), + (1.,math.nextafter(.25,0.)),(1.,math.nextafter(.25,math.inf))]: + assert predicate(tt,dd)==fixed(tt,dd) + assert count(tt,dd)==count_fixed(tt,dd) +print('extreme/binade cases: 5/5') +rng=random.Random(41) +for _ in range(10000): + dd=math.ldexp(rng.uniform(.5,1),rng.randint(-900,900)) + mm=rng.randint(2,2**32) + tt=mm*dd + for _ in range(rng.randint(0,2)): + tt=math.nextafter(tt,rng.choice([0.,math.inf])) + assert predicate(tt,dd)==fixed(tt,dd) +print('seed41 cell-vs-fixed: 10000/10000') +assert count(float(2**32-1),1.)==2**32-1 +assert count(float(2**32),1.)==2**32 +print('uint32 INC limit: N=4294967295 accept; N=4294967296 reject') +``` + +
+ ## 8. 반력 열률, 에너지와 안정성 반력은 하중 제거 전 full residual에서 구한다.