Кодовые уравнения

Чтобы симулировать систему, создайте функцию, которая возвращает вектора-столбца производных состояний, заданных состояний и временных значений. Эти две переменные $\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 =0.01$ и $\beta =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<t<15$. Используйте начальное условие $\mathit{x}\left(0\right)=\mathit{y}\left(0\right)=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. The 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')

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