back to table · edit · history · where entries came from · files · download
4212 bytes, as of the version from 2026-10-30 06:50 (current). Recorded here, not run.
"""Axial inertia of Lane–Emden polytropes -- numberdb.org/T616.
Adapted from the live T613 solver.py attachment. Adds the density moment
J=integral xi^4 theta^n dxi and H=integral xi^2 theta dxi. The identity
J=mu*xi_1^2-6H provides a separate integral check at the surface.
Only numberdb.sage initialises Sage. No symbolic power-series machinery is used.
All divisions involving computed quantities are in a named RealField.
"""
import numberdb.sage as numberdb
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from math import comb
def horner(a,h):
y=a[-1]
for v in a[-2::-1]: y=y*h+v
return y
def center(n,F,order,x):
"""theta=sum a[k]*xi**(2k); b=theta**n, with a[0]=b[0]=1."""
n=F(n);a=[F(1)];b=[F(1)]
for k in range(order):
a.append(-b[k]/((2*k+2)*(2*k+3)))
j=k+1
b.append(sum(((n+1)*i-j)*a[i]*b[j-i] for i in range(1,j+1))/j)
theta=horner(a,x*x)
q=-x*horner([2*k*a[k] for k in range(len(a))],x*x)
J=x**5*horner([b[k]/(2*k+5) for k in range(len(b))],x*x)
H=x**3*horner([a[k]/(2*k+3) for k in range(len(a))],x*x)
return theta,q,J,H
def surface_series(t,x0,q0,j0,h0,m,order):
"""x'=-2t x²/q, q'=-2t**m x⁴/q; derivatives with respect to t."""
F=x0.parent();zero=F(0)
x=[x0];q=[q0];js=[j0];hs=[h0];xx=[];inv=[];c=[];d=[];e=[]
power=[F(comb(m,j))*t**(m-j) for j in range(m+1)]
cube=[t**3,3*t*t,3*t,F(1)]
for k in range(order):
xx.append(sum((x[j]*x[k-j] for j in range(k+1)),zero))
inv.append(F(1)/q0 if k==0 else -sum(q[j]*inv[k-j] for j in range(1,k+1))/q0)
c.append(sum(xx[j]*inv[k-j] for j in range(k+1)))
d.append(sum(xx[j]*c[k-j] for j in range(k+1)))
e.append(sum(xx[j]*d[k-j] for j in range(k+1)))
js.append(-2*sum(power[j]*e[k-j] for j in range(min(m,k)+1))/(k+1))
hs.append(-2*sum(cube[j]*d[k-j] for j in range(min(3,k)+1))/(k+1))
x.append(-2*(t*c[k]+(c[k-1] if k else zero))/(k+1))
q.append(-2*sum(power[j]*d[k-j] for j in range(min(m,k)+1))/(k+1))
return x,q,js,hs
def inertia(n,working=150,order=160,divisor=12,start='1/2',trace=False):
"""Return a fixed-precision approximation; there is no certified ODE bound.
Tail-based steps use the final eight Taylor coefficients. Their predicted
contribution is <10**(-(working-20)); the cap also changes with divisor.
Precision, order, initial matching radius and step cap are varied in checks.
"""
n=QQ(n);m=int(2*n+1)
if n<0 or n>=5 or 2*n not in ZZ:
raise ValueError('This solver supports integer/half-integer 0 <= n < 5')
F=RealField(numberdb.bits(working));x=F(QQ(start))
theta,q,J,H=center(n,F,2*working,x)
t=theta.sqrt();count=0;tol=F(10)**(-(working-20));max_tail=F(0)
while t>0:
a,b,c,d=surface_series(t,x,q,J,H,m,order)
# Center singularity at t=1; possible distant complex singularities
# are additionally detected by the tail estimate, not assumed absent.
cap=min(F(1)/divisor,(1-t)/3,t)
for arr,scale in [(a,max(F(1),abs(x))),(b,max(F(1),abs(q))),(c,max(F(1),abs(J))),(d,max(F(1),abs(H)))]:
for k in range(order-7,order+1):
if arr[k]:
cap=min(cap,(tol*scale/abs(arr[k]))**(F(1)/k)/2)
h=-cap
xn=horner(a,h);qn=horner(b,h)
tail=max(abs(a[-1]*h**order)/max(F(1),abs(xn)),abs(b[-1]*h**order)/max(F(1),abs(qn)))
max_tail=max(max_tail,tail)
assert xn>x and qn>0 and qn>=q*(1-F(10)**(-working+10))
J,H=horner(c,h),horner(d,h)
x,q=xn,qn;t+=h;count+=1
if count>3000:raise RuntimeError('too many continuation steps')
meta={'steps':count,'last_term_max':str(max_tail),'working':working,'order':order,'divisor':divisor,'start':start}
value=2*J/(3*q*x*x)
alternate=F(2)/3-4*H/(q*x*x)
discrepancy=abs(value-alternate)/abs(value)
assert discrepancy<F(10)**(-working+30), ('integral identity',n,discrepancy)
meta.update(radius=str(x),mass=str(q),density_moment=str(J),enthalpy_moment=str(H),identity_relative_error=str(discrepancy))
return (value,meta) if trace else value