Creating a counter inside the ode45 function

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
end

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.

Answer 1

[Star Strider]

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.

.

[Chad]

The issue with using the interp1 function is that I have no control over the "real t" inside the funciton, even if I specify the time-steps, it is only generated in the output. I cannot map the tspan vector with a counter vector if I do not know the "real" elements on tspan (those generated by automatically by ode45). Am I missing something?

[Jan]

What a pity that you can find this bad example of using numerical software in the Matlab's documentation: ODE with Time-Dependent Terms. Matlab's ODE integrators are designed to handle differentiable functions only. The linear interpolation is not smooth, such that the behaviour of the step size controller is undefined.

I've seen too many scientists without experiences in numerical math, who evaluatre their ODEs with ODE45, but cannot answer, why they have chosen this integrator. They do not know, what stiffness means and omit a variation of the parameters and initial values to estimate the sensitivity of the process.

Of course, ODE45 does reply a final value, but considering this as "result" is like using a drilling machine to hit a nail into the wood.

[Star Strider]

@Chad — The MATLAB ODE solvers pass the ‘t’ value as an an argument to the ODE function you are integrating, so that value of ‘t’ is the time step.

I have no idea what you want ‘A’ to do, however it is entirely possible to create it (as a function of ‘t’) as an external function to the ODE function, and then have it evaluate and return whatever values you want to your ODE function. (It would be necessary to pass it as an additional argument to your ODE function in order for it to be defined there and for the ODE function to use its results. I have done this previously, and it works.) This would work with essentially all the MATLAB ODE solvers (at least the ones I am familiar with, such as ode45, ode15s and similar functions, and I assume it would work with all of them).

.

[Paul]

@Jan

In response to this comment ...

Can you clarify what you mean by the "behaviour of the step size controller is undefined?"

More generally, suppose I have a system model that includes components that can only be modeled based on measured data. For example, I have a table of displacement vs force from measurements of a non-linear spring and I need to model that spring in my odefun as F = k(x). So I need to have some way to convert that table into a function that is, at a miminum, continuous in x. Are you suggesting that linear interpolation is never appropriate to use in this situation with ode45? If not linear interpolation, then what are the minimum requirements for k(x)? Also have the first derivative continuous? Second derivative? Higher derivatives? I'm not arguing for one approach or another, just curious as to your opinon on the matter.

I agree that the linked example in your comment is not helpful, and might even be counter productive. To the uninitiated, it could appear that the doc is saying to always use that interp1 approach for time-dependent terms, even for a case where the functional form of those terms is known explicitly (as in that example).

[Jan]

ODE45 computes the steps with 2 methods, one with order 4 and one with order 5. The difference between the 2 solutions is used to determine a good size for the next step. If the difference is small, the size is increased, if it is large, the size iis reduced, and if the size exceeds the limits defined by the tolerances, the current step is repeated with the half step size.

This design is based on the assumption, that the function to be integrated is smooth. For each time step from t to t+dt the function is evaluated a certain number of times to get the results for the two orders. If there is a discontinuity in the trajectory, it depends on several factors, if the controller detects them and reduces the step size. If it does, the number of steps might explode, e.g. you get 1e7 steps instead of one. Because the evaluation of the function to be integrated has a limited accuracy only, this can amplify the numerical noise exponentially in the worst case.

If you are less lucky, the problem is hidden by cancellation errors and the integration does not react to the discontinuity. Then the result of the concerned interval is wrong and this can influence the final value of the trajectory massively.

The rule is to use numerical algorithms only inside their specified design. For the integrator of Matlab this means, that the function to be integrated must be smooth. Otherwise ODE45 is a very expensive random number generator.

My professor for numerical maths told us:

• Inside a function to be integrated, all jumps are forbidden. Therefore these command are not allowed in general: min, max, round, ceil, floor, if, abs, ... For other commands you need to be extremely careful, if they reach a discontinuity: tan, 1/x, ... No linear interpolations, no usage of dates (jumps at 24:00), no integer types.

"convert that table into a function that is, at a miminum, continuous in x"

To by a valid input for ODE45, it must be continuously differentiable. A cubic or spline interpolation solves this, but of course this can have other side effects. E.g. a spline interpolation can reply outputs, which are outside the range of inputs.

Creating a model, which needs to approximate data from a table, is a really hard job. Remember to vary the initial values and the parameters of the tables by a "small" value and check, if the claculated trajectory reacts "strongly". This would be a signal, that the model is instable and not "valid".

The validation of a model can be extremly complicated. Integrators used for reliable scientific work should contain an internal numerical differentiation, which computes the sensitivity in each time step (changes of the result to small variations of the initial values). The sensitivity should affec the step size control also. A clean solution is to detect each discontinuity and stop+restart the integration.

Choosing the tolerances is a serious job also. A lot of inexperienced users set the tolerance such, that they do not have to wait "too long". This can be a wrong choice. A tiny tolerance can cause very small time steps and in consequence the accumulated rounding errors can dominate the trajectory. With a large tolerance the steps are big and the discretization error gets important. Actually you need a numerical analysis of the model to decide, which tolerances lead to the most accurate result.

Numerical maths is a science. Ask a professional to apply numerical software reliably.

[Paul]

@Jan

Thanks for the comprehensive response. I agree with you in spirit, but I'm not as adamant that the function to be integrated, f(x,t), must have continuous first derivatives wrt to both of its arguments. Surely there are many problems where the solver will produce acceptable results (with proper settings for step size, tolerances, etc.) when f(x,t) includes continuous, piecewise linear equations as would be the case when using interp1().

[Jan]

Please define "acceptable". A user can ignore numerical instabilities and be satisfied with an inaccurate result. But this is not the scientific way of working.

For a piecewise linear function, a step of the integration can contain points before and after the change of the slope. The stepsize controller was not designed or tested to handle this correctly. If ignoring the fact, that the tool is not designed to solve the problem accurately, means "accepting" the final value as a "result", that you are right.

My professor for numerical maths was asked to find the cause of problems concerning the validation of a modelling of a car in a crash test simulation. The software was used by a famous German manufacturer. He found discontiuities in the model and added some code for the internal numerical differentiation to control the sensitivity, such that similar problems will be detected at least.

Do you want the design of the car, you are driving in, be optimized by software, which is driven outside the specifications? Well, I admit, a rhetorical question. Sorry.

As mentioned already, I've seen to many scientific publications based on integrations of ODEs without understanding, how such a computation has to be solved reliably. If I simulate the motion of a pen standing on its point, I could accpet the result that it will stay in this position forever - even without a discontinuity. But this is not science. Of course ignoring the discontinuities of ODEs can reply similar trajectories as if ODE45 is restarted at each jump. But exceptions do not proove, that the method is applicable in general.

[Walter Roberson]

"I'm not as adamant that the function to be integrated, f(x,t), must have continuous first derivatives wrt to both of its arguments."

In the following, the failure for the first case is not explicitly pointed out if the 1e16 are replaced with 1e15. The failure for the third case is not explicitly pointed out if delta is changed to 1e-14.

The Dt and D variables are used in the third equation to create the breakpoints for interp1 purposes: the slope is 1 from 1 to just less than 5, then interp1 is allowed to interpolate over a short distance to reach 1e16, then it is held constant at that until just before 10, then permitted to drop over a short span just before 10. This is approximating piecewise constant 1 until time 5 then 1e16 until time 10 then 1 again.

Certainly times [0 5 10], value [1 1e16 1] could have been used, but I think you will agree that would represent a very different function.

These results suggest that a guideline might be that ode45 will specifically complain for areas with slope about +/-1e31 (such as change of 1e16 over time 1e-15).

This does not mean that it gets correct results for smaller slopes: only that until then it adapts to produce something

tspan1 = [0 5]; tspan2 = [5 10]; tspan_full = [0 10];
y0 = 0;
delta = 1e-15;
Dt = [0, 5-delta,    5, 10-delta, 10];
D  = [1,       1, 1e16,     1e16,  1];

odefun1 = @(t,y) 1;
odefun2 = @(t,y) 1e16;
odefun = @(t,y) (t<=5).*odefun1(t,y) + (t>5).*odefun2(t,y);
odefun3 = @(t,y) interp1(Dt, D, t);

[t1, y1] = ode45(odefun, tspan_full, y0);

Warning: Failure at t=5.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-14) at time t.

plot(t1, y1)

[t2a, y2a] = ode45(odefun1, tspan1, y0);
y02 = y2a(end,:);
[t2b, y2b] = ode45(odefun2, tspan2, y02);
t2 = [t2a;t2b]; y2 = [y2a; y2b];
plot(t2, y2)

[t3, y3] = ode45(odefun3, tspan_full, y0);

Warning: Failure at t=5.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-14) at time t.

plot(t3, y3)

[Walter Roberson]

format long g
tspan1 = [0 5]; tspan2 = [5 10]; tspan_full = [0 10];
y0 = 0;
delta = 1e-13;
S1 = 1;
S2 = 1e15;
Dt = [0, 5-delta,    5, 10-delta, 10];
D  = [S1,     S1,   S2,       S2, S1];

odefun1 = @(t,y) S1;
odefun2 = @(t,y) S2;
odefun = @(t,y) (t<=5).*odefun1(t,y) + (t>5).*odefun2(t,y);
odefun3 = @(t,y) interp1(Dt, D, t);

[t1, y1] = ode45(odefun, tspan_full, y0);
plot(t1, y1)

approx_area1 = trapz(t1, y1)

approx_area1 = 
                  1.25e+16

[t2a, y2a] = ode45(odefun1, tspan1, y0);
y02 = y2a(end,:);
[t2b, y2b] = ode45(odefun2, tspan2, y02);
t2 = [t2a;t2b]; y2 = [y2a; y2b];
plot(t2, y2)

approx_area2 = trapz(t2, y2)

approx_area2 = 
                  1.25e+16

[t3, y3] = ode45(odefun3, tspan_full, y0);
plot(t3, y3)

approx_area3 = trapz(t3, y3)

approx_area3 = 
      1.24996476089707e+16

syms t
y(t) = piecewise(t <= 5, S1*t, t > 5 & t < 10, S2*t-S2*5+5)

y(t) = 

int(y(t), t, 0, 10)

ans = 

double(ans)

ans = 
                  1.25e+16

syms t
y(t) = piecewise(t <= 5, S1*t, t > 5 & t < 10, S2*t)

y(t) = 

int(y(t), t, 0, 10)

ans = 

double(ans)

ans = 
                  3.75e+16

The last two symbolic cases show the huge difference in output depending on how you handle the boundary condition The first of the symbolic cases shows that if you say that the boundary must be treated as continuous, then the result is only about 1/3 of the result for the case where the boundary is treated as discontinuous

[Paul]

@Jan

Of course there is not a single definition of "acceptable." Acceptance of the result, or validation of the model, is problem dependent and depends on several factors, one of which is the intended use. Much more rigor needs to be applied in using models to design and implement safety critical systems, for example, as opposed to develping a model to simulate a mass attached to a spring for an introdcutory physics class.

"Of course ignoring the discontinuities of ODEs can reply similar trajectories as if ODE45 is restarted at each jump. But exceptions do not proove, that the method is applicable in general."

I never said, nor did I mean to imply, that the method [of ignoring discontinuities] is applicable in general.

@Jan and @Walter Roberson

Perhaps I wasn't clear in the point I was trying to make, so I'm going to amend my statement with an additional sentence:

I agree with you in spirit, but I'm not as adamant that the function to be integrated, f(x,t), must have continuous first derivatives wrt to both of its arguments. Surely there are many problems where the solver will produce acceptable results (with proper settings for step size, tolerances, etc.) when f(x,t) includes continuous, piecewise linear equations as would be the case when using interp1(). Similarly, there are many problems where integrating piecewise linear equations will produce unacceptable results without precise step control (or other mitigations) in the solver that accounts for the discontinuties.

Walter's example of integrating through an extremely large change in the table data as in the comments above is an excelent example of the latter type of problem.

Are there no examples of the former?

[Jan]

@Paul: Thanks for your explanations.

The problem is not trivial. There are no "smooth" function at all in numerical maths: all values are quantized due to the limited precision of floating point values in the IEEE754 representation. If you have the output of a function, you cannot prove numerically, that it is differentiable or not. Each kink might be perfectly smooth in a microscopic scale which cannot be represented with 15 digits.

So what does "smooth" mean in numerical math, e.g. in a numerical integration? The algorithm of the integrator is devlopped based on a mathematical theory and real numbers. All underlying maths requires, that the functions are differentiable. As soon as the algorithm is implemented, it does not work on real numbers, but in quantized IEEE754 doubles. Two generations of computer scientist struggled with the problem, if this is valid. They tested thier theories exhaustively and found out, that the results of e.g. Runge-Kutta-integrators are good approximations of closed solutions of the ODEs (if there are any).

I do not know any publication concerning Runge-Kutta integrators working on non-smooth function without further methods, to detect the discontinuities and to properly stop and restart the integration and the step size controller. The integration might or might not reply an acceptable output.

Maybe there is such a publication already and I just did not read it yet. It could explain, that the integration works fine, if the n.th derivative does not have jumps larger than 1e-Q times the minimal step size, when all components of the trajectory are normalized to mean=1.0 and std=1.0, and the 2-norm is used to control the error with a specific weighting method.

Did you see the code to determine the step size for the next step in ODE45:

  % If there were no failures compute a new h.
  if nofailed
    % Note that absh may shrink by 0.8, and that err may be 0.
    temp = 1.25*(err/rtol)^pow;
    if temp > 0.2
      absh = absh / temp;
    else
      absh = 5.0*absh;
    end
  end  

This is voodoo already. Why 0.8? Where do these limits come from? Some integrations deliver more accurate results, if ODE45 does not divide the step size by 2, but by exp(1) in case of a rejected step. Why? Nobody has examined it yet, because for some other functions sqrt(2) or pi is better.

Step size controller require some heuristic parameters. The same applies for finding a "good" step size for the quotient of differences to approximate the derivative. An educated guess requires the size of the 2nd derivatives and the problem gets worse.

Matlab offers a set of different integrators. All of them are very basic and better for educational purpose than for the application on real scientific problems. Nevertheless, many scientists with deeper knowledge in numerical maths use them to produce outcomes for publications. There will be many examples, where their results are meaningful. But some results are just junk.

Looking for examples and counter examples is not useful here, because this does not allow to draw conclusions:

  0.1 + 0.1             == 0.2  % TRUE
  0.1 + 0.1 + 0.1 + 0.1 == 0.4  % TRUE
  0.1 + 0.1 + 0.1       == 0.3  % FALSE?!?

The rule for comparing floating point numbers is, that a "useful" limit is required and that the size of this limits depends on the problem you want to solve. This rule is ugly and you will find so many counter examples. But this is not trustworthy.

The ODE integrators shipped with Matlab are designed to handle smooth functions only. The development of the algorithm was based on this fact. You will find many non smooth functions, which are treated with an "acceptable" accuracy also. But this does not allow to estimate, if this holds true for a new problem also.

[Paul]

Is there TMW staff that participates in this forum with expertise in TMW ODE solvers? If there is, it would be interesting to call them with an @ and hear their response to this statement:

"Matlab offers a set of different integrators. All of them are very basic and better for educational purpose than for the application on real scientific problems." I assume that "integrator" in this context means "ODE solver."

I take no position on this statement myself.

FWIW, on the Simulink side, the solvers and models can be configured to deal with many types of discontinuities as we are discussing here. But it also offers a library of blocks for table lookups, and one option in those blocks is linear interpolation where the solver won't attempt to detect the breakpoints (at least as I understand how it works).

And I should add that the Matlab ODE solvers also offer event detection as well, so you can find and solve up to a discontinuity and then restart the solver. I don't normally use the ODE solvers in Matlab so don't have an opinion on how well this works.

[Walter Roberson]

Is there TMW staff that participates in this forum with expertise in TMW ODE solvers?

No, or at least not that I can recall.

I should add that the Matlab ODE solvers also offer event detection as well

Yes, but most people do not use them, and there is no facility to automatically restart. The user has to know to program a loop of calls.

[Paul]

It is too bad that the user has to program a loop, especially when the Simulink solvers make this transparent to the user.

Within the loop, the funcion odextend() looks like it would be helpful.

[Jan]

@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.

Answer 2

[Jan]

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;
end
counter = counter + 1;
a = A(:,counter);
xdot = a*x
end

There is no way to estimate in advance, how many columns A need.

[Walter Roberson]

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.

Answer 3

[Walter Roberson]

function xdot = x(t,x,A)
a = A(:,counter);
xdot = a*x
end

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.