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;
                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;
                     p= (2.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;

The plot I have obtained for $a$ vs $t$ is <span class=$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

<span class=$\rho$ vs $t$ plot" />

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

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


         void fcn(double t, double *y, double *f){
                double p,hub;
                 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]+ 

                         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] + 
                    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] - 
                   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] + 
                    err[i] =  h*((71.0/57600.0)*k1[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 = sqrt(errabs);
               if( errabs < eps){
                    t+= h;
                    for(i = 0; i < N; i++){
              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;      

             fprintf(fpw ,"%.25lf %.25lf %.25lf %.25lf \n", 
             if(iteration > 30000) break;

             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

      #define each ODE

      return [dAdt,dBdt,dCdt]

      #define initialcondition

      #declare a time vector


      #plot the results
  • 2
    $\begingroup$ 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. $\endgroup$
    – Richard
    Aug 6, 2021 at 20:28
  • $\begingroup$ @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. $\endgroup$
    – Dori
    Aug 6, 2021 at 22:10
  • $\begingroup$ 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. $\endgroup$
    – Richard
    Aug 6, 2021 at 22:16
  • $\begingroup$ The indentation of the code is painful $\endgroup$ Aug 6, 2021 at 22:35
  • $\begingroup$ 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/… $\endgroup$
    – vibe
    Aug 7, 2021 at 6:31


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.