back to table · edit · history · where entries came from · files · download
3678 bytes, as of the version from 2026-10-30 04:30 (current). Recorded here, not run.
"""Independent radial-coordinate mass and pressure integrals for T617.
Imported by independent_checks.py. Adapted from T614/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 in the calling checks.
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)}