pb résolution numérique ( supprimé à la demande du posteur )

[fabio123]

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 :

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.
-----

[physeb]

Bonjour,

tu n'ai pas obliger d'avoir un facteur d"échelle qui duminue avec k=1.

Tu as remplit ton univers avec uniquement de la matière et cela en faible quantité (la bonne quantité certe et justement elle est faible devant ). C'est alors normal d'obtenir un facteur d'échelle qui croit.

Je ne sais pas si tu veux mettre toutes les quantités d'energies mais là il faut que tu ajoute d'autres '' avec leur équation d'état associée (Eq 4).

Pour finir, tu peux relier la valeur de k à la somme des en comparant avec . Mais en fait il y a une difficulté que je n'ai jamais calculé par moi même: tu utilises k=-1,0,1
En fait k peut prendre n'importe quelle valeur sauf si tu renormalise tes équations par la valeur absolue de k. Tout ça pour te dire que je ne sais pas trop ce que ça va donner à la fin.

Bon courage!

[fabio123]

Bonjour,

j'ai continué à chercher la recontractation future du facteur d'échelle avec k=0 car d'après ce que vous m'avez dit, je ne suis pas obligé de prendre le paramètre de courbure égal à 1 (k=1). J'ai donc essayé de jouer sur le paramètre et . Les résultats sont visibles sur la figure en pièce jointe ("scale_factor.png"). J'ai pris les valeurs actuelles pour obtenir la ligne bleue ( et ) avec k=0.
sont liés par la relation :

.

la ligne rouge est obtenue avec toujours avec k=0.
Enfin, la ligne verte pour avec k=0 (.)
On pouvait s'attendre dans ce dernier cas à ce qu'il y ait un big crunch mais il n'apparait pas pour cette ligne verte.
J'ai essayé aussi avec k=1 sans succès. Est ce que le choix de "k" se fait en fonction de la condition ? c'est à dire on prend k=-1 ou k=0 pour et k=1 pour .

Voici le code matlab :

function friedmann

al=9.461*10^(15);
H0=71000/(3*10^(22));


t_begin=380000*(3600*24*365);

t_today=13.7*10^(9)*(3600*24*365);

init=[ 0.001; H0*10^(3/2)];
t_final=50*10^(9)*(3600*24*365);

disp(t_begin);
disp(t_today);

[t1,y1]=ode23tb(@syseq1,[t_begin t_final],init);
[t2,y2]=ode45(@syseq2,[t_begin t_final],init);
[t3,y3]=ode45(@syseq3,[t_begin t_final],init);

figure(1);
t1=t1/(3600*24*365);
t2=t2/(3600*24*365);
t3=t3/(3600*24*365);
plot(t1,y1(:,1),'b',t2,y2(:,1),'r',t3,y3(:,1),'g'); 
end




function dydt1 = syseq1(t,y) 
H0=71000/(3*10^(22));
p=0;
k=0;
cv=3*10^(8);
G=6.673*10^(-11);
rho_crit=3*H0^(2)/(8*pi*G);
rho_0=0.3*rho_crit;
Lambda_red=(k*cv^(2)/(H0^(2)))+1-(rho_0/rho_crit)
Lambda=Lambda_red*3*H0^(2)/cv^(2);

 dydt1=zeros(2,1);
 dydt1=[ y(2) ;                        
       (-(4*pi*G)/(3*cv^(2))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda*cv^(2)/3)*(y(2))*((8*pi*G*(rho_0/y(1)^(3))+Lambda*cv^(2)/3-k*cv^(2)/(y(1)^(2))))^(-1/2);          
       ];
 
   
end


function dydt2 = syseq2(t,y) 
H0=71000/(3*10^(22));
p=0;
k=0;
cv=3*10^(8);
G=6.673*10^(-11);
rho_crit=3*H0^(2)/(8*pi*G);
rho_0=rho_crit;
Lambda_red=(k*cv^(2)/(H0^(2)))+1-(rho_0/rho_crit);
Lambda=Lambda_red*3*H0^(2)/cv^(2);

 dydt2=zeros(2,1);
 dydt2=[ y(2) ;                        
       (-(4*pi*G)/(3*cv^(2))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda*cv^(2)/3)*(y(2))*((8*pi*G*(rho_0/y(1)^(3))+Lambda*cv^(2)/3-k*cv^(2)/(y(1)^(2))))^(-1/2);          
       ];
 
   
end

function dydt3 = syseq3(t,y) 
H0=71000/(3*10^(22));
p=0;
k=0;
cv=3*10^(8);
G=6.673*10^(-11);
rho_crit=3*H0^(2)/(8*pi*G);
rho_0=2*rho_crit;
Lambda_red=((k*cv^(2)/(H0^(2)))+1-(rho_0/rho_crit))
Lambda=Lambda_red*3*H0^(2)/cv^(2);

 dydt3=zeros(2,1);
 dydt3=[ y(2) ;                        
       (-(4*pi*G)/(3*cv^(2))*((rho_0/y(1)^(3))*cv^(2)+3*p)+Lambda*cv^(2)/3)*(y(2))*((8*pi*G*(rho_0/y(1)^(3))+Lambda*cv^(2)/3-k*cv^(2)/(y(1)^(2))))^(-1/2);          
       ];
 
   
end

Je ne prends en compte que la matière car celle-ci est prédominante aux âges considérés ici.

Merci pour votre aide.