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;
}
For radiation
void fcn(double t, double *y, double *f){
/* a=y[0],adot=y[1],rho=y[2] */
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 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.
