ode45

Resolver ecuaciones diferenciales no rígidas; método de orden intermedio

Sintaxis

[t,y] = ode45(odefun,tspan,y0)

[t,y] = ode45(odefun,tspan,y0,options)

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)

sol = ode45(___)

Descripción

[t,y] = ode45(odefun,tspan,y0), donde tspan = [t0 tf], integra el sistema de ecuaciones diferenciales $$y\text{'}=f\left(t,y\right)$$ de t0 a tf con condiciones iniciales y0. Cada fila del arreglo de solución y se corresponde con un valor devuelto en el vector columna t.

Todos los solvers de ODE de MATLAB^® pueden resolver sistemas de ecuaciones con el formato $$y\text{'}=f\left(t,y\right)$$, o problemas que incluyen una matriz de masa, $$M\left(t,y\right)y\text{'}=f\left(t,y\right)$$. Todos los solvers utilizan sintaxis similares. El solver ode23s solo puede resolver problemas con una matriz de masa si la matriz de masa es constante. ode15s y ode23t pueden resolver problemas con una matriz de masa que es singular, conocidos como ecuaciones algebraicas diferenciales (DAE). Especifique la matriz de masa utilizando la opción Mass de odeset.

ode45 es un versátil solver de ODE y es el primer solver que debería probar para solucionar la mayoría de problemas. Sin embargo, si el problema es rígido o requiere mucha precisión, existen otros solvers de ODE que podrían ser más adecuados para resolverlo. Para obtener más información, consulte Elegir un solver de ODE.

[t,y] = ode45(odefun,tspan,y0,options) también utiliza la configuración de integración definida por options, que es un argumento creado utilizando la función odeset. Por ejemplo, utilice las opciones AbsTol y RelTol para especificar tolerancias a errores absolutas y relativas, o la opción Mass para proporcionar una matriz de masa.

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) también encuentra dónde son cero las funciones de (t,y), denominadas funciones de evento. En la salida, te es el tiempo del evento, ye es la solución en el momento del evento e ie es el índice del evento activado.

Para cada función de evento, especifique si la integración debe terminar en un cero y si la dirección del cruce por cero tiene importancia. Hágalo estableciendo la propiedad 'Events' en una función, como myEventFcn o @myEventFcn, y creando una función correspondiente: [value,isterminal,direction] = myEventFcn(t,y). Para obtener más información, consulte Ubicación de eventos de EDO.

sol = ode45(___) devuelve una estructura que puede usar con deval para evaluar la solución en cualquier punto del intervalo [t0 tf]. Puede utilizar cualquiera de las combinaciones de argumentos de entrada de las sintaxis anteriores.

Ejemplos

ODE con un solo componente de solución

Las ODE sencillas que tienen un solo componente de solución pueden especificarse como función anónima en la llamada al solver. La función anónima debe aceptar dos entradas (t,y), aunque una de las entradas no se utilice en la función.

Resuelva la ODE

$$y^{\prime }=2t.$$

Especifique un intervalo de tiempo de [0 5] y la condición inicial y0 = 0.

tspan = [0 5];
y0 = 0;
[t,y] = ode45(@(t,y) 2*t, tspan, y0);

Represente la solución.

plot(t,y,'-o')

Resolver ecuación no rígida

La ecuación de van der Pol es una ODE de segundo orden

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0,$$

donde $$\mu >0$$ es un parámetro de escalar. Vuelva a escribir esta ecuación como un sistema de ODE de primer orden haciendo la sustitución $$y_1^{\prime }=y_2$$. El sistema resultante de las ODE de primer orden es

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\mu (1-y_1^2)y_2-y_1.\end{array}$$

El archivo de función vdp1.m representa la ecuación de van der Pol utilizando $$\mu =1$$. Las variables $$y_1$$ y $$y_2$$ son las entradas y(1) e y(2) de un vector de dos elementos dydt.

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

Resuelva la ODE utilizando la función ode45 en el intervalo de tiempo [0 20] con los valores iniciales [2 0]. La salida resultante es un vector columna de puntos de tiempo t y un arreglo de solución y. Cada fila en y se corresponde con un tiempo devuelto en la fila correspondiente de t. La primera columna de y se corresponde con $$y_1$$ y la segunda se corresponde con $$y_2$$.

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

Represente las soluciones para $$y_1$$ y $$y_2$$ con respecto a t.

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

Pasar parámetros extra a una función ODE

ode45 funciona solo con funciones que utilizan dos argumentos de entrada, t e y. Sin embargo, puede pasar parámetros extra definiéndolos fuera de la función y pasándolos cuando especifique el identificador de función.

Resuelva la ODE

$$y^{\prime \prime }=\frac{A}{B}ty.$$

Si vuelve a escribir la ecuación como un sistema de primer orden, obtendrá

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\frac{A}{B}t{y}_1.\end{array}$$

odefcn, una función local incluida al final de este ejemplo, representa este sistema de ecuaciones como una función que acepta cuatro argumentos de entrada: t, y, A y B.

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

Resuelva la ODE utilizando ode45. Especifique el identificador de función para que pase los valores predefinidos para A y B a odefcn.

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);

Represente los resultados.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

Resolver ODE con varias condiciones iniciales

Para sistemas de ODE sencillos con una ecuación, puede especificar y0 como un vector que contenga varias condiciones iniciales. Esta técnica crea un sistema de ecuaciones independientes a través de una expansión escalar, una para cada valor inicial, y ode45 resuelve el sistema para generar resultados para cada valor inicial.

Cree una función anónima para representar la ecuación $\mathit{f}\left(\mathit{t},\mathit{y}\right)=-2\mathit{y}+2\cos \left(\mathit{t}\right)\sin \left(2\mathit{t}\right)$. La función debe aceptar dos entradas para t e y.

yprime = @(t,y) -2*y + 2*cos(t).*sin(2*t);

Cree un vector de condiciones iniciales distintas en el rango $\left[-5,5\right]$.

y0 = -5:5; 

Resuelva la ecuación para cada condición inicial en el intervalo de tiempo $\left[0,3\right]$ utilizando ode45.

tspan = [0 3];
[t,y] = ode45(yprime,tspan,y0);

Represente los resultados.

plot(t,y)
grid on
xlabel('t')
ylabel('y')
title('Solutions of y'' = -2y + 2 cos(t) sin(2t), y(0) = -5,-4,...,4,5','interpreter','latex')

Esta técnica es útil para resolver ODE sencillas con varias condiciones iniciales. Sin embargo, la técnica también tiene algunos inconvenientes:

• No puede resolver sistemas de ecuaciones con varias condiciones iniciales. La técnica solo funciona cuando se resuelve una ecuación con varias condiciones iniciales.

• La unidad de tiempo elegida por el solver en cada salto se basa en la ecuación del sistema que necesita realizar los saltos más cortos. Esto significa que el solver puede realizar saltos cortos para satisfacer la ecuación para una condición inicial, pero las otras ecuaciones, si se resolviesen de forma independiente, utilizarían distintos tamaños de salto. A pesar de esto, en general es más rápido resolver para varias condiciones iniciales a la vez que resolver las ecuaciones de forma independiente utilizando un bucle for.

Para obtener más información sobre esta técnica, consulte Solve System of ODEs with Multiple Initial Conditions.

ODE con términos que dependen del tiempo

Considere la siguiente ODE con parámetros que dependen del tiempo

La condición inicial es . La función f(t) la define el vector de n por 1 f evaluado en los momentos ft. La función g(t) la define el vector de m por 1 g evaluado en los momentos gt.

Cree los vectores f y g.

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

Escriba una función llamada myode que interpole f y g para obtener el valor de los términos que dependen del tiempo en el momento especificado. Guarde la función en la carpeta actual para ejecutar el resto del ejemplo.

La función myode acepta argumentos de entrada extra para evaluar la ODE en cada unidad de tiempo, pero ode45 solo utiliza los primeros dos argumentos de entrada t e y.

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t

Resuelva la ecuación en el intervalo de tiempo [1 5] utilizando ode45. Especifique la función utilizando un identificador de función de forma que ode45 utilice los dos primeros argumentos de entrada de myode. Además, relaje los umbrales de error utilizando odeset.

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

Represente la solución, y, como una función de los puntos de tiempo, t.

plot(t,y)

Evaluar y ampliar la estructura de una solución

La ecuación de van der Pol es una ODE de segundo orden

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0.$$

Resuelva la ecuación de van der Pol con $$\mu =1$$ utilizando ode45. La función vdp1.m está incluida en MATLAB® y codifica las ecuaciones. Especifique una sola salida para devolver una estructura que contiene información sobre la solución, como el solver y los puntos de evaluación.

tspan = [0 20];
y0 = [2 0];
sol = ode45(@vdp1,tspan,y0)

sol = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [0 1.0048e-04 6.0285e-04 0.0031 0.0157 0.0785 0.2844 0.5407 0.8788 1.4032 1.8905 2.3778 2.7795 3.1285 3.4093 3.6657 3.9275 4.2944 4.9013 5.3506 5.7998 6.2075 6.5387 6.7519 6.9652 7.2247 7.5719 8.1226 8.6122 9.1017 9.5054 ... ] (1x60 double)
          y: [2x60 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

Utilice linspace para generar 250 puntos en el intervalo [0 20]. Evalúe la solución en estos puntos utilizando deval.

x = linspace(0,20,250);
y = deval(sol,x);

Represente el primer componente de la solución.

plot(x,y(1,:))

Amplíe la solución a $$t_f=35$$ utilizando odextend y añada el resultado a la gráfica original.

sol_new = odextend(sol,@vdp1,35);
x = linspace(20,35,350);
y = deval(sol_new,x);
hold on
plot(x,y(1,:),'r')

Argumentos de entrada

odefun — Funciones que resolver

identificador de función

Funciones que resolver, especificadas como identificador de función que define las funciones que se desea integrar.

La función dydt = odefun(t,y), para un escalar t y un vector columna y, debe devolver un vector columna dydt de tipo de datos single o double que se corresponda con $$f\left(t,y\right)$$. odefun debe aceptar tanto argumentos de entrada t como y, aunque uno de los argumentos no se utilice en la función.

Por ejemplo, para resolver $$y\text{'}=5y-3$$, utilice la función:

function dydt = odefun(t,y)
dydt = 5*y-3;
end

Para un sistema de ecuaciones, la salida de odefun es un vector. Cada elemento del vector es el valor calculado del lado derecho de una ecuación. Por ejemplo, considere el sistema de dos ecuaciones

$$\begin{array}{l}y{\text{'}}_1=y_1+2y_2\\ y{\text{'}}_2=3y_1+2y_2\end{array}$$

Una función que calcula el valor del lado derecho de cada ecuación en cada unidad de tiempo es

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);
end

Para obtener información sobre cómo proporcionar parámetros adicionales a la función odefun, consulte Parametrizar funciones.

Ejemplo: @myFcn

Tipos de datos: function_handle

tspan — Intervalo de integración

vector

Intervalo de integración, especificado como vector. Como mínimo, tspan debe ser un vector de dos elementos [t0 tf] que especifique el momento inicial y final. Para obtener soluciones en momentos específicos entre t0 y tf, utilice un vector más largo con el formato [t0,t1,t2,...,tf]. Los elementos de tspan deben ser todos crecientes o todos decrecientes.

El solver impone las condiciones iniciales proporcionadas por y0 en el momento inicial tspan(1) y, después, integra desde tspan(1) hasta tspan(end):

• Si tspan tiene dos elementos [t0 tf], el solver devuelve la solución evaluada en cada salto de integración interno dentro del intervalo.

• Si tspan tiene más de dos elementos [t0,t1,t2,...,tf], el solver devuelve la solución evaluada en los puntos proporcionados. Sin embargo, el solver no salta con precisión a cada punto especificado en tspan. En su lugar, el solver utiliza sus propios saltos internos para calcular la solución y, a continuación, la evalúa en los puntos solicitados en tspan. Las soluciones generadas en los puntos especificados son del mismo orden de precisión que las soluciones calculadas en cada salto interno.

  El hecho de especificar varios puntos intermedios no afecta mucho a la eficiencia del cálculo, pero puede afectar a la administración de la memoria en el caso de sistemas grandes.

El solver utiliza los valores de tspan para calcular valores adecuados para InitialStep y MaxStep:

• Si tspan contiene varios puntos intermedios [t0,t1,t2,...,tf], los puntos especificados dan una noción de la escala para el problema, que puede afectar al valor de InitialStep utilizado por el solver. Por lo tanto, es posible que la solución obtenida por el solver sea distinta en función de si especifica tspan como un vector de dos elementos o como un vector con puntos intermedios.

• Los valores iniciales y finales en tspan se utilizan para calcular el tamaño máximo del salto MaxStep. Por lo tanto, si cambia los valores iniciales o finales en tspan, es posible que el solver utilice una secuencia de saltos distinta, lo que podría cambiar la solución.

Ejemplo: [1 10]

Ejemplo: [1 3 5 7 9 10]

Tipos de datos: single | double

y0 — Condiciones iniciales

vector

Condiciones iniciales, especificadas como vector. y0 debe tener la misma longitud que el vector de salida de odefun, de forma que y0 contenga una condición inicial para cada ecuación definida en odefun.

Tipos de datos: single | double

options — Estructura de opciones

arreglo de estructura

Estructura de opciones, especificada como arreglo de estructuras. Use la función odeset para crear o modificar la estructura de options. Consulte Resumen de opciones de ODE para obtener una lista de las opciones compatibles con cada solver.

Ejemplo: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) especifica una tolerancia a errores relativa de 1e-5, activa la visualización de las estadísticas del solver y especifica la función de salida @odeplot para representar la solución mientras se calcula.

Tipos de datos: struct

Argumentos de salida

t — Puntos de evaluación

vector columna

Puntos de evaluación, devueltos como vector columna.

• Si tspan contiene dos elementos [t0 tf], t contiene los puntos de evaluación internos utilizados para realizar la integración.

• Si tspan contiene más de dos elementos, t es lo mismo que tspan.

y — Soluciones

arreglo

Soluciones, devueltas como arreglo. Cada fila en y se corresponde con la solución en el valor devuelto en la fila correspondiente de t.

te — Momento de eventos

vector columna

Momento de eventos, devuelto como vector columna. Los momentos de evento en te se corresponden con las soluciones devueltas en ye, e ie especifica qué evento ha ocurrido.

ye — Solución en momento de eventos

arreglo

Solución en momento de eventos, devuelta como arreglo. Los momentos de evento en te se corresponden con las soluciones devueltas en ye, e ie especifica qué evento ha ocurrido.

ie — Índice de función de evento activado

vector columna

Índice de función de evento activado, devuelto como vector columna. Los momentos de evento en te se corresponden con las soluciones devueltas en ye, e ie especifica qué evento ha ocurrido.

sol — Estructura para evaluación

arreglo de estructura

Estructura de evaluación, especificada como arreglo de estructuras. Utilice esta estructura con la función deval para evaluar la solución en cualquier punto del intervalo [t0 tf]. El arreglo de estructuras sol siempre incluye estos campos:

Campo de estructura	Descripción
sol.x	Vector fila de los saltos elegidos por el solver.
sol.y	Soluciones. Cada columna sol.y(:,i) contiene la solución en el momento sol.x(i).
sol.solver	Nombre del solver.

Además, si especifica la opción Events de odeset y se detectan eventos, sol también incluye estos campos:

Campo de estructura	Descripción
sol.xe	Puntos donde ocurrieron los eventos. sol.xe(end) contiene el punto exacto de un evento terminal, en su caso.
sol.ye	Soluciones que se corresponden con eventos en sol.xe.
sol.ie	Índices en el vector devueltos por la función especificada en la opción Events. Los valores indican qué evento detectó el solver.

Algoritmos

ode45 se basa en una fórmula explícita de Runge-Kutta (4,5), el par de Dormand-Prince. Se trata de un solver de salto único; cuando se calcula y(tn), solo necesita la solución en el punto de tiempo inmediatamente anterior, y(tn-1) [1], [2].

Referencias

[1] Dormand, J. R. and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comp. Appl. Math., Vol. 6, 1980, pp. 19–26.

[2] Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.

Capacidades ampliadas

Generación de código C/C++

Genere código C y C++ mediante MATLAB® Coder™.

Notas y limitaciones de uso:

• Todos los argumentos de la opción odeset deben ser constantes.

• La generación de código no es compatible con una matriz de masa constante en la estructura de opciones. Proporcione una matriz de masa como función.

• Debe proporcionar al menos los dos argumentos de salida T e Y.

• Los tipos de entrada deben ser homogéneos; todos dobles o todos simples.

• Debe estar activada la compatibilidad con los tamaños variables. La generación de código requiere la asignación de memoria dinámica cuando tspan tiene dos elementos o utiliza funciones de evento.

Historial de versiones

Introducido antes de R2006a

Consulte también

ode23 | ode78 | ode89 | ode113 | ode15s | odeset | odeget | deval | odextend

Temas

• Elegir un solver de ODE
• Resumen de opciones de ODE
• Resolver ecuaciones diferenciales ordinarias (EDO) no rígidas
• Funciones anónimas
• Troubleshoot Common ODE Problems