ode45

비경직성(Nonstiff) 미분방정식 풀기 — 중간 차수 방법

구문

[t,y] = ode45(odefun,tspan,y0)

[t,y] = ode45(odefun,tspan,y0,options)

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)

sol = ode45(___)

설명

[t,y] = ode45(odefun,tspan,y0)입니다. 여기서 tspan = [t0 tf]는 t0에서 tf까지의 구간에서 초기 조건 y0을 사용하여 연립미분방정식 $$y\text{'}=f\left(t,y\right)$$를 적분합니다. 해 배열 y의 각 행은 열 벡터 t에 반환된 값에 대응합니다.

모든 MATLAB^® ODE 솔버는 $$y\text{'}=f\left(t,y\right)$$ 형식의 연립방정식이나 질량 행렬이 있는 문제 $$M\left(t,y\right)y\text{'}=f\left(t,y\right)$$를 풀 수 있습니다. 솔버는 모두 유사한 구문을 사용합니다. ode23s 솔버는 상수 질량 행렬을 갖는 문제만 풀 수 있습니다. ode15s와 ode23t는 특이 질량 행렬을 포함하는 문제(즉, 미분대수 방정식(DAE))를 풀 수 있습니다. odeset의 Mass 옵션을 사용하여 질량 행렬을 지정합니다.

ode45는 다용도 ODE 솔버이며 대부분의 문제에서 최우선적으로 시도해 봐야 할 솔버입니다. 하지만, 문제가 경직성 문제이거나 높은 정확성을 요하는 경우 다른 ODE 솔버가 더 적합할 수 있습니다. 자세한 내용은 ODE 솔버 선택하기 항목을 참조하십시오.

[t,y] = ode45(odefun,tspan,y0,options)는 odeset 함수를 사용하여 생성된 인수인 options로 정의된 적분 설정도 사용합니다. 예를 들어, 절대 허용오차와 상대 허용오차를 지정하려면 AbsTol 옵션과 RelTol 옵션을 사용하고 질량 행렬을 제공하려면 Mass 옵션을 사용하십시오.

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)는 이벤트 함수라고 하는 (t,y)의 함수가 0인 위치를 추가로 찾습니다. 출력값에서 te는 이벤트 발생 시간이고, ye는 이벤트 발생 시 계산된 해이며, ie는 트리거된 이벤트의 인덱스입니다.

각 이벤트 함수에 대해, 0에서 적분을 종료할지 여부와 영점교차의 방향을 고려할지 여부를 지정합니다. 이를 수행하려면 'Events' 속성을 함수(예: myEventFcn 또는 @myEventFcn)로 설정하고 대응 함수 [value,isterminal,direction] = myEventFcn(t,y)를 생성합니다. 자세한 내용은 ODE 이벤트 위치 항목을 참조하십시오.

sol = ode45(___)는 deval과 함께 사용하여 구간 [t0 tf] 내의 임의의 점에서 해를 계산할 수 있는 구조체를 반환합니다. 위에 열거된 구문에 나와 있는 입력 인수를 원하는 대로 조합하여 사용할 수 있습니다.

예제

단일 해 성분을 갖는 ODE

단일 해 성분을 갖는 단순한 ODE는 솔버 호출 시 익명 함수로 지정할 수 있습니다. 익명 함수는 입력값 중 하나가 함수에 사용되지 않는 경우에도 두 입력값 (t,y)를 모두 받아야 합니다.

다음 ODE를 풉니다.

$$y^{\prime }=2t.$$

시간 구간 [0 5]와 초기 조건 y0 = 0을 지정합니다.

tspan = [0 5];
y0 = 0;
[t,y] = ode45(@(t,y) 2*t, tspan, y0);

해를 플로팅합니다.

plot(t,y,'-o')

비경직성(Nonstiff) 방정식 풀기

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

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0,$$

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

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

함수 파일 vdp1.m은 $$\mu =1$$를 사용하여 반데르폴 방정식을 나타냅니다. 변수 $$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) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

ODE 함수에 추가 파라미터 전달하기

ode45는 두 개의 입력 인수 t와 y를 사용하는 함수에만 동작합니다. 하지만, 함수 외부에 추가 파라미터를 정의하고 함수 핸들을 지정할 때 이를 전달하는 방식으로 추가 파라미터를 전달할 수 있습니다.

다음 ODE를 풉니다.

$$y^{\prime \prime }=\frac{A}{B}ty.$$

방정식을 1계 시스템으로 다시 작성하면 다음 결과가 생성됩니다.

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\frac{A}{B}t{y}_1.\end{array}$$

이 예시의 끝에 포함된 로컬 함수 odefcn은 이 연립방정식을 4개의 입력 인수 t, y, A, B를 받는 함수로 나타냅니다.

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

ode45를 사용하여 ODE를 풉니다. A와 B의 미리 정의된 값을 odefcn에 전달하도록 함수 핸들을 지정합니다.

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);

결과를 플로팅합니다.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

여러 초기 조건을 사용하여 ODE 풀기

방정식이 한 개인 단순한 연립 ODE의 경우 y0을 여러 초기 조건을 포함하는 벡터로 지정할 수 있습니다. 이 기법은 스칼라 확장을 통해 각 초기값에 대해 하나의 독립적인 연립방정식을 만들고, ode45는 이 연립방정식을 풀어서 각 초기값에 대한 결과를 생성합니다.

방정식 $\mathit{f}\left(\mathit{t},\mathit{y}\right)=-2\mathit{y}+2\cos \left(\mathit{t}\right)\sin \left(2\mathit{t}\right)$를 나타내는 익명 함수를 만듭니다. 이 함수는 t와 y에 대해 두 개의 입력값을 받아야 합니다.

yprime = @(t,y) -2*y + 2*cos(t).*sin(2*t);

범위 $\left[-5,5\right]$에 있는 여러 초기 조건으로 구성된 벡터를 만듭니다.

y0 = -5:5; 

ode45를 사용하여 시간 구간 $\left[0,3\right]$에서 각 초기 조건의 방정식을 풉니다.

tspan = [0 3];
[t,y] = ode45(yprime,tspan,y0);

결과를 플로팅합니다.

plot(t,y)
grid on
xlabel('t')
ylabel('y')
title('Solutions of y'' = -2y + 2 cos(t) sin(2t), y(0) = -5,-4,...,4,5','interpreter','latex')

이 기법은 여러 초기 조건을 가진 단순한 ODE를 풀려는 경우 유용합니다. 그러나 이 기법에는 장단점이 있습니다.

• 여러 초기 조건을 가진 여러 연립방정식은 풀 수 없습니다. 이 기법은 여러 초기 조건을 가진 한 개의 방정식을 풀려는 경우에만 사용할 수 있습니다.

• 솔버가 각 스텝에서 선택하는 시간 스텝은 연립방정식에서 가장 작은 스텝을 취해야 하는 방정식을 기반으로 합니다. 즉, 솔버는 한 초기 조건의 방정식을 충족하기 위해 작은 스텝을 취할 수 있지만, 나머지 방정식은 독립적으로 풀 경우 다른 스텝 크기를 사용합니다. 그럼에도 불구하고, 한 번에 여러 초기 조건의 방정식을 푸는 것이 for 루프를 사용하여 방정식을 개별적으로 푸는 것보다 일반적으로 더 빠릅니다.

이 기법에 대한 자세한 내용은 여러 초기 조건을 사용하여 연립 ODE 풀기 항목을 참조하십시오.

시간 종속적인 항을 갖는 ODE

다음과 같은 시간 종속적인 파라미터를 갖는 ODE가 있다고 가정해 보겠습니다.

초기 조건은 입니다. 함수 f(t)는 시간 ft에서 계산된 n×1 벡터 f에 의해 정의됩니다. 함수 g(t)는 시간 gt에서 계산된 m×1 벡터 g에 의해 정의됩니다.

벡터 f와 g를 만듭니다.

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

f와 g를 보간하는 myode라는 함수를 작성하여 지정된 시간에서의 시간 종속적인 항의 값을 구합니다. 예제의 나머지 부분을 실행하려면 이 함수를 현재 폴더에 저장하십시오.

myode 함수는 각 시간 스텝에서 ODE를 계산하기 위해 추가 입력 인수를 받지만, ode45는 처음 두 개의 입력 인수 t와 y만을 사용합니다.

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t

ode45를 사용하여 시간 구간 [1 5]에 대한 방정식을 풉니다. ode45가 myode의 처음 두 개 입력 인수만 사용하도록 함수 핸들을 사용하여 함수를 지정합니다. 또한, odeset을 사용하여 오차 임계값을 완화합니다.

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

해 y를 시간 지점 t로 구성된 함수로 플로팅합니다.

plot(t,y)

해 구조체를 계산하고 확장하기

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

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0.$$

ode45를 사용하여 $$\mu =1$$인 반데르폴 방정식을 풉니다. 함수 vdp1.m은 MATLAB®과 함께 제공되는 함수로, 이 방정식을 담고 있습니다. 솔버 및 계산 지점 등 해에 대한 정보를 포함하는 구조체를 반환하는 단일 출력값을 지정합니다.

tspan = [0 20];
y0 = [2 0];
sol = ode45(@vdp1,tspan,y0)

sol = struct with fields:
     solver: 'ode45'
    extdata: [1×1 struct]
          x: [0 1.0048e-04 6.0285e-04 0.0031 0.0157 0.0785 0.2844 0.5407 0.8788 1.4032 1.8905 2.3778 2.7795 3.1285 3.4093 3.6657 3.9275 4.2944 4.9013 5.3506 5.7998 6.2075 6.5387 6.7519 6.9652 7.2247 7.5719 8.1226 8.6122 9.1017 9.5054 … ] (1×60 double)
          y: [2×60 double]
      stats: [1×1 struct]
      idata: [1×1 struct]

linspace를 사용하여 구간 [0 20] 내에 250개의 지점을 생성합니다. deval을 사용하여 이들 지점에서 해를 계산합니다.

x = linspace(0,20,250);
y = deval(sol,x);

해의 첫 번째 성분을 플로팅합니다.

plot(x,y(1,:))

odextend를 사용하여 해를 $$t_f=35$$로 확장하고 그 결과를 원래 플롯에 추가합니다.

sol_new = odextend(sol,@vdp1,35);
x = linspace(20,35,350);
y = deval(sol_new,x);
hold on
plot(x,y(1,:),'r')

입력 인수

odefun — 풀려는 함수

함수 핸들

풀려는 함수로, 적분할 함수를 정의하는 함수 핸들로 지정됩니다.

스칼라 t와 열 벡터 y에 대한 함수 dydt = odefun(t,y)는 $$f\left(t,y\right)$$에 대응하는 single이나 double 데이터형의 열 벡터 dydt를 반환해야 합니다. odefun은 입력 인수 t와 y 모두를 받아야만 하는데, 이는 이 중 하나가 함수에 사용되지 않더라도 마찬가지입니다.

예를 들어, $$y\text{'}=5y-3$$을 풀려면 다음 함수를 사용하십시오.

function dydt = odefun(t,y)
dydt = 5*y-3;
end

연립방정식에 대한 odefun의 출력값은 벡터입니다. 벡터의 각 요소는 한 개 방정식 우변의 계산된 값입니다. 예를 들어, 다음과 같이 두 개로 이루어진 연립방정식이 있다고 가정하겠습니다.

$$\begin{array}{l}y{\text{'}}_1=y_1+2y_2\\ y{\text{'}}_2=3y_1+2y_2\end{array}$$

각 시간 스텝에서 각 방정식의 우변 값을 계산하는 함수는 다음과 같습니다.

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);
end

함수 odefun에 추가 파라미터를 제공하는 방법에 대한 자세한 내용은 함수를 파라미터화하기 항목을 참조하십시오.

예: @myFcn

데이터형: function_handle

tspan — 적분 구간

벡터

적분 구간으로, 벡터로 지정됩니다. 최소한, tspan은 초기 시간과 최종 시간을 지정하는, 요소를 2개 가진 벡터 [t0 tf]여야 합니다. t0과 tf 사이 특정 시간에서의 해를 구하려면 [t0,t1,t2,...,tf] 형식의 더 긴 벡터를 사용해야 합니다. tspan의 요소는 모두 증가하거나 모두 감소해야 합니다.

솔버는 초기 시간 tspan(1)에 y0으로 주어진 초기 조건을 적용하고 tspan(1)에서 tspan(end)까지의 구간에 대해 적분을 계산합니다.

• tspan이 두 개의 요소 [t0 tf]를 갖는 경우, 솔버는 구간 내 각 내부 적분 스텝에서 계산된 해를 반환합니다.

• tspan이 세 개 이상의 요소 [t0,t1,t2,...,tf]를 갖는 경우 솔버는 지정된 지점에서 계산된 해를 반환합니다. 하지만, 솔버가 tspan에 지정된 각 지점으로 정확하게 이동하는 것은 아닙니다. 대신에 솔버는 자체 내부 스텝을 사용하여 해를 계산한 다음에 tspan의 요청된 지점에서 해를 계산합니다. 지정된 지점에서 구해진 해는 각 내부 스텝에서 계산된 해와 동일한 정도의 정확성을 가집니다.

  중간 지점을 여러 개 지정해도 계산의 효율성에는 거의 영향을 미치지 않지만 대규모 방정식 시스템의 경우에는 메모리 관리에 영향을 미칠 수 있습니다.

tspan의 값은 다음과 같이 솔버가 InitialStep과 MaxStep에 적합한 값을 계산하는 데 사용됩니다.

• tspan에 여러 개의 중간 지점 [t0,t1,t2,...,tf]가 있는 경우 지정된 지점은 문제의 크기를 나타내며, 이는 솔버에서 사용되는 InitialStep의 값에 영향을 미칠 수 있습니다. 따라서, 솔버에서 구한 해는 tspan을 요소를 2개 가진 벡터, 또는 중간 지점을 갖는 벡터로 지정하는지 여부에 따라 달라질 수 있습니다.

• tspan의 초기값과 최종 값은 최대 스텝 크기 MaxStep을 계산하는 데 사용됩니다. 따라서, tspan에서 초기값이나 최종 값을 변경하면 솔버가 다른 스텝 시퀀스를 사용하게 되어 해가 바뀔 수 있습니다.

예: [1 10]

예: [1 3 5 7 9 10]

데이터형: single | double

y0 — 초기 조건

벡터

초기 조건으로, 벡터로 지정됩니다. y0이 odefun에 정의된 각 방정식에 대한 초기 조건을 포함하도록 y0은 odefun의 벡터 출력값과 길이가 동일해야 합니다.

데이터형: single | double

options — options 구조체

구조체형 배열

options 구조체로, 구조체형 배열로 지정됩니다. odeset 함수를 사용하여 options 구조체를 만들거나 수정할 수 있습니다. 각 솔버와 호환되는 옵션 목록은 ODE 옵션 요약 항목을 참조하십시오.

예: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot)은 1e-5의 상대 허용오차를 지정하고, 솔버 통계량 표시를 설정하고, 해가 계산될 때 해를 플로팅하도록 출력 함수 @odeplot을 지정합니다.

데이터형: struct

출력 인수

t — 계산 지점

열 벡터

계산 지점으로, 열 벡터로 반환됩니다.

• tspan이 두 개의 요소 [t0 tf]를 포함하는 경우 t는 적분이 수행된 내부 계산 지점을 포함합니다.

• tspan이 세 개 이상의 요소를 포함하는 경우 t는 tspan과 동일합니다.

y — 해

배열

해로, 배열로 반환됩니다. y의 각 행은 대응하는 t의 행에 반환된 값에서 계산된 해에 해당합니다.

te — 이벤트 시간

열 벡터

이벤트 시간으로, 열 벡터로 반환됩니다. te의 이벤트 시간은 ye로 반환되는 해에 대응되며, ie는 발생한 이벤트를 지정합니다.

ye — 이벤트 발생 시 계산된 해

배열

이벤트 발생 시 계산된 해로, 배열로 반환됩니다. te의 이벤트 시간은 ye로 반환되는 해에 대응되며, ie는 발생한 이벤트를 지정합니다.

ie — 트리거된 이벤트 함수의 인덱스

열 벡터

트리거된 이벤트 함수의 인덱스로, 열 벡터로 반환됩니다. te의 이벤트 시간은 ye로 반환되는 해에 대응되며, ie는 발생한 이벤트를 지정합니다.

sol — 계산할 구조체

구조체형 배열

계산할 구조체로, 구조체형 배열로 반환됩니다. 이 구조체를 deval 함수와 함께 사용하여 구간 [t0 tf] 내에 있는 임의의 지점에서 해를 계산합니다. sol 구조체형 배열은 항상 다음 필드를 포함합니다.

구조체 필드	설명
sol.x	솔버에서 선택한 스텝으로 구성된 행 벡터입니다.
sol.y	해입니다. 각 열 sol.y(:,i)에는 시간 sol.x(i)에서 계산된 해가 포함됩니다.
sol.solver	솔버 이름입니다.

또한, odeset의 Events 옵션을 지정한 상태에서 이벤트가 발견되면 sol은 다음 필드도 포함합니다.

구조체 필드	설명
sol.xe	이벤트가 발생한 지점입니다. sol.xe(end)는 종료 이벤트(있을 경우)의 정확한 지점을 포함합니다.
sol.ye	sol.xe의 이벤트에 대응하는 해입니다.
sol.ie	Events 옵션에 지정된 함수가 반환하는 벡터에 대한 인덱스. 이러한 값은 솔버가 감지한 이벤트를 나타냅니다.

알고리즘

ode45는 명시적 룽게-쿠타(Explicit Runge-Kutta) (4,5) 공식인 도르만드-프린스(Dormand-Prince) 쌍을 기반으로 합니다. 이는 단일 스텝(Single-step) 솔버입니다. 즉, y(tn)을 계산할 경우 이 솔버는 바로 이전 시간 점 y(tn-1)에서의 해만 필요로 합니다[1], [2].

참고 문헌

[1] Dormand, J. R. and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comp. Appl. Math., Vol. 6, 1980, pp. 19–26.

[2] Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.

확장 기능

C/C++ 코드 생성

MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

사용법 관련 참고 및 제한 사항:

• 모든 odeset 옵션 인수는 상수여야 합니다.

• 코드 생성 시 options 구조체의 상수 질량 행렬은 지원되지 않습니다. 질량 행렬은 함수로 제공하십시오.

• 최소 2개의 출력 인수 T와 Y를 제공해야 합니다.

• 입력값 유형은 동일해야 합니다. 즉, 모두 double형이거나 모두 single형이어야 합니다.

• 가변 크기 지원이 활성화되어 있어야 합니다. tspan 에 두 요소가 있거나 이벤트 함수를 사용하는 경우 코드 생성 시 동적 메모리 할당이 필요합니다.

버전 내역

R2006a 이전에 개발됨

참고 항목

ode23 | ode78 | ode89 | ode113 | ode15s | odeset | odeget | deval | odextend

도움말 항목

• ODE 솔버 선택하기
• ODE 옵션 요약
• 비경직성(Nonstiff) ODE 풀기
• 익명 함수
• 일반적인 ODE 문제 해결하기