В этом примере показано, как решить дифференциальное уравнение, представляющее модель хищника/добычи, использующую оба ode23
и ode45
. Эти функции являются для числового решения размера шага переменной использования обыкновенных дифференциальных уравнений методами интегрирования Рунге-Кутта. ode23
использует простую 2-ю и 3-ю пару порядка формул для средней точности и ode45
использует 4-ю и 5-ю пару порядка в более высокой точности.
Считайте пару обыкновенных дифференциальных уравнений первого порядка известной как уравнения Лотки - Вольтерры или модель добычи хищника:
Переменные и измерьте размеры добычи и популяций хищников, соответственно. Квадратичный перекрестный термин составляет взаимодействия между разновидностями. Популяция жертв увеличивается, когда никакие хищники не присутствуют, и уменьшения популяции хищников, когда добыча недостаточна.
Чтобы симулировать систему, создайте функцию, которая возвращает вектор-столбец производных состояния, учитывая состояние и временные стоимости. Эти две переменные и может быть представлен в MATLAB как первые два значения в векторном y
. Точно так же производные являются первыми двумя значениями в векторном yp
. Функция должна принять значения для t
и y
и возвратите значения, произведенные уравнениями в yp
.
yp(1) = (1 - alpha*y(2))*y(1)
yp(2) = (-1 + beta*y(1))*y(2)
В этом примере уравнения содержатся в файле под названием lotka.m
. Этот файл использует значения параметров и .
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
на интервале . Используйте начальное условие так, чтобы популяции хищников и добычи были равны.
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')
Результаты показывают, что решение дифференциальных уравнений с помощью различных численных методов может произвести немного отличающиеся ответы.