1

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.

Easyquestionsonly
  • 103
  • 6

1 Answers1

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.

Community
  • 1
  • 1
horchler
  • 18,384
  • 4
  • 37
  • 73
  • Thank you, I feared it might be that I have to call two ODE functions. I tried to avoid it though, there is no easier way? Im trying not to make an already complicated script more complicated. – Easyquestionsonly Feb 11 '15 at 13:03
  • The constants im passing is necessary because I needed to make the constants editable in an excel sheet, import as variables in matlab workspace and finally turn them into structs to pass into a function. The whole reason was to avoid having to use matlab in editing the script and to avoid globals (those constants had to be known by at least two functions) which ate a lot of memory. See here: http://stackoverflow.com/questions/27218341/what-is-the-best-way-to-pass-lots-of-constants-and-workspace-variables-to-an-fso – Easyquestionsonly Feb 11 '15 at 13:06
  • 1
    @Easyquestionsonly: There's no other way unless you want inaccurate and possibly completely wrong results. Yes, you [should not use globals](http://mathworks.com/matlabcentral/answers/51946-systematic-do-not-use-global-don-t-use-eval). The link you gave is for `fsolve` –the objective function has only one required input. The ODE solvers, like `ode15s`, have two. If `c` is a fixed parameter or a set of them, then they should be passed as the third argument. See [this answer](http://stackoverflow.com/a/16681767/2278029) of mine where `n` is passed as a parameter to the integration function `f`. – horchler Feb 11 '15 at 18:20
  • I managed to find another way to do it. It has the general form where the function which ODE is loading has a load c.mat command at the beginning. this c. struct file contains my variables. Once i pass a critical value in the function i save c.mat, else i do nothing. When the function is called in the same loop again, the first thing it sees is load c.mat with the changed value. Like this i have my one event saved in the same ODE loop. – Easyquestionsonly Feb 12 '15 at 16:19
  • 1
    @Easyquestionsonly: I don't understand that because you've not provided any code, but if it works and you can validate that your results are accurate and correct: excellent. The key thing is that you don't introduce discontinuities, steps, `if` statements, etc. into your integration function, the function that describes your ODEs and is called by `ode15s` (`f1` in the case of your question). These will cause your system to become much stiffer, reduce accuracy, and possibly cause completely wrong results. – horchler Feb 12 '15 at 16:41
  • Thank you for pointing this out. The main point is to introduce some kind of perturbation (actually to react to a perturbation), in a quite dynamic environment. Artifacts can be easily recognized in my case but maybe not in yours. – Easyquestionsonly Feb 12 '15 at 18:47