MATLABオイラー法

Question

オイラー法を使ってフリードマン方程式を近似しようとしていますが、私のコードはうまくいかないようです。私はapでもslopeでも値を得ていません。

N=12; 
    rho0=1; 
    h=0.001; 
    j=N*h; 
    t=0:h:j; 
    a0=1; 
    G=6.67*10^-11; 

    rhop=-3*(sqrt((8*pi*G)/3)); 
    rho=zeros(size(N+1)); 
    rho(1)=rho0; 
    a=zeros(size(N+1)); 
    ap=zeros(size(N+1)); 
    slope=zeros(size(N+1)); 
    a(1)=1; 


    for i= 1:(N); 
    slope(i)=rhop*((rho(i))^(3/2)); 
    ap(i)=sqrt((8*pi*G)/3)*((rho(i))^(1/2))*a(i); 
    rho(i+1)=rho(i)+slope(i)*h; 
    a(i+1)=a(i)+h*ap(i); 

    end 

    plot(t,a) 
    plot(t,rho) 

私は間違っていますか？

Answer 1

あなたの数学をチェックしましたか？あなたのコードを見ると、あなたのループは間違ってコーディングされているようです（免責事項：私は天体物理学者ではありません）。

また、あなたが必要としているのはこのcode snippetです。これは、あなたが過度に単純化しようとしている可能性があると私に信じさせる、あなたがやっていることよりもはるかに複雑です。

function friedmann 

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

%%%% Time arrays 
t_begin(1)=13.7; 
t_begin(2)=13.7; 
t_final(1)=0; 
t_final(2)=40; 

%%%% Solving in two parts : backwards and forwards 
for i=1:2 

%%%% Initial Omega parameters  
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; 

%%%% Initial scale factor and derivated scale factor 
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 for solver 
options = odeset('Events',@events,'RelTol',10^(-10),'AbsTol',10^(-10)); 

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

figure(1); 
hold on; 
%%%% Plot solutions 
plot(t1,y1(:,1),'b',t2,y2(:,1),'r',t3,y3(:,1),'k',t4,y4(:,1),'g'); 
box on; 
grid on; 
xlabel('Cosmic Time in Gyr'); 
ylabel('Relative size of the universe'); 
ylim([ 0 4]); 
%%%% Set label for each curve 
posX = 5; 
posX1 = posX; 
posX2 = posX; 
posX3 = posX; 
posX4 = posX; 
posY1 = 3.4; 
posY2 = 3.0; 
posY3 = 2.6; 
posY4 = 2.2; 
strMax = ['\Omega_{m} = 0.3 , \Omega_{\Lambda} = 0']; 
text(posX1, posY1, strMax, 'FontSize', 15, 'Color', 'blue', 'VerticalAlignment', 'top'); 
strMax = ['\Omega_{m} = 0.3 , \Omega_{\Lambda} = 0.7']; 
text(posX2, posY2, strMax, 'FontSize', 15, 'Color', 'red', 'VerticalAlignment', 'top'); 
strMax = ['\Omega_{m} = 5 , \Omega_{\Lambda} = 0']; 
text(posX3, posY3, strMax, 'FontSize', 15, 'Color', 'black', 'VerticalAlignment', 'top'); 
strMax = ['\Omega_{m} = 1 , \Omega_{\Lambda} = 0']; 
text(posX4, posY4, strMax, 'FontSize', 15, 'Color', 'green', 'VerticalAlignment', 'top'); 
title('Relative size of the universe vs Cosmic Time'); 

end 
end 

%%%% Differential system function 
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 

%%%% Integration stopping function 
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