back to table · edit · history · where entries came from · files · download
4002 bytes, as of the version from 2026-10-30 06:27 (current). Recorded here, not run.
"""Taylor continuation of Lane–Emden backgrounds and their tidal fields.
Adapted from the live T618 solver.py (archived and hashed in this run).
Use named Sage fields after numberdb.sage. No certified ODE error bound.
Y=xi H'/H; H is regular as xi**ell at the center. The surface coordinate
is t=sqrt(theta), so the tidal density derivative has the integer power
2*n-1 for each positive selected index, including n=1/2.
"""
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):
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)
c=[F(1)]
for j in range(1,order+1):
c.append(sum((n*i-j)*a[i]*c[j-i] for i in range(1,j+1))/j)
ys=[]
for ell in (2,3):
v=[F(1)]
for k in range(order):
v.append(-n*sum(c[j]*v[k-j] for j in range(k+1))/((2*k+2)*(2*k+2*ell+3)))
ys.append(F(ell)+horner([2*k*v[k] for k in range(len(v))],x*x)/horner(v,x*x))
return theta,q,ys
def surface_series(t,x0,q0,y0,n,order):
F=x0.parent();z=F(0);m=int(2*n+1)
x=[x0];q=[q0];ys=[[v] for v in y0]
xx=[];inv=[];c=[];d=[];e=[];f=[]
yy=[[],[]];g=[[],[]];eg=[[],[]]
power=[F(comb(m,j))*t**(m-j) for j in range(m+1)]
tide=[F(comb(m-2,j))*t**(m-2-j) for j in range(m-1)] if n else []
nn=F(n)
for k in range(order):
xx.append(sum((x[j]*x[k-j] for j in range(k+1)),z))
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(x[j]*inv[k-j] for j in range(k+1)))
f.append(sum(x[j]*c[k-j] for j in range(k+1)))
for i,ell in enumerate((2,3)):
y=ys[i]
yy[i].append(sum(y[j]*y[k-j] for j in range(k+1)))
g[i].append((F(ell*(ell+1)) if k==0 else z)-y[k]-yy[i][k])
eg[i].append(sum(e[j]*g[i][k-j] for j in range(k+1)))
forcing=2*nn*sum((tide[j]*f[k-j] for j in range(min(m-2,k)+1)),z) if n else z
y.append((-2*(t*eg[i][k]+(eg[i][k-1] if k else z))+forcing)/(k+1))
x.append(-2*(t*c[k]+(c[k-1] if k else z))/(k+1))
q.append(-2*sum(power[j]*d[k-j] for j in range(min(m,k)+1))/(k+1))
return [x,q]+ys
def integrate(n,working=150,order=160,divisor=12,start='1/2'):
n=QQ(n)
if n<0 or n>=5 or 2*n not in ZZ:raise ValueError('Need integer or half-integer 0 <= n < 5')
F=RealField(numberdb.bits(working));x=F(QQ(start))
theta,q,ys=center(n,F,2*working,x)
state=[x,q]+ys;t=theta.sqrt();count=0;tol=F(10)**(-(working-20));max_tail=F(0)
while t>0:
arrays=surface_series(t,state[0],state[1],state[2:],n,order)
cap=min(F(1)/divisor,(1-t)/3,t)
for arr,v in zip(arrays,state):
scale=max(F(1),abs(v))
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;new=[horner(a,h) for a in arrays]
max_tail=max(max_tail,max(abs(a[-1]*h**order)/max(F(1),abs(v)) for a,v in zip(arrays,new)))
assert new[0]>state[0] and new[1]>0
state=new;t+=h;count+=1
if count>3000:raise RuntimeError('Too many steps')
x,q=state[:2];ys=state[2:]
if n==0:ys=[y-3 for y in ys] # surface density jump, Y(R+)=Y(R-)-3
loves=[(ell-y)/(2*(ell+1+y)) for ell,y in zip((2,3),ys)]
return loves,{'radius':str(x),'mass':str(q),'surface_Y':list(map(str,ys)),
'steps':count,'last_term_max':str(max_tail),'working':working,
'order':order,'divisor':divisor,'start':start}