back to table · edit · history · where entries came from · files · download
5254 bytes, as of the version from 2026-10-30 02:01 (current). Recorded here, not run.
"""Independent radial-coordinate mass and virial checks of T614.
Run with SageMath: sage -python independent_checks.py
Uses theta and -xi² theta' as functions of xi, and stops while theta>0.
Stop at theta<1e-30, add the leading mass tail, and accumulate
J = integral xi^2 theta^(n+1) dxi. The independent Pohozaev/virial
identity is J=(n+1)*mu^2/((5-n)*xi_1). Errors from endpoint approximation
are heuristic. Exact n=0,1 are controls; n=5 uses a separate quadrature.
This is build validation, not the campaign'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*z*val**(nn+1)
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=[]
yp=[sum(ys[i]*power[k-i] for i in range(k+1)) for k in range(degree)]
for k in range(degree):
jcoeff.append((r*r*yp[k]+(2*r*yp[k-1] if k else 0)+(yp[k-2] if k>=2 else 0))/(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**4*y**(n+2)/(q*(n+2))
virial=(n+1)*mass**2/((5-n)*root)
err=abs(J-virial)/abs(virial)
assert err<F(10)**-40, ('virial',index,err)
return mass,{'radius':str(root),'steps':steps,'terminal_theta':str(y),'dps':dps,'degree':degree,'mass_tail':str(tail),'pressure_integral':str(J),'virial_relative_error':str(err)}
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=2*F(6).sqrt() if n=='0' else F.pi()
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']['mass']))/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)
import mpmath as mp
with mp.workdps(140):
# x=sqrt(3)*u/sqrt(1-u^2) maps the infinite n=5 integral to [0,1].
# Substitution reduces the integrand to 3*sqrt(3)*u^2.
# Run quadrature with its own nonzero known control first.
assert abs(mp.quad(lambda u:u*u,[0,1])-mp.mpf(1)/3)<mp.mpf('1e-135')
mu=mp.quad(lambda u:3*mp.sqrt(3)*u*u,[0,1])
err=abs(mu-mp.sqrt(3))/mp.sqrt(3)
assert err<mp.mpf('1e-135')
assert abs(mu*(1+mp.mpf('1e-25'))-mp.sqrt(3))>mp.mpf('1e-40')
output['5']={'value':str(mu),'quadrature_relative_error':str(err),'working_dps':140}
print('DATA '+json.dumps(output),flush=True)
if __name__=='__main__':main()