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

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

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

Таким образом, правильным решением будет обработка событий 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