Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab-Friedmann Equations, help?

  1. Aug 14, 2008 #1
    I have been asked to create a programme on Matlab to solve the Friedmann equations. However being quite computer illiterate I am not really sure where to begin. I have been reading some books on Matlab, but still not sure how to go about this. Could anyone offer any tips?
     
  2. jcsd
  3. Aug 17, 2011 #2
    here is a numerical solving program with matlab "ode" solvers,
    most useful for ΛCDM models where there is not analytic solution :


    Code (Text):

    function friedmann

    H0=71000/(3*10^(22))*(3600*24*365*10^(9));

    t_begin(1)=13.7;
    t_begin(2)=13.7;
    t_final(1)=0;
    t_final(2)=40;
    for j=1:2
       
    Omega1_m=0.3;
    Omega1_red=0;
    Omega1_k=1-Omega1_m-Omega1_red;
    Omega2_m=0.3;
    Omega2_red=0.7;
    Omega2_k=1-Omega2_m-Omega2_red;
    Omega3_m=5;
    Omega3_red=0;
    Omega3_k=1-Omega3_m-Omega3_red;
    Omega4_m=1;
    Omega4_red=0;
    Omega4_k=1-Omega4_m-Omega4_red;

    a1_0=1;
    a1_prim_0=H0*(Omega1_m/a1_0+Omega1_k)^(1/2);
    init1=[ a1_0 ; a1_prim_0 ];
    a2_0=1;
    a2_prim_0=H0*(Omega2_m/a2_0+Omega2_k+Omega2_red*a2_0)^(1/2);
    init2=[ a2_0 ; a2_prim_0 ];
    a3_0=1;
    a3_prim_0=H0*(Omega3_m/a3_0+Omega3_k)^(1/2);
    init3=[ a3_0 ; a3_prim_0 ];
    a4_0=1;
    a4_prim_0=H0*(Omega4_m/a4_0+Omega4_k)^(1/2);
    init4=[ a4_0 ; a4_prim_0 ];
    options = odeset('Events',@events,'RelTol',10^(-10),'AbsTol',10^(-10));


    [t1,y1]=ode45(@(t1,y1) syseq(t1,y1,Omega1_m,Omega1_red,Omega1_k),[t_begin(j) t_final(j)],init1,options);
    [t2,y2]=ode45(@(t2,y2) syseq(t2,y2,Omega2_m,Omega2_red,Omega2_k),[t_begin(j) t_final(j)],init2,options);
    [t3,y3]=ode45(@(t3,y3) syseq(t3,y3,Omega3_m,Omega3_red,Omega3_k),[t_begin(j) t_final(j)],init3,options);
    [t4,y4]=ode45(@(t4,y4) syseq(t4,y4,Omega4_m,Omega4_red,Omega4_k),[t_begin(j) t_final(j)],init4,options);

    figure(1);
    hold on;
    plot(t1,y1(:,1),'b',t2,y2(:,1),'r',t3,y3(:,1),'k',t4,y4(:,1),'g');

    xlabel('Time in Gyr');
    ylabel('Relative size of the universe');
    ylim([ 0 4]);

    end
    end

    function dydt = syseq(t,y,Omega_m,Omega_red,Omega_k)

    H0=71000/(3*10^(22))*(3600*24*365*10^(9));

    dydt=[ y(2) ;                        
     (-(H0^(2)/2)*(Omega_m/y(1)^(3)-2*Omega_red))*Omega_m*(1/(H0^(2))*y(2)^2-Omega_k-Omega_red*y(1)^(2))^(-1);                
           ];
     
       
    end


    function [value,isterminal,direction] = events(t,y)
    % Locate the time when height passes through zero in a
    % decreasing direction and stop integration.
    value = y(1);     % Detect height = 0
    isterminal = 1;   % Stop the integration
    direction = -1;   % Negative direction only
    end

     
    Regards
     
    Last edited by a moderator: Oct 6, 2011
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?