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

라이브 스크립트 열기

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

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

$\begin{array}{l} \frac{dx}{dt} = x - α xy \\ \frac{dy}{dt} = - y + β xy . \end{array}$

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

방정식 코딩하기

시스템을 시뮬레이션하기 위해, 상태 값과 시간 값이 제공될 경우 상태 도함수로 구성된 열 벡터를 반환하는 함수를 만듭니다. MATLAB®에서 두 변수 $x$ 와 $y$ 는 벡터 y의 처음 2개 값으로 나타낼 수 있습니다. 마찬가지로 도함수 역시 벡터 yp의 처음 2개 값입니다. 함수는 t와 y에 대한 값을 받은 후 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;

시스템 시뮬레이션하기

문제와 초기 조건을 정의하기 위한 ode 객체를 만듭니다. 초기 포식자 개체군과 피식자 개체군이 동일하도록 초기 조건 $x (0) = y (0) = 20$ 을 사용합니다. ode23 솔버를 사용하도록 지정합니다. 그런 다음 solve 메서드를 사용하여 시간 간격 $0 < t < 15$ 에 대해 시스템을 시뮬레이션합니다.

y0 = [20; 20];   
F = ode(ODEFcn=@lotka,InitialValue=y0,Solver="ode23")

F = 
  ode with properties:

   Problem definition
               ODEFcn: @lotka
          InitialTime: 0
         InitialValue: [2x1 double]

   Solver properties
    AbsoluteTolerance: 1.0000e-06
    RelativeTolerance: 1.0000e-03
               Solver: ode23

  Show all properties

t0 = 0;
tfinal = 15;
S = solve(F,t0,tfinal);
t = S.Time;
y = S.Solution;

결과 플로팅하기

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

plot(t,y,"-o")
title("Predator/Prey Populations Over Time")
xlabel("t")
ylabel("Population")
legend("Prey","Predators",Location="north")

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

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

여러 솔버의 결과 비교하기

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

F.Solver = "ode45";
S2 = solve(F,t0,tfinal);
t2 = S2.Time;
y2 = S2.Solution;

plot(y(1,:),y(2,:),"o",y2(1,:),y2(2,:),".");
title("Phase Plane Plot")
legend("ode23","ode45")

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

참고 항목

ode45 | ode23