Main Content

비경직성(Nonstiff) ODE 풀기

스크립트 열기

이 페이지에는 ode45를 사용한 비경직성 상미분 방정식 풀이에 대한 두 가지 예제가 포함되어 있습니다. MATLAB®에는 비경직성 ODE에 대한 여러 개의 솔버가 있습니다.

  • ode45

  • ode23

  • ode78

  • ode89

  • ode113

ode45는 대부분의 비경직성 문제에서 가장 잘 동작합니다. 그러나, 조금 덜 엄격한 허용오차가 용인되는 문제나 반경직성이 있는 문제에는 ode23이 권장됩니다. 마찬가지로, 엄격한 허용오차를 가지는 문제에서나 ODE 함수를 계산하는 데 시간이 많이 걸리는 경우 ode113이 ode45보다 더 효율적일 수 있습니다. ode78 및 ode89는 안정성에 있어 정확도가 중요한 역할을 하는 긴 구간에서의 적분에 뛰어난 고차수 솔버입니다.

비경직성 솔버로 문제를 풀 때 시간이 오래 걸리거나 적분에 지속적으로 실패할 경우에는 문제가 경직성 문제일 수 있습니다. 자세한 내용은 경직성(Stiff) ODE 풀기 항목을 참조하십시오.

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

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

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

여기서 $\mu > 0$는 스칼라 파라미터입니다. 대입 $y'_1 = y_2$를 수행하여 이 방정식을 연립 1계 ODE로 재작성합니다. 결과로 생성되는 1계 연립 ODE는 다음과 같습니다.

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

연립 ODE는 ODE 솔버가 사용할 수 있는 함수 파일로 코딩되어야 합니다. ODE 함수의 일반 함수 시그니처는 다음과 같습니다.

  dydt = odefun(t,y)

즉, 이 함수는 계산에 t를 사용하지 않는 경우에도 t와 y를 모두 입력값으로 받아야 합니다.

함수 파일 vdp1.m은 $\mu = 1$를 사용하여 반데르폴 방정식을 코딩합니다. 변수 $y_1$와 $y_2$는 y(1)과 y(2)로 표현되며, 요소를 2개 가진 열 벡터 dydt에는 $y'_1$ 및 $y'_2$에 대한 표현식이 포함되어 있습니다.

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

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

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

초기값 [2 0]을 사용하고 시간 구간 [0 20]에 대해 ode45 함수를 사용하여 ODE를 풉니다. 출력값은 시간 지점 t의 열 벡터와 해 배열 y입니다. y의 각 행은 t의 해당 행에 반환된 시간에 대응됩니다. y의 첫 번째 열은 $y_1$에 대응되고, 두 번째 열은 $y_2$에 대응됩니다.

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

t에 대해 $y_1$과 $y_2$의 해를 플로팅합니다.

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) using ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

vdpode 함수는 동일한 문제를 풀지만, $\mu$에 대한 사용자 지정 값을 받습니다. 반데르폴 방정식은 $\mu$가 증가함에 따라 경직성(Stiff)이 됩니다. 예를 들어, 값 $\mu = 1000$을 사용하는 경우 ode15s와 같은 경직성 솔버를 사용하여 시스템을 풀어야 합니다.

예제: 비경직성 오일러(Euler) 방정식

외부 힘이 작용하지 않는 강체에 대한 오일러 방정식은 비경직성 문제용으로 제작된 ODE 솔버에 대한 표준 테스트 문제입니다.

방정식은 다음과 같습니다.

$$ \begin{array}{cl} y'_1 &= y_2y_3 \\ y'_2 &= -y_1y_3 \\ y'_3 &=
-0.51y_1y_2. \end{array}$$

함수 파일 rigidode는 이 1계 연립방정식을 정의하고 초기값 $y_1$, $y_2$, $y_3$에 대응하는 초기 조건 벡터 [0; 1; 1]을 사용하여 시간 구간 [0 12]에 대해 이 연립방정식을 풉니다. 로컬 함수 f(t,y)는 이 연립방정식을 인코딩합니다.

rigidode는 출력 인수 없이 ode45를 호출하므로 솔버는 디폴트 출력 함수 odeplot을 사용하여 각 스텝이 완료된 후 해에 해당하는 점을 자동으로 플로팅합니다.

function rigidode
%RIGIDODE  Euler equations of a rigid body without external forces.
%   A standard test problem for non-stiff solvers proposed by Krogh.  The
%   analytical solutions are Jacobian elliptic functions, accessible in
%   MATLAB.  The interval here is about 1.5 periods; it is that for which
%   solutions are plotted on p. 243 of Shampine and Gordon.
%
%   L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary
%   Differential Equations, W.H. Freeman & Co., 1975.
%
%   See also ODE45, ODE23, ODE113, FUNCTION_HANDLE.

%   Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94
%   Copyright 1984-2014 The MathWorks, Inc.

tspan = [0 12];
y0 = [0; 1; 1];

% solve the problem using ODE45
figure;
ode45(@f,tspan,y0);

% --------------------------------------------------------------------------

function dydt = f(t,y)
dydt = [    y(2)*y(3)
   -y(1)*y(3)
   -0.51*y(1)*y(2) ];

rigidode 함수를 호출하여 비경직성 오일러 방정식을 풉니다.

rigidode
title('Solution of Rigid Body w/o External Forces using ODE45')
legend('y_1','y_2','y_3','Location','Best')

참고 항목

| | | |

관련 항목

외부 웹사이트