Siguiente Anterior

Solución numérica de ecuaciones diferenciales mediante ode45

MATLAB dispone de varias funciones para resolver mediante procedimientos numéricos ecuaciones diferenciales: ode23, ode45, ode113, etc, (véase en el sistema de ayuda para qué tipos de problemas es más adecuado cada uno de los procedimientos). Eligiremos ode45 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)

En la mayor parte de los ejemplos, utilizaremos los tres primeros parámetros: llamaremos a 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.

Vamos a volver a resolver los problemas planteados en este capítulo mediante la función MATLAB ode45.

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.

dq dt = V ε R q RC

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 desintagración radioactiva. A-->B-->C donde C es un elemento estable.

{ dx dt =ax dy dt =axby

En la matriz x que devuelve la función ode45, 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 radiactivo_1 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

Alternativamente, vamos a definir las funciones f (t,x,y) y g(t,x,y) en un fichero .M y le pasamos los valores de los parámetros a y b.

function z=func_radioactivo(t,x,a,b)
    z=[-a*x(1);a*x(1)-b*x(2)]; % x(1) es x, x(2) es y
end

Elaboramos el script radioactivo_2 para establecer los valores de los parámetros a y b, las condiciones iniciales y llamar a la función que realiza la integración numérica ode45. El primer parámetro de ode45 es el handler (manejador de la función) a integrar que se especifica del siguiente modo @nombre_funcion.

[t,x]=ode45(@func_radioactivo,tspan,x0);

Ahora bien, func_radioactivo precisa de los valores de los parámetros a y b. Hay dos formas de hacerlo. La más sencilla es definir una función anónima fg en términos de func_radioactivo. En el problema 3 "Sistemas de ecuaciones de Lorentz" describimos el segundo procedimiento.

fg=@(t,x) func_radioactivo_1(t,x,a,b);

Véase la misma situación en la llamada a función fzero al final de la página Raíces de una ecuación

a=input('parámetro a: ');
b=input('parámetro b: ');
%condiciones iniciales
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) func_radioactivo(t,x,a,b);
[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 radioactivo_2

>> radioactivo_2
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

d x d t = v

d v d t = f ( t , x , v )

En este sistema de dos ecuaciones diferenciales de primer orden x(1) 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.

Como ejemplo, estudiamos las oscilaciones amortiguadas, que hemos descrito en la página anterior.

d 2 x d t 2 +2γ dx dt + ω 0 2 x=0{ dx dt =v dv dt =2γv ω 0 2 x

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)). Véase la página Vectores y matrices

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 difrenciales de primer orden

d x d t = v x

d v x d t = f ( t , x , v x , y , v y )

d y d t = v y

d v y d t = g ( t , x , v x , y , v y )

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.

Como ejemplo, estudiamos el movimiento de un planeta alrededor del Sol o de un satélite artificial alrededor de la Tierra.

d 2 x d t 2 =4 π 2 x ( x 2 + y 2 ) 3/2 { dx dt = v x d v x dt =4 π 2 x ( x 2 + y 2 ) 3/2 d 2 y d t 2 =4 π 2 y ( x 2 + y 2 ) 3/2 { dy dt = v y d v y dt =4 π 2 y ( x 2 + y 2 ) 3/2

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

Opciones de ode45

Imaginemos que estudiamos el movimiento de caída de un cuerpo, no sabemos cuanto tiempo tardará en llegar al suelo, desconocemos el valor del elemento tf en el vector tspan. Sin embargo, queremos detener el proceso de integración numérica de la ecuación diferencial que describe el movimiento cuando la posición del móvil sea cero. La función MATLAB ode45 dispone de un parámetro adicional options donde podemos indicarlo, pero es bastante lioso e intentaremos explicarlo mediante ejemplos.

Volvemos a resolver la ecuación diferencial que describe las oscilaciones amortiguadas y detendremos el proceso de integración cuando el móvil alcance la posición máxima, su velocidad es nula.

Supongamos que el oscilador amortiguado estudiado anteriormente, de frecuencia natural ω0=2, constante de amortiguamiento γ=0.25, parte de la posición x0=2.5 con velocidad nula, queremos detener el proceso de integración cuando el móvil alcance la posición máxima, cuando su velocidad es nula, tal como se muestra en la figura, con lo que se completa un periodo.

Los pasos a seguir son los siguientes:

1.-Definimos la función cuyo nombre es opcion_ode45

function [detect,stopin,direction]=opcion_ode45(t,x)
    detect=[x(1) x(2)]; %[posición, velocidad]
    stopin=[0 1]; %  1 indice que detiene la integración cuando la velocidad se hace cero 
    direction=[0 -1]; % 1 crece, -1 decrece, 0 no importa
end

2.-Creamos la estructura opts con la llamada a la función odeset

opts=odeset('events',@opcion_ode45);

Cuando se utiliza options la función ode45 devuelve los tiempos te en los cuales ocurren los 'eventos' y los correspondientes valores de las variables dependientes (posición, velocidad) en el vector xe. Finalmente, ie es un índice que es útil cuando se vigilan varios eventos.

3.-Le pasamos opts a la función ode45 en su cuarto parámetro

[t,x,te,xe,ie]=ode45(odefunc,tspan,x0,opts);

Escribimos el script oscilador_2 para resolver la ecuación diferencial de segundo orden y detener el proceso de integración de acuerdo con lo estipulado en el parámetro opts.

w0=2;
g=0.25;
x0=[2.5 0];
tf=10;
f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; 
tspan=[0 10];
opts=odeset('events',@opcion_ode45);
[t,x,te,xe,ie]=ode45(f,tspan,x0,opts);
te,xe,ie
plot(t,x(:,1),'r')
grid on
xlabel('t')
ylabel('x');
title('oscilador amortiguado')

Cuando corremos el script oscilador_2 en la ventana de comandos se imprime los siguientes datos relativos a los eventos.

Tiempo, te Posición x(1) Velocidad x(2) Indice ie
0.0000 2.5000 -0.0000 2
2.4378 0.0000 2.7173 1
3.1662 1.1322 -0.0000 2

La columna de tiempos nos porporciona el periodo de la oscilación, te=3.1662.

Se sugiere al lector, cambiar en la función opcion_ode45

direction=[0 1];

Vamos ahora a marcar en la representación gráfica de la oscilación amortiguada, las posiciones de máxima amplitud x(2)=0 y cuando pasa por el origen x(1)=0 sin detener el proceso de integración numérica.

Definimos una nueva versión de la función opcion1_ode45

function [detect,stopin,direction]=opcion1_ode45(t,x)
    detect=[x(1) x(2)]; %[posición, velocidad]
    stopin=[0 0];  
    direction=[0 0]; 
end

Creamos la estructura opts mediante odeset y se la pasamos al procedimiento de integración ode45.

w0=2;
g=0.25;
x0=[2.5 0];
tf=10;
f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; 
tspan=[0 7];
opts=odeset('events',@opcion1_ode45);
[t,x,te,xe,ie]=ode45(f,tspan,x0,opts);

hold on
plot(t,x(:,1),'r')
plot(te(ie==1),xe(ie==1),'o','markersize',6,'markerfacecolor','k')
plot(te(ie==2),xe(ie==2),'o','markersize',6,'markerfacecolor','b')
grid on
xlabel('t')
ylabel('x');
title('oscilador amortiguado')
hold off

Si solamente estamos interesados en los máximos definimos una nueva versión de la función opcion2_ode45

function [detect,stopin,direction]=opcion2_ode45(t,x)
    detect=x(2); 
    stopin=0;
    direction=-1; 
end

Modificamos el script oscilador_4

w0=2;
g=0.25;
x0=[2.5 0];
tf=10;
f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; 
tspan=[0 7];
opts=odeset('events',@opcion1_ode45);
[t,x,te,xe,ie]=ode45(f,tspan,x0,opts);
te,xe,ie
hold on
plot(t,x(:,1),'r')
plot(te(ie==1),xe(ie==1),'o','markersize',6,'markerfacecolor','b')
grid on
xlabel('t')
ylabel('x');
title('oscilador amortiguado')
hold off

Corremos el script oscilador_4 en la ventana de comandos y observamos los resultados

Paso de parámetros a la función

Como hemos visto, a ode45 se le pasa la función (handle) a integrar en su primer argumento. Si la función contiene parámetros como la frecuencia angular ω0, no hay problema si la función se define como anónima, tal como se ha mostrado en los ejemplos previos. Si la función se define en un fichero entonces a la función se le pasan los valores de los parámetros en el quinto argumento params de ode45. Los pasos se explican en el siguiente ejemplo:

El sistema de ecuaciones de Lorentz es un sistema de tres ecuaciones diferenciales de primer orden

dx dt =σx+σy dy dt =ρxyxz dz dt =βz+xy

donde σ=10, β=8/3 y ρ=28

Vamos a resolver el sistema de tres ecuaciones diferenciales con las condiciones iniciales siguientes: en el instante t=0, x0=-8, y0=8 z0=27.

  1. Escribir una función denominada lorentz(t,x,p) como fichero.M que contenga las tres ecuaciones y dependa de los tres parámetros σ=p(1), β=p(2) y ρ=p(3). Observar que la variable x se guarda en el primer elemento x(1), la variable y en el segundo x(2) y la variable z en el tercer x(3) elemento del vector x.
  2. Escribir un script denominado lorentz_script que llame a la función MATLAB ode45, para resolver el sistema de ecuaciones diferenciales para las condiciones iniciales especificadas.
  3. Pasar los parámetros σ, β y ρ como elementos de un vector p a la función ode45 para que a su vez se los pase a la función lorentz.
  4. Dibujar el atractor de Lorentz de z en función de x hasta el instante tf=20 en una primera ventana gráfica.
  5. Dibujar x, y y z en función del tiempo en tres zonas de una segunda ventana gráfica.
  6. Examinar el comportamiento del sistema para otras condiciones iniciales, t=0, x0=1, y0=2 z0=3.

Definimos la función lorentz como fichero .M

function fg=lorentz(t,x,p)    
%x(1) es x, x(2) es y y x(3) es z
% p(1) es sigma, p(2) es beta y p(3) es rho
    fg=[-p(1)*x(1)+p(1)*x(2); p(3)*x(1)-x(2)-x(1)*x(3); -p(2)*x(3)+x(1)*x(2)]; 
end

Elaboramos el script lorentz_script

x0=[-8 8 27]; %valores iniciales
tspan=[0 20]; 
p=[10 8/3 28]; %parámetros
%no pasamos nada [] en el parámetro options de ode45
[t,x]=ode45(@lorentz,tspan,x0,[],p);
figure
plot(x(:,1),x(:,3),'r')
xlabel('x')
ylabel('z');
title('Atractor de Lorentz')

figure
subplot(3,1,1)
plot(t,x(:,1))
ylabel('x');
subplot(3,1,2)
plot(t,x(:,2))
ylabel('y');
subplot(3,1,3)
plot(t,x(:,3))
ylabel('z');
xlabel('t')

En la ventana de comandos corremos el script lorentz_script

>> lorentz_script

Aparecen dos ventanas gráficas, la primera con el atractor de Lorentz, la representación z(x) y la segunda ventana dividida en tres regiones que muestran x(t), y(t) y z(t)

 

Siguiente Anterior