Main Content

포식자-피식자 방정식 풀기

이 예제에서는 ode23ode45를 모두 사용하여 포식자/피식자 모델을 나타내는 미분 방정식을 푸는 방법을 보여줍니다. ode23과 ode45는 가변 스텝 크기 룽게-쿠타 적분 방법을 사용하여 상미분 방정식의 수치 해를 구하기 위한 함수입니다. ode23은 중간 정도의 정확도를 얻기 위해 단순하게 식의 2차와 3차 항을 사용하고, ode45는 더 높은 정확도를 얻기 위해 4차와 5차 항을 사용합니다.

로트카-볼테라 방정식 또는 포식자-피식자 모델로 알려진 1계 상미분 방정식 쌍을 사용해 보겠습니다.

dxdt=x-αxydydt=-y+βxy.

변수 xy는 각각 피식자 개체군과 포식자 개체군의 크기를 측정합니다. 2차 교차 항은 종 간의 상호 작용을 설명합니다. 포식자가 없으면 피식자 개체군이 증가하고 피식자가 부족하면 포식자 개체군이 감소합니다.

방정식 코딩하기

시스템을 시뮬레이션하기 위해, 상태 값과 시간 값이 제공될 경우 상태 도함수로 구성된 열 벡터를 반환하는 함수를 만듭니다. MATLAB®에서 두 변수 xy는 벡터 y의 처음 2개 값으로 나타낼 수 있습니다. 마찬가지로 도함수 역시 벡터 yp의 처음 2개 값입니다. 함수는 ty에 대한 값을 받은 후 yp의 방정식에서 생성한 값을 반환해야 합니다.

yp(1) = (1 - alpha*y(2))*y(1)

yp(2) = (-1 + beta*y(1))*y(2)

이 예제에서 방정식은 lotka.m이라는 파일에 포함되어 있습니다. 이 파일은 파라미터 값 α=0.01β=0.02를 사용합니다.

type lotka
function yp = lotka(t,y)
%LOTKA  Lotka-Volterra predator-prey model.

%   Copyright 1984-2014 The MathWorks, Inc.

yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

시스템 시뮬레이션하기

ode23을 사용하여 구간 0<t<15에 대해 lotka에 정의된 미분 방정식을 풉니다. 포식자 개체군과 피식자 개체군이 동일하도록 초기 조건 x(0)=y(0)=20을 사용합니다.

t0 = 0;
tfinal = 15;
y0 = [20; 20];   
[t,y] = ode23(@lotka,[t0 tfinal],y0);

결과 플로팅하기

시간에 대해 결과로 나타나는 개체군을 플로팅합니다.

plot(t,y)
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
legend('Prey','Predators','Location','North')

Figure contains an axes object. The axes object with title Predator/Prey Populations Over Time contains 2 objects of type line. These objects represent Prey, Predators.

이제 개체군 간을 비교해 플로팅합니다. 결과로 나타나는 위상 평면 플롯에서는 개체군 간 순환 관계가 훨씬 명확해집니다.

plot(y(:,1),y(:,2))
title('Phase Plane Plot')
xlabel('Prey Population')
ylabel('Predator Population')

Figure contains an axes object. The axes object with title Phase Plane Plot contains an object of type line.

여러 솔버의 결과 비교하기

ode23 대신 ode45를 사용하여 다시 시스템을 풉니다. ode45 솔버는 각 스텝마다 시간이 더 오래 걸릴 뿐만 아니라 더 많은 스텝을 진행합니다. 그렇지만 ode45의 출력값은 매끄럽게 표시됩니다. 이는 기본적으로 솔버가 연속적 확장 식을 사용하여 각 스텝 범위 안에서 시간 지점을 4개의 균일한 간격으로 나눠 출력값을 생성하기 때문입니다. ('Refine' 옵션을 사용하면 지점 개수를 조정할 수 있습니다.) 비교를 위해 2개의 해를 모두 플로팅합니다.

[T,Y] = ode45(@lotka,[t0 tfinal],y0);

plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'-');
title('Phase Plane Plot')
legend('ode23','ode45')

Figure contains an axes object. The axes object with title Phase Plane Plot contains 2 objects of type line. These objects represent ode23, ode45.

이 결과는 서로 다른 수치 계산법을 사용하여 미분 방정식을 풀어도 생성되는 답은 약간의 차이만 있음을 보여줍니다.

참고 항목

|

관련 항목