# Numerical solution to background equations of cosmology

I am trying to solve the background equations of cosmology numerically using Runge-Kutta Dormand prince method with simplified assumption $$8\pi G=1$$ and $$c=1$$. The equations are $$\ddot a = -\frac{1}{6}(\rho + 3p)a$$ and $$\dot{\rho}= -3 \frac{\dot a}{a}(\rho + p)$$ with equation of state for matter $$(p=0)$$ and radiation $$(p=\frac{1}{3}\rho)$$.

I should get $$a\propto t^{2/3}$$ for matter and $$a\propto t^{1/2}$$ for radiation. But I am getting same variation for scale factor $$a$$ in both cases. I am an amateur student. I am unable to figure out what wrong I did. Please help me. I am sharing the equation I have given to my code.

For matter,

     void fcn(double t, double *y, double *f){
double p,hub;
p=0;
hub=y[1]/y[0];
f[0] = y[1];
f[1] = -(y[0]*y[2]+3*p*y[0])/6.0;
f[2] = -3*hub*y[2]-3*hub*p;

}


       void fcn(double t, double *y, double *f){
double p,hub;
hub=y[1]/y[0];
p= (1.0/3.0)*y[2];
f[0] = y[1];
f[1] = -(y[0]*y[2]+3*p*y[0])/6.0;
f[2] = -3*hub*y[2]-3*hub*p;
}


I am also attaching the plot I have obtained for $$a$$ vs $$t$$ $a$ vs $$t$$" />

One was just trying to see what is the type of functional dependence we get, not putting the adjact initial conditions according to cosmology. Can this be because of wrong initial conditions? One more thing is that variation of $$\rho$$ with time should be same in both the cases ($$\rho \propto t^{-2}$$) but I am getting different result for this also.

My plot for $$\rho$$ vs $$t$$ is as below

$\rho$ vs $$t$$ plot" />

The code I used to solve this is

        #include<stdio.h>
#include<math.h>
#include<stdbool.h>
#include  <stdlib.h>

#define N 3

double ddopri5(void fcn(double, double *, double *),
double *y);
double alpha;
void fcn(double t, double *y, double *f);
double eps;

int main(void){
double y[N];
//eps = 1.e-9;
printf("Enter epsilon:\n");
scanf("%lg", &eps);
y[0]=1.0;
y[1]=12.00;
y[2]=30.00000567845;
ddopri5(fcn, y);

}

void fcn(double t, double *y, double *f){
double p,hub;
p=0;
hub=y[1]/y[0];

f[0] = y[1];
f[1] = -(y[0]*y[2]+3*p*y[0])/6.0;
f[2] = -3*hub*y[2]-3*hub*p;

}

double ddopri5(void fcn(double, double *, double *),
double *y){
double t, h, a, b, tw, chi;
double w[N], k1[N], k2[N], k3[N], k4[N], k5[N],
k6[N],k7[N], err[N], dy[N];
int i;
double errabs;
int iteration;
FILE *fpw;
fpw=fopen("case1.dat", "w");

iteration = 0;
//eps = 1.e-9;
h = 1.0e-6;
a = 0.0;
b = 30;
t = a;
while(t < b -eps){
fcn(t, y, k1);
tw = t+ (1.0/5.0)*h;
for(i = 0; i < N; i++){
w[i] = y[i] + h*(1.0/5.0)*k1[i];
}
fcn(tw, w, k2);
tw = t+ (3.0/10.0)*h;
for(i = 0; i < N; i++){
w[i] = y[i] + h*((3.0/40.0)*k1[i] +
(9.0/40.0)*k2[i]);
}

fcn(tw, w, k3);
tw = t+ (4.0/5.0)*h;
for(i = 0; i < N; i++){
w[i] = y[i] + h*((44.0/45.0)*k1[i] -
(56.0/15.0)*k2[i] +
(32.0/9.0)*k3[i]);
}
fcn(tw, w, k4);
tw = t+ (8.0/9.0)*h;
for(i = 0; i < N; i++){
w[i] = y[i] + h*((19372.0/6561.0)*k1[i] -
(25360.0/2187.0)*k2[i] +
(64448.0/6561.0)*k3[i]
-(212.0/729.0)*k4[i]);
}
fcn(tw, w, k5);
tw = t + h;
for(i = 0; i < N; i++){
w[i] = y[i] + h*((9017.0/3168.0)*k1[i] -
(355.0/33.0)*k2[i] +
(46732.0/5247.0)*k3[i] +
(49.0/176.0)*k4[i] -
(5103.0/18656.0)*k5[i]) ;
}
fcn(tw, w, k6);
tw = t + h;
for(i = 0; i < N; i++){
w[i] = y[i] + h*((35.0/384.0)*k1[i] +
(500.0/1113.0)*k3[i] +
(125.0/192.0)*k4[i] -
(2187.0/6784.0)*k5[i] +
(11.0/84.0)*k6[i]);
}
fcn(tw, w, k7);
errabs = 0;
for(i = 0; i < N; i++){
dy[i] = h*((35.0/384.0)*k1[i] +
(500.0/1113.0)*k3[i] +
(125.0/192.0)*k4[i] -
(2187.0/6784.0)*k5[i] +
(11.0/84.0)*k6[i]);

err[i] =  h*((71.0/57600.0)*k1[i]  -
(71.0/16695.0)*k3[i]
+ (71.0/1920.0)*k4[i] -
(17253.0/339200.0)*k5[i] +
(22.0/525.0)*k6[i] -
(1.0/40.0)*k7[i]);
errabs+=err[i]*err[i];
}

errabs = sqrt(errabs);

if( errabs < eps){
t+= h;

for(i = 0; i < N; i++){
y[i]+=dy[i];
}
}
chi=errabs/eps;
chi = pow(chi, (1.0/6.0));
if(chi > 10)    chi = 10;
if(chi < 0.1)   chi = 0.1;
h*=  0.95/chi;
if( t + h > b ) h = b - t;

iteration++;
fprintf(fpw ,"%.25lf %.25lf %.25lf %.25lf \n",
t,y[0],y[1],y[2]);

if(iteration > 30000) break;
}

fclose(fpw);
return 0;
}


The initial conditions doesnot have and physical significance.

• maybe consider scicomp.stackexchange.com Aug 6 '21 at 18:43
• What solver are you using? Aug 6 '21 at 20:36
• @AndersSandberg , I am using a C code for Dormand Prince method which was doing excellent to solve other coupled differential equations.
– Dori
Aug 6 '21 at 21:52

A curious choice of equations. You may like to try instead this form:

$$\frac{H^{2}}{H_{0}^{2}} = \Omega_{0,m}a^{-3} + \Omega_{0,r}a^{-4} + \Omega_{0,k}a^{-2} + \Omega_{0,\Lambda}$$, making sure that the $$\Omega$$ factors sum up to 1.

Alternatively, try the following system:

$$H^{2} = \frac{\rho}{3}$$

$$\dot{\rho} = -3H\rho(1+w)$$

for $$p = w\rho$$

These will prevent the need to consider the second derivative.

You should confirm that the following holds:

$$d_{t}\ln{\rho}=\frac{\dot{\rho}}{\rho} = -3(1+w)\frac{\dot{a}}{a} = -3(1+w) d_{t}\ln(a)$$ so that

$$\rho \propto a^{-3(1+w)}$$

Lastly, I suggest doing a simple Euler method solver and starting your solution at a=1 (you can integrate towards a=0, towards the past, too), as there may be problems with your functions' interface with the solver, and there may be problems with starting at a=0.

• Thank you for your answer . I will definately try this. Actually, during a project work I was assigned to solve the pair of equations I mention in the question numerically.
– Dori
Aug 6 '21 at 21:55