Output a variable each time step using ODE45

Hello,

I am trying to output a value from an ODE45 function at each time step, filling up an array, following this I need to read the array each time step as well. The idea is to alter a value in the ODE45 function using the condition of the previous time steps.

To write the value to be output I was attempting to use the Persistent variable and write it out to the main work space but not having much luck, is this the right approach I should be taking? I have included the start of the function below showing my so far fruitless attempts!

Thanks in advance for any assistance!

function y = eight_notches(t,P,T,~,~,x,ve,tfinal)

persistent h

i=t*1000;

h(i) = x;

assignin('base', 'h', h)

[Torsten]

For ODE45, there does not exist something you call "previous time step" because the solver adaptively changes its time step size.

Maybe you could explain what you ultimatively try to do and why you need variables from the previous time step.

[Steven Taggart]

Hi Torsten,

It is a safety relief valve model I have setup, after it is fully open the flow into the model is shut off, this involves having a condition in the ODE45 function that changes a variable to 0 after a set amount of time steps. To find the number of time steps I have to run the model several times to find the point at which it fully opens. This can take some time and isn't really productive. So instead of changing the variable to 0 after a set number of time steps I want to check whither the valve is fully open then switch the value to 0. So I need to store the current and previous values of the valve opening height to check wither it currently is fully open or has been in previous time steps.

What id like to do is fill an array in the main work space with the current valve opening height at each time step then check the full array to decide whither the in flow value should be set to 0.

I hope that clears it up slightly! :)

Cheers

Answer 1

[Jan]

The idea of a "previous time step" in the execution of ODE45 is not meaningful. Remember that the integrator adjusts the stepsize to run the calculationsw inside the given tolerances. If this tolerance is violated, a step is rejected. Such rejected steps would store the wanted parameter for a time in the future with a value which was rejected due to an bad accuracy.

In addition ODE45 evaluates the function to be integrated multiple times per step. Therefore a "previous step" is not clearly defined. Maybe you mean the "previously accepted time step". This is important. But even then the approach is not correct: ODE45 (and the other integrators of Matlab) can handle differentiable functions only. Using any parameter, which modifies the output in a non-smooth way, leads to unexpected behavior. If you are lucky, the integration stops with an error, because the the required step size is below the limits. With less luck, you get a final value, which is dominated by rounding errors, and then ODE45 is nothing but a weak random number generator. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047.

Filling the base workspace by assignin is a bad idea also. It suffers from the same problem as global variables. Please search in the forum or the net for "global variables problem".

The suggestion of Torsten is the only reliable and numerically correct way: Use event functions to stop the integration, when a certain parameter or value is reached. Then restart the integration with the changed parameters.

I've seen too many poor applications of integrators. Neither smoothness nor stiffness are considered and the results are not controlled by a variation of the initial values and parameters. Because an integration has the physical meaning of a simulation, a "final value" is not a "result", but you need an estimation of the accuracy in addition. You can e.g. simulate a pencil standing on the tip and find the trajectory, that it does not move in any direction. Of course this does not have a physical meaning. A trustworthy simulation needs to vary the initial position by e.g. +/- 1e-8 degree. Then the completely different final positions demonstrate that the system is instable.

Unfortunately Matlab offers the barbone integrators ODExy without mentioning these pitfalls in the documentation. It is even worse: The docs of ODE45 contain a non-smooth example:

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t);  % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t);  % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g;     % Evaluate ODE at time t

Bad. The linear interpolation of interp1 is not smooth. The best idea is to avoid any of these commands in the function to be integrated: if, max, round/fix/ceil/floor, abs, ... Even a tan is dangerous, except if you control the argument and stop with an error if it gets near to pi*(2n+1)/2.

But in spite of this fundamental rules of numerical maths, a brute integration of non-smooth functions without controlling the sensitivity occurs in many PhD theses and publications.

[Steven Taggart]

Hi Jan,

Thanks for the explanation, I will attempt as you have said, this is being used for a PhD thesis as well so I would prefer it to be correct!

Cheers

Steven

Answer 2

[Torsten]

Use the "event" facility of ODE45. It will help you to exactly locate the time when the valve is open.

Take a look at the "ballode" example.

Best wishes

Torsten.

[Steven Taggart]

Hi Torsten,

Thanks I will try that out!

Cheers

Steven.

[Jan]

+1. This is the only valid approach to handle discontinuities. I've elaborated this only a bit.