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

라이브 스크립트 열기

이 예제에서는 가변 스텝 크기 룽게-쿠타 적분 방법을 사용하여 포식자/피식자 모델을 나타내는 미분 방정식을 푸는 방법을 보여줍니다. 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]
         EquationType: standard

   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")

Figure contains an axes object. The axes object with title Predator/Prey Populations Over Time, xlabel t, ylabel Population 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, xlabel Prey Population, ylabel Predator Population contains an object of type line.

여러 솔버의 결과 비교하기

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")

Figure contains an axes object. The axes object with title Phase Plane Plot contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent ode23, ode45.

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

참고 항목

ode45 | ode23