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

Question

我正在使用ode15s来模拟/求解一组ODE。我想实现一个功能，在模拟过程中达到给定条件时，模型中的一个数字会在固定的时间内以程序方式改变(例如，指标常数(，然后返回。

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

dx/dt=αx-βx*y

dy/dt=(delta+指示剂(xy-γy+ε指示剂

指示器以0开头。比方说，当x达到10时，我想将指示器切换到1，持续10个时间单位，然后将其翻转回0。

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

非常感谢您的建议！

Answer 1

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

因为您无法预测ODE在什么时间点以什么顺序函数进行评估。LutzL

因此，准确的解决方案是处理ODE事件。下面给出了修改后的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循环中完成，在那里进行集成。

---

(旧帖(以下解决方案可能会导致不准确的结果，应优先处理第一部分中解释的事件。

您可以将问题封装在一个类中，该类能够保存一个状态(即其属性(。类应该有一个方法，该方法用作可变步长积分器的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