ode15s

Resolver ecuaciones diferenciales y DAE rígidas; método de orden variable

Sintaxis

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

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

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

sol = ode15s(___)

Descripción

[t,y] = ode15s(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.

[t,y] = ode15s(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] = ode15s(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 = ode15s(___) 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 }=-10t.$$

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

tspan = [0 2];
y0 = 1;
[t,y] = ode15s(@(t,y) -10*t, tspan, y0);

Represente la solución.

plot(t,y,'-o')

Resolver ODE rígidas

Un ejemplo de un sistema rígido de ecuaciones son las ecuaciones de van der Pol en oscilación de relajación. El ciclo límite tiene regiones en las que los componentes de la solución cambian lentamente y el problema es bastante rígido, alternadas con regiones de cambios muy bruscos donde no es rígido.

El sistema de ecuaciones es:

Las condiciones iniciales son y . La función vdp1000 está incluida en MATLAB® y codifica las ecuaciones.

function dydt = vdp1000(t,y) %#ok<INUSD>
%VDP1000  Evaluate the van der Pol ODEs for mu = 1000.
%
%   See also ODE15S, ODE23S, ODE23T, ODE23TB.

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

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

Resolver este sistema usando ode45 con las tolerancias a errores relativas y absolutas (1e-3 y 1e-6, respectivamente) es extremadamente lento, por lo que se requiere varios minutos para resolver y representar la solución. ode45 requiere millones de unidades de tiempo para completar la integración, debido a las áreas de rigidez donde tiene dificultades para cumplir las tolerancias.

Esta es una gráfica de la solución obtenida por ode45, que tarda mucho tiempo en calcularse. Observe la enorme cantidad de unidades de tiempo necesarias para pasar por áreas de rigidez.

Resuelva el sistema rígido usando el solver ode15s y, después, represente la primera columna de la solución y con respecto a los puntos de tiempo t. El solver ode15s pasa por áreas rígidas con muchos menos pasos que ode45.

[t,y] = ode15s(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

Pasar parámetros extra a una función ODE

ode15s 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

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

odefcn.m 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);

Resuelva la ODE utilizando ode15s. 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] = ode15s(@(t,y) odefcn(t,y,A,B), tspan, y0);

Represente los resultados.

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

Comparar solvers ODE rígidos

El solver ode15s es una buena primera elección para la mayoría de problemas rígidos. Sin embargo, es posible que el resto de solvers rígidos sea más eficiente para determinados tipos de problemas. Este ejemplo resuelve una ecuación de prueba rígida usando los cuatro solvers ODE rígidos.

Considere la ecuación de prueba

$$y^{\prime }=-\lambda y.$$

La ecuación se vuelve cada vez más rígida a medida que la magnitud de $$\lambda$$ aumenta. Utilice $$\lambda =1\times 10^9$$ y la condición inicial $$y(0)=1$$ en el intervalo de tiempo [0 0.5]. Con estos valores, el problema es lo suficientemente rígido como para que ode45 y ode23 tengan dificultades para integrar la ecuación. Además, use odeset para pasar la jacobiana constante $$J=\frac{\partial f}{\partial y}=-\lambda$$ y activar la visualización de las estadísticas del solver.

lambda = 1e9;
y0 = 1;
tspan = [0 0.5];
opts = odeset('Jacobian',-lambda,'Stats','on');

Resuelva la ecuación con ode15s, ode23s, ode23t y ode23tb. Haga subgráficas para compararlas.

subplot(2,2,1)
tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc

104 successful steps
1 failed attempts
212 function evaluations
0 partial derivatives
21 LU decompositions
210 solutions of linear systems
Elapsed time is 1.499378 seconds.

title('ode15s')
subplot(2,2,2)
tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc

63 successful steps
0 failed attempts
191 function evaluations
0 partial derivatives
63 LU decompositions
189 solutions of linear systems
Elapsed time is 0.460420 seconds.

title('ode23s')
subplot(2,2,3)
tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc

95 successful steps
0 failed attempts
125 function evaluations
0 partial derivatives
28 LU decompositions
123 solutions of linear systems
Elapsed time is 0.237496 seconds.

title('ode23t')
subplot(2,2,4)
tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc

71 successful steps
0 failed attempts
167 function evaluations
0 partial derivatives
23 LU decompositions
236 solutions of linear systems
Elapsed time is 0.231518 seconds.

title('ode23tb')

Todos los solvers rígidos funcionan bien, pero ode23s completa la integración con el menor número de pasos y se ejecuta de la forma más rápida para este problema concreto. Dado que se especifica la jacobiana constante, ninguno de los solvers necesita calcular derivadas parciales para calcular la solución. Especificar la jacobiana beneficia principalmente a ode23s, dado que normalmente evalúa la jacobiana en cada paso.

Para problemas de rigidez generales, el rendimiento de los solvers rígidos varía en función del formato del problema y las opciones especificadas. Proporcionar la matriz jacobiana o el patrón de dispersión siempre mejora la eficiencia del solver para problemas rígidos. Pero, dado que los solvers rígidos utilizan la jacobiana de manera diferente, la mejora puede variar de manera significativa. En la práctica, si un sistema de ecuaciones es muy grande o necesita resolverse muchas veces, vale la pena investigar el rendimiento de los diferentes solvers para minimizar el tiempo de ejecución.

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 \left(1-y_1^2\right)y_1^{\prime }+y_1=0.$$

Resuelva la ecuación de van der Pol con $$\mu =1000$$ utilizando ode15s. La función vdp1000.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 3000];
y0 = [2 0];
sol = ode15s(@vdp1000,tspan,y0)

sol = struct with fields:
     solver: 'ode15s'
    extdata: [1×1 struct]
          x: [0 1.4606e-05 2.9212e-05 4.3818e-05 1.1010e-04 1.7639e-04 2.4267e-04 3.0896e-04 4.5006e-04 5.9116e-04 7.3226e-04 8.7336e-04 0.0010 0.0012 0.0013 0.0015 0.0017 0.0018 0.0021 0.0024 0.0027 0.0030 0.0033 0.0044 0.0055 0.0066 … ] (1×592 double)
          y: [2×592 double]
      stats: [1×1 struct]
      idata: [1×1 struct]

Utilice linspace para generar 2500 puntos en el intervalo [0 3000]. Evalúe el primer componente de la solución en estos puntos utilizando deval.

x = linspace(0,3000,2500);
y = deval(sol,x,1);

Represente la solución.

plot(x,y)

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

tf = 4000;
sol_new = odextend(sol,@vdp1000,tf);
x = linspace(3000,tf,350);
y = deval(sol_new,x,1);
hold on
plot(x,y,'r')

Resolver un problema de Robertson como ecuaciones algebraicas diferenciales (DAE) semiexplícitas

Este ejemplo reformula un sistema de ODE como un sistema de ecuaciones algebraicas diferenciales (DAE). El problema de Robertson que se encuentra en hb1ode.m es un problema de prueba clásico para programas que resuelven ODE rígidas. El sistema de ecuaciones es

hb1ode resuelve este sistema de ODE hasta el estado estacionario con las condiciones iniciales , y . Pero las ecuaciones también cumplen una ley de conservación lineal,

En términos de la solución y las condiciones iniciales, la ley de conservación es

El sistema de ecuaciones se puede reescribir como un sistema de DAE usando la ley de conservación para determinar el estado de . Esto reformula el problema como el sistema de DAE

El índice diferencial de este sistema es 1, dado que solo se requiere una única derivada de para convertir esto en un sistema de ODE. Por lo tanto, no se requieren más transformaciones antes de resolver el sistema.

La función robertsdae codifica este sistema de DAE. Guarde robertsdae.m en la carpeta actual para ejecutar el ejemplo.

function out = robertsdae(t,y)
out = [-0.04*y(1) + 1e4*y(2).*y(3)
   0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2
   y(1) + y(2) + y(3) - 1 ];

El código de ejemplo completo para esta formulación del problema de Robertson está disponible en hb1dae.m.

Resuelva el sistema de DAE usando ode15s. Las condiciones iniciales consistentes para y0 son obvias de acuerdo con la ley de conservación. Utilice odeset para establecer las opciones:

• Utilice una matriz de masa constante para representar el lado izquierdo del sistema de ecuaciones.

• Establezca la tolerancia a errores relativa en 1e-4.

• Utilice una tolerancia a errores absoluta de 1e-10 para el segundo componente de la solución, dado que la escala varía drásticamente del resto de componentes.

• Deje la opción 'MassSingular' en su valor predeterminado 'maybe' para probar la detección automática de un DAE.

y0 = [1; 0; 0];
tspan = [0 4*logspace(-6,6)];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);
[t,y] = ode15s(@robertsdae,tspan,y0,options);

Represente la solución.

y(:,2) = 1e4*y(:,2);
semilogx(t,y);
ylabel('1e4 * y(:,2)');
title('Robertson DAE problem with a Conservation Law, solved by ODE15S');

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

ode15s es un solver de paso variable y orden variable (VSVO) basado en fórmulas de diferenciación numérica (NDF) de órdenes 1 a 5. Opcionalmente, puede usar las fórmulas de diferenciación regresiva (BDF, también conocidas como el método de Gear) que son generalmente menos eficientes. Al igual que ode113, ode15s es un solver de pasos múltiples. Utilice ode15s si ode45 falla o es muy ineficiente y sospecha que el problema es rígido, o cuando resuelva una ecuación diferencial algebraica (DAE) [1], [2].

Referencias

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

[2] Shampine, L. F., M. W. Reichelt, and J.A. Kierzenka, “Solving Index-1 DAEs in MATLAB and Simulink,” SIAM Review, Vol. 41, 1999, pp. 538–552.

Capacidades ampliadas

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

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

Notas y limitaciones de uso:

• La compatibilidad con los tamaños variables y la asignación de memoria dinámica deben estar activadas.

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

• Debe proporcionar puntos de evaluación y soluciones como argumentos de salida.

Generación de código de GPU

Genere código CUDA® para GPU NVIDIA® mediante GPU Coder™.

Consulte las notas y limitaciones de uso en la sección de generación de código C/C++. Las mismas notas y limitaciones de uso son aplicables a la generación de código GPU.

Historial de versiones

Introducido antes de R2006a

Consulte también

ode45 | ode23s | ode23t | ode23tb | odeset | odeget | deval

Temas

• Elegir un solver de ODE
• Resumen de opciones de ODE
• Solve Stiff ODEs
• Solve Differential Algebraic Equations (DAEs)
• Funciones anónimas
• Troubleshoot Common ODE Problems