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

이 예제에서는 여러 초기 조건 세트로 연립상미분방정식을 푸는 두 가지 기법을 비교합니다. 기법은 다음과 같습니다.

• for 루프를 사용하여 각 초기 조건 세트에 대해 하나씩 여러 개의 시뮬레이션을 수행합니다. 이 기법은 사용하기 간단하지만, 큰 연립방정식의 경우에는 최상의 성능을 제공하지 않습니다.

• ODE 함수를 벡터화하여 모든 초기 조건 세트에 대해 동시에 연립방정식을 풉니다. 이 기법은 큰 연립방정식을 더 빨리 풀 수 있는 방법이지만, 입력 형태를 적절하게 변경하도록 ODE 함수를 재작성해야 합니다.

이러한 기법을 설명하는 데 사용되는 방정식은 잘 알려진 로트카-볼테라 방정식으로, 포식자 개체군과 피식자 개체군을 나타내는 1차 비선형 미분 방정식입니다.

문제 설명

로트카-볼테라 방정식은 생물학적 시스템에서 포식자 개체군과 피식자 개체군을 나타내는 두 개의 1차 비선형 연립 ODE입니다. 시간이 경과하면 포식자 개체군과 피식자 개체군은 방정식에 따라 바뀝니다.

$$\begin{array}{l}\frac{\mathrm{dx}}{\mathrm{dt}}=\alpha \mathit{x}-\beta \mathrm{xy},\\ \frac{\mathrm{dy}}{\mathrm{dt}}=\delta \mathrm{xy}-\gamma \mathit{y}.\end{array}$$

이 방정식의 변수는 다음과 같습니다.

• $\mathit{x}$는 피식자 개체군 크기입니다.

• $\mathit{y}$는 포식자 개체군 크기입니다.

• $\mathit{t}$는 시간입니다.

• $\alpha$, $\beta$, $\delta$, $\gamma$는 두 종 간의 상호 작용을 나타내는 상수 파라미터입니다. 이 예제에서는 파라미터 값 $\alpha =\gamma =1$, $\beta =0.01$, $\delta =0.02$를 사용합니다.

이 문제의 경우 $\mathit{x}$ 및 $\mathit{y}$에 대한 초기값은 초기 개체군 크기입니다. 방정식을 풀면 종이 상호 작용하는 동안 시간이 지남에 따라 개체군이 어떻게 변경되는지에 대한 정보를 알 수 있습니다.

하나의 초기 조건으로 방정식 풀기

MATLAB®에서 로트카-볼테라 방정식을 풀려면 방정식을 인코딩하는 함수를 작성하고 적분에 사용할 시간 구간을 지정하고 초기 조건을 지정하십시오. 그런 다음 ODE 솔버 중 하나(예: ode45)를 사용하여 시간 경과에 따라 시스템을 시뮬레이션할 수 있습니다.

방정식을 인코딩하는 함수는 다음과 같습니다.

function dpdt = lotkaODE(t,p)
% LOTKA Lotka-Volterra predator-prey model
delta = 0.02;
beta = 0.01;

dpdt = [p(1) .* (1 - beta*p(2));
        p(2) .* (-1 + delta*p(1))];
end

(이 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.)

시스템에 두 개의 방정식이 있으므로 dpdt는 각 방정식에 대해 하나의 요소가 있는 벡터입니다. 또한 해 벡터 p는 각 해 성분에 대해 하나의 요소를 갖습니다. p(1)은 원래 방정식에서 $\mathit{x}$를 나타내고 p(2)는 원래 방정식에서 $\mathit{y}$를 나타냅니다.

다음으로 적분의 시간 구간을 $\left[0,15\right]$로 지정하고 $\mathit{x}$ 및 $\mathit{y}$에 대한 초기 개체군 크기를 50으로 설정합니다.

t0 = 0;
tfinal = 15;
p0 = [50; 50];

ODE 함수, 시간 길이 및 초기 조건을 지정하여 ode45로 연립방정식을 풉니다. 결과로 발생하는 개체군을 시간에 대해 플로팅합니다.

[t,p] = ode45(@lotkaODE,[t0 tfinal],p0);
plot(t,p)
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
legend('Prey','Predators')

해는 주기성을 보이므로 위상 플롯에서 해 서로에 대해 플로팅합니다.

plot(p(:,1),p(:,2))
title('Phase Plot of Predator/Prey Populations')
xlabel('Prey')
ylabel('Predators')

결과 플롯은 주어진 초기 개체군 크기에 대한 해를 보여줍니다. 서로 다른 개체군 크기에 대한 방정식을 풀려면 p0의 값을 변경하고 시뮬레이션을 다시 실행하십시오. 하지만 이 방법은 한 번에 하나의 초기 조건에 대해서만 방정식을 풉니다. 다음 두 섹션에서 다양한 초기 조건에 대해 해를 구하는 기법을 설명합니다.

방법 1: for- 루프를 사용하여 여러 초기 조건 계산하기

초기 조건이 여러 개인 연립 ODE를 푸는 가장 간단한 방법은 for 루프를 사용하는 것입니다. 이 기법은 단일 초기 조건 기법과 동일한 ODE 함수를 사용하지만, for 루프가 풀이 과정을 자동화합니다.

예를 들어, $\mathit{x}$에 대한 초기 개체군 크기 상수를 50으로 고정하고 for 루프를 사용하여 $\mathit{y}$에 대한 초기 개체군 크기를 10에서 400 사이로 다르게 적용할 수 있습니다. y0에 대한 개체군 크기의 벡터를 만든 다음, 값을 순회하여 각 초기 조건 세트에 대해 방정식을 풉니다. 모든 반복의 결과가 포함된 위상 플롯을 플로팅합니다.

y0 = 10:10:400;
for k = 1:length(y0)
    [t,p] = ode45(@lotkaODE,[t0 tfinal],[50 y0(k)]);
    plot(p(:,1),p(:,2))
    hold on
end
title('Phase Plot of Predator/Prey Populations')
xlabel('Prey')
ylabel('Predators')
hold off

위상 플롯은 서로 다른 초기 조건 세트에 대해 계산된 모든 해를 보여줍니다.

방법 2: 벡터화된 ODE 함수를 사용하여 여러 초기 조건 계산하기

여러 초기 조건에 대한 연립 ODE를 푸는 또 다른 방법은 모든 방정식이 동시에 풀리도록 ODE 함수를 다시 작성하는 것입니다. 이 방법을 수행하는 단계는 다음과 같습니다.

• 모든 초기 조건을 ode45에 행렬로 제공합니다. 행렬의 크기는 s×n이며 여기서 s는 해 성분의 개수이고 n은 풀려는 초기 조건의 개수입니다. 그러면 행렬의 각 열은 시스템에 대한 하나의 완전한 초기 조건 세트를 나타냅니다.

• ODE 함수는 초기 조건의 개수인 n에 대한 추가 입력 파라미터를 받아야 합니다.

• ODE 함수 내부에서 솔버가 해 성분 p를 열 벡터로 전달합니다. ODE 함수는 벡터 형태를 크기가 s×n인 행렬로 변경해야 합니다. 그러면 행렬의 각 행에 각 변수에 대한 모든 초기 조건이 포함됩니다.

• 표현식이 해 성분에 대한 벡터를 받을 수 있도록 ODE 함수는 벡터화된 형식으로 방정식을 풀어야 합니다. 즉, f(t,[y1 y2 y3 ...])은 [f(t,y1) f(t,y2) f(t,y3) ...]을 반환해야 합니다.

• 마지막으로, ODE 솔버가 각 함수 호출에서 벡터를 다시 받을 수 있도록 ODE 함수는 출력의 형태를 다시 벡터로 변경해야 합니다.

이러한 단계를 수행하면 ODE 함수가 벡터를 행렬로 형태 변경하고 모든 초기 조건에 대해 각 해 성분을 푸는 동안 ODE 솔버는 해 성분에 대한 벡터를 사용하여 연립방정식을 풀 수 있습니다. 그 결과 한 번의 시뮬레이션으로 모든 초기 조건에 대해 연립방정식을 풀 수 있습니다.

로트카-볼테라 연립방정식에 이 방법을 구현하려면 먼저 초기 조건의 개수 n을 찾은 다음 초기 조건으로 구성된 행렬을 만드십시오.

n = length(y0);
p0_all = [50*ones(n,1) y0(:)]';

다음으로 n을 입력으로 받을 수 있도록 ODE 함수를 다시 작성합니다. n을 사용하여 해 벡터의 형태를 행렬로 변경한 다음, 벡터화된 연립방정식을 풀고 출력의 형태를 다시 벡터로 변경합니다. 이러한 작업을 수행하는 수정된 ODE 함수는 다음과 같습니다.

function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

ode45를 사용하여 모든 초기 조건에 대한 연립방정식을 풉니다. ode45는 ODE 함수가 두 개의 입력을 받기를 요구하므로 n의 값을 작업 공간에서 lotkasystem에 전달하기 위해 익명 함수를 사용합니다.

[t,p] = ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all);

출력 벡터를 크기가 (numTimeSteps*s)×n인 행렬로 형태 변경합니다. 출력 p(:,k)의 각 열에는 하나의 초기 조건 세트에 대한 해가 포함되어 있습니다. 해 성분으로 구성된 위상 플롯을 플로팅합니다.

p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:end,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off

해당 결과는 for 루프 기법으로 구한 결과와 유사합니다. 하지만 벡터화된 풀이 기법 중 몇 가지 속성에 유의해야 합니다.

• 계산된 해는 단일 초기 입력에서 계산된 해와 약간 다를 수 있습니다. 이 차이는 ODE 솔버가 시간 스텝의 크기를 계산하기 위해 전체 시스템에 노름 검사를 적용하기 때문에 발생하며, 따라서 해의 시간 스텝 동작이 약간 다릅니다. 일반적으로 시간 스텝을 변경하면 해의 정확성에는 영향을 미치지 않지만 해가 계산되는 시점에 영향을 미칩니다.

• 시스템의 숫자형 야코비 행렬을 자동으로 계산하는 경직성(Stiff) ODE 솔버(ode15s ode23s, ode23t, ode23tb)의 경우, odeset의 JPattern 옵션을 사용하여 야코비 행렬의 블록 대각 희소성 패턴을 지정하면 계산의 효율성을 높일 수 있습니다. 야코비 행렬의 블록 대각 형식은 재작성된 ODE 함수에서 수행된 입력 형태 변경의 결과로 만들어집니다.

시간 측정 결과 비교하기

timeit을 사용하여 위에서 설명한 방법들의 시간을 각각 측정합니다. 하나의 초기 조건 세트를 가진 방정식 풀이 시간은 다른 방법들의 확장 성능을 확인하기 위한 기준치로 포함되어 있습니다.

% Time one IC
baseline = timeit(@() ode45(@lotkaODE,[t0 tfinal],p0),2);

% Time for-loop
for k = 1:length(y0)
    loop_timing(k) = timeit(@() ode45(@lotkaODE,[t0 tfinal],[50 y0(k)]),2);
end
loop_timing = sum(loop_timing);

% Time vectorized fcn
vectorized_timing = timeit(@() ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all),2);

시간 측정 결과가 포함된 표를 만듭니다. 모든 결과에 1e3을 곱하여 시간을 밀리초로 표현합니다. 풀이당 시간(풀려는 초기 조건 개수로 각 시간을 나눈 값)이 있는 열을 포함시킵니다.

TimingTable = table(1e3.*[baseline; loop_timing; vectorized_timing], 1e3.*[baseline; loop_timing/n; vectorized_timing/n],...
    'VariableNames',{'TotalTime (ms)','TimePerSolution (ms)'},'RowNames',{'One IC','Multi ICs: For-loop', 'Mult ICs: Vectorized'})

TimingTable=3×2 table
                            TotalTime (ms)    TimePerSolution (ms)
                            ______________    ____________________

    One IC                     0.40718              0.40718       
    Multi ICs: For-loop         11.967              0.29916       
    Mult ICs: Vectorized        1.6168              0.04042       

TimePerSolution 열에서 벡터화된 기법이 세 가지 방법 중 가장 빠른 방법임을 알 수 있습니다.

로컬 함수(Local Function)

다음은 ode45가 해를 계산하기 위해 호출하는 로컬 함수입니다.

function dpdt = lotkaODE(t,p)
% LOTKA Lotka-Volterra predator-prey model
delta = 0.02;
beta = 0.01;

dpdt = [p(1) .* (1 - beta*p(2));
        p(2) .* (-1 + delta*p(1))];
end
%------------------------------------------------------------------
function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

참고 항목

ode45 | odeset

관련 항목

• ODE 솔버 선택하기
• 포식자-피식자 방정식 풀기