Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

Supplementary material for the paper "Cosmology from the two-dimensional renormalization group acting as the Ricci flow" (one pdf file showing calculations and three SageMath notebooks performing numerical calculations).

Project: cosmology
Views: 637
Kernel: SageMath (stable)

8.5 Numerical calculation of the invariant cSc_S

This SageMath notebook uses the optional TIDES package which provides an arbitrary precision ode solver.

To add the TIDES package to SageMath (on a unix machine) run in a shell     $ sage -i tides as instructed at http://doc.sagemath.org/html/en/reference/misc/sage/misc/package.html

Definitions

decimal_precision = 120 # number of decimal digits in floating point numbers # binary_precision = ceil(decimal_precision*log(10,2)) R = RealField(binary_precision); RealNumber = R print "decimal precision =", decimal_precision, "binary precision =",binary_precision # solver_step_size = 0.001 # step_size in T for the trajectories returned by the ode solver solver_tolrel = 1e-80 # the relative tolerance for the solver solver_tolabs = 1e-80 # the absolute tolerance for the solver #solver_tolrel = 1e-60 # the relative tolerance for the solver #solver_tolabs = 1e-60 # the absolute tolerance for the solver # sqrt3=R(sqrt(3)) nu = (sqrt3-1)/2 lambdaplus = 2+sqrt3 lambdaminus = 2-sqrt3 betaplus=1+sqrt3/3 betaminus=1-sqrt3/3 # ode_function(T,hminus,hplus)=[-1-hminus^2+lambdaminus*(hplus*hminus+1),\ -1-hplus^2+lambdaplus*(hplus*hminus+1)] def calc_traj(point_i,T_f): [T_i,hminus_i,hplus_i] = point_i # integrate backwards or forwards? if T_f < T_i: signed_solver_step_size = -solver_step_size else: signed_solver_step_size = solver_step_size traj=desolve_tides_mpfr( ode_function, [hminus_i,hplus_i], T_i, T_f, signed_solver_step_size, solver_tolrel, solver_tolabs, decimal_precision) # prune runaway points traj_pruned = [point for point in traj if point[1] != NaN and point[2] != NaN] hh_list = [point[1:] for point in traj_pruned] return(traj_pruned,hh_list,len(traj)-len(traj_pruned)) # # large T expansion # Ry.<y> = R[] f_poly = 0.0*y v_poly = 0.0*y f_tilde_coeffs = [1.0] v_tilde_coeffs = [] def calc_fv_polys(): global f_poly, v_poly f_tilde_coeffs = [1.0] v_tilde_coeffs = [] for m_tilde in range(N_tilde_coeffs+1): ff_sum=sum(f_tilde_coeffs[mp]*f_tilde_coeffs[m_tilde-mp]/binomial(m_tilde,mp) for mp in range(m_tilde+1)) v_tilde_coeffs.append(-2*f_tilde_coeffs[m_tilde] + ff_sum/(2*m_tilde+1)) fv_sum=sum(f_tilde_coeffs[mp]*v_tilde_coeffs[m_tilde-mp]/binomial(m_tilde,mp) for mp in range(m_tilde+1)) f_tilde_coeffs.append(-(m_tilde+1/2)*f_tilde_coeffs[m_tilde]/(m_tilde+1)+(ff_sum+fv_sum)/(m_tilde+1)) f_poly = 0.0*y v_poly = 0.0*y for mm in range(N_tilde_coeffs): f_poly += f_tilde_coeffs[mm]*factorial(mm) *y^(2*mm+1) v_poly += v_tilde_coeffs[mm]*factorial(mm) *y^(2*mm+1) # # find point on S using large T expansion # def s_point(T_val): y_val = 1.0/T_val f_T_val = f_poly(y=y_val) v_T_val = -T_val + v_poly(y=y_val) hminus_val = f_T_val + betaminus * v_T_val hplus_val = f_T_val + betaplus * v_T_val return([T_val,hminus_val,hplus_val]) # return({'f_T':f_T_val,'v_T':v_T_val,'hminus':hminus_val,'hplus':hplus_val}) # # solve ode starting from point on S at large T # def separatrix(T,T_f): point_i=s_point(T) (s_traj,s_hh_list,number_pruned) = calc_traj(point_i,T_f) print number_pruned, " points pruned", len(s_traj), " points" return(s_traj,s_hh_list)
decimal precision = 120 binary precision = 399

Integrate the ode from large TiT_i to somewhat larger TfT_f to see the instability of the separatrix SS

Calculate the large TT expansion of SS to order 100

N_tilde_coeffs =100 calc_fv_polys()

Choose Ti=20T_i=20 (arbitrarily). The function separatrix(TiT_i,TfT_f) uses the large TT expansion at T=TiT=T_i to calculate the initial condition, then integrates the ode to TfT_f.

# T_i = 20, T_f = 25.4 (sep_traj,sep_hh_list)=separatrix(20,25.25) list_plot(sep_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
0 points pruned 5251 points
Image in a Jupyter notebook

Points are pruned from the trajectory when the ode runs away. The value of TfT_f is chosen by trial and error to capture the instability.

try higher order 120 in the large TT expansion

N_tilde_coeffs =120 calc_fv_polys()
(sep_traj,sep_hh_list)=separatrix(20,25.75) list_plot(sep_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
0 points pruned 5751 points
Image in a Jupyter notebook

The instability shows up very slightly later in TT, suggesting that the higher order expansion might be slightly more accurate at placing the initial point on the separatrix.

Now integrate backwards in TT to find the separatrix at large h−h_{-}. Choose the earliest TfT_f such that the solver does not run away (i.e., no points are pruned).

(sep_traj_back,sep_back_hh_list) = separatrix(20,-0.6220) list_plot(sep_back_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
1 points pruned 20622 points
Image in a Jupyter notebook
(sep_traj_back,sep_back_hh_list) = separatrix(20,-0.621) list_plot(sep_back_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
0 points pruned 20622 points
Image in a Jupyter notebook

Get the earliest point on the trajectory.

back_point = sep_traj_back[-1]; back_point
[-0.62099999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999998, 11979.738443413583082144782481069482787725314129673935381110229234360188556339628695461441469072904293690294961463237222, -0.000083474276564673047602103434215648075436230636068420002091783368642468601744304745701367323280449836251819145030300315359]

Integrate foward from that earliest point until just before runaway.

(sep_traj_forward,sep_forward_hh_list,number_pruned) =\ calc_traj(back_point,12.82) list_plot(sep_forward_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
Image in a Jupyter notebook

Plot the points on the trajectory with h−<10.0h_{-} < 10.0, displaying the instability.

pruned=[point for point in sep_forward_hh_list if point[0]< 10.0] list_plot(pruned,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
Image in a Jupyter notebook

Perturb the initial point very slightly to display the instability in both directions.

T_back = back_point[0] hminus_back = back_point[1] hplus_back = back_point[2]
point_mod=[T_back,hminus_back,hplus_back - 1.31211060e-88 ] (sep_traj_mod,sep_mod_hh_list,number_pruned) =\ calc_traj(point_mod,13.4) pruned=[point for point in sep_mod_hh_list if point[0]< 10.0] list_plot(pruned,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
Image in a Jupyter notebook
point_mod_2=[T_back,hminus_back,hplus_back - 1.31211061e-88] (sep_traj_mod_2,sep_mod_2_hh_list,number_pruned_2) =\ calc_traj(point_mod_2,13.4) pruned=[point for point in sep_mod_2_hh_list if point[0]< 10.0] list_plot(pruned,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
Image in a Jupyter notebook

The difference in the initial h+h_{+} values is O(10−96)O(10^{-96}). This is much smaller than the tolerance parameters of the ode solver.

Define a function to calculate the double expansion eS(x,y)e^S(x,y) to O(xNyN)O(x^N y^N) and a function to use eS(x,y)e^S(x,y) to calculate cc.

Rxy.<xx,yy> = R[] def calc_expS(N): j_size = N+1 k_size = N+1 zklist = [0.0 for k in range(k_size)] Sjk = [[0.0 for k in range(k_size)] for j in range(j_size)] Sjk[0][0] = 1.0 if j_size > 1: Sjk[1][0] = 1.0 for k in range(1,k_size): for j in range(1,j_size): Sjk[j][k] = ((2-j-(1+nu)*k)*Sjk[j-1][k] \ + (lambdaminus*j+(1-nu)*(k-1))*Sjk[j][k-1]) \ /(j+(2+nu)*k) expS = 0.0 for j in range(j_size): for k in range(k_size): expS+= Sjk[j][k]*xx^j*yy^k return(expS) def calc_c(point,expS): hminus_val = point[0] hplus_val = point[1] x_val = 1.0/hminus_val^2 y_val = (hplus_val*hminus_val+1.0)*x_val expS_val = expS(xx=x_val,yy=y_val) c_val = expS_val*y_val*(x_val^(-2.0-nu)) return(c_val)

For one of the two trajectories bracketing the separatrix, calculate cc at the earliest point for several values of NN.

cvals={} for N in range(10): expS = calc_expS(N) cvals[N] = calc_c(point_mod[1:],expS)
cvals
{0: 0.28445788665466583381471805147876193750348369531743315490387388017381802005436613900703825930517289597954613816749432817, 1: 0.28445788863675554418877978914259637544544616653289475791447457114119762673701745040770625196759798492977030445020689204, 2: 0.28445788863675554418877978914259637544053557077520960232297491369849384079444308496582665196457907472588138752667272458, 3: 0.28445788863675554418877978914259637544053557079029671698615895337795517756130003494036406043263051121782074883638954280, 4: 0.28445788863675554418877978914259637544053557079029671693057357826578273290504777101635027895311635213621996382967154836, 5: 0.28445788863675554418877978914259637544053557079029671693057357849535473125067261485784686777059223949814030804694403173, 6: 0.28445788863675554418877978914259637544053557079029671693057357849535473022464798120993511983711278571617736358154813162, 7: 0.28445788863675554418877978914259637544053557079029671693057357849535473022464798606926241540811483036562653646955135020, 8: 0.28445788863675554418877978914259637544053557079029671693057357849535473022464798606926239134773811471851460068024544137, 9: 0.28445788863675554418877978914259637544053557079029671693057357849535473022464798606926239134773823811940895094995129935}

Check that the expansion improves with NN for small NN.

NN=0 logdevs = [[N,log(abs(cvals[N]-cvals[9]),10)] for N in cvals if N >=NN and N < 9] list_plot(logdevs)
Image in a Jupyter notebook

Evaluate the deviation in cc along the early part of the separatrix, i.e., small TT, large h−h_{-} to check the accuracy of the expansion of es(x,y)e^s(x,y).

hminus_min = 30 expS = calc_expS(9) Tc_list=[[point[0],calc_c(point[1:],expS)-cvals[9]] for point in sep_traj_mod if point[1]>hminus_min] list_plot(Tc_list,axes_labels=[r"$T$",r"$\Delta c$"])
Image in a Jupyter notebook

check numerical precision

nu*(1+nu)-0.5
0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
lambdaplus*lambdaminus-1
1.5490367659397273091272122533145407390280876768195163461951051765377665392390743859450728937108058387632230170030217870e-120

check large T expansion of separatrix

N_tilde_coeffs = 10 calc_fv_polys() test_ode = y*(0.5*y^2*f_poly.derivative(y)-f_poly*v_poly-f_poly^2-1)+f_poly test_ode
5.9109714707410239459590521448978416357858977397635351659452142136058193323590324881683916140189492021067141957989139353e10*y^39 - 1.2371689146222568392662362600849397377664170507382153409043180489598123902677680903280798725102520982241349403752203940e10*y^37 + 2.0918528906674824373038519446503936596269825857957406873269514622352053620857438680864347427268892151141347116033286798e9*y^35 - 3.4192511731384801012290318940091859706049161508458237596236575505396044057852560930905279090897117784985942695370246068e8*y^33 + 5.7801148446586037557755313033172040022643511902463346288739693965031963479205623926760244535217772375883483574423198136e7*y^31 - 1.0582397791770538048771035564399633674067094158777555783081640148113461085373497061819114286825017841070985975671496604e7*y^29 + 2.2082062131082423016774882521181040588803694045790884526378724057795686447147031380724858115814498367519576002969360312e6*y^27 - 570053.89406978727818174708041765976050801091928654064619561929597054758230339573988144527471369319219570431166180577056*y^25 + 218938.04630259429567173498079193690771655290796280525505548848219249713180618876230454194973335963065188088530758932253*y^23 - 2.7048765579483540656291006431062453471417056851230520698307813153751529101669157691566655152089325758793545908391846767e6*y^21 + 2.0303534698525193786192196446644348374588950797608764652848482569915813583114395791459259432446274289837316728462007167e-115*y^19 + 2.5379418373156492232740245558305435468236188497010955816060603212394766978892994739324074290557842862296645910577508958e-116*y^17 - 6.3448545932891230581850613895763588670590471242527389540151508030986917447232486848310185726394607155741614776443772396e-117*y^15 + 3.9655341208057019113656633684852242919119044526579618462594692519366823404520304280193866078996629472338509235277357748e-118*y^13

check accuracy of large T expansion

N_tilde_coeffs =100 calc_fv_polys() [sep_traj_list_1,sep_hh_list_1] = separatrix(100,101) #list_plot(sep_hh_list_1) diff_list = [[point[0],point[1]-s_point(point[0])[1]] for point in sep_traj_list_1] list_plot(diff_list)
0 points pruned 1001 points
Image in a Jupyter notebook

Try the large TT expansion for larger TiT_i

(sep_traj,sep_hh_list)=separatrix(30,34.25) list_plot(sep_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
0 points pruned 4251 points
Image in a Jupyter notebook
(sep_traj,sep_hh_list)=separatrix(50,52.65) list_plot(sep_hh_list,axes_labels=[r"$h_{-}$",r"$h_{+}$"])
0 points pruned 2651 points
Image in a Jupyter notebook