ODE solver (45) giving two solutions (i.e. two plots) for an ODE! How do I plot only the one I need?

This is giving me two solutions (see the code and attached plot) but only one is relevent. Kindly let me know how can I choose and plot only the relevent one?!

Here is the function:

function [DY_Dt] = ScFld(N, Y)
%%constants
Mp=1.22*10^19; %GeV
G=1/Mp^2; %Gravitational constant units Gev^-2 
p= pi;
L = 0.001;
V0= 2.6975*10^-47; %GeV^4
k=(8*p*G)^0.5;

%variables 
phi = Y(1);
p = Y(2);
rho_m = Y(3);
rho_r=Y(4);

V=V0*exp(-L*k*phi);
H2= (k^2*(V+rho_m+rho_r))/(3- (k^2/2)*p^2);
dH_H2= -(k^2/2)*(p^2+((rho_m+(4/3)*rho_r)/H2));

DY_Dt = [p ; -p*(3+dH_H2)+(1/H2)*(k*L*V0*exp(-L*k*phi)) ; -3*rho_m ; -4*rho_r];   
end

And Here the code for solving and plotting:

O_m0=0.3111;
O_r0=10^-4;
H0=67660/(3*10^(22)*1.5192*10^24);
Mp=1.22*10^19; 
G=1/Mp^2;
p= pi;
L = 0.001;
V0= 2.6975*10^-47;
k=(8*p*G)^0.5;
rho_c=(3*H0^2)/k^2; 
rho_m0= O_m0*rho_c;
rho_r0=O_r0*rho_c;

Ic1=0.001 ;   
Ic2=0;  
Ic3= rho_m0; 
Ic4= rho_r0; 
Ic = [Ic1 Ic2 Ic3 Ic4];

%Independent variable range
N_f=-log(2001);
N_i=-log(1);
N_=linspace(N_i, N_f, 300);

options = odeset('RelTol',10^(-20),'AbsTol',10^(-20));
 % Solving differential Equations 
 [N, Y] = ode45(@ScFld, N_, Ic, options);
 
 V=V0*exp(-L*k.*Y(:,1));
 H2= (k^2*(V+Y(:,3)+Y(:,4)))/(3-((k^2/2).*Y(:,2).^2));
rho_phi=((H2.*Y(:,2).^2)/2)+V;

%plotting
figure(1);
 hold on;
 plot(N, log(rho_phi), 'b') 
 box on;
 grid on;

[Jan]

Hint: 2.6975*10^-47 is a multiplication and an expensive power operation. 2.6975e-47 is a cheap constant.

It is meaningless to reduce the tolerances of the integration to such a tiny value.

[M Sinha]

May I know how will that help if I change the way of writing the number?

Answer 1

[Star Strider]

Use element-wise division

H2= (k^2*(V+Y(:,3)+Y(:,4)))./(3-((k^2/2).*Y(:,2).^2))

↑ ← HERE

then ‘H2’ returns a column vector (rather than a matrix) and only onee solution.

% And Here the code for solving and plotting:
O_m0=0.3111;
O_r0=10^-4;
H0=67660/(3*10^(22)*1.5192*10^24);
Mp=1.22*10^19; 
G=1/Mp^2;
p= pi;
L = 0.001;
V0= 2.6975*10^-47;
k=(8*p*G)^0.5;
rho_c=(3*H0^2)/k^2; 
rho_m0= O_m0*rho_c;
rho_r0=O_r0*rho_c;
Ic1=0.001 ;   
Ic2=0;  
Ic3= rho_m0; 
Ic4= rho_r0; 
Ic = [Ic1 Ic2 Ic3 Ic4];
%Independent variable range
N_f=-log(2001);
N_i=-log(1);
N_=linspace(N_i, N_f, 300);
options = odeset('RelTol',10^(-20),'AbsTol',10^(-20));
 % Solving differential Equations 
 [N, Y] = ode45(@ScFld, N_, Ic, options)

Warning: RelTol has been increased to 2.22045e-14.

N = 300×1

0 -0.0254 -0.0508 -0.0763 -0.1017 -0.1271 -0.1525 -0.1780 -0.2034 -0.2288

Y = 300×4

1.0e+19 * 0.0000 0 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000

 
 V=V0*exp(-L*k.*Y(:,1));
 H2= (k^2*(V+Y(:,3)+Y(:,4)))./(3-((k^2/2).*Y(:,2).^2));

H2 = 300×1

1.0e+-71 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

rho_phi=((H2.*Y(:,2).^2)/2)+V;

rho_phi = 300×1

1.0e+-33 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

%plotting
figure(1);
 hold on;
 plot(N, log(rho_phi), 'b') 
 box on;
 grid on;

function [DY_Dt] = ScFld(N, Y)
%%constants
Mp=1.22*10^19; %GeV
G=1/Mp^2; %Gravitational constant units Gev^-2 
p= pi;
L = 0.001;
V0= 2.6975*10^-47; %GeV^4
k=(8*p*G)^0.5;
%variables 
phi = Y(1);
p = Y(2);
rho_m = Y(3);
rho_r=Y(4);
V=V0*exp(-L*k*phi);
H2= (k^2*(V+rho_m+rho_r))/(3- (k^2/2)*p^2);
dH_H2= -(k^2/2)*(p^2+((rho_m+(4/3)*rho_r)/H2));
DY_Dt = [p ; -p*(3+dH_H2)+(1/H2)*(k*L*V0*exp(-L*k*phi)) ; -3*rho_m ; -4*rho_r];   
end

.

[M Sinha]

Thank you so much Mr. Star Strider! I completely overlooked the part where I should have put '.' before '/'. It worked!

[Star Strider]

As always, my pleasure!

No worries! That is the most frequent problem I see here with respect to vectorising operations. You are definitely not the first, and I am certain will not be the last to encounter that problem.

.

[M Sinha]

Haha, I suppose so.