在特定条件下使用标志更改ODE计算中的常数 [英] Change a constant in ODE calculations under particular conditions with a flag

问题描述

我有一个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的答案，在这种情况下您必须小心，以免这样的步骤可能会导致不准确，并且您必须小心检查代码是否符合预期.

解决方案

要准确地执行此操作，您应该使用ode15s并继续积分.这看似粗糙，但它是如何准确完成这种事情的.

您可以在此处看到基本事件检测的示例.在命令窗口中键入ballode以查看其代码.您可以在我在这里的答案中看到该演示的稍微复杂一些的版本.这是一个示例，该示例使用事件在指定时间(而不是您指定状态值的情况)准确地更改ODE. /p>

注意:我感到奇怪的是，您将所谓的常量"(c6)作为第二个参数传递给

UPDATE I: PROBLEM SOLVED

Here a first own solution, it is not polished and requires a structure file approach. Which means, the constants which should suddenly changed upon an event, must be struct files which are handed in the ODE function into the function that should be evaluated (look above for the ODE syntax). The function accepts the inputs like this:

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

Please look into Horchlers kind answer where you have to take care that such a step may introduce inaccuracies and you have to be careful to check if your code does what it is supposed to do.

解决方案

To do this accurately, you should use event detection within the ODE solver. I can't give you a specific answer because you've only provided the ode15s call it in your question, but you'll need to write an events function and then specify it via odeset. Something like this:

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

You'll want to use the events to integrate right to the point where the "acidicity reaches a critical point" and stop the integration. Then call ode15s again with the new value and continue the integration. This may seem crude, but it how this sort of thing can be done accurately.

You can see an example of basic event detection here. Type ballode in your command window to see the code for this. You can see a slightly more complex version of this demo in my answer here. Here's an example of using events to accurately change an ODE at specified times (rather than your case of specified state values).

Note: I find it strange that you're passing what you call "constants", c, as the second argument to ode15s. This function has strict input argument requirements: the first is the independent variable (often time), and the second is the array of state variables (same as your initial condition vector). Also if f1 only takes two arguments, @(x,c)f1(x,c) is superfluous – it's sufficient to pass in @f1.

这篇关于在特定条件下使用标志更改ODE计算中的常数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！