Решение жестких ОДУ

Открыть скрипт

Эта страница содержит два примера решения жестких обыкновенных дифференциальных уравнений с помощью ode15s. MATLAB ® имеет четыре решателя, предназначенных для жестких ОДУ.

ode15s
ode23s
ode23t
ode23tb

Для большинства жестких задач ode15s работает лучше всего. Однако ode23s, ode23t, и ode23tb может быть более эффективным, если проблема позволяет допускать грубую ошибку.

Что такое жесткая ОДУ?

Для некоторых задач ОДУ размер шага, взятый решателем, принуждается к необоснованно малому уровню по сравнению с интервалом интегрирования, даже в области, где кривая решения гладкая. Эти размеры шага могут быть настолько маленькими, что для прохождения короткого временного интервала могут потребоваться миллионы оценок. Это может привести к тому, что решатель не справится с интегрированием, но даже если он преуспеет, это займет очень много времени.

Уравнения, которые вызывают такое поведение в решателях ОДУ, считаются жесткими. Проблема, которую ставят жесткие ОДУ, заключается в том, что явные решатели (такие как ode45) несостоятельно медленны в достижении решения. Вот почему ode45 классифицируется как нежесткий решатель наряду с ode23 и ode113.

Решатели, которые разработаны для жестких ОДУ, известных как жесткие решатели, обычно выполняют больше работы на шаг. Окупаемость заключается в том, что они способны делать гораздо большие шаги и имеют улучшенную числовую стабильность по сравнению с нежесткими решателями.

Опции решателя

Для жестких задач установка матрицы Якобия с помощью odeset особенно важно. Жесткие решатели используют матрицу Якобиана $\partial f_i / \partial y_j$ , чтобы оценить локальное поведение ОДУ по мере начала интегрирования, поэтому предоставление якобиевой матрицы (или, для больших разреженных систем, её шаблона разреженности) очень важно для эффективности и надежности. Используйте Jacobian, JPattern, или Vectorized опции odeset для указания сведений о якобиане. Если вы не поставляете якобиан, то решатель оценивает его численно, используя конечные различия.

См. odeset полный список других опций решателя.

Пример: Жесткое уравнение Ван дер Поля

Уравнение Ван дер Поля является ОДУ второго порядка

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

где $\mu > 0$ является скалярным параметром. Когда, $\mu = 1$ полученная система ОДУ является нежесткой и легко решается с помощью ode45. Однако, если вы увеличиваетесь $\mu$ до 1000, то решение резко меняется и демонстрирует колебания на гораздо более длительной временной шкале. Аппроксимация решения задачи начального значения становится более трудной. Поскольку эта конкретная задача является жесткой, решатель, предназначенный для нежестких задач, таких как ode45, слишком неэффективно, чтобы быть практичным. Используйте жесткий решатель, такой как ode15s вместо этого для этой задачи.

Переписать уравнение Ван дер Поля как систему ОДУ первого порядка путем замены. $y'_1 = y_2$ Получившаяся система ОДУ первого порядка

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

The vdp1000 функция оценивает уравнение ван дер Пол используя. $\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)];

Используйте ode15s функция, чтобы решить задачу с вектором начальных условий [2; 0], в течение временного интервала [0 3000]. По причинам масштабирования постройте график только первого компонента решения.

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

The vdpode функция также решает ту же задачу, но принимает заданное пользователем значение для. $\mu$ Уравнения становятся все более жесткими по мере $\mu$ увеличения.

Пример: Разреженная Брюссельская система

Классическая брюссельская система уравнений потенциально большая, жесткая и разреженная. Система Брюсселя моделирует диффузию в химической реакции и представлена системой уравнений с включением, $u$ ,, $v$ $u'$ и. $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}$

Файл функции brussode решает этот набор уравнений на временном интервале [0,10] с. $\alpha = 1/50$ Начальные условия

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

где $x_j = i/N+1$ для. $i=1,...,N$ Поэтому $2N$ в системе существуют уравнения, но Якобиан является $\partial f / \partial y$ перемычкой с постоянной шириной 5, если уравнения упорядочены как. По $u_1,v_1,u_2,v_2,...$ $N$ мере увеличения проблема становится все более жесткой, и якобиан становится все более скудным.

Вызов функции brussode(N), for, $N \ge 2$ задает значение для N в системе уравнений, соответствующих количеству точек сетки. По умолчанию brussode использует. $N = 20$

brussode содержит несколько подфункций:

Вложенная функция f(t,y) кодирует систему уравнений для задачи Брюсселатора, возвращая вектор.
Локальная функция jpattern(N) возвращает разреженную матрицу 1s и 0s, показывающую местоположения ненули в якобиане. Эта матрица присвоена JPattern поле структуры опций. Решатель ОДУ использует этот шаблон разреженности, чтобы сгенерировать якобиан численно как разреженную матрицу. Предоставление этого шаблона разреженности в задаче значительно сокращает количество вычислений функции, необходимых для генерации 2N-by-2N якобиана, с 2N вычислений до всего 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
% ---------------------------------------------------------------------------