back to table · edit · history · where entries came from · files · download
5645 bytes, as of the version from 2026-10-30 06:50 (current). Recorded here, not run.
"""Independent radial-coordinate axial inertia checks -- numberdb.org/T616.
Run: sage -python independent_checks.py, beside numerical-evidence.json.
Adapted from T614's live radial-coordinate check, with its pressure moment
replaced by the axial density moment. No production solver is imported.
The integration uses xi, theta and q, with exact rational central coefficients,
and stops at positive theta < 1e-30. A leading surface-tail correction is
included. Endpoint and Taylor errors are heuristic, not rigorous bounds.
Analytic n=0,1 controls run first. The Yip-Leung Appendix C.4 approximant
provides a separate literature check over its claimed 0<=n<=3 range.
This build check does not replace the controller's independent critique.
"""
import numberdb.sage as numberdb
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from fractions import Fraction
from pathlib import Path
import time,json
def evalpoly(coeff,x):
ans=x.parent()(0)
for c in reversed(coeff):ans=ans*x+c
return ans
def initial(index,F):
# Exact rational central coefficients, computed independently.
a=[Fraction(1)];power=[Fraction(1)];n=Fraction(index)
for k in range(200):
a.append(-power[k]/((2*k+2)*(2*k+3)))
j=k+1
power.append(sum((n*i-(j-i))*a[i]*power[j-i] for i in range(1,j+1))/j)
r=F(1)/4;u=r*r
af=[F(QQ(c.numerator)/QQ(c.denominator)) for c in a]
y=evalpoly(af,u)
q=-r*evalpoly([2*k*v for k,v in enumerate(af)],u)
return r,y,q
def radial(index,dps=70,degree=80):
F=RealField(numberdb.bits(dps));n=F(QQ(index));r,y,q=initial(index,F)
tol=F(10)**(-dps+20);stop=F(10)**-30;steps=0
# Integrate the central series to initialize J independently.
# J is initialized by high-precision quadrature of the analytic local IVP.
import mpmath as mp
with mp.workdps(dps):
nn=mp.mpf(str(n));rr=mp.mpf(str(r))
aa=[mp.mpf(1)];bb=[mp.mpf(1)]
for k in range(200):
aa.append(-bb[k]/((2*k+2)*(2*k+3)))
j=k+1
bb.append(sum(((nn+1)*i-j)*aa[i]*bb[j-i] for i in range(1,j+1))/j)
def central(z):
val=mp.polyval(list(reversed(aa)),z*z)
return z**4*val**nn
J=F(str(mp.quad(central,[0,rr])))
while y>stop:
ys=[y];qs=[q];power=[y**n]
recip=[F((-1)**k*(k+1))/r**(k+2) for k in range(degree)]
for k in range(degree):
ys.append(-sum(qs[i]*recip[k-i] for i in range(k+1))/(k+1))
rhs=r*r*power[k]
if k>=1:rhs+=2*r*power[k-1]
if k>=2:rhs+=power[k-2]
qs.append(rhs/(k+1))
j=k+1
power.append(sum(((n+1)*i-j)*ys[i]*power[j-i] for i in range(1,j+1))/(j*y))
jcoeff=[]
r4=[r**4,4*r**3,6*r*r,4*r,F(1)]
for k in range(degree):
jcoeff.append(sum(r4[i]*power[k-i] for i in range(min(4,k)+1))/(k+1))
h=min(r/3, y*r*r/(4*q),F(1)/8)
for coeff,scale in [(ys,abs(y)),(qs,abs(q))]:
for k in range(degree-5,degree+1):
if coeff[k]:h=min(h,(tol*scale/abs(coeff[k]))**(F(1)/k)/2)
new_y=evalpoly(ys,h);new_q=evalpoly(qs,h)
assert 0<new_y<y and new_q>0
J+=h*evalpoly(jcoeff,h)
r+=h;y=new_y;q=new_q;steps+=1
if steps>6000:raise RuntimeError('too many radial steps')
# Exact zero of theta_ext(rho)=y+q(1/rho-1/r).
root=r/(1-r*y/q)
tail=r**4*y**(n+1)/(q*(n+1))
mass=q+tail
J+=r**6*y**(n+1)/(q*(n+1))
value=2*J/(3*mass*root*root)
return value,{'radius':str(root),'mass':str(mass),'steps':steps,'terminal_theta':str(y),'dps':dps,'degree':degree,'mass_tail':str(tail),'density_moment':str(J)}
def main():
evidence_path=Path(__file__).with_name('numerical-evidence.json')
expected=json.loads(evidence_path.read_text())['values']
output={}
for n in ['0','1','1/2','3/2','2','5/2','3','7/2','4','9/2']:
started=time.time();value,meta=radial(n);F=value.parent()
if n in ['0','1']:
control=F(2)/5 if n=='0' else F(2)/3-4/F.pi()**2
err=abs(value-control)/control
assert err<F(10)**-40,(n,err)
meta['control_relative_error']=str(err)
assert abs(value*(1+F(10)**-25)-control)/control>F(10)**-40
error=abs(value-F(expected[n]['210']['inertia']))/value
assert error<F(10)**-40,(n,error)
meta['relative_difference_from_surface_solver']=str(error)
meta['value']=str(value);meta['seconds']=time.time()-started;output[n]=meta
print('CHECK',n,json.dumps(meta),flush=True)
# Published approximation, with its documented finite accuracy, not a
# high-precision oracle. The reported 0.21 percent is rounded: the
# actual n=3 error is 0.2141340818 percent, which rounds to 0.21.
# Record that discrepancy explicitly; use 0.215 percent for the check.
F=RealField(256)
for n in ['0','1/2','1','3/2','2','5/2','3']:
x=F(QQ(n));b=[F(24)/5,F('3.452950061777'),F('0.08417922198458'),F('0.01816712663676'),F('519.3645962035'),F('-155.4882028232')]
pade=F(2)/3*(5-x)**3*(b[0]+x*b[1]+x*x*b[2]+x**3*b[3])/(1000+x*b[4]+x*x*b[5])
val=F(output[n]['value']);err=abs(pade-val)/val
assert err<F('0.00215'),(n,err)
output[n]['literature_pade_relative_error']=str(err)
output[n]['literature_rounded_percent_limit']='0.21 (interpreted as rounded, not a strict bound)'
output[n]['test_percent_tolerance']='0.215'
if n=='3':
assert F('0.205')<100*err<F('0.215')
print('DATA '+json.dumps(output),flush=True)
if __name__=='__main__':main()