경직성(Stiff) ODE 풀기

이 페이지에는 ode15s를 사용한 경직성 상미분 방정식 풀이에 대한 두 가지 예제가 포함되어 있습니다. MATLAB®에는 경직성 ODE용으로 설계된 네 가지 솔버가 있습니다.

ode15s
ode23s
ode23t
ode23tb

ode15s는 대부분의 경직성 문제에서 가장 잘 동작합니다. 그러나, 조금 덜 엄격한 허용오차가 용인되는 문제의 경우 ode23s, ode23t, ode23tb가 더 효율적일 수 있습니다.

경직성(Stiff) ODE란?

일부 ODE 문제의 경우, 해 곡선이 매끄럽게 나타나는 영역에서도 솔버가 취하는 스텝 크기가 적분 구간과 비교하여 합리적이지 않은 작은 수준까지 강제로 줄어듭니다. 이러한 스텝 크기는 짧은 시간 구간을 지날 때에도 수백만 번의 계산이 필요할 정도로 매우 작을 수 있습니다. 이 경우 솔버가 적분에 실패할 수 있으며, 성공하더라도 매우 긴 시간이 걸리게 됩니다.

ODE 솔버에서 이러한 동작을 초래하는 방정식을 경직성(Stiff) 방정식이라고 합니다. 경직성 ODE의 문제는 양함수 솔버(예: ode45)로 해를 구할 때 참을 수 없을 정도로 느리다는 점입니다. 이 때문에 ode45는 ode23, ode78, ode89 및 ode113과 함께 비경직성 솔버로 분류됩니다.

경직성 ODE용으로 설계된 이른바 경직성 솔버는 일반적으로 스텝당 더 많은 작업을 수행합니다. 대신 더 큰 스텝을 취할 수 있고 비경직성 솔버에 비해 향상된 수치적 안정성을 가집니다.

솔버 옵션

경직성(Stiff) 문제를 풀려는 경우 odeset을 사용하여 야코비 행렬을 지정하는 것이 특히 중요합니다. 경직성 솔버는 야코비 행렬 $\partial f_i / \partial y_j$ 를 사용하여 적분이 진행됨에 따른 ODE의 로컬 동작을 추정하므로, 야코비 행렬(또는 큰 희소 시스템의 경우 그 희소성 패턴)을 제공하는 것이 효율성과 안정성 측면에서 매우 중요합니다. odeset의 Jacobian, JPattern 또는 Vectorized 옵션을 사용하여 야코비 행렬 관련 정보를 지정합니다. 야코비 행렬을 제공하지 않으면 솔버가 유한 차분을 사용하여 수치적으로 이에 대한 추정값을 계산합니다.

다른 솔버 옵션에 대한 전체 목록은 odeset 항목을 참조하십시오.

예제: 경직성 반데르폴 방정식(Stiff van der Pol equation)

반데르폴 방정식은 2계 ODE입니다.

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

여기서 $\mu > 0$ 는 스칼라 파라미터입니다. $\mu = 1$ 이면 결과로 생성되는 연립 ODE는 비경직성이고 ode45를 사용하여 쉽게 풀 수 있습니다. 그러나, $\mu$ 를 1000으로 증가시키면 해가 급격히 변하고 훨씬 더 긴 시간 스케일에서 진동이 나타납니다. 초기값 문제의 해에 대한 근삿값을 구하는 것은 더욱 어려워집니다. 이러한 특정 문제는 경직성 문제이기 때문에 ode45와 같이 비경직성 문제용으로 제작된 솔버는 너무 비효율적이라 실용적이지 못합니다. 이 문제에는 ode15s와 같은 경직성 솔버를 대신 사용하십시오.

대입을 수행하여 반데르폴 방정식을 1계 연립 ODE로 바꿉니다. 결과로 생성되는 1계 연립 ODE는 다음과 같습니다.

$
\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 함수를 사용하여 시간 구간 [0 3000]에 대해 초기 조건 벡터 [2; 0]을 적용해 문제를 풉니다. 스케일링 때문에 해의 첫 번째 성분만 플로팅합니다.

[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$ 가 증가함에 따라 방정식도 점점 더 경직성을 띠게 됩니다.

예제: 희소 브뤼셀레이터 시스템(Sparse Brusselator System)

전형적인 브뤼셀레이터 연립방정식은 크기가 크고, 경직성(Stiff)을 가지며 희소 형식일 수 있습니다. 브뤼셀레이터 시스템은 화학 반응의 확산을 모델링하고 , , , 가 포함된 연립방정식으로 표현됩니다.

$\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는 $\alpha = 1/50$ 을 사용하여 시간 구간 [0,10]에 대해 이 방정식 세트를 풉니다. 초기 조건은 다음과 같습니다.

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

여기서 에 대해 입니다. 따라서 개의 연립방정식이 있지만, 방정식이 순서로 정렬된 경우 야코비 행렬 $\partial f / \partial y$ 는 너비가 5인 일정한 띠가 있는 행렬입니다. 이 증가하면 문제의 경직성이 증가하며, 야코비 행렬도 더 희소화됩니다.

$N \ge 2$ 에 대한 함수 호출 brussode(N)은 연립방정식에서 N의 값을 그리드 점의 개수에 해당하는 값으로 지정합니다. 기본적으로, brussode는 을 사용합니다.

brussode는 몇 가지 부 함수를 포함합니다.

중첩 함수 f(t,y)는 브뤼셀레이터 문제(Brusselator Problem)에 대한 연립방정식을 인코딩하며 벡터를 반환합니다.
로컬 함수 jpattern(N)은 1과 0으로 구성된 희소 행렬을 반환하여 야코비 행렬에서 0이 아닌 요소의 위치를 보여줍니다. 이 행렬은 options 구조체의 JPattern 필드에 할당됩니다. ODE 솔버는 이 희소성 패턴을 사용하여 야코비 행렬을 수치적으로 희소 행렬로 생성합니다. 문제에 이 희소성 패턴을 제공하면 2N×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
% ---------------------------------------------------------------------------