通过在 Matlab 中计算 ODE 函数中的 ODE 值来检测稳态答案

【问题标题】：Detecting steady-state by calculating ODE values in ODE function in Matlab通过在 Matlab 中计算 ODE 函数中的 ODE 值来检测稳态
【发布时间】：2016-08-09 05:29:51
【问题描述】：

我有一个要求解的 ODE 方程组，但有一个棘手的部分是，当系统达到稳态时，我想更改一个（或多个）参数的值。例如，考虑以下情况：

function dydt = diff(t,x,params)
   F = params(1);
   G = params(2);
   dydt = zeros(2,1);
   dydt(1) = F*x(1) - G*x(1)*x(2);
   dydt(2) = (F-G)*x(2);
end

例如，我希望我的代码能够在系统达到稳态时，将 F 的值更改为 10，将 G 的值更改为 2。我正在考虑使用例如

来检测 dydt(1) 和 dydt(2) 的值

if norm(dydt)<1
   F = 10;
   G = 2;
end

如何在 Matlab 中对 ODE 表达式执行此操作？如果我将这个 if 条件放在 ODE 表达式之前，则 dydt 的值将始终为零。但是如果我把这个 If 条件放在 ODE 表达式之后，If 条件将无法更正 ODE 表达式。

谢谢！

【问题讨论】：

执行此操作的正确方法是使用event location 根据您的稳态条件 (example here) 停止集成。然后更改您的参数并开始新的集成。避免在 ODE 函数中添加 if 语句来更改参数。
我考虑了事件位置，但这是否意味着我必须停止集成，然后重新开始一个新的集成过程？如果是这样，我是否必须在找到稳态之前将初始值更新为第一个积分步骤的最后一点？

标签： matlab if-statement conditional-statements ode

【解决方案1】：

假设 ODE 的参数是固定的，并且不依赖于系统的状态。您要做的是模拟分段连续 ODE。这类问题通常用event location 解决（假设存在解决方案）。您需要在关键点停止模拟，更改参数，然后重新开始新的模拟，初始条件与上一次结束时相同。

以下是您的 ODE 函数的一个小示例。我不知道你的初始条件，初始参数值，或者其他设置，所以这里只是为了演示这个方案：

function eventsdemo
params = [-1.5 1];
tspan = [0 10];
x0 = [1;1];
opts = odeset('Events',@events);
[t1,x1] = ode45(@(t,x)f(t,x,params),tspan,x0,opts); % Simulate with events

% Change parameters, set initial conditions based on end of previous run
params = [1.5 1];
x0 = x1(end,:);
tspan = [t1(end) 10];
[t2,x2] = ode45(@(t,x)f(t,x,params),tspan,x0); % Simulate again

% Concatenate results, removing duplicate points
t = [t1;t2(2:end)];
x = [x1;x2(2:end,:)];

figure;
plot(t,x);
hold on;
plot(t2(1),x2(1,:),'k*'); % Plot event location


function dxdt = f(t,x,params) %#ok<INUSL>
F = params(1);
G = params(2);
dxdt = [F*x(1) - G*x(1)*x(2);
        (F-G)*x(2)];


function [value,isterminal,direction] = events(t,x) %#ok<INUSL>
value = norm(x)-1e-3; % Don't try to detect exact zero for asymptotics
isterminal = true;
direction = -1;

在这种情况下，解会渐近地逼近 (0,0) 处的稳定不动点。重要的是不要试图准确地检测 (0,0)，因为解决方案可能永远不会达到那个点。您可以改为使用较小的容差。但是，根据您的系统，选择此容差可能会影响您更改参数后的行为。您还可以考虑将此问题重新表述为边值问题。我不知道你想用这个系统做什么，所以我不能说太多（这可能与这个网站无关）。

【讨论】：