Описание задачи

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

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

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

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

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

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

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

Затем перепишите функцию ODE так, чтобы она приняла 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 требует, чтобы функция ODE приняла два входов, использовала анонимную функцию, чтобы пройти в значении n из рабочей области в lotkasystem.

[t,p] = ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all);

Измените форму вектора выхода на матрицу с (numTimeSteps*s) размера-by- 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      

The 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