使用 ode15s 将状态变量设置为给定值

Question

我正在使用 ode15s 来 simulate/solve 一组常微分方程。我想实现一个功能，在模拟过程中达到给定条件后，模型中的数字会以编程方式更改（例如，指标常量）一段固定的时间，然后恢复原状。

这可能是，例如使用 Lotka-Volterra 方程：

dx/dt = alphax - betax*y

dy/dt = (delta+指标)xy - gammay + epsilon指标

指标从 0 开始。假设当 x 达到 10 时，我想将指标切换为 1 10 个时间单位，然后将其翻转回 0。

这可以通过使用全局变量以肮脏的方式完成，但是，这是我想避免的事情（不可能并行化 + 一般避免使用全局变量）。使用 ode15s 时是否有一种巧妙的替代方法（即，我不知道时间步长）？

非常感谢您的任何建议！

Answer 1

编辑： 正如 LutzL 正确指出的那样，在不处理事件的情况下用非平滑状态包装 ODE 可能会导致不准确的结果

as you can not predict at what time points in what order the ODE function is evaluated. LutzL

所以准确的解决办法是处理ODE events。下面给出了修改后的 Lotka-Volterra 方程的示例，如果 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 循环中完成，在该循环中进行积分。

(Old post)以下解决方案可能会导致结果不准确，处理事件时，应该优先采用第一部分的解释。

您可以将您的问题包装在一个 class 中，它能够保持状态（即它的属性）。 class 应该有一个方法，用作可变步长积分器的 odefun。另请参阅 here 了解如何在 MATLAB 中编写 classes。

下面的示例演示了如何实现您提供的示例：

% 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

这个class可以按如下方式使用：

% 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

使用 ode15s 将状态变量设置为给定值

setting state variable to a given value using ode15s

matlab

ode