Ecuaciones Diferenciales en Matlab PDF

Preview

Summary

Practica de laboratorio, acerca del desarrollo y solución de problemas de ecuaciones diferenciales utilizando Matlab....

Description

    Tema: ECUACIONES DIFERENCIALES EN MATLAB

Resolución exacta de ecuaciones diferenciales ordinarias Para resolver de forma exacta una o varias ecuaciones diferenciales, Matlab dispone de la orden dsolve. Por defecto, la variable independiente es t, pero se puede usar cualquier otra variable si se incluye como último argumento:

La letra D se utiliza para representar la derivación con respecto a la variable independiente, es decir, u' se escribe Du; las derivadas orden superior u'', u''',... se escriben D2u, D3u,. . . Cuando se resuelve un sistema de ecuaciones diferenciales, Matlab proporciona las funciones solución en orden léxico-gráfico. Ejemplos 1. Obtener la solución general de una ecuación diferencial:

>> u = dsolve('Du=u/2') u= C1*exp(1/2*t) Esto significa que la solución de la ecuación diferencial u' = u/2 es: u(t) = C1 e

t/2

;

(C1 € R)

2. Resolver el problema de valores iniciales:

Basta escribir >> u = dsolve('Du = u/2','u(0) = 1/4') Se obtiene u = 1/4*exp(1/2*t)

   

Es decir, la solución es

La solución se puede representar gráficamente usando ezplot; por ejemplo, en el intervalo [0,3]. >> ezplot(u,[0 3])

3. Resolver el problema:

>> u = dsolve('Du = r*u','u(0) = u0') y se obtiene u = u0*exp(r*t)

Ejercicios Propuestos Calcula la solución general de la ecuación logística u'= (a – b*u ) u

Calcula la solución del problema anterior para t >=0, teniendo como condición inicial u(0)= u0 Para los siguientes valores de u0: 10, 20, 50

   

RESOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES ORDINARIAS Cuando se resuelve de manera exacta un problema de valor inicial del tipo

Se obtiene como solución una función u(.) que puede ser evaluada en cualquier instante t del intervalo [TI ; TF ]. En cambio, cuando el problema se resuelve mediante un método numérico, sólo es posible obtener una colección finita {u1,u2, …, un } de aproximaciones de la solución en ciertos instantes t1 < t2 < …, tN, con TI = t1 y TF = tN, es decir,

Matlab dispone de una amplia variedad de métodos numéricos para resolver problemas de valor inicial. Destacamos dos de ellos:

En ambos casos, se obtiene como respuesta un vector columna t que contiene los instantes en los que se ha calculado la solución aproximada, y el vector u - o una matriz, en el caso de un sistemaque contiene tales aproximaciones. Los instantes contenidos en t son elegidos automáticamente por el método numérico; si se desean instantes predeterminados, basta reemplazar el argumento [TI TF] por [T1 T2 , … TN]. La solución aproximada se puede representar gráficamente con la orden plot. MATLAB dispone de varias funciones para resolver mediante procedimientos numéricos ecuaciones diferenciales. Elegiremos ode45 y ode23tb para resolver la mayor parte de los problemas.

   

La función ode45 Su sintaxis es la siguiente

[t,x]=ode45(odefun,tspan,x0, options, params) x es una matriz donde cada columna corresponde a las variables dependientes y t es el vector tiempo. odefun es el nombre de la función, tspan especifica el intervalo de tiempo, un vector de dos números tspan=[ti,tf], tiempo inicial y final. Para obtener valores de las variables dependientes en instantes concretos0 t, t 1 , t2, ... tn . se escribe tspan=[t0,t1....tn]; x0 es un vector que contiene los valores iniciales. options es una estructura que se crea con la función odeset, params son parámetros que queremos pasar a la función odefun

En la mayor parte de los ejemplos, se emplean los tres primeros parámetros: la función ode45 y le pasaremos la función odefunc, los instantes inicial y final en el vector tspan y las condiciones iniciales en el vector x0. La función ode23 [x,y] = ode23('función',a,b,inicial)

Esta instrucción regresa un conjunto de coordenadas "x" y "y" que representan a la función y=f(x), los valores se calculan a través de métodos Runge-Kuta de segundo y tercer orden. El nombre "función", define una función que representa a una ecuación diferencial ordinaria, ODE23 proporciona los valores de la ecuación diferencial y'=g(x,y). Los valores "a" y "b" especifican los extremos del intervalo en el cual se desea evaluar a la función y=f(x). El valor inicial y = f(a) especifica el valor de la función en el extremo izquierdo del intervalo [a,b].

EJEMPLOS DE ECUACIONES DIFERENCIALES DE PRIMER ORDEN EN MATLAB

1. Simular la solución de la ecuación diferencial:

    Para las condiciones iniciales (to = 0, yo = 1.0) y determinar el valor de la ordenada y para el correspondiente tiempo t = 0.2. Se utiliza el programa ode23tb, con un incremento del tiempo de ฀t = 0.01, el cual es utilizado para la determinación de la precisión en la obtención de los cálculos.

Editor- C:\ Users\ Ecuaciones\ \metodosnumericos.m

Documents\

MATLAB\

Ecuac.

%El programa resuelve la ecuación diferencial y´= 1- t + 4y, para y(o) = 1 function Ecuaciones tspan = 0 : 0.01: 0.2; y0 = [1]; ode23tb(@f,tspan,y0); function y1 = f(t,y); y1 = [1-t+4.*y]; [t ,y, y1] RESULTADOS:

Dif.

De

1er

orden

   

PRIMERA CLASE DE PROGRAMA Resolver la ecuación diferencial:

Editor- C:\ Users\ Ecuaciones\ Documents\ MATLAB\ Ecuac. Dif. De 1er orden \Ejercicio1.m %El programa resuelve la ecuación diferencial y´= 20 Cos(60t)-5y, yo=0 function Ejercicio1 tspan = 0 : 0.001: 0.5; y0 = [0]; ode23tb(@f,tspan,y0); function y1 = f(t,y); y1 = [20.*cos(60.*t)-5.*y]; [t ,y, y1] Respuestas, primeros valores:

   

Ejercicio Propuesto: Resolver la Ecuación diferencial:

SEGUNDA CLASE DE PROGRAMA

Resolver la Ecuación diferencial

Editor- C:\ Users\ Ecuaciones\ Documents\ MATLAB\ Ecuac. Dif. De 1er orden \Ejercicio2.m %El programa resuelve la ecuación diferencial y´= Cos(5t),para y(0) = 2.0 function Ejercicio2 tspan = 0 : 0.01: 3; y0 =[2]; [T,Y] = ode23tb(@f,tspan,y0); t = 0 : 0.01 : 3; Z = [cos(5.*t)]; [T ,Y,Z.'] plot(T,Z),title('Derivada de la función'),xlabel('Tiempo t'), ylabel('derivada y´'),pause; plot(T,Y(:,1),'-')%Gráfica de y title('Solución de la ecuación diferencial'); %Título del gráfico legend('y');%Leyenda del gráfico xlabel('Tiempo t');%Nombre del eje x ylabel('Función y');%Nombre del eje y hold on function y1 = f(t,y); y1 = [cos(5.*t)];

    RESULTADOS GRÁFICOS

   

Resolver la Ecuación diferencial

Editor- C:\ Users\ Ecuaciones\ Documents\ MATLAB\ Ecuac. Dif. De 1er orden \modelo3.m %El programa resuelve la ecuación diferencial dy/dt= 0.819(100y -y2),para y(0)= 10 function modelo3 tspan = 0 : 0.01: 10; y0 =[10]; [T,Y,Z] = ode23tb(@f,tspan,y0); [T ,Y] [Z.'] plot(T,Y,'-')%Gráfica de y title('Solución de la ecuación diferencial'); %Título del gráfico legend('y');%Leyenda del gráfico xlabel('Tiempo t');%Nombre del eje x ylabel('Función y');%Nombre del eje y hold on function dy = f(t,y) dy = 0.008109.*(100.*y(1) -y(1)^2); Z = dy

    RESULTADOS GRÁFICOS:

Ejercicio Propuesto: Resolver la ecuación diferencial:

   

SIMULACIONES DE ECUACIONES DIFERENCIALES DE SEGUNDO ORDEN EN MATLAB Resolver la ecuación diferencial:

Programa de simulación: Editor- C:\ Users\ Ecuaciones\ Documents\ MATLAB\ Ecuac. Dif. de 2do orden \Ejemplo.m %Ejemplo.m soluciona la ecuac dif d2y/dt2 + 6dy/dt + 9y = 0, %para y(0)=1,y´(0)= 0.5 ,,,,dy(2) = d2y/dt2, dy(1)= dy/dt, function Ejemplo tspan = [0:0.01:3];%Vector de la variable independiente y0 = [1, 0.5]; %Condiciones iniciales [T,Y] = ode45(@f,[0:0.01:3],[1 0.5]); %esta sentencia se puede cambiar ¨ode45(@funcion,tspan,y0) [T,Y]% imprime los valores plot(T,Y(:,1),'-',T,Y(:,2), '-')%Gráfica de y y y´ Esto para que queden de distinto color title('Solución'); %Título del gráfico legend('y','y´');%Leyenda del gráfico xlabel('Tiempo t');%Nombre del eje x ylabel('y y y´');%Nombre del eje y hold on function dy = f(t,y)%Función que depende de y , t dy = zeros(2,1); %Vector columna con el número de ecuaciones a resolve dy(1) = y(2); % Ecuación donde decimos que y(2)es la derivada de y(1) dy(2) = -(6)*y(2)-(9)*y(1);%Ecuación diferencial a resolver RESULTADOS GRÁFICOS

por

   

SOLUCIÓN ANALÍTICA:

   

EJERCICIO PROPUESTO Resolver la ecuación diferencial y comparar el resultado numérico obtenido con el programa con el obtenido analíticamente.

   

EJEMPLOS DE APLICACIONES EN ECUACIONES DIFERENCIALES

ECUACIÓN DIFERENCIAL DE PRIMER ORDEN Elaboramos el script titulado carga _1 para integrar la ecuación diferencial de primer orden que describe la carga de un condensador. V0=10; R=input('Resistencia R: '); C=input('Capacidad C: '); tf=input('tiempo final, tf: '); f=@(t,x) V0/R-x/(R*C); tspan=[0 tf]; x0=0; [t,x]=ode45(f,tspan,x0); plot(t,x,'r') xlabel('t') ylabel('q'); title('carga del condensador')

En la ventana de comandos corremos el script carga_1

>> carga_1 Resistencia R: 2 Capacidad C: 0.8 tiempo final, tf: 10

SISTEMA DE DOS ECUACIONES DIFERENCIALES DE PRIMER ORDEN Elaboramos el script titulado radiactivo_1 para integrar el sistema de dos ecuaciones diferenciales de primer orden que describe la serie de desintegración radioactiva. A-->B-->C donde C es un elemento estable. En la matriz x que devuelve la función ode 45, x(1) representará los sucesivos valores de la variable x y x(2) representará a la variable y. El mismo criterio se empleará para determinar el vector x0 de las condiciones iniciales. La definición de las funciones f (t,x,y) y g(t,x,y) aparecen en un vector columna, separadas por ; (punto y coma). fg=@(t,x) [-a*x(1);a*x(1)-b*x(2)]; % x(1) es x, x(2) es y

    El script será el siguiente: a=input('parámetro a: '); b=input('parámetro b: '); %condiciones iniciales en el vector x0 x0=zeros(1,2); x0(1)=input('valor inicial de x: '); x0(2)=input('valor inicial de y: '); tf=input('tiempo final, tf: '); tspan=[0 tf]; fg=@(t,x) [-a*x(1);a*x(1)-b*x(2)]; [t,x]=ode45(fg,tspan,x0); plot(t,x) xlabel('t') ylabel('x,y'); title('dx/dt=-ax, dy/dt=ax-by') En la ventana de comandos corremos el script radiactivo_1 >> radioactivo_1 parámetro a: 0.1 parámetro b: 0.2 valor inicial de x: 100 valor inicial de y: 0 tiempo final, tf: 20

ECUACIÓN DIFERENCIAL DE SEGUNDO ORDEN Una vez que se ha entendido como resolver un sistema de dos ecuaciones diferenciales de primer orden es posible entender la resolución de cualquier ecuación diferencial o sistema. Podemos definir las funciones de forma anónima o explícitamente en un fichero .M

En este sistema de dos ecuaciones diferenciales de primer orden x(1) se representará los sucesivos valores de la variable x y x(2) representará a la variable v. El mismo criterio se empleará para determinar el vector x0 de las condiciones iniciales. Las funciones a integrar v, y f (t,x,v) aparecen en un vector columna, separadas por ; (punto y coma) f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; % x(1) es x, x(2) es v

    Elaboramos el script oscilador_1 para resolver la ecuación de segundo grado que describe las oscilaciones amortiguadas w0=input('frecuencia angular, w0: '); g=input('rozamiento, gamma: '); %condiciones iniciales x0=zeros(1,2); x0(1)=input('posición inicial, x0: '); x0(2)=input('velocidad inicial, v0: '); tf=input('tiempo final, tf: '); f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; tspan=[0 tf]; [t,x]=ode45(f,tspan,x0); plot(t,x(:,1),'r') grid on xlabel('t') ylabel('x'); title('oscilador amortiguado')

Si en el comando plot ponemos plot(t,x), se representa la posición x(1) y la velocidad x(2) en función del tiempo (en dos colores asignados por MATLAB). Si solamente queremos representar la posición x(1) en función del tiempo t, se escribe plot(t,x(:,1)). En la ventana de comandos corremos el script oscilador_1 >> oscilador_1 frecuencia angular w0: 2 rozamiento, gamma: 0.5 posición inicial, x0: 0 velocidad inicial,v0: 10 tiempo final, tf: 10

Sistema de dos ecuaciones diferenciales de segundo orden En este caso tenemos un sistema de cuatro ecuaciones diferenciales de primer orden

    En este sistema x(1) representará los sucesivos valores de la variable x y x(2) representará a la variable vx, x(3) a la variable y y x(4) a la variable vy . El mismo criterio se empleará para determinar el vector x0 de las condiciones iniciales. Elaboramos el script orbita_1 para resolver el sistema de dos ecuaciones de segundo grado que describe el movimiento de un cuerpo celeste. %condiciones iniciales x0=zeros(1,4); x0(1)=input('posición inicial x: '); x0(2)=input('velocidad incial x: '); x0(3)=0; x0(4)=input('velocidad incial y: '); tf=input('tiempo final, tf: '); tspan=[0 tf]; fg=@(t,x)[x(2);-4*pi*pi*x(1)/(sqrt(x(1)*x(1)+x(3)*x(3)))^3; x(4); -4*pi*pi*x(3)/(sqrt(x(1)*x(1)+x(3)*x(3)))^3]; [t,x]=ode45(fg,tspan,x0); plot(x(:,1),x(:,3),'r') xlabel('x') ylabel('y'); title('Órbita de un planeta') En Figure Window representamos la trayectoria, es decir, los puntos de abscisas x(1) que guardan los valores x y las ordenadas x(3) que guardan los valores y, en función del tiempo t, se escribe plot(t,x(:,1),x(:,3)). En la ventana de comandos corremos el script orbita_1

>> orbita_1 posición inicial x: 1 velocidad incial x: 0 velocidad incial y: 6.27 tiempo final, tf: 1...