Solucin de Ecuaciones Diferenciales Ordinarias

Ing. Martn Rivera Toledo

1. Problemas con valores iniciales

Las ecuaciones diferenciales ordinarias no slo es una de las partes ms bellas de las
matemticas, sino que tambin es una herramienta esencial para modelar muchas situaciones fsicas: el
flujo de un reactivo en un reactor qumico de tipo tubular, el cambio de temperatura en una barra inmersa
en un fluido a menor temperatura, la respuesta dinmica y de control de un proceso, ... etc.
La importancia de las ecuaciones diferenciales ordinarias en matemticas aplicadas se debe al
hecho que la mayor parte de las leyes cientficas se expresan ms fcilmente en trminos de razones de
cambio, por ejemplo, el enfriamiento de una aleta:

dT
5/4
= 0.27(T 60)
dt
con su condicin inicial t = 0

T = T0

Muchos sistemas fsicos presentan

ordenes de magnitud muy diferentes entre sus
valores propios, esto provoca un comportamiento
conocido como rigidez (stiff), por ejemplo:
Cdigo ejemplo 1

y
y' = [A] y con y = 1
y2
y la condicin inicial

t =0

2
500.5 499.5
y = y [A] =

1
499.5 500.5

FIGURA 1

cuya solucin est dada por

1. 5
0.5 1000t
y = e t +
dado que los valores propios son 1 = 1000 y 2 = 1
e
1. 5
0. 5

Estructura para resolver ODES con MATLAB

Tipos de rutinas

Problemas

1. La profundidad h del agua al vaciarse un tanque cilndrico vertical por un agujero en su fondo est
descrita por

dh
A
= 0
dt
Aw

2 gh con g = 32 .2

ft
s2

en donde Aw y Ao son las reas transversales del tanque y del agujero, respectivamente. Resuelva la
ecuacin con una profundidad inicial del agua de 20 ft, Aw = 50 ft2 y Ao = 0.5 ft2. En qu momento queda
vaco el tanque?
2. Resuelva la ecuacin diferencial ordinaria de segundo orden

d2x
+ 4x = 5sen( 2t ) + 3 cos( 2t ) con las condiciones iniciales x( 0) = 1 y x' (0) = 1
dt 2
3. Resuelva el sistema de ecuaciones diferenciales ordinarias siguiente

dx
+ y=t
dt

4x +

dy
= 0 con las condiciones iniciales x(0 ) = 1 y y(0 ) = 2
dt

Ejemplo 2
El modelo matemtico de un CSTR descrito por Shacham, Brauner y Cutlip (1994) Chem. Eng. Educ.
est dado por las expresiones siguientes

V dCa/dt = Fo (Cao - Ca) V k Ca

Cp V dT/dt = Cp Fo (To - T) V k Ca U A (T - Tj)
j Cpj Vj dTj/dt = j Cpj Fj (Tjo - Tj) + U A (T - Tj)
donde

PARMETRO
Fo
Cao
V
Fj
R
Vj

k = exp [ -E / (R T)]

VALOR
40 ft3/h
0.50 mol/ft3
48 ft3
49.9 ft3/h
1.99 BTU/mol R
3.85 ft3
7.08*1010 1/h
30000 BTU/mol

A
Tjo
To
Cp
Cpj

150 BTU/ hft 2 F

-30000 BTU/mol
250 ft2
530 R
530 R
0.75 BTU/lb R
1.0 BTU/lb R
50 lb/ft3
62.3 lb/ft 3

Reporte
a. Los tres estados estacionarios
b. La simulacin dinmica usando
los valores nominales de CA, T y
Tj en el estado estacionario
c. La simulacin dinmica usando
los valores nominales de CA, T
5 % y Tj en el estado estacionario
d. La simulacin dinmica usando
un valor de Fj 5% los valores

nominales de CA, T y Tj en el
estado estacionario
e. La simulacin dinmica usando
un valor de F0 5% los valores
nominales de CA, T y Tj en el
estado estacionario

Solucin
a. Estados estacionarios: Cdigo
i. Ca = 0.0591 mol/ft3 T = 651.0596 R

Tj = 641.7920 R

ii. Ca = 0.4717 mol/ft3 T = 537.1114 R Tj = 536.5670 R

iii. Ca = 0.2437 mol/ft3 T = 600.1893R

Tj = 594.8160 R

La simulacin dinmica, partiendo de las condiciones iniciales es la siguiente Cdigo

b. La simulacin dinmica usando los valores nominales de CA, T y Tj en el estado estacionario

Ca = 0.2437 mol/ft3 T = 600.1893R Tj = 594.8160 R es

d. La simulacin dinmica usando un valor de Fj 5% los valores nominales de CA, T y Tj en el estado

estacionario

e. La simulacin dinmica usando un valor de F 5% los valores nominales de CA, T y Tj en el

estado estacionario

Problemas de Tarea [Para entregar el jueves 1 4 de Julio de 2005]

Problema 1
Para la ecuacin diferencial siguiente

dy
= t2 t
dt

Se tiene la solucin analtica dada por

y(1) = 0

t3 t2 1
y= +
3 2 6

Trace las grficas para

las respuestas de las soluciones analtica y numrica,

Problema 2
El modelo matemtico de un circuito elctrico se da por la ecuacin

0.5

Con

Q=0

d 2Q
dt 2

+6

dQ
+ 50Q = 24 sen 10 t
dt

dQ
e I=
en t = 0
dt

Exprselo como un sistema de dos ecuaciones de primer orden y trace su respuesta

2. Solucin de sistemas de ecuaciones diferenciales y algebraicas [DAES]

La forma ms general de un DAE es F (t , y, y ') = 0 en donde F puede ser singular, un caso
especial de una representacin semi-implcita de un DAE o un ODE con restricciones es

x ' = f (t , x, z )

0 = g (t , x, z)
Para la rutinas de solucin de ecuaciones diferenciales ordinarias y algebraicas en MATLAB, se considera
una representacin con la forma siguiente

[M (t , y )] d y =
dt

f (t , y )

en donde

[M (t , y )]

es la matriz de masa,

la cual tiene valores de 1 en la diagonal principal si le corresponde una ecuacin diferencial, de lo contrario
tendr un valor de cero para una ecuacin algebraica, por ejemplo, para el sistema DAE siguiente:
Ejemplo 1
Trace la respuesta del sistema algebro-diferencial siguiente

dy2
=0
dt

dy1
= y3
dt

0 = y1 y 2 + y 3 (1 y 2 ) t

con las condiciones iniciales t = 0 y1 = 1, y2 = 0 & y3 = 1 tome en cuenta que su correspondiente

matriz de masa ser

1 0 0

[M (t, y)] = 0 1 0
0 0 0

Solucin del problema DAE 1 Cdigo

Problemas de tarea
Problema 3. Reactor de flujo continuo no isotrmico
La reaccin orgnica elemental irreversible en fase lquida

A + B > C
se llevar a cabo en un reactor de flujo continuo, las especies A & B se alimentan al reactor a una temperatura
de 50C y 20C, respectivamante, con un flujo de 10 mol/min cada una

Reporte:
A. Las concentraciones de los componentes A, B, C y la temperatura a la salid a para CSTR con un
volumen de 100 dm3 para los casos:
a. Con enfriamiento a travs de una interfase metlica con un coeficiente UA = 1000 cal / min K y una
temperatura ambiente de 50C.
b. Adiabtico

Datos:
k = 0.01 dm3 / mol min a 300 K con
Vo = 50 dm3

Ea = 10 kcal/mol

Q = 1 dm3 /min

CAo = CBo = 10 mol / dm3

CpA = 15 cal/ mol K

CpB = 15 cal/ mol K

CpC = 30 cal/ mol K

HAo = -20 kcal/ mol

HBo = -15 kcal/ mol

HCo = -41 kcal/ mol

Solucin
1. Representacin esquemtica y modelo

rA = kC AC B

rB = r A

rC = rA

A. CSTR [Cdigo MATLAB]

dC A C A0 C A
=
+ rA
dt

dC B C B 0 C B
=
+ rB
dt

dC C C C 0 C c
=
+ rC
dt

con

V
Q

UA(Ta T ) 1
+ [C A0 Cp A (TA0 T ) + CB 0 CpB (TB 0 T )] + ( Hrxn)( rA )
dT
V

=
dt
C ACp A + C B CpB + CC CpC
CI

t =0

C A = C A0

C B = CB 0

CC = CC 0

T = T0

Problema 4. Flujo por gravedad [Cdigo]

En cunto tiempo se vaciar el tanque?

Formule sus modelos y resulvalos a travs del uso de las herramientas
matemticas de MATLAB
Caractersticas del sistema:
L1: 34 m
L2: 100 m
D: 1 acero comercial ced. 40
ZA: 37 m
Kentrada: 0.5
Ksalida: 1.0
L/D codo:30
: 1.5*10-4
Dint: 1.049 in
Dext: 1.315 in

De la ecuacin diferencial de balance de masa se tiene

dm
=m
& A m& B y tomando en cuenta que se pueden expresar los flujos en masa
dt
2
D A Z
en funcin de los flujos volumtricos m& A = Q A , m& B = Q B & m =
4
dZ 4 (Q A QB )
=
sta se simplifica como
y su condicin inicial t = 0 Z = ZA
2
dt
D A

Considerado una condicin quasi-estacionaria para el balance de energa mecnica entre los puntos A &
B la ecuacin de Bernoulli ser

( Z Z B ) g = vB 2
gc

1
+
f
+ K ent + K sal con la ecuacin de Colebrook-White

2 gc
D

1
2.51
= 2.0 log
+
f
3.7D Re f
Re =

DvB

para calcular el factor de friccin y el nmero de Reynolds como

y finalmente el flujo volumtrico

Q=

D 2v B
4

Bibliografa

1.

Ascher, U.M. & Petzold, L.R. (1998), Computer Methods for Ordinary Differential Equations
and Diferencia-Algebraic Equations, Ed. SIAM

2.

Gerald,C.F. & Wheatley, P.O. (1999), Anlisis Numrico con Aplicaciones 6a Edicin
Prentice-Hall Mxico

3.

Zill, D. (1997) ecuaciones Diferenciales con Aplicaciones de Modelado, 6 Edicin ,

International Thomson Editor

4.

Seinfield, J. H. , Lapidus, L. & Hwang, M. (1970) Review of Numerical Integration

Techniques for Stiff Ordinary Differential Equations, Ind. Eng. Chem. Fundam. 9(2) pp.266275

5.

Flores, T. (2002), [Link]

Cdigo ejemplo 1 regresar

% inicio del archivo odestiff.m
% archivo para resolver un par de ecuaciones
diferenciales
% ordinarias de primer orden con valores iniciales
% y' = [A] y con ci t=0 y=y0
% donde A=[-500.5 499.5;
%
499.5 -500.5];
%
y = [y1 y2]' & y0 = [2 1]'
%
% inicio de la funcion odestiff
%-------------------------------function odestiff
%------------------------------ -clc; clear all; format compact;
tspan = linspace(0,0.1);
y0=[2 1];

y2a = 1.5*exp(-tspan)-0.5*exp(-1000*tspan);

[t,y]=ode45(@edostiff1,tspan,y0)
y1=y(:,1);y2=y(:,2);

dydt = A*y;
% fin de la funcion de edostiff
% fin del archivo odestiff.m

y1a = 1.5*exp(-tspan)+0.5*exp(-1000*tspan);

figure(1),plot(t,y1,t,y2,'--'),xlabel('t'),ylabel('y')
legend('y_{1}','y_{2}')
figure(2),plot(tspan,y1a,tspan,y2a,'-'),xlabel('t'),ylabel('y')
legend('y_{1}','y_{2}')
% fin de la funcion odestiff
% inicio de la funcion de edostiff1
%-------------------------------function dydt = edostiff1(t,y)
%-------------------------------A=[-500.5 499.5;
499.5 -500.5];

Cdigo ejemplo 2 regresar

% inicio del archivo pshacham.m
clc; clear all; format compact; format long;
global F0 Ca0 V Fj R Vj alfa E U lambda A Tj0 T0 Cp
Cpj rho rhoj
% parametros del modelo
F0=40;
% ft3/h
Ca0=0.50*1.1;
%mol/ft3
V=48;
%ft3
Fj=49.9;
%ft3/h
R=1.99;
%BTU/mol-R
Vj=3.85;
%ft3
alfa=7.08e10; % 1/h
E=3e4;
%BTU/mol
U=150;
%BTU/ h-ft2-R
lambda=-3e4; %BTU/mol
A=250;
%ft2
Tj0=530;
%R
T0=530;
%R
Cp=0.75;
%BTU/bl-R
Cpj=1.0;
%BTU/bl-R
rho=50;
%lb/ft3
rhoj=62.3;
%lb/ft3
x0=[Ca0,T0,Tj0];
[x,fval,flag]=fsolve('eshacham',x0)
Ca=x(1)
T=x(2)

Tj=x(3)
% fin del archivo pshacham.m
% inicio del archivo eshacham.m
function [f]=eshacham(x)
global F0 Ca0 V Fj R Vj alfa E U lambda A Tj0 T0 Cp
Cpj rho rhoj
Ca=x(1); T=x(2); Tj=x(3); % reasignacion de variables
RT=R*T;
k=alfa*exp( -E/RT);
VkCa=V*k*Ca;
rhoCpF0=rho*Cp*F0;
lambdaVkCa=lambda*VkCa;
UA=U*A;
rhoCpV=rho*Cp*V;
rhojCpjVj=rhoj*Cpj*Vj;
rhojCpjFj=rhoj*Cpj*Fj;
% expresiones del balance de masa y energia
f(1)=(F0*(Ca0-Ca)-VkCa)/V;
f(2)=(rhoCpF0*(T0-T)-lambdaVkCa-UA*(T Tj))/rhoCpV;
f(3)=(rhojCpjFj*(Tj0-Tj)+UA*(T-Tj))/rhojCpjVj;
% fin del archivo eshacham.m

Cdigo para respuesta dinmica del CSTR regresar

% inicio del archivo pshacham.m
clc; clear all; format compact; format long;
global F0 Ca0 V Fj R Vj alfa E U lambda A Tj0 T0
Cp Cpj rho rhoj
% parametros del modelo
F0=40;
% ft3/h
Ca0=0.50;
%mol/ft3
V=48;
%ft3
Fj=49.9;
%ft3/h
R=1.99;
%BTU/mol-R
Vj=3.85;
%ft3
alfa=7.08e10; % 1/h
E=3e4;
%BTU/mol
U=150;
%BTU/ h-ft2-R
lambda= -3e4; %BTU/mol
A=250;
%ft2
Tj0=530;
%R
T0=530;
%R
Cp=0.75;
%BTU/bl-R

Cpj=1.0;
%BTU/bl-R
rho=50;
%lb/ft3
rhoj=62.3;
%lb/ft3
Ca10 = 0.2437; % mol/ft3
T10 = 600.1893;%R
Tj10 = 594.8160;%R
%tr=linspace(0,10); % h tiempo de operacion
tr=[0,10];
x0=[Ca10,T10,Tj10];
[t,x]=ode15s('shacham',tr,x0);
Ca=x(:,1); T=x(:,2); Tj=x(:,3);
figure(1)
plot(t,Ca,'o-')
grid
xlabel('Tiempo [ h ]')
ylabel('Ca [ lbmol/ft3 ]')
title('Concentracion vs Tiempo')
figure(2)
plot(t,T,'o-',t,Tj,'*-')

grid
legend('T','Tj')
xlabel('Tiempo [ h ]')
ylabel('T [ R ]')
title('Temperatura vs Tiempo')
% fin del archivo pshacham.m

VkCa=V*k*Ca;
rhoCpF0=rho*Cp*F0;
lambdaVkCa=lambda*VkCa;
UA=U*A;
rhoCpV=rho*Cp*V;
rhojCpjVj=rhoj*Cpj*Vj;
rhojCpjFj=rhoj*Cpj*Fj;

% inicio del archivo shacham.m

function [dfdt]=shacham(t,x)
global F0 Ca0 V Fj R Vj alfa E U lambda A Tj0 T0
Cp Cpj rho rhoj
Ca=x(1); T=x(2); Tj=x(3); % reasignacion de
variables

% expresiones del balance de masa y energia

dCadt=(F0*(Ca0-Ca)-VkCa)/V;
dTdt=(rhoCpF0*(T0-T)-lambdaVkCa-UA*(T Tj))/rhoCpV;
dTjdt=(rhojCpjFj*(Tj0-Tj)+UA*(T-Tj))/rhojCpjVj;

RT=R*T;
k=alfa*exp(-E/RT);

dfdt=[dCadt dTdt dTjdt]';

% fin del archivo shacham.m

Cdigo para resolver el DAE regresar

% inicio del archivo daesex1.m
% ---------------------function daesex1
% ---------------------clc; clear all; format compact;
% declaracion de la matriz de masa
mass_matrix = [1 0 0; 0 1 0 ; 0 0 0];
% condiciones iniciales
y0 = [1; 0 ; 1];
% intervalo de integracion
tspan = [0 10];
% declaracion de los parametros para odeset
options = odeset ('Mass',mass_matrix,'RelTol',1e4,'AbsTol',...
[1e-6 1e-10 1e-6]);
% integracion del conjunto DAE
[t,y] = ode15s (@daes1,tspan,y0,options);
% trazo de los resultados
figure(1), plot(t,y), xlabel('tiempo'),
ylabel('y'), legend('y_{1}','y_{2}','y_{3}',3)
title ('solucion del DAE')
% ----------------------------function dydt_g = daes1 (t,y)
% ----------------------------%reasignacion de variables
y1 = y(1); y2 = y(2); y3 = y(3);
% declaracion de las ecuaciones diferenciales
% y algebraicas
dy1dt = y3;
dy2dt = 0;
g = y1*y2 - y3*(1-y2)-t;
dydt_g = [dy1dt dy2dt g]';
% fin del archivo daesex1.m

Cdigo MATLAB caso CSTR dinmico regresar

% inicio del archivo pcstrni.m
% problema 9.5b fogler pag 575, 3a ed
% septiembre 2004 martin rivera toledo uia
clc; clear all; format compact;
global Cpa Cpb Cpc ha0 hb0 hc0 Tr R kTr Tk Ea...
Ca0 Cb0 Cc0 Ta0 Tb0 Q UA Ta V tau
%-----------------------------------------------------% datos del problema
%-----------------------------------------------------Cpa = 15; % cal/ mol K
Cpb = 15; % cal/ mol K
Cpc = 30; % cal/ mol K
ha0 = -20e3;% % cal/ mol
hb0 = -15e3;% cal/ mol
hc0 = -41e3;% cal/ mol
Tr = 273; % K
R = 1.982; % cal / mol K
Ea = 10e3; % cal/ mol
Fa0 = 10; % mol/min
Fb0 = 10; % mol/min
Fc0 = 0; % mol/min
Ta0 = 20 + 273.15; % K
Tb0 = 50 + 273.15; % K
Q = 1; % dm3/min
%UA = 1000; % cal / min K
UA=0;
V = 100; % dm3
tau =V/Q; % min tiempo de residencia
Ta = 50 + 273.15; % temperatura ambiente
Tk=300;%K temperatura del coeficiente cinetico

kTr=0.01;% dm3 / mol min

%-----------------------------------------------------------% condiciones iniciales del sistema
%-----------------------------------------------------------Ca0 = Fa0/Q;% mol
Cb0 = Fb0/Q;% mol
Cc0 = Fc0/Q;% mol
T0 = 25 + 273.15; % K
x0 = [Ca0 Cb0 Cc0 T0];
%-----------------------------------------------------% calculo de la operacion en regimen dinamico
%-----------------------------------------------------t0 = linspace(0,10);
[t xr] = ode15s('cstrni',t0,x0);
Ca=xr(:,1); Cb=xr(:,2); Cc=xr(:,3);
T=xr(:,4); % reasignacion de variables
%-----------------------------------------------------% respuesta grafica
%-----------------------------------------------------figure(1), plot(t,Ca,t,Cb,t,Cc)
xlabel ('tiempo [min]'),ylabel ('Ci [ mol i /dm^{3}]')
legend ('C_{A}','C_{B}','C_{C}')
figure(2), plot(t,T)
xlabel ('tiempo [min]'),ylabel ('Temperatura [ K ]')
% fin del archivo pcstrni.m
% inicio del archivo cstrni.m
% problema 9.5b fogler pag 575, 3a ed
% septiembre 2004 martin rivera toledo uia
%-----------------------------------------------------function [dfdt]=cstrni(t,x)
%-----------------------------------------------------global Cpa Cpb Cpc ha0 hb0 hc0 Tr R kTr Tk Ea...
Ca0 Cb0 Cc0 Ta0 Tb0 Q UA Ta V tau
Ca=x(1 ); Cb=x(2); Cc=x(3); T=x(4); % reasignacion de
variables
%
Qr = UA*(Ta-T); % remocion de calor
CiCpi = Ca*Cpa+Cb*Cpb+Cc*Cpc;
dhr0 = hc0-(ha0+hb0);% calor de reaccion a Tr
dcp = Cpc -(Cpa+Cpb); % diferencia de capacidad
calorifica
dhr = dhr0+dcp*(T -Tr); % calor de reaccion a la
temperatura de reaccion
k = kTr*exp(Ea/R*(1/Tk-1/T)); % coeficiente cinetico
a la temperatura de reaccion
%-----------------------------------------------------% expresiones de rapidez
%-----------------------------------------------------ra = -k*Ca*Cb; % rapidez de reaccion para
componente A
rb = ra; % rapidez de reaccion para componente B

rc = -ra; % rapidez de reaccion para componente C

%-----------------------------------------------------% ecuaciones del modelo
%-----------------------------------------------------dCadt = (Ca0-Ca)/tau + ra;
% balance de masa
para especie A
dCbdt = (Cb0-Cb)/tau + rb;
% balance de masa
para especie B
dCcdt = (Cc0-Cc)/tau + rc;
% balance de masa
para especie C
dTdt = ( Qr/V + (Ca0*Cpa*(Ta0 -T) + Cb0*Cpb*(Tb0
-T))/tau...
+ (-dhr)*(-ra) ) / CiCpi; % balance de energia
%dTdt = 0;
%-----------------------------------------------------% vector resultante de ecuaciones diferenciales
%-----------------------------------------------------dfdt=[dCadt dCbdt dCcdt dTdt]';
% fin del archivo cstrni.m

Codigo flujo por gravedad [regresar]

% inicio del archivo Eetanque.m
%---------------------function eetanque
%---------------------% programa principal para la respuesta dinamica del
modelo
% de flujo por gravedad
% elaborado por martin rivera
%
clc;clear all; format compact;
%declaracion de variables globales A B C ... etc
global Da Db L h Qa Qb E visc ro g gc
%---------------------------------------------------%
declaracion de los parametros del sistema
%---------------------------------------------------Da = 1.0; %m diametro del tanque
Db = 0.05; %m diametro del tanque
L = 7;
%m longitud del tubo de descarga
E = 4.921e-4;%m rugosidad absoluta
visc = 0.001;%kg /m s viscosidad del fluido
ro = 998.2; %kg/m3 densidad del fluido
g = 9.81; %m/s2 aceleracion de la gravedad
gc = 9.81; %kgm/s2/kgf factor conversion gravitatorio
a absoluto
h = 3;
%m altura del nivel del liquido del tanque
Qa = 0; %m3/s flujo volumetrico de alimentacion
%----------------------------------------------------

%
definicion de l a matriz de masa
%---------------------------------------------------mass_matrix=eye(3);
mass_matrix(2,2)=0;mass_matrix(3,3)=0;
abstolv(1:3)=1e-6; % vector de tolerancia absoluta para
las variables
options = odeset ('Mass',mass_matrix,'RelTol',1e4,'AbsT ol',abstolv);
%---------------------------------------------------%
condiciones iniciales
%---------------------------------------------------h0=3; Q0=0.01; f0=0.03;
x0=[h0, Q0, f0]; % condiciones iniciales
tr=linspace(0,250); % definicion del
integracion

tiempo

de

%---------------------------------------------------[t,x]=ode15s(@ode_tanque,tr,x0,options);
%---------------------------------------------------h=x(:,1); Q=x(:,2); f=x(:,3);
figure(1), plot(t,h)
xlabel ('tiempo [ s ]'), ylabel ('altura nivel [m]')
title('nivel tanque vs tiempo')
figure(2), plot(t,Q)
xlabel ('tiempo [ s ]'), ylabel ('Q [m^{3}/s]')
title('flujo volumetrico vs tiempo')
figure(3), plot(t,f)
xlabel ('tiempo [ s ]'), ylabel ('factor de friccion ')
title('factor de friccion vs tiempo ')
%fin del archivo dtanque.m

% inicio del archivo ode_tanque.m

%---------------------------------------------------function [dfdt_g]= ode_tanque(t,x)
%---------------------------------------------------%declaracion de variables globales A B C ... etc
global Da Db L h Qa E visc ro g gc
% reasignacion de variables originales var(i)=x(i);
h = x(1); % altura del nivel del liquido
Qb=x(2); % flujo columetrico
fb=x(3); % factor de friccion
% calculo del numero de Reynolds
Re=4*Qb*ro/(pi*Db*visc);
Er=E/Db; % rugosidad relativa
% expresiones del modelo f(i) = expresion;
dhdt = 4*(Qa - Qb)/(pi*Da^2); % balance de masa
g1
=
(h
+
L)*g/gc
8*Qb^2*(1+fb*L/Db)/(pi^2*gc*Db^4);% balance
energia mecanica

de

%g2
=
1/fb^0.5+2*log10(E/(3.7*Db)+2.51*pi*Db*visc/(4*Qb*r
o*fb^0.5)); % colebrook-white
g2=10(3.7*fb^0.5*Re/(Er*fb^0.5*Re+3.7*2.51))^(2.0*fb^0.5);
dfdt_g = [dhdt g1 g2]';
% fin del archivo ode_tanque.m