back to table · edit · history · where entries came from · files · download
2645 bytes, as of the version from 2026-10-30 06:27 (current). Recorded here, not run.
"""Independent radial-coordinate DOP853 checks of T619.
No production solver imports. Integrate the modified potential y=xi H'/H-
xi^3 theta^n/q (Yip and Leung eq. 24), which handles n=0's density jump.
Double precision is used only for this independent check, never for table data.
Two tighter tolerances and initial radii are compared. Controls run first.
"""
import json,math,sys
import numpy as np
from scipy.integrate import solve_ivp
INDICES=('0','1','1/2','3/2','2','5/2','3','7/2','4','9/2')
def radial(n,rtol=3e-13,start=1e-4):
# Independent central expansion through xi^4 for theta and q/xi^3.
x=start
theta=1-x*x/6+n*x**4/120
q=x**3*(1/3-n*x*x/30+n*(8*n-5)*x**4/2520)
# Modified logarithmic derivative: y=ell-3 + [n/5-n/(2ell+3)]xi² + O(xi⁴).
ys=[ell-3+(n/5-n/(2*ell+3))*x*x for ell in (2,3)]
def fun(x,v):
th,q,y2,y3=v
rho=max(th,0.)**n if n else 1.
s=x**3*rho/q
return [-q/x**2,x*x*rho]+[(ell*(ell+1)-y*(y+1)-2*s*(y+2))/x for ell,y in zip((2,3),(y2,y3))]
def event(x,v):return v[0]
event.terminal=True;event.direction=-1
sol=solve_ivp(fun,(x,1000),[theta,q]+ys,method='DOP853',events=event,rtol=rtol,atol=rtol/100,max_step=.15)
assert sol.success and len(sol.t_events[0])==1
radius=sol.t_events[0][0];th,q,*ys=sol.y_events[0][0]
return [(ell-y)/(2*(ell+1+y)) for ell,y in zip((2,3),ys)],radius,q,len(sol.t)
def main(path):
source=json.load(open(path));out={'method':'radial modified-potential DOP853, independent of production Taylor source','models':{},'controls':{}}
for nstr in INDICES:
from fractions import Fraction
n=float(Fraction(nstr))
a=radial(n,1e-11,2e-4);b=radial(n,3e-13,1e-4)
expected=list(map(float,source['models'][nstr]['runs']['210']['values']))
errors=[abs(x-y)/abs(y) for x,y in zip(b[0],expected)]
assert max(errors)<2e-8,(nstr,errors)
if nstr in ('0','1'):
anchors=[.75,.375] if nstr=='0' else [(15-math.pi**2)/(2*math.pi**2),(105-10*math.pi**2)/(6*math.pi**2)]
assert max(abs(x-y) for x,y in zip(b[0],anchors))<1e-10
assert not all(abs(x-(y+.001))/abs(y)<2e-8 for x,y in zip(b[0],anchors))
out['controls'][nstr]={'analytic_passed':True,'displaced_value_rejected':True}
out['models'][nstr]={'coarser':a,'finer':b,'relative_errors':errors}
out['worst_relative_error']=max(max(v['relative_errors']) for v in out['models'].values())
print(json.dumps(out,indent=2))
if __name__=='__main__':main(sys.argv[1] if len(sys.argv)>1 else __file__.replace('radial_check.py','numerical-evidence.json'))