matlab - Change a constant in ODE calculations under particular conditions with a flag

Question

I have an ODE for calculating how acidicity changes. Everything is working just fine, only I would like to change a constant whenever acidicity reaches a critical point. It is supposed to be some kind of irreversible effect I wish to simulate.

My constants are coming from a structure file (c) I load once in the ODE function.

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

The main problem I have here is not telling Matlab to change the constant but remember if it happened already during the simulation once. so Matlab should take the irreversibly changed constant rather than the one I supply at the beginning.

I would like to write a flag that is saved while the ODE is running and an if condition, "if flag exists, change constant". How to do that?

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.

Answer 1

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.