установка переменной состояния в заданное значение с помощью ode15s - PullRequest
Новая
       19

установка переменной состояния в заданное значение с помощью ode15s

0 голосов
/ 19 ноября 2018

Я использую ode15s для моделирования / решения набора ODE.Я хотел бы реализовать функцию, при которой после достижения заданного условия во время моделирования число в модели изменяется программно (например, константа индикатора) в течение фиксированного периода времени, а затем возвращается обратно.

Это может быть, например, использование уравнений Лотки-Вольтерра:

dx / dt = альфа x - бета x * y

dy / dt = (дельта + индикатор) x y - гамма y + epsilon индикатор

индикатор начинается с 0. Предположим, что когда x достигнет 10, я бы хотел переключить индикатор на 1 для 10единицы времени, а затем переверните его обратно на 0.

Это может быть сделано грязным способом с использованием глобальных переменных, однако я бы хотел этого избежать (невозможно распараллелить + общее избегание глобальных переменных).переменные).Есть ли альтернативный способ при использовании ode15s (то есть, я не знаю шаг по времени)?

Большое спасибо за любые предложения!

1 Ответ

0 голосов
/ 19 ноября 2018

Редактировать: Как правильно отмечено LutzL , перенос ODE с негладким состоянием без обработки событий может привести к неточным результатам

, поскольку выне может предсказать, в какие моменты времени и в каком порядке оценивается функция ODE. LutzL

Таким образом, точное решение состоит в том, чтобы иметь дело с событиями ODE .Пример для модифицированных уравнений Лотки-Вольтерры приведен ниже, где событие срабатывает, если x получает> 10 и индикатор будет включен на 10 секунд: 

% parameters and times:
params = ones(5,1); % [alpha, ..., epsilon]
z_start = [2, 1];
t_start = 0;
t_end = 30;

options = odeset('Events',@LotkaVolterraModEvents); % set further options here, too.

% wrap ODE function with indicator on and off
LVModODE_indicatorOn = @(t,z)LotkaVolterraModODE(t,z,1, params);
LVModODE_indicatorOff = @(t,z)LotkaVolterraModODE(t,z,0, params);

% storage for simulation values:
data.t = t_start;
data.z = z_start;
data.teout = [];
data.zeout = zeros(0,2);
data.ieout = [];

% temporary loop variables:
z_0 = z_start;
t_0 = t_start;
isIndicatorActive = false;

while data.t(end) < t_end % until the end time is reached
    if isIndicatorActive
        % integrate for 10 seconds, if indicator is active
        active_seconds = 10;
        [t, z, te,ze,ie] = ode15s(LVModODE_indicatorOn, [t_0 t_0+active_seconds], z_0, options);
    else
        % integrate until end or event, if indicator is not active.
        [t, z, te,ze,ie] = ode15s(LVModODE_indicatorOff, [t_0 t_end], z_0, options);
        isIndicatorActive = true;
    end

    %append data to storage
    t_len = length(t);
    data.t = [data.t; t(2:end)];
    data.z = [data.z; z(2:end,:)];
    data.teout = [data.teout; te];
    data.zeout = [data.zeout; ze];
    data.ieout = [data.ieout; ie];

    % reinitialize start values for next iteration of loop
    t_0 = t(end);
    z_0 = z(end, :);

    % set the length of the last instegration
    options = odeset(options,'InitialStep',t(end) - t(end-1));
end

%% plot your results:
figure;
plot(data.t, data.z(:,1), data.t, data.z(:,2));
hold all
plot(data.teout, data.zeout(:,1), 'ok');
legend('x','y', 'Events in x')

%% Function definitions for integration and events:
function z_dot = LotkaVolterraModODE(t, z, indicator, params)
    x = z(1); y= z(2);
    % state equations: modified Lotka-Volterra system
             z_dot =  [params(1)*x - params(2)*y;
                       (params(4) + indicator)*x*y - params(3)*y + params(5)*indicator];
end

function [value, isTerminal, direction] = LotkaVolterraModEvents(t,z)
    x = z(1);
    value = x-10; % event on rising edge when x passes 10
    isTerminal = 1; %stop integration -> has to be reinitialized from outer logic
    direction = 1; % only event on rising edge (i.e. x(t_e-)<10 and x(t_e+)>10)
end

Основная работа выполненав цикле while, где происходит интеграция.

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

Вы можете заключить вашу задачу в класс, который может содержать состояние (то есть его свойства).Класс должен иметь метод, который используется как odefun для интегратора с переменным шагом.См. Также здесь о том, как писать классы в MATLAB.

Пример ниже демонстрирует, как этого можно достичь в приведенном вами примере: 

% file: MyLotkaVolterra.m
classdef MyLotkaVolterra < handle
    properties(SetAccess=private)
        %define, if the modified equation is active
        indicator; 

        % define the start time, where the condition turned active.
        tStart; 

        % ode parameters [alpha, ..., epsilon]
        params;
    end

    methods
        function self = MyLotkaVolterra(alpha, beta, gamma, delta, epsilon)
            self.indicator = 0;
            self.tStart = 0;
            self.params = [alpha, beta, gamma, delta, epsilon];
        end

        % ODE funciton for the state z = [x;y] and time t
        function z_dot = odefun(self, t, z)
             x = z(1); y= z(2);
             if (x>=10 && ~self.indicator)
                 self.indicator = 1;
                 self.tStart = t;
             end

             %condition to turn indicator off:
             if (self.indicator && t - self.tStart >= 10)
                self.indicator = false; 
             end

             % state equations: modified Lotka-Volterra system
             z_dot =  [self.params(1)*x - self.params(2)*y;
                       (self.params(4) + self.indicator)*x*y - ...
                        self.params(3)*y + self.params(5)*self.indicator];
        end
    end
end

Этот классможно использовать следующим образом: 

% your ode using code:
% 1. create an object (`lvObj`) from the class with parameters alpha = ... = epsilon = 1
lvObj = MyLotkaVolterra(1, 1, 1, 1, 1);
% 2. pass its `odefun`method to the integrator (exaple call with ode15s)
[t,y] = ode15s(@lvObj.odefun, [0,5], [9;1]); % 5 seconds
...