Creating a counter inside the ode45 function

19 views (last 30 days)
Chad on 8 Sep 2021
Commented: Jan on 13 Sep 2021
I would like to set up a simple counter inside the ode45 function for the purpose of changing a variable which affects the ode. A simple example:
function xdot = x(t,x,A)
a = A(:,counter);
xdot = a*x
My question is how can I create a variable "counter" where it starts at 1 and increases by 1 with each time-step until ode45 is done solving? My input "t" is usually non-integer, but does have specifice time-steps (for example: t=[0.01:0.01:1]). I would like to avoid using "t" to create "counter" if possible.

Answers (3)

Star Strider
Star Strider on 8 Sep 2021
See the ode45 documentation section on ODE with Time-Dependent Terms since that is most likely what you actually want to do. (It is possible to create ‘A’ as an external function of ‘t’ and pass it as an argument to the ODE function if necessary.) The interp1 function has a number of interpolation methods, including 'nearest', 'previous', and 'next' so it is not always necessary to do continuous interpolation if the intent is to use ‘step-wise’ interpolation.
Jan on 13 Sep 2021
@Paul: Some members of the Mathworks team visited my university. During the discussions we asked them for methods for parameter estimation. For this job a set of parameters of an ODE or DAE must be adjusted to match the trajectory to a wanted output or measurement. The Mathworks team has explained, that single shooting methods are not the best choice for these problems, because it is too unlikly that the initial values are suitable to produce a fair start point for a Newton method. Therefore they suggested multiple shooting approachs.
We stopped the discussion soon with a smile. They forgot, that my professor has developped the multiple shooting method together with Deuflhard and Bulirsch.
Numerical maths is not a core skill of MathWorks. Matlab's ODE solvers are bullet-proof, well tested and state of the art for solving toy problems and homework question. Some real scientific problems can be solves also, but this is not the general case.
The treatment of discontinuities might be better in Simulink solvers. I did not examine this.

Sign in to comment.

Jan on 8 Sep 2021
Edited: Jan on 9 Sep 2021
This is not useful. ODE45 has a stepsize controller, which rejectes steps if they do not match the tolerances. This means, that the function to be integrated can be evaluated several times for the same value ot t. In addition each step requires mutliple evaluations of the function to be integrated. Some calls are used for the order 4 and some in addition to process the order 5 for the stepsize control. This means, that a ODE78 integrator would compute a completely different result, because the function evaluations happen for other intermediate times.
In cosequence, the result depends more on the chosen integrator than on the actual mathematical defintion of the problem.
Even using interp1 would create a non-smooth function. In consequence the stepsize control might react unexpectedly: Either ODE45 stops with an error, or the step size is reduced until the rounding errors dominate the jump. Then you get a final value, but this can be massively influenced by the accumulated errors due to a huge number of steps.
Of course you can use a persistent variable to implement such a counter, but this will not be meaningful:
function xdot = x(t,x,A)
persistent counter % WARNING: I recommend to avoid this!
if isempty(counter)
counter = 0;
counter = counter + 1;
a = A(:,counter);
xdot = a*x
There is no way to estimate in advance, how many columns A need.
  1 Comment
Walter Roberson
Walter Roberson on 9 Sep 2021
And remember that it is normal for the variable step ode*() routines to take steps that assessed as failing; are you sure you would want to record the data for those steps?
It is part of the process for the ode*() routines to evaluate at several "nearby" points in order to estimate the curvature to figure out which point will be selected. Are you sure you want to record the data for those nearby points?
Remember that at the time of the call to your ode function, you cannot tell whether that location is going to be rejected, or if the point being evaluated is one of the "nearby" points being used to estimate slope.
Furthermore... the points that are output are not necessarily points that have been passed though the ode function! The ode*() routines estimate the value at the projected location, and do not pass the location through the ode function to get the more precise value.
Basically.... the information you would gather this way would be mostly useless except for something like doing a surface plot.

Sign in to comment.

Walter Roberson
Walter Roberson on 11 Sep 2021
function xdot = x(t,x,A)
a = A(:,counter);
xdot = a*x
My question is how can I create a variable "counter" where it starts at 1 and increases by 1 with each time-step until ode45 is done solving?
Make the current value of A to use (so, a) a shared variable. Also keep a record of the last time boundary that was successfuly passed.
Use an event function. The return value of the event function should be positive if the current input t is greater than the next time boundary after the last successful one. When you detect that the current input t is equal to the next time boundary (to within tolerance), update the shared a and update the record of last successful time boundary to be the current time, and return 0.
event functions are run only after a successful step, and they are run by a section of code that is devoted to finding the zero crossing. The assumption is that if the integration succeeded within tolerance to a location, that with the assumption of continuity, there should be a point between the old and the new where the event function returns 0 and which should also be within integration tolerance.
The assumption of course can fail if the function is discontinuous in one of the first two derivatives, since you might have stepped over a discontinuous area and thought you had succeeded in integration when you really had not succeeded.
So the above is how you can update a according to time.
Is it a good idea? No. But this answers the question you asked.


Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!