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

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

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

  • ode15s

  • ode23s

  • ode23t

  • ode23tb

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

Что такое Жесткое ОДУ?

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

Уравнения, которые вызывают это поведение в решателях ОДУ, как говорят, жестки. Проблема, что жесткое положение ОДУ то, что явные решатели (такие как ode45) являются ненадежно медленными в достижении решения. Это то, почему ode45 классифицируется как нежесткий решатель наряду с ode23, ode78, ode89, и 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$ до 1 000, затем решение изменяется существенно и показывает колебание на намного более длительном масштабе времени. Аппроксимация решения задачи с начальными значениями становится более трудной. Поскольку эта конкретная проблема жестка, решатель, предназначенный для нежестких проблем, такова как 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$ увеличения.

Пример: Разреженная система Brusselator

Классическая система уравнений Brusselator является потенциально большой, жесткой, и разреженной. Диффузия системных моделей Brusselator в химической реакции, и представлена включением системы уравнений$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) кодирует систему уравнений для проблемы Brusselator, возвращая вектор.

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

Решите систему Brusselator для$N=20$ путем выполнения функционального brussode.

brussode

Решите систему для$N=50$ путем определения входа к brussode.

brussode(50)

Смотрите также

| | |

Похожие темы

Документация MATLAB

Поддержка

Для просмотра документации необходимо авторизоваться на сайте