Решите систему ОДУ с несколькими начальными условиями

Этот пример сравнивает два метода, чтобы решить систему обыкновенных дифференциальных уравнений с несколькими наборами начальных условий. Методы:

  • Используйте for- цикл, чтобы выполнить несколько симуляций, один для каждого набора начальных условий. Этот метод прост в использовании, но не предлагает лучшую эффективность для больших систем.

  • Векторизуйте функцию ОДУ, чтобы решить систему уравнений для всех наборов начальных условий одновременно. Этот метод является более быстрым методом для больших систем, но требует перезаписи функции ОДУ так, чтобы это изменило входные параметры правильно.

Уравнения, используемые, чтобы продемонстрировать эти методы, являются известными уравнениями Lotka-Volterra, которые являются нелинейными дифференциальными уравнениями первого порядка, которые описывают популяции хищников и добычи.

Описание проблемы

Уравнения Lotka-Volterra являются системой двух первых порядков, нелинейные ОДУ, которые описывают популяции хищников и добычи в биологической системе. В зависимости от времени популяции хищников и добычи изменяются согласно уравнениям

dxdt=αx-βxy,dydt=δxy-γy.

Переменные в этих уравнениях

  • x численность населения добычи

  • y численность населения хищников

  • t время

  • α, β, δ, и γ постоянные параметры, которые описывают взаимодействия между двумя разновидностями. Этот пример использует значения параметров α=γ=1, β=0.01, и δ=0.02.

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

Решите уравнения с одним начальным условием

Чтобы решить уравнения Lotka-Volterra в MATLAB, запишите функцию, которая кодирует уравнения, задайте временной интервал для интегрирования и задайте начальные условия. Затем можно использовать один из решателей ОДУ, таких как ode45, симулировать систему в зависимости от времени.

Функция, которая кодирует уравнения,

function dpdt = lotkaODE(t,p)
% LOTKA Lotka-Volterra predator-prey model
delta = 0.02;
beta = 0.01;

dpdt = [p(1) .* (1 - beta*p(2));
        p(2) .* (-1 + delta*p(1))];
end

(Эта функция включена как локальная функция в конце примера.)

С тех пор существует два уравнения в системе, dpdt вектор с одним элементом для каждого уравнения. Кроме того, вектор решения p имеет один элемент для каждого компонента решения: p(1) представляет x в исходных уравнениях и p(2) представляет y в исходных уравнениях.

Затем задайте временной интервал для интегрирования как [0,15] и набор размеры начальной генеральной совокупности для x и y к 50.

t0 = 0;
tfinal = 15;
p0 = [50; 50];

Решите систему с ode45 путем определения функции ОДУ, отрезка времени и начальных условий. Постройте получившиеся популяции по сравнению со временем.

[t,p] = ode45(@lotkaODE,[t0 tfinal],p0);
plot(t,p)
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
legend('Prey','Predators')

Поскольку решения показывают периодичность, строят решения друг против друга в графике фазы.

plot(p(:,1),p(:,2))
title('Phase Plot of Predator/Prey Populations')
xlabel('Prey')
ylabel('Predators')

Получившиеся графики показывают решение для данных размеров начальной генеральной совокупности. Чтобы решить уравнения для различных размеров начальной генеральной совокупности, измените значения в p0 и повторно выполненный симуляция. Однако этот метод только решает уравнения для одного начального условия за один раз. Следующие два раздела описывают методы, чтобы решить для многих различных начальных условий.

Метод 1: вычислите несколько начальных условий с for-цикл

Самый простой способ решить систему ОДУ для нескольких начальных условий с for- цикл. Этот метод использует ту же функцию ОДУ в качестве одного начального метода условия, но for- цикл автоматизирует процесс решения.

Например, можно содержать размер начальной генеральной совокупности для x постоянный в 50, и использование for- цикл, чтобы варьироваться размер начальной генеральной совокупности для y между 10 и 400. Создайте вектор из численностей населения для y0, и затем цикл по значениям, чтобы решить уравнения для каждого набора начальных условий. Постройте график фазы с результатами всех итераций.

y0 = 10:10:400;
for k = 1:length(y0)
    [t,p] = ode45(@lotkaODE,[t0 tfinal],[50 y0(k)]);
    plot(p(:,1),p(:,2))
    hold on
end
title('Phase Plot of Predator/Prey Populations')
xlabel('Prey')
ylabel('Predators')
hold off

График фазы показывает все вычисленные решения для различных наборов начальных условий.

Метод 2: вычислите несколько начальных условий с векторизованной функцией ОДУ

Другой метод, чтобы решить систему ОДУ для нескольких начальных условий должен переписать функцию ОДУ так, чтобы все уравнения были решены одновременно. Шаги, чтобы сделать это:

  • Обеспечьте все начальные условия к ode45 как матрица. Размером матрицы является s- n, где s количество компонентов решения и n количество начальных условий, решаемых для. Каждый столбец в матрице затем представляет один полный набор начальных условий для системы.

  • Функция ОДУ должна принять дополнительный входной параметр для n, количество начальных условий.

  • В функции ОДУ решатель передает компоненты решения p как вектор-столбец. Функция ОДУ должна изменить вектор в матрицу с размером s- n. Каждая строка матрицы затем содержит все начальные условия для каждой переменной.

  • Функция ОДУ должна решить уравнения в векторизованном формате, так, чтобы выражение приняло векторы для компонентов решения. Другими словами, f(t,[y1 y2 y3 ...]) должен возвратить [f(t,y1) f(t,y2) f(t,y3) ...].

  • Наконец, функция ОДУ должна изменить свой выход назад в вектор так, чтобы решатель ОДУ получил вектор назад от каждого вызова функции.

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

Чтобы реализовать этот метод для системы Лотки - Вольтерры, запустите путем нахождения количества начальных условий n, и затем сформируйте матрицу начальных условий.

n = length(y0);
p0_all = [50*ones(n,1) y0(:)]';

Затем перепишите функцию ОДУ так, чтобы она приняла n как вход. Используйте n чтобы изменить вектор решения в матрицу, затем решите векторизованную систему и измените выход назад в вектор. Модифицированная функция ОДУ, которая выполняет эти задачи,

function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

Решите систему уравнений для всех начальных условий с помощью ode45. Начиная с ode45 требует, чтобы функция ОДУ приняла два входных параметров, использовала анонимную функцию, чтобы передать в значении n от рабочей области до lotkasystem.

[t,p] = ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all);

Измените выходной вектор в матрицу с размером (numTimeSteps*s)- n. Каждый столбец выхода p(:,k) содержит решения для одного набора начальных условий. Постройте график фазы компонентов решения.

p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:end,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off

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

  • Расчетные решения могут немного отличаться, чем вычисленные из одного начального входа. Различие возникает, потому что решатель ОДУ применяет проверки нормы к целой системе, чтобы вычислить размер временных шагов, таким образом, продвигающееся во время поведение решения немного отличается. Изменение во временных шагах обычно не влияет на точность решения, а скорее в каких временах решение оценено.

  • Для жестких решателей ОДУ (ode15sode23sode23tode23tb) это автоматически оценивает числовой якобиан системы, задавая шаблон разреженности диагонали блока якобиана с помощью JPattern опция odeset может повысить эффективность вычисления. Форма диагонали блока якобиана является результатом входного изменения, выполняемого в переписанной функции ОДУ.

Сравните результаты синхронизации

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

% Time one IC
baseline = timeit(@() ode45(@lotkaODE,[t0 tfinal],p0),2);

% Time for-loop
for k = 1:length(y0)
    loop_timing(k) = timeit(@() ode45(@lotkaODE,[t0 tfinal],[50 y0(k)]),2);
end
loop_timing = sum(loop_timing);

% Time vectorized fcn
vectorized_timing = timeit(@() ode45(@(t,p) lotkasystem(t,p,n),[t0 tfinal],p0_all),2);

Составьте таблицу результатами синхронизации. Умножьте все результаты 1e3, чтобы описать времена в миллисекундах. Включайте столбец со временем на решение, которое делится каждый раз на количество начальных условий, решаемых для.

TimingTable = table(1e3.*[baseline; loop_timing; vectorized_timing], 1e3.*[baseline; loop_timing/n; vectorized_timing/n],...
    'VariableNames',{'TotalTime (ms)','TimePerSolution (ms)'},'RowNames',{'One IC','Multi ICs: For-loop', 'Mult ICs: Vectorized'})
TimingTable=3×2 table
                            TotalTime (ms)    TimePerSolution (ms)
                            ______________    ____________________

    One IC                      1.4634                1.4634      
    Multi ICs: For-loop         46.113                1.1528      
    Mult ICs: Vectorized        2.5454              0.063636      

TimePerSolution столбец показывает, что векторизованный метод является самым быстрым из этих трех методов.

Локальные функции

Перечисленный здесь локальные функции что ode45 вызовы, чтобы вычислить решения.

function dpdt = lotkaODE(t,p)
% LOTKA Lotka-Volterra predator-prey model
delta = 0.02;
beta = 0.01;

dpdt = [p(1) .* (1 - beta*p(2));
        p(2) .* (-1 + delta*p(1))];
end
%------------------------------------------------------------------
function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

Смотрите также

|

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте