使用标志在特定条件下更改 ODE 计算中的常数

Question

我有一个用于计算酸度如何变化的 ODE。一切都很好，只是我想在酸度达到临界点时改变一个常数。这应该是我想模拟的某种不可逆的效果。

我的常量来自一个结构文件 (c) 我在 ODE 函数中加载过一次。

[Time,Results] = ode15s(@(x, c) f1(x, c),[0 c.length],x0,options);

我在这里遇到的主要问题不是告诉 Matlab 更改常数，而是记住如果它已经在模拟期间发生过一次。所以 Matlab 应该采用不可逆变化的常数，而不是我一开始提供的那个。

我想编写一个在 ODE 运行时保存的标志和一个if 条件，“如果标志存在，则更改常量”。该怎么做？

更新一： 问题解决了

这是第一个自己的解决方案，它没有经过打磨，需要结构文件的方法。这意味着，应该在事件中突然改变的常量必须是结构文件，这些文件在 ODE 函数中传递到应该评估的函数中（查看上面的 ODE 语法）。该函数接受如下输入：

function [OUTPUT] = f1(t, x, c)

% So here, the constants all start with c. followed by the variable name. (because they are structs instead of globals)

%% write a flag when that event happens

if c.someODEevent <= 999 && exist ('flag.txt')  == 0
    dlmwrite ('flag.txt',1);
end

%% next iteration will either write the flag again or not. more importantly, if it exists, the constant will change because of this.

if exist ('flag.txt')  == 2
        c.changingconstant = c.changingconstant/2;
end



end

请查看 Horchlers 的友好回答，您必须注意这样的步骤可能会引入不准确之处，并且您必须小心检查您的代码是否完成了它应该做的事情。

Answer 1

要准确地做到这一点，您应该在 ODE 求解器中使用 event detection。我不能给你一个具体的答案，因为你只在你的问题中提供了ode15s 调用它，但你需要编写一个events 函数，然后通过odeset 指定它。像这样的：

function acidity_main
% load c data
...
x0 = ...
options = odeset('Events',@events); % add any other options too

% integrate until critical value and stop
[Time1,Results1] = ode15s(@(x,c)f1(x,c),[0 c.length],x0,options);

x0 = Results(end,:); % set new initial conditions
% pass new parameters -it's not clear how you're passing parameters
% if desired, change options to turn off events for faster integration
[Time2,Results2] = ode15s(@(x,c)f1(x,c),[0 c.length],x0,options);

% append outputs, last of 1 is same as first of 2
Time = [Time1;Time2(2:end)];
Results = [Results1;Results2(2:end,:)];
...

function y=f1(x,c)
% your integration function
...

function [value,isterminal,direction] = events(x,c)
value = ... % crossing condition, evaluates to zero at event condition
isterminal = 1; % stop integration when event detected
direction = ... % see documentation

您需要使用事件来整合到“酸度达到临界点”的点并停止整合。然后使用新值再次调用ode15s 并继续集成。这可能看起来很粗糙，但它是如何准确地完成这种事情的。

您可以查看基本事件检测示例here。在命令窗口中输入 ballode 以查看此代码。您可以在my answer here 中看到这个演示的稍微复杂的版本。这是an example 使用事件在指定时间准确更改 ODE（而不是指定状态值的情况）。

注意：我觉得奇怪的是，您将所谓的“常量”c 作为第二个参数传递给ode15s。该函数有严格的输入参数要求：第一个是自变量（通常是时间），第二个是状态变量数组（与您的初始条件向量相同）。此外，如果f1 只接受两个参数，@(x,c)f1(x,c) 是多余的——传入@f1 就足够了。