Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

html

Inicio MATLAB Procedimientos numéricos

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.

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 concretos
t0, t1, 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, que explicaremos al final de
esta página ya que es un asunto bastante complicado.
params son parámetros que queremos pasar a la función odefun

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 un script para integrar la ecuación diferencial de primer orden que describe la
carga de un condensador.

1 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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

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

Sistema de dos ecuaciones diferenciales de primer orden

Elaboramos un script 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.

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 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')

2 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

En la ventana de comandos corremos el script

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 un script 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.

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

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

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

3 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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.

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 un script 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

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

4 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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.

Elaboramos un script 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(x(:,1),x(:,3)).

En la ventana de comandos corremos el script

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

5 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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

6 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

function [detect,stopin,direction]=opcion_ode45(t,x)
detect=[x(1) x(2)]; %[posición, velocidad]
%1 indice que detiene la integración cuando la velocidad se hace cero
stopin=[0 1];
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. 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 un script 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 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

Cuando parte de la posición inicial x(1)=2.5 la velocidad es cero x(2)=0, detecta

velocidad (índice ie=2).
Cuando pasa por el origen x(1)=0 detecta posición (índice ie=1), pero no se detiene ya
que en stopin se ha puesto un cero.
Cuando la posición es x(1)=1.1322 detecta velocidad nula x(2)=0, (índice ie=2) y la
integración numérica se detiene ya que en stopin se ha puesto un uno y la velocidad
decrece en direction se ha puesto un -1.

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

7 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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];

8 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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

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

9 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

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

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');

10 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

xlabel('t')

En la ventana de comandos corremos el 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)

Ejemplos en el curso de Física

ode45
Movimiento sobre una superficie semicircular cóncava

Movimiento sobre una superficie cóncava en forma de cicloide

Movimiento sobre una superficie parabólica

Un péndulo simple situado sobre una plataforma móvil

Oscilador de masa variable

Movimiento del extremo de una cadena bajo la acción de un fuerza constante

Estudio del movimiento de una cadena con una máquina de Atwood

Rozamiento en el bucle

Movimiento de dos bloques atados por una cuerda

Movimiento de una partícula atada a una cuerda que se enrolla en un cilindro horizontal.

Solución numérica de las ecuaciones del movimiento

Dos partículas unidas por una cuerda.

El problema de Euler de los tres cuerpos

Viaje de la Tierra a la Luna

Caída de una varilla inclinada

Caída de un lápiz en posición vertical

Movimiento bajo la acción de fuerzas centrales

Simulación de los giros del patinador de hielo

Fuerza de rozamiento proporcional al cuadrado de la velocidad

Descenso de un paracaidista en una atmósfera no uniforme

Una caja que puede volcar

El péndulo simple

Oscilaciones amortiguadas

11 de 12 27/01/2017 04:57 p. m.
Solución numérica de ecuaciones diferenciales mediante <em>ode</em>45 [Link]

Oscilador amortiguado por una fuerza de módulo constante.

El oscilador forzado. El estado estacionario

Vibración de una molécula diató[Link] potencial de Lennard-Jones

El oscilador de “Atwood”

Oscilador no lineal

Oscilaciones forzadas en un sistema formado por partículas unidas por muelles. Fuerza
periódica

Aproximación al equilibrio de dos gases contenidos en un recinto adiabático y separados por

un émbolo (II)

El motor electrostático de Franklin

Movimiento en campos eléctrico y magnético cruzados (II)

Oscilaciones de un cilindro que rueda sobre un plano inclinado con un imán en su interior

Oscilaciones de un imán colgado de un muelle

Caída de un imán

Movimiento de un imán en un tubo metálico vertical

Oscilaciones de un imán colgado de un muelle amortiguadas por una bobina

Osciladores magnéticamente acoplados

Oscilaciones de un imán colgado de un muelle amortiguadas por una placa metálica

Péndulo accionado por las fuerzas de marea

ode15s
Vaciado de un depósito cerrado

Cohete propulsado por agua

Grado en Ingeniería de Energías Renovables

Angel Franco García, Copyright © 2016

12 de 12 27/01/2017 04:57 p. m.