Matlab - Sistema de 3 ecuaciones diferenciales

E = 0.1

rhoa = 1030 % kg/m3

h = 1

syms w(t) b(t) rho(t)

ec1 = 'Dw = (b*(rho*w^2+h(rhoa-rho)))/(w*(b*rho+2*E*h*rhoa))'

ec2 = 'Db = (w^2*(b+2*E*h))*(b*rho+2*E*h*rhoa)/(b*(rho*w^2+h*(rhoa-rho)))'

ec3 = 'Drho = (b*rho+2*E*h*rhoa)/(b+2*E*h)'

[w,b,rho] = dsolve(ec1,ec2,ec3)

%C1 = 'w(0) = 0.1' % m/s

%C2 = 'b(0) = 1' % m

%C3 = 'rho(0) = 100' % kg/m3

CI = 'w(0) = 0.1, b(0) = 1, rho(0) = 100'

[w,b,rho] = dsolve(ec1,ec2,ec3,CI)

function ode45caballero13

y0=[0.1 1 2];

global E rhoa h

E = 0.1;

rhoa = 1030; % kg/m3

h = 1;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);

 [t,y]=ode45(@ode45caballero,[0 0.3],[y0(1) y0(2) y0(3)],options  );

figure(gcf)

plot(t,y(:,1), t,y(:,2),t,y(:,3));

axis([0 0.3,0 10])

legend('w(t)','b(t)','rho(t)')

grid on

[w,b,rho]=dsolve1

end

 function dy = ode45caballero(t,y)

 global E rhoa h

dy = zeros(3,1);    % a column vector

dy(1) = y(2)*(y(3)*y(1)^2+h*(rhoa-y(3)))/(y(1)*(y(2)*y(3)+2*E*h*rhoa));

dy(2) = -(y(2)^2*(y(2)+2*E*h))*(y(2)*y(3)+2*E*h*rhoa)/(y(2)*(y(3)*y(1)^2+h*(rhoa-y(3))));

dy(3) = (y(2)*y(3)+2*E*h*rhoa)/(y(2)+2*E*h);

 end

 function [w,b,rho]=dsolve1

E = 0.1;

rhoa = 1030; % kg/m3

h = 1;

syms t w(t) b(t) rho(t)

 ec1 = diff(w)==b*(rho*w^2+h*(rhoa-rho))/(w*(b*rho+2*E*h*rhoa))

  ec2 = diff(b) ==(w^2*(b+2*E*h))*(b*rho+2*E*h*rhoa)/(b*(rho*w^2+h*(rhoa-rho)))

 ec3 = diff(rho)== (b*rho+2*E*h*rhoa)/(b+2*E*h)

 [w,b,rho] = dsolve(ec1,ec2,ec3, w(0) == 0.1,  b(0) == 1,rho(0) == 100 )

 end

>> ode45caballero13

ec1(t) =

diff(w(t), t) == (b(t)*(rho(t)*w(t)^2 - rho(t) + 1030))/(w(t)*(b(t)*rho(t) + 206))

ec2(t) =

diff(b(t), t) == (w(t)^2*(b(t) + 1/5)*(b(t)*rho(t) + 206))/(b(t)*(rho(t)*w(t)^2 - rho(t) + 1030))

ec3(t) =

diff(rho(t), t) == (b(t)*rho(t) + 206)/(b(t) + 1/5)

Warning: Explicit solution could not be found.

> In dsolve (line 201)

  In ode45caballero13>dsolve1 (line 33)

  In ode45caballero13 (line 14)

w =

[ empty sym ]

b =

     []

rho =

     []

w =

[ empty sym ]

b =

     []

rho =

     []