Source: https://issuu.com/didaktica/docs/uso-de-matlab-en-la-solucio__n-de-a
Timestamp: 2018-02-22 09:30:44+00:00

Document:
Uso de Matlab en la solución de algunos sistemas dinámicos Héctor Escobosa by Didáktica - issuu
AUTOR: HÉCTOR ISAIAS ESCOBOSA IRIBE
PONENCIA: USO DE MATLAB EN LA SOLUCIÓN DE ALGUNOS SISTEMAS DINÁMICOS (INVESTIGACIÓN TEÓRICA EN DESARROLLO). Presentada en el CIIDET de la Ciudad de Querétaro, en el mes de Noviembre de 2010.
INSTITUCIÓN DE ADSCRIPCIÓN: INSTITUTO TECNOLÓGICO DE CULIACÁN.
E-MAIL: hescobosa@hotmail.com ÁREA
ENTORNOS VIRTUALES DE APRENDIZAJE: APLICACIÓN DE TECNOLOGIAS DE LA INFORMACIÓN EN LA EDUCACIÓN.
RESUMEN: Una de las aplicaciones más importantes de las Ecuaciones Diferenciales (Ordinarias y Parciales), se tiene en el estudio de Sistemas Dinámicos, los cuales representan situaciones ó problemas diversos que se presentan en la vida real. Dichos Sistemas Dinámicos pueden resolverse utilizando diferentes métodos analíticos, como los de Euler, Taylor y Cauchy, entre otros, para obtener Soluciones que satisfagan las condiciones iniciales y/ó de frontera, de las Ecuaciones Diferenciales que modelan dichos Sistemas. También pueden obtenerse Soluciones Numéricas para tales Sistemas, usando algoritmos como los de Euler y Runge-Kutta, entre otros. Este trabajo resulta más sencillo, utilizando las Tecnologías de la Información disponibles en la actualidad. Así, tenemos que usando programas (Software), tales como MATLAB, MATHCAD, MATHEMATICA, MAPLE, entre otros, podemos realizar estudios detallados sobre el comportamiento de las Soluciones de Sistemas Dinámicos de interés. En este trabajo, primero presento la forma en que podemos obtener las soluciones de Ecuaciones Diferenciales Ordinarias (EDO), mediante instrucciones ó comandos de MATLAB, para diferentes condiciones iniciales. En segundo término, presento el estudio de cinco Sistemas Dinámicos representativos de ese amplio, interesante e importante tema, por sus aplicaciones en Ciencias e Ingenierías, resueltos mediante el uso de MATLAB. La finalidad del presente trabajo, es contribuir en el estudio de los Sistemas Dinámicos que se modelan mediante EDO y ampliar su difusión en las Instituciones de Educación Superior, como Tecnológicos y Universidades.
Gráfica movimiento de Punto A a un Punto B, en Espiral R3
Punto A 0.5
-1 Eje x
1.-) Uso de Software MATLAB, en la Resolución de EDO. En esta parte, veremos diferentes instrucciones del Software Matlab para determinar soluciones analíticas y numéricas de Ecuaciones Diferenciales Ordinarias (EDO) de orden n ³ 1, mediante la aplicación de algunos principios ó métodos básicos, como el de Euler, el de Cauchy y los de Runge-Kutta (ode23 y ode45). Utilizaremos dichas instrucciones ó comandos de MATLAB, en diferentes modelos matemáticos que se utilizan con frecuencia en los cursos de Matemáticas IV , V y que se pueden usar en otras materias de Ingeniería, tales como: Circuitos Eléctricos, Vibraciones Mecánicas, Cinética Química, Operaciones Unitarias, Sistemas de Control, etc.
EDO y su resolución analítica con MATLAB. Las EDO no siempre son fáciles de resolver; incluso, en algunas ocasiones, es dificil obtener la Solución Analítica de los modelos matemáticos que representan. Es en estos casos, en los que el Software Matlab, nos proporciona una herramienta poderosa, que nos ayuda a encontrar Soluciones de dichos modelos EDO y obtener las gráficas de su comportamiento. La función dsolve determina las soluciones simbólicas ó analíticas de EDO. Cuando se trabaja con Ecuaciones Diferenciales, en MATLAB, se necesita una forma de incluir las diferenciales en una expresión. Por lo tanto, la sintáxis de dsolve es un poco diferente a otras funciones de Matlab. Las ecuaciones son especificadas usando la letra D para indicar diferenciación; D2, D3, etc, para indicar derivada de 2º y 3er orden respectivamente, y así sucesivamente. Cualquier letra seguida de D# es una variable dependiente. Por ejemplo, la ecuación diferencial: d2Y + Y
d t2 es representada, en Matlab, por la expresión simbólica: D2Y + Y = 0 La variable independiente en la EDO, se puede especificar libremente por el usuario; sin embargo, por default es la variable tiempo ( t ) en Matlab. Por ejemplo, la solución general de la EDO de primer orden siguiente: dy = 1 + y2
se puede resolver fácilmente con Matlab, mediante:
>> dsolve ('DY = 1 + Y^2 ')
% determina la solución general de la EDO
% la respuesta de Matlab tan( t + C1 )
% C1 es una constante de integración
Resolviendo la misma EDO para la condición inicial Y(0) = 1, se tendrá: 3
>>Y= dsolve( ' DY = 1 + Y^2 ' , ' Y(0) = 1 ' ) % almacena respuesta de Matlab en Y. Y=
Tan( t + (Pi / 4 ) )
Ahora, vamos a resolver una EDO de 2ยบ orden, para las condiciones iniciales dadas:
DY(0) = 0 = Cos(2 t ) - Y ;
>> Y = dsolve('D2Y = Cos(2*t ) - Y', 'DY(0) = 0', 'Y(0) = 1' ) % almacena la respuesta de Matlab en Y. Y= (-2/3)*Cos( t )^2 + 1/3 + (4/3)*Cos( t )
>> Y = simple(Y)
% instrucciรณn para que Matlab simplifique la expresiรณn, % un poco mรกs, si es posible.
Y= (-1/3)*Cos(2* t ) + (4/3)*Cos( t )
Veremos otro ejemplo, resolviendo la siguiente EDO: D2Y - 2 D Y - 3 Y = 0
>>Y = dsolve( ' D2Y - 2*DY - 3*Y = 0 ' ) Y=
% almacena la respuesta de Matlab en Y
C1* exp(- t ) + C2* exp(3*t )
Aplicando las condiciones iniciales: Y(0) = 0 e Y(1) = 1, tendremos: >>Y = dsolve( ' D2Y - 2*DY - 3*Y = 0 ', 'Y(0) = 0', 'Y(1) = 1' ) 1
1 * exp(- t ) -
Y= (exp(-1) - exp(3)
* exp(3*t ) exp(-1) - exp(3)
Para simplificar al mรกximo, usamos 4
>> Y = simple(Y) Y=
% simplifica al máximo la expresión para Y
- ( exp(- t ) - exp(3* t ) ) / ( exp(3) - exp(-1) )
Para mejorar la presentación, usamos la instrucción pretty: >> Pretty(Y) exp(-t ) - exp(3*t ) exp(3) - exp(-1)
Ahora, se construye la gráfica de la solución de à D2Y - 2 D Y - 3 Y = 0.
1/(exp(t) (1/exp(1) - exp(3))) - exp(3 t)/(1/exp(1) - exp(3)) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
>>ezplot(Y, [0 , 1] )
% instrucción de Matlab, para gráficar la solución
2.-) Aplicación de EDO en algunos Sistemas Dinámicos. A continuación, presento aplicaciones de EDO en cinco DinámicosàProblemas interesantes e importantes en la realidad.
APLICACIÓN #1-) En la elaboración de una gran cantidad de productos de limpieza, tanto de uso industrial como para el hogar, se utilizan diversos componentes, sobresaliendo entre ellos, los derivados clorados del benceno; es por ello que la reacción química para lograr obtener dichos derivados, es de gran importancia. La Reacción Química para la Cloración del Benceno se presenta a continuación: K1 C6H6 + CL2
C6H5CL + HCL K2
C6H5CL + CL2
C6H4CL + HCL K3
C6H4CL + CL2
C6H3CL + HCL
Se está considerando, para éste sistema dinámico, que la temperatura permanece constante durante todo el tiempo que dura la reacción. Aquí, se establece que: Q X Y Z W
Concentración del Cloro en Molgr. Concentración del –Benceno en Molgr. Concentración del Monocloro-Benceno " Concentración del Dicloro-Benceno " Concentración del Tricloro-Benceno "
Las Ecuaciones Diferenciales Ordinarias que modelan a éste Sistema Dinámico son: a)
dX/dt = - K1 X Q
dY/dt = (K1 X - K2 Y) Q
dZ/dt = (K2 Y - K3 Z) Q
dW/dt = K3 Z Q
Tomaremos como condiciones iniciales, t = 0 X(0) = 4 molgr. y Q =1 molgr (base de cálculo). De la ecuación (a), tenemos: x dX
dt Ln X – Ln 4 = - K1 t dt 4 X 0 realizando la integración, sustituyendo límites y simplificando obtenemos: t e) X(t) = 4 Î- K1 = - K1 X Q
= - K1
La ecuación(e), nos dará la Concentración del Benceno en cualquier instante de tiempo t>0. 6
Ahora, sustituyendo la ecuación (e) en la (b), tenemos: dY
= K1 ( 4 Î- K1 t ) - K2 Y
+ K2 Y = 4 K1 Î- K1
ecuación diferencial lineal, cuya solución es:
dt 4 K1
- K2 t
(Î- ( K2
- K1 ) t
(K2 - K1)dt
(K2 - K1) Integrando, sustituyendo la condición inicial Y(0) = 0, y simplificando, obtenemos: 4 K1 Y(t) =
[ Î- K1
- Î- K2
... ec(f)
(K2-K1) Esta ecuación( f ), nos dará la concentración del monocloro-benceno en cualquier instante de tiempo t > 0. Ahora, sustituyendo la ecuación ( f ) en la ( c), obtenemos la ecuación diferencial que nos dará la concentración del dicloro-benceno en cualquier instante de tiempo t > 0. dZ
4 K1 K2 [Î- K1
Î- K2 t
] - K3 Z
(K2 - K1)
la cual es también, una ecuación diferencial lineal de primer orden, cuya solución estará dada por: 4 K1K2
Î- K2
- K1 t
Z(t) = (K2 - K1) 4 K1K2 +
( K3 - K1) Î- K3
( K3-K2 )
...ec(g)
(K3-K1)(K3-K2)
La ecuación (g), nos dará la concentración del dicloro-benceno en cualquier instante t > 0. Procediendo de manera análoga, obtenemos la ecuación-solución para W(t), que nos dara la concentración del tricloro-benceno en cualquier instante de tiempo t > 0. Sustituyendo la ecuación (g) en la ecuación (d), integrando y simplificando tenemos que: 7
4K2 K3
- 4 K1 K3
W(t) = (K2-K1)(K1-K3) - 4 K1 K2
(K2-K1)(K2-K3) Î- K3
+ C(k)
-4 K2K3 (K3-K2) + 4 K1K3 (K3-K1) + 4 K1K2 (K2 -K1) en donde, C(k) = ( K2-K1 ) ( K3-K1 ) ( K3-K2 ) Así, tenemos resuelto analíticamente el modelo de cinética química que representa al Sistema Dinámico para la CLORACION del BENCENO. La gráfica típica para esta solución analítica, se presenta a continuación:
X(0) = 4 Molgr.
W(t) Z(t) Y(t)
Comportamiento de cada uno de los componentes que se tienen en la reacción química para la cloración del benceno.
A continuación, presento un Programa en MATLAB, para analizar el comportamiento de la Cinética de una Reacción Química, del tipo: k1 A
C (Reversible)
% Programa de MATLAB para Resolver el Sistema de EDO que modelan % la Cinetica de la Reaccion Quimica A
% y graficar los Perfiles de Concentración de cada componente clc, clear all 8
% function Xpunto=dx(t,X) % Programa de MATLAB para construir una funcion % global K % Xpunto(1)=-K(1)*X(1)+K(2)*X(2); % Xpunto(2)=K(1)*X(1)-K(2)*X(2)-K(3)*X(2)+K(4)*X(3); % Xpunto(3)=K(3)*X(2)-K(4)*X(3); % Xpunto=Xpunto'; % fin del programa - Funcion de Matlab echo on pause global K % Datos de Entrada fprintf('\ Introduce las Constantes de la Cinética de Reacción \') K(1)=input(' A --> B , K1= ' ) ; K(2)=input(' B --> A , K2= ' ); K(3)=input(' B --> C , K3= ' ) ; K(4)=input(' C --> B , K4= ' ); % fprintf('\ Introduce las Condiciones ó Concentaciones iniciales \') X0(1)=input(' Concentración inicial de A = ' ); X0(2)=input(' Concentración inicial de B = ' ); X0(3)=input(' Concentración inicial de C = ' ); tmax=input('\ Dame el tiempo máximo \ tmax= '); % tmax= el valor que se introduzca; K=[1,0,2,3] ; X0=[2 0 0] ; t=[0:0.01:tmax]; [t,X]=ode45('dx',[0,tmax],X0); plot(t,X),xlabel('Tiempo (seg)'),ylabel('Concentraciones') title('Grafica-Cinetica de Reaccion A->B->C'),Legend('XA ','XB','XC') text(1,1.4,'Programa elaborado por Hector Isaias Escobosa') ; echo off % fin del programa Cinetica_Quimica 9
% Al Ejecutar este Programa, se obtiene la siguiente Grรกfica.
Grafica-Cinetica de Reaccion A->B->C 2 XA XB XC
Aplicación # 2 -) Modelo para un PÉNDULO SIMPLE. Este es uno de los problemas de aplicación más representativos en Física. Consideremos un péndulo simple de masa ( m ) que cuelga de un punto fijo por medio de una vara de longitud (L), como se muestra en la siguiente figura, y que se mueve en un ambiente viscoso que amortigua la velocidad, con un coeficiente de amortiguamiento “b” ( en kg/mt/seg ). El ángulo (Ө) es el que forma la vara con la vertical. La velocidad de la masa ( m ) estará dada en función de ( L ) ; de ( dӨ / dt ) , y de la fuerza tangencial FT , mediante la Ө L
ecuación de movimiento siguiente : modelo matemático : FT = - m g Sen(Ө) – b L
…ec(1) ,
en donde g es la aceleración de la gravedad. De la segunda Ley de Newton, sabemos que: ---ec(2)
Igualando las ecuaciones (1) y (2), se obtiene la Ecuación Diferencial Ordinaria (EDO) que describe el movimiento de un Péndulo simple: . . . ec(*)
Si hacemos X(1) = Ө , y X(2) = q =
= velocidad angular , entonces la ec(*) que es una
EDO de 2° orden, se puede reducir a un Sistema de Ecuaciones Diferenciales de 1er orden: X(1) = X(2) X(2) = ( - b / m ) X(2) + ( - w / mL ) Sen[X(1)] , siendo w = mg
% Estas ecuaciones se resuelven mediante el siguiente Programa MATLAB. % primero creamos una Función de MATLAB function xpunto=péndulo(t,x); w=2;
% peso del péndulo ;
% longitud de la vara
b=0.08; % coeficiente de amortiguamiento; g=9.8; m=w/g;
% gravedad terrestre
% masa del péndulo
xpunto=[x(2) - (b/m)*x(2) – w/(m*L)*sin( x(1) ) ];
% finaliza la Función 11
% Programa de MATLAB à estudia comportamiento de un Péndulo simple % para los valores dados a los parámetros del sistema. to=0;
% tiempo inicial
% tiempo final
xo=[pi/2 , 0];
% condiciones iniciales asignadas
[t , x]=ode45(‘pendulo’,[to , tf], xo); Plot(t , x),grid,xlabel(‘posicion’),ylabel(‘velocidad’) Legend(‘Angulo’,’Velocidad angular’) % Fin de Sbrutina- function xpunto=péndulo(t,x) % Programa para resolver una EDO de 2° Orden, que MODELA % el COMPORTAMIENTO de un PENDULO SIMPLE: % mL (d2Y/dt2) + bL(dY/dt) + mg Sen(Y)=0 % Y representa al Angulo(Tetha) y (dY/dt)--> representa la Velocidad angular % Para valores pequeños de Y (Tetha), Sen(Y) à Y % el Modelo se puede expresar como sigue: (haciendo Y = y) % (d2y/dt2 ) + (b/m)*(dy/dt) + (g/L)*y =0
; y = f(t)
clc,clear all, echo on % Determinar la Solucion de la EDO: (d2y/dt2) + (b/m)*(dy/dt) + (g/L)*y =0 % Para la Condicion inicial Dy(0)=0, y(o)=pi/2 % se asignaran los siguientes valores para los parametros: a=0.4;
% a= (b/m) = 0.4;
k=49/3;
% k= (g/L)=49/3;
pause % oprima cualquier tecla para resolver la Ecuación Diferencial G=dsolve('D2y+0.4*Dy+(49/3)*y=0','Dy(0)=0','y(0)=pi/3') 12
G=simple(G) % La Solución que se obtiene (al ejecutar el programa) es: G = (pi*(1222*cos((3666^(1/2)*t)/15) + 3666^(1/2)*sin((3666^(1/2)*t)/15)))/(3666*exp(t/5)) pause pretty(G) % Ahora se construira la grafica para el intervalo [-6, 6] % Oprima cualquier Tecla para que se construya dicha grafica pause ezplot(G,[0,6]),grid,xlabel('Tiempo t (seg)'),ylabel('Posicion Pendulo') Text(-4,3,'Grafica del comportamiento de un Péndulo Simple') % La GRAFICA que se obtiene al ejecutar este Programa es la siguiente: (p (1222 cos((36661/2 t)/15) + 36661/2 sin((36661/2 t)/15)))/(3666 exp(t/5)) 3
Gráfica del comportamiento de un Péndulo Simple
Posición Péndulo
2 1 0 -1 -2 -3 -6
0 Tiempo t (seg)
(p cos((5 71/2 t)/4))/(2 exp(t/4)) + (71/2 sin((5 71/2 t)/4) (p + 8))/(70 exp(t/4)) 6
Gráfica para la Sol. EDO-Péndulo, con Dy(0)=1, y(0)=pi/2
Posición y(t)
% Variando valores de parámetros, tendremos otro caso, continuación. clear,clc % Ahora se determinara´ la Solucion , para cuando (b/m)=2 y (g/L)=2 % obteniendose la EDO à D2y+2*Dy+2*y=0 tiene Raices Complejas % para las Condiciones iniciales --> DY(0)=0 , y Y(0)=1 pause
% Oprima cualquier tecla para continuar
F=dsolve('D2y+2*Dy+2*y=0','Dy(0)=0','y(0)=1') % Ahora se determina la Solución F = cos(t)/exp(t) + sin(t)/exp(t) Y=simple(F) Y = (cos(t) + sin(t))/exp(t) % Ahora se construirá la grafica para el intervalo [0,10] ezplot(F,[0,5]),grid,xlabel('Tiempo t (seg)'),ylabel('Posicion y(t) '), text(-4,200,'Grafica Solución de la EDO: D2y+2*Dy+2*y=0 ; para Dy(0)=0 , y(0)=1') 14
(cos(t) + sin(t))/exp(t)
Gráfica Solución de la EDO: D2y + 2*Dy + 2*y = 0 ; para Dy(0)=0 , y(0)=1
_______________________________________________________________ % Ahora se determina la Solución , para cuando (b/m) = 4 y (g/L) = 5 % obteniendose la EDO --> D2y+4*Dy+5*y=0 % para las Condiciones iniciales --> DY(0)=0 , y Y(0)=1 pause %
F=dsolve('D2y+4*Dy+5*y=0','Dy(0)=1','y(0)=2')
F = Y=dsolve('D2y+4*Dy+5*y=0','Dy(0)=1','y(0)=2') F = (2*cos(t))/exp(2*t) + (5*sin(t))/exp(2*t) % Ahora se construira´ la grafica para el intervalo [0,5] Pause
% oprima cualquier Tecla para que se construya dicha grafica
ezplot(F,[-1,6]),grid,xlabel('Tiempo t (seg)'),ylabel('Posicion y(t)'), text(1,1,’Solución del Modelo: D2y+4*Dy+5*y = 0, Dy(0)=1 , y(0) = 2’) 15
(2 cos(t))/exp(2 t) + (5 sin(t))/exp(2 t) 2
Solución del Modelo: D2y+4*Dy+5*y=0 , para Dy(0)=1 , y(0)=2
2 3 Tiempo t (seg)
% Cambiando el Rango de impresión de à [0 , 5] y el titulo de la Gráfica, se obtiene: Grafica Comportamiento Pendulo, realizo: Hector Isaias 1
Posicion Pendulo
2.5 3 Tiempo(seg)
APLICACIÓN #3-) Estudio de un Problema de Selección de un Producto Agrícola Perecedero. Dentro de una compañía empacadora y procesadora de alimentos perecederos, se presentó el siguiente problema: ¿Cómo separar en grandes cantidades los Tomates verdes de los ya maduros, de manera rápida, económica, confiable y automática? Los cultivadores de la región han encontrado que les resulta más económico cosechar sus "Huertas", por medio de máquinas que arrancan la planta completa (hojas, tallo, raíz, tomates, etc), para luego separar por medio de Tolvas los Tomates verdes y semi-maduros del resto de la planta. Existen, por supuesto, muchas ideas acerca de la manera en que puede llevarse a cabo la separación de los Tomates, como por ejemplo: · · · · · · · · · · · ·
Por el Color del Tomate. Aprovechando la propiedad de Emisividad-Reflectividad de la superficie del Tomate. Por la Dureza relativa del Tomate. Por el Tamaño del Tomate. Por el Peso del Tomate. Midiendo la Conductividad eléctrica del Tomate. Por el Volumen específico del Tomate. Por su Flotabilidad en el agua ú otro medio. Usando "Rayos X" para determinar el tamaño de las semillas. Midiendo la Conductividad Térmica del Tomate. Por su Capacidad Calorífica específica. Usando Números aleatorios (Randoms Numbers).
Establecimiento del Problema: Se quiere hacer la Evaluación de una idea que al parecer puede tener éxito; la idea consiste en separar los tomates verdes de los ya maduros o semi-maduros, aprovechando su propiedad de "Flotabilidad en el agua".
Una prueba experimental sencilla, demuestra que los tomates verdes flotan en el agua, en cambio los maduros y semi-maduros tienden a depositarse en el fondo de un tanque con agua; por ello, se requiere determinar el tiempo que tarda un tomate verde en emerger a la superficie del tanque cuando se le suelta desde una profundidad conocida. La siguiente figura, ilustra el presente problema.
Figura que ilustra el problema de separación de los Tomates verdes de los maduros. Tanque con Agua
Establecer (Definición) el Problema: Se requiere determinar el tiempo que tarda en llegar a la superficie de un tanque con agua, un tomate de aproximadamente 134 gramos de peso, con un volumen de 142 cm3, si se suelta desde una profundidad de 60 cm en un gran tanque conteniendo agua común. Problema - Encontrar una Relación entre: PT
à Peso específico del Tomate.
à Peso específico del Agua. Ho
à Diametro del Tomate.
à Profundidad.
Formulación del Modelo Matemático: Este es un Problema relacionado con el flujo lento de un fluido incompresible, moviendose alrededor de un objeto sumergido en agua, que tiene un eje o un plano de simetría paralelo a la velocidad de aproximación del medio líquido. Marco Teórico - se aplican: · Ecuaciones para la cantidad de movimiento · Coeficiente de Resistencia 18
Suposiciones: .- Se considera que la forma geométrica del Tomate es una esfera de radio uniforme. .- Se considera que la densidad del Tomate es constante, se desprecian los efectos de la Presión y de la Temperatura. Aplicación de Principios Fisicos: Haciendo un diagrama de cuerpo libre de la esfera y aplicando la segunda ley de Newton, se tiene: åF= m a
Do FW. En donde: à FB = PW VT g
Fuerza de Empuje del Agua.
Volumen del Tomate.
Aceleración local de la gravedad.
Peso del Tomate. à FW = PT VT g
Densidad del Tomate.
Fuerza de arrastre del fluido. à FD = CD PW V2 A /2
Coeficiente de arrastre del fluido.
Velocidad del Tomate.
Sección transversal de una esfera.
Diametro de la esfera (Tomate).
A = p Do / 4
Sustituyendo estas expresiones en la ec(2), tenemos:
PW VT g - PT VT g - CD PW V2 A /2 = PT VT a
De esta ecuación (3), se obtiene: 19
...ec(4) 2PT VT
Pero como a = ( d2Y / dt2 ) , y haciendo K1 = [(PW /PT) - 1] g ; y K2 = (pPW CD)/4Do Entonces, la ecuación (4) queda como sigue:
( dY2 / dt2 ) = K1 - K2 V2 ...ec(5) La ecuación (5), es una EDO que representa o modela el Problema analizado. Verificación de Resultados: Para llevar a cabo una verificación del modelo resultante, se requiere hacer alguna suposición adicional, por ejemplo:
1º Caso: Suponer que la fuerza de arrastre es poco significativa, en cuyo caso a partir de la ec(5) se determina el valor de la aceleración:
(dY2/dt2) = a = [ (PW /PT) - 1 ] g = [(1000/961.2) - 1 ] x 9.81 = 0.396 Mts/seg2. Si la aceleración es constante y Ho = 0.6 Mts, entonces, el tiempo que tarda el Tomate en subir a la superficie del tanque es:
t = [2 Ho / a ]1/2 = [2(0.60)/0.396]1/2 = 1.74 seg. Sin embargo, falta determinar si la fuerza de arrastre es verdaderamente poco significativa; para ello, es necesario evaluar su magnitud. Con el valor de tiempo ya determinado, calculamos la velocidad promedio del Tomate:
Vprom = ( Ho / t ) = (0.6 / 1.74 ) = 0.3448 mts/seg. Ahora, se evalua el Número de Reynolds:
RE = D V r / m = (0.3448)(0.065)(756.58) / (6.222x10-5) = 35800
Con este valor se determina el coeficiente de arrastre en la esfera: CD = 0.50
De la curva tomada de C.E. Lapple "Dust and Mist Collections, Chemical's Enginnering Handbook". J.H. Perry. Mc-Graw-Hill, New York. Su forma aproximada es la siguiente:
103 Factor de 102
Fricción ó
Coeficiente de arrastre 10
3.6x104.
Número de Reynolds RE = D V r / m
Por lo tanto, el valor promedio para la fuerza de arrastre es:
FD = CD (PW V2 A)/2 = (0.5){756.58*0.3448* p* (0.0646)2 }/(2*4) = 0.0098 Kg. Ahora, la Fuerza Neta debida al empuje del agua y el Peso del Tomate es:
FNeta = g (PW - PT) VT = (9.81){(1000 - 961.2)/9.80665}(142x10-6) = 0.00551 Kg.
Al comparar la fuerza promedio de arrastre y la fuerza neta debida al empuje del agua y el peso del Tomate, se puede observar que la Fuerza de arrastre de ninguna manera puede ser despreciada, por lo tanto se tendrá que hacer una suposición distinta.
2º Caso: La experiencia demuestra que una partícula que se desplaza con movimiento, ya sea hacia arriba o hacia abajo dentro de un fluido, alcanza su velocidad final muy rápidamente y luego se sigue moviendo con velocidad constante; ello significa, que la mayor parte del trayecto lo efectúa con velocidad constante (es decir, con aceleración nula a = 0), luego entonces, la ecuación de movimiento para este caso; la ecuación (4), queda como: (PW - PT ) g
[ (PW /PT) - 1 ] g = ( pPW V2 CD)/4Do PT.
Suponiendo un valor para el coeficiente de arrastre CD = 0.45, se obtiene una velocidad de V=0.27 mts/seg. El Número de Reynolds para este valor de velocidad es de 2.8x104, de la gráfica del coeficiente de arrastre Vs Número de Reynolds se observa que corresponde a 0.45, luego entonces éste es el valor correcto. Con este valor de velocidad se puede determinar el tiempo necesario para que el Tomate llegue a la superficie del tanque; esto es,
t = [ Ho / V ] = (0.60/0.27) = 2.22 seg. En la práctica, se há comprobado de manera experimental, que dicho valor de tiempo es correcto. _________________________________________________________________________
Ahora veremos la Respuesta que se obtiene al analizar este Problema usando Matlab. Para ello, consideraremos que k1 = 0.396 y k2 = 0.5 (valor aproximado) Quedando el Modelo Matemático, como sigue: D2y + 0.5 Dy = 0.396 , para obtener Solución General, hacemos >>Y= dsolve(‘D2y+0.5*Dy=0.396’) % la respuesta de Matlab es: Y = C1*exp(-t/2) + 99*t/125 + C2 . Ahora, para las condiciones iniciales siguientes: Dy(0) = 0.02 e y(0) = 0.6 , tendremos lo siguiente:
>>Y= dsolve(‘D2y+0.5*Dy=0.396’,’Dy(0)=0.02’,’y(0)=0.6’) % la respuesta de Matlab es: 22
Y = (193 / 125)*exp(-t/2) + (99 / 125)*t + (- 118 / 125) .
Gráfica para el Modelo: D2y+0.5*Dy=0.396 , con Dy(0)=0.02 , y(0)=0.6 3.5 3 2.5
Valor de F(t)
Se probaron diferentes condiciones iniciales Por ello se obtienen diversas gráficas Libro: Métodos Numéricos con MATLAB; Prof. Escobosa
1 Tiempo (seg)
Se observa que, aproximadamente, en t = 2.22 seg. à F(t) = 0
Aplicación #4-) Modelo General para un CIRCUITO ELÉCTRICO Elemental. ( una Ecuación Integro-Diferencial ) . . .ec(1) En donde: R
R = Resistencia ( se mide en Ohms )
L = Inductancia ( se mide en Henrryios ) C = Capacitancia ( se mide en Faradios )
E(t) = Fuerza Electromotríz (en Volts) I = I(t) à Corriente Eléctrica ( Amperios)
L = Inductancia (magnitud de Campo magnético inducido por una corriente eléctrica I ) presente en el circuito eléctrico básico ( en Henryos )
R= parámetro que representa la magnitud de una Resistencia presente en el Circuito ( en Ohms = Ω. ) C = Parámetro que representa la Conductancia de un Condensador presente en el Circuito (se mide en faradios. )
I = I(t) variable que representa la magnitud de la Corriente eléctrica que fluye a través del Circuito (en Amperios.) en cualquier instante de tiempo t > 0. E( t ) = Fuerza Electromotriz de la Fuente de Energía del Circuito ( en Voltios )
t = variable tiempo ( en segundos ) Nota: como
ʃ I dt
, de donde se obtiene el siguiente modelo.
Modelo General para un CIRCUITO ELÉCTRICO una Ecuación Diferencial de 2° orden, en donde q= q(t) es la carga eléctrica, variable dependiente del tiempo t > 0. . . . ec(2) para cuando se desea encontrar 24
q ( la carga eléctrica en función del tiempo t > 0 ) Modelo General para un CIRCUITO ELÉCTRICO una Ecuación Diferencial de 2° orden, en donde I = I(t) es la corriente eléctrica, variable dependiente del tiempo t > 0. para cuando se desea encontrar I ( la corriente eléctrica que fluye a través del circuito, en función del tiempo t > 0 )
% Enseguida veremos el Programa MATLAB para determinar la Solucion del Modelo ec(2) % Matematico ( EDO de 2° Orden ) que representa el comportamiento % de un Circuito Electrico LRC} % L(dI/dt)+ RI + (1/C) q(t) = E , o tambien: % I'(t) = (-R/L)*I(t)+ (-1/L*C)q(t)+ E/L clc,clear, clf, echo on, hold on % Oprima cualquier Tecla para iniciar a correr el programa pause % Primero introduciremos valores para los Parametros R=100; % Resistencia en Ohms C=0.001; % Capacitancia en faradios L=0.20; % Inductancia en Henrryos E=1.0; % Voltaje en Voltios h=0.00025; % Tolerancia en error de calculo n=1; % valor de inicio de numero de iteracciones t(1)=0; Y(:,1)=[0;0] % Y representa la corriente electrica I(t)=Y M=[0,1;-1/(L*C),-R/L] ; S=[0 ; E/L]; while n<101 % inicia ciclo iterativo Y(:,n+1)=Y(:,n)+h*(M*Y(:,n)+S); t(n+1)=n*h; n=n+1; end % fin del ciclo iterativo % Ahora se construye la Grafica que ilustra el comportamiento de % de la Carga y de la Corriente Eléctrica que fluye a través del circuito pause % Oprima cualquier Tecla para que se construya dicha grafica plot(t,100*Y(1,:),t,Y(2,:),'--') text(t(30),119*(Y(1,30)),'Carga Eléctrica à q*100') text(t(28),Y(2,30),' comportamiento de à I = Corriente Electrica') xlabel('tiempo ( t )'),ylabel('I(t) en Ampers , q(t) en Coulombs') % title(‘Gráfica Solución Modelo para Circuito Eléctrico’) 25
% A continuación se presenta la Grafica que se obtiene con Matlab xlabel('tiempo ( t )'),ylabel('I(t) en Ampers , q(t) en Coulombs') % Fin del Programa
Gráfica Solución Modelo para Circuito Eléctrico 0.025
comportamiento de --> I = Corriente Electrica Carga Eléctrica -> q*100
Aplicación #5-) Modelo General para un SISTEMA MASA-RESORTE. (Vibraciones Mecánicas Elementales )
m = masa del objeto unido al Resorte ( en Kg. ) b= representa Amortiguación en el Sistema (Kg / seg. )
K = parámetro que representa la Elasticidad del Resorte ( en Newton / mt. )
X = variable desplazamiento (posición) del Sistema en cualquier instante de tiempo t > 0 . F( t ) = Fuerza externa que actúa sobre el Sistema ( en Newton )
t = variable tiempo ( en segundos )
Aquí, analizaremos varios casos, a través de Programas Matlab. % Primero veremos un Programa, para resolver una EDO de 2° Orden % que represente el caso más sencillo, cuando se tiene un Sistema No-Amortiguado y además % No-Forzado; esto es, cuando b = 0 y F(t) = 0. clc, clear Pause % En Sistema Masa Resorte, la Posición se representara por à y = y(t) % se trata de determinar la Solucion General de la EDO à D2y + 9y = 0 % Un Oscilador Armonico Simple, para las condiciones iniciales % Dy(0)= 0
Y=dsolve('D2y+9*y=0','Dy(0)=0','y(0)=1') Y=cos(3*t)+sin(3*t) % Solución para este caso % Ahora se construira´ la grafica para la Solucion pause % Oprima cualquier Tecla para que se construya dicha grafica % t=-2*pi:pi/50:2*pi; Y=cos(3*t)+sin(3*t); plot(t,Y),grid,xlabel('eje x'),ylabel('eje y') title('Grafica Solución del Modelo D2y + 9*y = 0 , para Dy(0)=0 , y(0)=1') Movimiento Armónico Simple Gráfica Solución del Modelo D2y + 9*y = 0 , para Dy(0)=0 , y(0)=1 1
% Segundo caso: veremos un Programa, para resolver una EDO de 2° Orden % que represente a un Sistema Amortiguado, y además Forzado % en donde masa = m=1 kg, coeficiente de Amortiguación = b=4 kg/seg, % Elasticidad = k=5 Newton/mt, y Fuerza Externa=F(t)=exp(-0.1*t) clc, clear 28
% Resolveremos la EDO à D2y + 4*Dy + 5*y = exp(-0.1*t) F=dsolve('D2y+4*Dy+5*y=exp(-0.2*t)','Dy(0)=0','y(0)=1') % un Sistema Amortiguado y Forzado m=1, b=4, k=5 F(t)=exp(-0.1*t) ezplot(F,[-1,4],'b') xlabel('Tiempo t (seg)'),ylabel('Posición y(t)') text(0,0.5,’Gráfica Solución de D2y + 4*Dy + 5*y = exp(-t)') plot(t,g1),grid,xlabel('eje x'),ylabel('eje y'),title('Sistema Masa-Resorte Amortiguado')
cos(t)/(2 exp(2 t)) +...+ (cos(t) (cos(t) - sin(t)))/(2 exp(t)) 1 0.8 0.6 Posición y(t)
Gráfica Solución de D2y + 4*Dy + 5*y = exp(-t) 0.4 0.2 0 -0.2 -0.4 -0.6 -1
1.5 2 Tiempo t (seg)
Como un ejemplo adicional, Determinaremos la Solucion de la EDO à D2y+2*Dy+2*y=0 % que representa a un Sistema Masa-Resorte Amortiguado, pero No-Forzado % para las condiciones iniciales à Dy(0)=1 , y(0)=2 clear all,clf,clc t=-2*pi:pi/50:2*pi; pause F=dsolve('D2y+2*Dy+2*y=0') % Aquí Matlab arroja el siguiente Resultado (Solución General)
F=(C1*Cos(t) + C2*Sin(t) )/exp(t) % Ahora, resolveremos para la condición inicial dada: Y=dsolve('D2y+2*Dy+2*y=0','Dy(0)=0','y(0)=1') % La Solución que arroja Matlab es: 29
Y=(Cos(t) +Sin(t) )/exp(t) % un Sistema Amortiguado, donde à m=1 kg , b=2 kg/seg , k = 2 Newton/mt. F(t) = 0 Pause % Ahora se construye la Gráfica de la Solución Posición del Sistema ezplot(Y,[-1,4],'b') t=-2*pi:pi/40:2*pi ; g=(sin(t)+cos(t))/exp(t); plot(t,g),xlabel('Tiempo (seg)'),ylabel('Posición y(t)') text(0,-0.4,'Gráfica Solución de D2y+2*Dy+2*y = 0 ') cos(t)/exp(t) + sin(t)/exp(t) 1 0.8 0.6
Gráfica Solución de D2y + 2*Dy + 2*y = 0
Otro Programa de aplicación-) Resolución de la EDO de 2º orden que representa a un Sistema Masa - Resorte % Aquí, se utilizan las Ecuaciones de Runge-Kutta ( de cuarto orden). % k1 = hf(yn,tn) % k2 = hf(yn+k1/2,tn+1/2) % k3 = hf(yn+k2/2,tn+1/2) k4 = hf(yn+k3,tn+1) % fórmula recursiva % yn+1 = 1/6(k1 + 2k2 + 2k3 + k4)*yn clc, clear, clg , echo on M = 0.5; k = 100; B = 0 ; a = B/M ; b = k/M ; n = 0 ; t = 0 ; h = 0.025; 30
Y(: , 1) = [1 ; 0 ] ; t(1) = 0 ; % Condiciones Iniciales Q = [1 , 1 ; -b , a ]; % inicio de iteracciones while t<=1 n = n+1; k1 = h*Q*Y(: , n) ; k2 = h*Q*(Y(: , n) + k1/2); k3 = h*Q*(Y(: , n) + k2/2) ; k4 = h*Q*(Y(: , n) + k3); Y(: , n+1) = Y(: , n) + (k1 + 2*k2 + 2*k3 + k4)/6; t(n+1) = n*h ; end pause; plot(t ,Y(1 , :) ,'b', t , Y(2 , :) ,'m' ); xlabel('tiempo (t)'); ylabel(' X e V '); % X = posición , V = velocidad axis([ 0, 1, -20 , 20 ] ); % fin del programa se obtiene la siguiente gráfica de Matlab 20 15 10
t iem po (t)
Documento elaborado por: Héctor Isaias Escobosa Iribe
Uso de MATLAB para resolver EDO de Orden Superior. Una EDO de orden superior puede escribirse como un Sistema de Ecuaciones Diferenciales de primer orden acopladas, utilizando cambios de variable apropiados. Por ejemplo, consideremos la siguiente Resolución
de las EDO de Lotka y Volterra que representan a un Sistema Presa-Depredador, a través de un Programa "Script", como veremos a continuación. % % Resolución numérica de las ODE = Ordinary differential equations de Lotka y Volterra.
echo on , clc , clf % % % % % %
ODE23 y ODE45 son funciones para la solución numérica de ODE Estas funciones utilizan el procedimiento de integración de Runge-Kutta-Fehlberg ODE23 utiliza ecuaciones de integración de Runge-Kutta de 2º and 3º orden ODE45 utiliza ecuaciones de 4º y 5º orden para una mayor exactitud. Este Ejemplo muestra como utilizar dichas funciones para Resolver ODE
pause clc
% para hacer una pausa; "Oprime una tecla para continuar. % para limpiar la pantalla de video % Considere el Par de ODE conocidas como las "Ecuaciones de Lotka-Volterra" % que representan el comportamiento (ó modelo) de un Sistema Presa-Depredador. % % y1' = (1 - alpha*y2)*y1 % y2' = (-1 + beta*y1)*y2 % % las funciones y1 e y2 miden el tamaño de la Población de Presas y de Depredadores, % respectivamente. El término 2°grado proviene de las interacciones de las dos especies. % Note que la Población de Presas aumenta cuando no está presente Depredadores; % la Población de Depredadores disminuye cuando no hay Presas % Para Simular este Sistema, hay que crear una función (function), en un Archivo.M % Por ejemplo LOTKA.M; ésta función, regresa las derivadas de las variables de estado, % sus valores para diferentes valores de t tiempo. Para este caso, el archivo será
function yp = lotka(t, y) % Lotka-Volterra à Modelo Presa-Depredador
yp = diag( [ 1 - 0.01*y(2) , -1 + 0.02*y(1) ]*y ; Pause % Para Simular el Sistema de ODE definidas en LOTKA.M , para el intervalo % 0 < t < 15, utilizaremos la instrucción ODE23:
t0 = 0 ; tfinal = 15; y0 = [20 20]'; % Define las condiciones iniciales. % [t,y] = ode23('lotka',t0,tfinal,y0);
pause tol = 1.e-3; % grado de exactitud trace = 1; tfinal = tfinal*(1 + eps); % Continuación del Programa (Script)
[t, y] = ode23('lotka', t0, tfinal, y0, tol, trace); plot(t,y), grid, 32
title('Resolver Sistema de Ecuaciones de Lotka-Volterra, usando ODE23 '); pause % Ahora simularemos el Sistema utilizando ODE45 % [T,Y] = ode45('lotka', t0, tfinal, y0 );
pause [T,Y] = ode45('lotka', t0, tfinal, y0, tol, trace); plot(T,Y), grid pause title('Grรกfica: Sistema de Ecuaciones de Lotka-Volterra '); xlabel('Eje x '); pause echo off end Grรกfica: Sistema de Ecuaciones de Lotka-Volterra 350 300 250 200 150 100 50 0 0
80 100 Eje x
Otro ejemplo lo tenemos en la EDO de Van-Der-Pool , cuyo modelo es:
D2y = ( 1 - y2 ) Dy - y Y" = g(x , y , y' ) = y' (1 - y2 ) - y Primero se definen dos nuevas funciones : U1( x ) = y' U2( x ) = y Luego obtenemos el siguiente sistema de Ecuaciones Diferenciales de primer orden acopladas: U1' = y" = g( x , U2 , U1) = U1( 1 - (U2)2 ) - U2 U2' = U1 Para resolver este Sistema de Ecuaciones, primero se define ( se construye ) una Funciรณn que calcule los valores de las ecuaciones, en un Archivo.M
Function Uprima = ecua2 (x , U ) Uprima(1) = U(1) * (1 - U2^2 ) - U(2); Uprima(2) = U(1); % Luego, para resolver el Sistema de Ecuaciones Diferenciales de primer orden, sobre el intervalo [0,20] % para las condiciones iniciales y' (0) = 1.0 , y(0) = 0, utilizamos las siguientes instrucciones de MATLAB:
Inicial = [ 0 10 ] ; [x , ynum] = ode23('ecua2' , 0 , 10 , inicial ) ; Subplot(2 1 1), Plot(x , ynum(: , 1 ) ), . . . Subplot(2 1 2), Plot(x , ynum(: , 2 ) ), . . . Title(' Grรกficas correspondientes a la Soluciรณn de la EDO de VanDerPoll'), grid 5 4.5 4 3.5 3
Grรกficas correspondientes a la Soluciรณn de la EDO de Van er oll D
5.3 5.2 5.1 5 4.9 0
Problema del Paracaidísta: Un paracaidísta de masa m kg salta desde un avión en t = 0. Suponiendo que la velocidad vertical inicial del paracaidísta es v(0) = 0, que la caida es vertical y que el arrastre aerodinámico está dado por Faire = k v2 , donde k es una constante y v es la velocidad vertical (positiva hacia abajo), determine: a-) El modelo matemático que nos permita obtener la velocidad vertical del paracaidísta, en cualquier instante de tiempo t ³ 0 . b-) Construya la gráfica (usando un programa en Matlab), velocidad contra tiempo, para cuando m = 70 kg , k = 0.27 kg/mt , en el intervalo de tiempo t Î[ 0 , 20 ] seg . a-) El modelo matemático para este problema del paracaidísta, se obtiene aplicando la segunda ley de Newton, para el equilibrio de fuerzas à m ( dv / dt ) = - Faire + mg . Sustituyendo Faire = - k v2 , en esta ecuación, obtenemos: El modelo para la velocidad del paracaidísta : dv / dt = (- k / m ) v2 + g , para v ( 0 ) = 0 Nota: El modelo para la posición del paracaidísta estará dado por: d2y / dt2 = (- k / m) (dy / dt)2 + g ; dado que v = dy / dt , para y( 0 ) = 0 ^ y ' ( 0 ) = Dy( 0 ) = 0
b-) Programa para resolver, numericamente, el problema del Paracaidísta, y obtener la gráfica de su comportamiento . % Resolución numérica del problema, usando Software MATLAB. % Autor: M.C. Héctor Isaias Escobosa Iribe
Clc , Echo on , Clear, clf, hold off Pause % Oprima una Tecla para iniciar la Resolución del Problema t = 0 ; n = 0 ; v = 0; k = 0.27 ; m = 70 ; g = 9.8 ; h = 0.1 ; % t_rec(1) = t ; v_rec(1) = v ; % While t < = 20 n=n+1; v = v + h * (( -k / m)*v^2 + g ); t=t+h; v_rec(n + 1) = v ; t_rec(n + 1 ) = t ; end pause % Se presenta la Solución Gráfica del Problema del Paracaidísta plot(t_rec , v_rec) , grid , xlabel('Tiempo t (seg)') , ylabel('valor de Vx10 (mt/seg)'), title(' Soluciones al problema del Paracaidísta') pause % Fin del Programa end 35
% Gráfica de la Solución al problema del Paracaidísta, usando MATLAB. Soluciones al Problema del Paracaidista 3.5
Valor de Vx10 mt/seg
2.5 3 Tiempo t (seg)
Gráfica de la Solución al problema del Paracaidísta 80 70 60 50 40 v (mt/seg) 30 20 10 0 0
Tiempo (seg) 36
Resolución Numérica de EDO usando MATLAB. MATLAB contiene dos Funciones para determinar soluciones numéricas de EDO, Ecuaciones Diferenciales Ordinarias: Ode23 y Ode45. Estas funciones utilizan el procedimiento de integración de Runge-KuttaFehlberg. La Función Ode23 utiliza ecuaciones de integración de Runge-Kutta de 2º and 3º orden La Función Ode45 utiliza ecuaciones de 4º y 5º orden para una mayor exactitud. Las dos Funciones tienen los mismos tipos de argumentos, por lo cual sólo estusiaremos la utilización de la Función Ode23. La forma más simple de la Función Ode23 requiere de 4 argumentos: Ø El primer argumento es el nombre ( 'entre comillas' ) de una Función de Matlab que devuelve el valor de la Ecuación Diferencial y' = Dy = g(x , y) para los valores asignados a x) e (y). Dicha Función debe creada como un Archivo.M ( M-File ) de Matlab. Ø El segundo y tercer argumentos representan los extremos del intervalo sobre el cual se evalurá la función y =f(x). Ø El cuarto argumento contiene la "Condición Inicial" ó el punto límite que se necesita para determinar una solución única de la EDO. El siguiente ejemplo muestra como utilizar dichas funciones para Resolver ODE. Por ejemplo, supongamos que se desea obtener la solución numérica de la siguiente EDO de primer orden: 2 (dy / dx ) = 3 X , para la condición: y (2 ) = 0.5; esto es, y = 0.5 cuando x = 2. 3
En este sencillo caso, la solución analítica es y = X - 7.5 Las instrucciones de Matlab que definen y resuelven la EDO anterior, asumiendo como entradas escalares para ( x ) e ( y ), son las siguientes: % Recuerde que primero debe crear la Función como un Archivo.M % Por ejemplo la Función G1.M % A continuación se construye el Programa de Matlab que resuelve la EDO clc % para limpiar la pantalla de Video
function dy = g1( x , y ) dy = 3*X^2 ; [x , ynum] = ode23( 'g1' , 2 , 4 , 0.5 ) ; yana = X..^3 - 7.5 Plot( x , ynum , x , yana , ' 0 ' ) , grid, xlabel('Eje x'), ylabel('Eje y') Nota-1-) La Variable ynum es la Solución numérica de la EDO ; en tanto que yana contiene la Solución analítica. Nota-2-) Todas estas instrucciones se pueden introducir (editando) en un Programa (Archivo.M) de Matlab, que le puede servir de base para resolver otras EDO. Nota-3-) Las Funciones Ode23 y Ode45 pueden ser utilizadas con dos Parámetros (ó argumentos) adicionales; un quinto argumento puede utilizarse para especificar la tolerancia que se relaciona con el tamaño del incremento utilizado por las funciones; las tolerancias predeterminadas son 0.001 para Ode23 y 0.000001 para Ode45. Un sexto parámetro puede utilizarse para que se impriman los resultados intermedios. 37
Como un segundo ejemplo, supongamos que se desea obtener la solución numérica de la siguiente EDO de primer orden: (dy / dx ) = 3 y + Î2X , para la condición: y (0 ) = 3; esto es, y = 3 cuando x = 0. En este otro sencillo caso, la solución analítica es y = 4 Î3 - Î2 Las instrucciones de Matlab que definen y resuelven la EDO anterior, asumiendo como entradas escalares para ( x ) e ( y ), son las siguientes: % Recuerde que primero debe crear la Función como un Archivo.M % Por ejemplo la Función G2 . M % A continuación se construye el Programa de Matlab que resuelve la EDO clc % para limpiar la pantalla de Video echo on % para que se desplieguen las instruccines en la pantalla al correr el programa pause % Oprima una Tecla para continuar t
function dy = g2( x , y ) dy = 3*y + exp(2 * x) ; [x , ynum] = ode23( 'g2' , 0 , 5 , 3 ) ; yana = 4*exp(3*x) - exp(2*x) ; Plot( x , ynum , x , yana , ' 0 ' ) , grid, xlabel('Eje Tiempo t'), ylabel('valor de Y') % 7
4 exp(3 t) - exp(2 t)
2 3 Eje Tiempo t
Documento elaborado por: M.C. Héctor Isaias Escobosa Iribe. Ponencia CIIDET ; Congreso de Nov, de 2010. 38
Uso de Matlab en la solución de algunos sistemas dinámicos Héctor Escobosa
Una de las aplicaciones más importantes de las Ecuaciones Diferenciales (Ordinarias y Parciales), se tiene en el estudio de Sistemas Dinám...

References: Resolución 
 resolución 
 Resolución 
 Resolución

 Resolución 
 Resolución 
 Resolución 

Resolución