# Numerical Soultion 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)$$

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 the cases. I am a beginner. I am unable to figure out where I have gone wrong . Please help me.

The equations I pass to my code are as follows

For matter

       void fcn(double t, double *y, double *f){

double p,hub;
p=0;
hub=y/y;
f = y;
f = -(y*y+3*p*y)/6.0;
f = -3*hub*y-3*hub*p;

}


           void fcn(double t, double *y, double *f){
double p,hub;
hub=y/y;
p= (2.0/3.0)*y;
f = y;
f = -(y*y+3*p*y)/6.0;
f = -3*hub*y-3*hub*p;
}


The plot I have obtained for $$a$$ vs $$t$$ is $a$ vs $$t$$ plot" />

I am just trying to see what is the type of functional dependence we get , not putting the initial conditions as per cosmology. Can this be because of wrong initial conditions?

One more thing is that variation of $$\rho$$ with time should be same in both 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" />

My code was just working excellent to solve other coupled differential equations.

         #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=1.0;
y=12.00;
y=30.00000567845;
ddopri5(fcn, y);

}

void fcn(double t, double *y, double *f){
double p,hub;
p=0;
hub=y/y;
f = y;
f = -(y*y+3*p*y)/6.0;
f = -3*hub*y-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] -n
(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,y,y);
if(iteration > 30000) break;
}

fclose(fpw);
return 0;
}


I also tried with Scipy . I am sharing my Scipy code below

      from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

def odes(x,t):
#assign each ode to a vector
A=x
B=x
C=x

#define each ODE
dBdt=-(1/6)*A*C
dCdt=-(3*C*B)/A

#define initialcondition
x0=[1,5,13]

#declare a time vector
t=np.linspace(0,30,1000)
x=odeint(odes,x0,t)

A=x[:,0]
B=x[:,1]
C=x[:,2]

#plot the results
plt.semilogy(t,A)
plt.semilogy(t,B)
plt.semilogy(t,C)
plt.show()

• What if the problem is your integrator, which you don't talk much about, and not your equations? Are you open to using Python for this problem? For such simple equations C is not a user-friendly approach. Aug 6 at 20:28
• @Richard I am using a C code for Dormand Prince method . When I was unable to get the desired result I tried with Python Scipy.integrator also. There also I was getting the same results.
– Dori
Aug 6 at 22:10
• if you are able to post your entire Scipy code instead of the C code, I think you'll find that you get better help faster because it will allow folks to reproduce your problem. Aug 6 at 22:16
• The indentation of the code is painful Aug 6 at 22:35
• If you want to use C, why not use the GSL library which has many well-tested ODE algorithms: gnu.org/software/gsl/doc/html/…
– vibe
Aug 7 at 6:31