Описание проблемы

Уравнения Lotka-Volterra являются системой двух первых порядков, нелинейные ОДУ, которые описывают популяции хищников и добычи в биологической системе. В зависимости от времени популяции хищников и добычи изменяются согласно уравнениям

Переменные в этих уравнениях

Для этой проблемы, начальных значений для $\mathit{x}$ и $\mathit{y}$ размеры начальной генеральной совокупности. Решение уравнений затем предоставляет информацию о том, как популяции изменяются в зависимости от времени, когда разновидности взаимодействуют.

Решите уравнения с одним начальным условием

Чтобы решить уравнения Lotka-Volterra в MATLAB, запишите функцию, которая кодирует уравнения, задайте временной интервал для интегрирования и задайте начальные условия. Затем можно использовать один из решателей ОДУ, таких как 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];

Решите систему с 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-цикл

Самый простой способ решить систему ОДУ для нескольких начальных условий с for- цикл. Этот метод использует ту же функцию ОДУ в качестве одного начального метода условия, но 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: вычислите несколько начальных условий с векторизованной функцией ОДУ

Другой метод, чтобы решить систему ОДУ для нескольких начальных условий должен переписать функцию ОДУ так, чтобы все уравнения были решены одновременно. Шаги, чтобы сделать это:

Если вы выполняете эти шаги, то решатель ОДУ может решить систему уравнений с помощью вектора для компонентов решения, в то время как функция ОДУ изменяет вектор в матрицу и решает каждый компонент решения для всех начальных условий. Результат состоит в том, что можно решить систему для всех начальных условий в одной симуляции.

Чтобы реализовать этот метод для системы Лотки - Вольтерры, запустите путем нахождения количества начальных условий n, и затем сформируйте матрицу начальных условий.

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

Затем перепишите функцию ОДУ так, чтобы она приняла n как вход. Используйте n чтобы изменить вектор решения в матрицу, затем решите векторизованную систему и измените выход назад в вектор. Модифицированная функция ОДУ, которая выполняет эти задачи,

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 требует, чтобы функция ОДУ приняла два входных параметров, использовала анонимную функцию, чтобы передать в значении 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- метод цикла. Однако существуют некоторые свойства векторизованного метода решения, который необходимо иметь в виду:

Сравните результаты синхронизации

Время каждый из предыдущих методов с помощью 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.69141               0.69141      
    Multi ICs: For-loop         24.206               0.60516      
    Mult ICs: Vectorized        1.8598              0.046496      

TimePerSolution столбец показывает, что векторизованный метод является самым быстрым из этих трех методов.

Локальные функции

Перечисленный здесь локальные функции что 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