В этом примере показано, как решить дифференциальное уравнение, представляющее модель хищника/добычи, использующую оба 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Используйте 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')
Результаты показывают, что решение дифференциальных уравнений с помощью различных численных методов может произвести немного отличающиеся ответы.