Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Resolver ODE rígidas

Esta página contiene dos ejemplos de resolución de ecuaciones diferenciales ordinarias con ode15s. MATLAB® tiene cuatro solvers diseñados para ODE rígidas.

  • ode15s

  • ode23s

  • ode23t

  • ode23tb

Para la mayoría de los problemas rígidos, ode15s funciona mejor. Sin embargo, ode23s, ode23t y ode23tb pueden resultar más eficientes si el problema permite una tolerancia a errores cruda.

¿Qué es una ODE rígida?

Para algunos problemas de ODE, el tamaño del paso dado por el solver se fuerza a un nivel injustificadamente pequeño en comparación con el intervalo de integración, incluso en una región donde la curva de solución es suave. Estos tamaños de los pasos pueden ser tan pequeños que recorrer un breve intervalo de tiempo puede requerir millones de evaluaciones. Esto puede hacer que el solver falle en la integración, pero incluso si tiene éxito tardará mucho en hacerlo.

A las ecuaciones que provocan este comportamiento en los solvers de ODE se las denomina rígidas. El problema que plantean las ODE rígidas es que los solvers explícitos (como ode45) resultan insosteniblemente lentos para lograr una solución. Esta es la razón por la que ode45 se clasifica como un solver no rígido junto con ode23 y ode113.

Los solvers diseñados para ODE rígidas, conocidos como solvers rígidos, suelen hacer más trabajo por cada paso. La ventaja es que pueden dar pasos mucho mayores, y tienen una estabilidad numérica mejorada en comparación con los solvers no rígidos.

Opciones de solver

Para problemas rígidos, es especialmente importante especificar la matriz jacobina con odeset. Los solvers rígidos utilizan la matriz jacobina $\partial f_i / \partial y_j$ para estimar el comportamiento local de las ODE a medida que avanza la integración, por lo que proporcionar la matriz jacobina (o, en el caso de grandes sistemas dispersos, el patrón de dispersión) es fundamental para la eficiencia y la fiabilidad. Utilice las opciones Jacobian, JPattern o Vectorized de odeset para especificar información sobre la jacobina. Si no proporciona la jacobina, entonces el solver la estima numéricamente con diferencias finitas.

Consulte odeset para obtener un listado completo de las opciones del solver.

Ejemplo: ecuación de van der Pol rígida

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

$$y''_1 - \mu \left( 1 - y_1^2\right) y'_1+y_1=0,$$

donde $\mu > 0$ es un parámetro escalar. Cuando $\mu = 1$, el sistema resultante de ODE es no rígido y se resuelve fácilmente con ode45. Sin embargo, si aumenta $\mu$ a 1000, entonces la solución cambia drásticamente y muestra oscilación en una escala de tiempo mucho mayor. Aproximarse a la solución del valor inicial se vuelve más difícil. Dado que este problema particular es rígido, un solver pensado para problemas no rígidos, como ode45, resulta demasiado ineficiente para ser práctico. En su lugar, utilice un solver rígido como ode15s para este problema.

Reescriba la ecuación van del Pol como un sistema de ODE de primer orden realizando la sustitución $y'_1 = y_2$. El sistema resultante de ODE de primer orden es

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

La función vdp1000 evalúa la ecuación de van der Pol con $\mu = 1000$.

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

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

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

Utilice la función ode15s para resolver el problema con un vector inicial de condiciones de [2; 0], durante un intervalo de tiempo de [0 3000]. Por razones de escala, represente solo el primer componente de la solución.

[t,y] = ode15s(@vdp1000,[0 3000],[2; 0]);
plot(t,y(:,1),'-o');
title('Solution of van der Pol Equation, \mu = 1000');
xlabel('Time t');
ylabel('Solution y_1');

La función vdpode también resuelve el mismo problema, pero acepta un valor especificado por el usuario para $\mu$. Las ecuaciones se vuelven más rígidas a medida que $\mu$ aumenta.

Ejemplo: sistema de Brusselator disperso

El clásico sistema de ecuaciones Brusselator es potencialmente grande, rígido y disperso. El sistema de Brusselator modela la difusión en una reacción química, y se representa con un sistema de ecuaciones que implica $u$, $v$, $u'$ y $v'$.

$$ \begin{array}{cl} u'_i &= 1+u_i^2v_i-4u_i+ \alpha \left( N + 1 \right)
^2 \left( u_{i-1}-2_i+u_{i+1} \right)\\ v'_i &= 3u_i-u_i^2v_i + \alpha
\left( N+1 \right) ^2 \left( v_{i-1} - 2v_i+v_{i+1} \right) \end{array}$$

El archivo de función brussode resuelve este conjunto de ecuaciones en el intervalo de tiempo [0,10] con $\alpha = 1/50$. Las condiciones iniciales son

$$\begin{array}{cl} u_j(0) &= 1+\sin(2 \pi x_j)\\ v_j(0) &=
3,\end{array}$$

donde $x_j = i/N+1$ para $i=1,...,N$. Por lo tanto, existen $2N$ ecuaciones en el sistema, pero la jacobina $\partial f / \partial y$ es una matriz de bandas con un ancho constante de 5 si las ecuaciones se ordenan como $u_1,v_1,u_2,v_2,...$. A medida que $N$ aumenta, el problema se vuelve más rígido y la jacobina se vuelve más dispersa.

La llamada de la función brussode(N), para $N \ge 2$, especifica un valor para N en el sistema de ecuaciones, que se corresponde con el número de puntos de la cuadrícula. De forma predeterminada, brussode utiliza $N = 20$.

brussode contiene algunas subfunciones:

  • La función anidada f(t,y) codifica el sistema de ecuaciones para el problema de Brusselator y devuelve un vector.

  • La función local jpattern(N) devuelve una matriz dispersa de unos y ceros que muestra las ubicaciones de los no ceros en la jacobina. Esta matriz se asigna al campo JPattern de la estructura de opciones. El solver de ODE utiliza este patrón de dispersión para generar la jacobina numéricamente como una matriz dispersa. Proporcionar este patrón de dispersión en el problema reduce de forma considerable el número de evaluaciones de función necesarias para generar la jacobina 2N por 2N, de 2N evaluaciones a solo 4.

function brussode(N)
%BRUSSODE  Stiff problem modelling a chemical reaction (the Brusselator).
%   The parameter N >= 2 is used to specify the number of grid points; the
%   resulting system consists of 2N equations. By default, N is 20.  The
%   problem becomes increasingly stiff and increasingly sparse as N is
%   increased.  The Jacobian for this problem is a sparse constant matrix
%   (banded with bandwidth 5).
%
%   The property 'JPattern' is used to provide the solver with a sparse
%   matrix of 1's and 0's showing the locations of nonzeros in the Jacobian
%   df/dy.  By default, the stiff solvers of the ODE Suite generate Jacobians
%   numerically as full matrices.  However, when a sparsity pattern is
%   provided, the solver uses it to generate the Jacobian numerically as a
%   sparse matrix.  Providing a sparsity pattern can significantly reduce the
%   number of function evaluations required to generate the Jacobian and can
%   accelerate integration.  For the BRUSSODE problem, only 4 evaluations of
%   the function are needed to compute the 2N x 2N Jacobian matrix.
%
%   Setting the 'Vectorized' property indicates the function f is
%   vectorized.
%
%   E. Hairer and G. Wanner, Solving Ordinary Differential Equations II,
%   Stiff and Differential-Algebraic Problems, Springer-Verlag, Berlin,
%   1991, pp. 5-8.
%
%   See also ODE15S, ODE23S, ODE23T, ODE23TB, ODESET, FUNCTION_HANDLE.

%   Mark W. Reichelt and Lawrence F. Shampine, 8-30-94
%   Copyright 1984-2014 The MathWorks, Inc.

% Problem parameter, shared with the nested function.
if nargin<1
   N = 20;
end

tspan = [0; 10];
y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)];

options = odeset('Vectorized','on','JPattern',jpattern(N));

[t,y] = ode15s(@f,tspan,y0,options);

u = y(:,1:2:end);
x = (1:N)/(N+1);
figure;
surf(x,t,u);
view(-40,30);
xlabel('space');
ylabel('time');
zlabel('solution u');
title(['The Brusselator for N = ' num2str(N)]);

% -------------------------------------------------------------------------
% Nested function -- N is provided by the outer function.
%

   function dydt = f(t,y)
      % Derivative function
      c = 0.02 * (N+1)^2;
      dydt = zeros(2*N,size(y,2));      % preallocate dy/dt
      
      % Evaluate the 2 components of the function at one edge of the grid
      % (with edge conditions).
      i = 1;
      dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + c*(1-2*y(i,:)+y(i+2,:));
      dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + c*(3-2*y(i+1,:)+y(i+3,:));
      
      % Evaluate the 2 components of the function at all interior grid points.
      i = 3:2:2*N-3;
      dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
         c*(y(i-2,:)-2*y(i,:)+y(i+2,:));
      dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
         c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:));
      
      % Evaluate the 2 components of the function at the other edge of the grid
      % (with edge conditions).
      i = 2*N-1;
      dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + c*(y(i-2,:)-2*y(i,:)+1);
      dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + c*(y(i-1,:)-2*y(i+1,:)+3);
   end
% -------------------------------------------------------------------------

end  % brussode

% ---------------------------------------------------------------------------
% Subfunction -- the sparsity pattern
%

function S = jpattern(N)
% Jacobian sparsity pattern
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end
% ---------------------------------------------------------------------------

Resuelva el sistema Brusselator para $N=20$ ejecutando la función brussode.

brussode

Resuelva el sistema para $N=50$ especificando una entrada para brussode.

brussode(50)

Consulte también

| | |

Temas relacionados