Bonjour,

j'essaie de résoudre les équations de Friedmann dans le cas général de manière numérique sous Matlab. Pour cela, j'utilise les solveurs "ode" de Matlab. Les équations s'écrivent :

eq n°1 :


eq n°2 :


On peut donc reécrire eq n° 1 sous la forme :

eq n°3 :

avec l'expression de la densité :

eq n°4 :


avec


On peut mettre cette équation sous la forme d'un système différentiel :

y1'=R'=y2
y2'=R"

Puis, on utilise le solveur ode45 qui est similaire à Runge-Kutta. Voici le source Matlab :

Code:
function friedmann

H0=71000/(3*10^(22)); % H0=71km/s/Mpc
init=[ 0.001; 1000*H0]; % R_begin mille fois plus petit qu aujourd'hui, R'_begin= mille fois l'expansion actuelle
t_begin=380000*(3600*24*365); % 380.000 années
t_final=20*10^(9)*(3600*24*365); %20 Milliard d'années

[t,y]=ode45(@syseq,[t_begin t_final],init);
figure(1);  
plot(t,y(:,1));  %plot du facteur d'échelle en fonction du temps
end

function dydt = syseq(t,y) 
H0=71000/(3*10^(22));
Lambda_red=0.7; % constante cosmologique réduite
p=0;
k=1; % paramètre de courbure, sphérique ici 
cv=3*10^(8);
G=6.673*10^(-11);
rho_0=5*10^(-27); % densité massique de la matière baryonique
Lambda=Lambda_red*3*H0^(2);


 dydt=zeros(2,1);
 dydt=[ y(2) ;                        
       (-(4*pi*G)/(3*cv^(2))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda/3)*(y(2))*((8*pi*G*(rho_0/y(1)^(3))+Lambda/3-k/(y(1)^(2))))^(-1/2);          
       ];
   
   
end
Le problème se situe au niveau des résultats. Pour k=0 ou k=-1, j'obtiens une expansion infinie mais pour k=1 (univers de type sphérique), je ne retrouve pas la recontraction future, le facteur d'échelle ne cesse de croître. La méthode numérique mise en place est-elle correcte ? d'où peut venir ce résultat ?

merci par avance de votre aide.