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

Чтобы симулировать систему, создайте функцию, которая возвращает вектор-столбец производных состояния, учитывая состояние и временные стоимости. Эти две переменные $\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. 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')

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