Integrating two ODE'S

1 view (last 30 days)
Ian DSouza
Ian DSouza on 28 Dec 2017
Commented: Ian DSouza on 30 Dec 2017
I'm trying to integrate these two equations: d/da(rho*a^3) = -3*w*rho*a^2; ((da/dt)/a)^2 + k/a^2 = 8*pi*G/3*rho
where a==a(t), rho==rho(a),w=w(rho).
For now, I am taking 'w' to b a constant. I found out d(rho)/dt = d(rho)/da * da/dt to get 'ode1' and expressed rhoas a function of time 't'. 'ode2' is just the second equation.
Code:
---------------------
syms a(t) rho(t) w(rho)
G=1;
w=1/3;
k=0;
ode1 = diff(rho) == -3*rho*(w+1)*sqrt(8*pi*G/3*rho);
ode2 = (diff(a)/a)^2 + k/a^2 == 8*pi*G/3*rho;
odes = [ode1;ode2];
Sol = dsolve(odes)
aSol(t) = Sol.a;
rhoSol(t) = Sol.rho
---------------------------- Output:
Sol =
struct with fields:
a: [1×1 sym]
rho: [1×1 sym]
aSol(t) =
C5
rhoSol(t) =
0
K>>
This is of course not the solution. Can someone please tell me where I'm going wrong. Thanks.
I am open to numerically integrating these as well
  2 Comments
David Goodmanson
David Goodmanson on 30 Dec 2017
Hello Ian, For w = constant, since rho = rho(a) only depends on time implicitly through a(t), it appears that the solution to the first equation is simply
rho = C*a^(-3(w+1))
and then a(t) can be calculated from the second equation.
Ian DSouza
Ian DSouza on 30 Dec 2017
Thanks David. Yes,that is the solution for w=constant. I decided to keep the ODE because I wanted to eventually vary 'w' as a function of rho. But I got it to work using ODE45 now.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!