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

Открыть сценарий

Эта страница содержит два примера решения жестких обыкновенных дифференциальных уравнений с использованием 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, то решение резко меняется и проявляет колебания в гораздо более длительных временных масштабах. Аппроксимация решения задачи начального значения становится более сложной. Поскольку эта конкретная проблема является жесткой, решатель, предназначенный для проблем nonstiff, таких как 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}
$

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

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), для, $N \ge 2$ задает значение для N в системе уравнений, соответствующих количеству точек сетки. По умолчанию brussode использует. $N = 20$

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

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