Уравнения кода

Чтобы смоделировать систему, создайте функцию, которая возвращает вектор столбца производных состояний с заданными значениями состояния и времени. Две переменные $\mathit{x}$ и $\mathit{y}$ могут быть представлены в MATLAB как первые два значения в векторе y. Аналогично, производные являются первыми двумя значениями в векторе. yp. Функция должна принимать значения для t и y и возвращает значения, полученные из уравнений в yp.

yp(1) = (1 - alpha*y(2))*y(1)

yp(2) = (-1 + beta*y(1))*y(2)

В этом примере уравнения содержатся в файле с именем lotka.m. В этом файле используются значения параметров $\alpha \mathrm{=}\mathrm{}$0,01 $и\mathrm{\beta }\mathrm{=}$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;

Моделирование системы

Использовать ode23 для решения дифференциального уравнения, определенного в lotka за интервал $0<\mathrm{t}$< 15. Используйте начальное условие $\mathit{}\mathrm{}x\mathit{}\left(0\right)\mathrm{=}$ y (0) = 20, чтобы популяции хищников и добычи были равны.

t0 = 0;
tfinal = 15;
y0 = [20; 20];   
[t,y] = ode23(@lotka,[t0 tfinal],y0);

Результаты графика

Постройте график результирующих популяций по времени.

plot(t,y)
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')

Сравнение результатов различных решателей

Повторное решение системы с помощью ode45, вместо ode23. ode45 решатель занимает больше времени для каждого шага, но он также занимает больше шагов. Тем не менее, выход ode45 является гладким, потому что по умолчанию решатель использует формулу непрерывного расширения для получения выходных данных в четырех равномерно разнесенных точках времени в диапазоне каждого выполненного шага. (Можно настроить количество точек с помощью 'Refine' вариант.) Постройте график обоих решений для сравнения.

[T,Y] = ode45(@lotka,[t0 tfinal],y0);

plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'-');
title('Phase Plane Plot')
legend('ode23','ode45')

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