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

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

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

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

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

Для решения уравнений Лотка-Вольтерра в 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[\mathrm{}\mathrm{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-loop автоматизирует процесс решения.

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

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

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

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

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

Время использования каждого из предыдущих методов 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