Is there any way to use multiple tolerances in ODE solvers?

I'm frequently using systems of ODEs to simulate the dynamics of ecological communities. The number of equations is large: from hundreds to thousands, and running times are sometimes desperately long. The diversity of species’ characteristics implies that some of them have dynamics which are orders of magnitude faster than others: like the difference between a bacteria and a tree. The ‘faster’ species require tighter tolerances which, in turn, result in smaller time steps, but that is usually unnecessarily small for the ‘slower’ species, which are essentially constants (or very slowly varying) at that time scale (e.g. hours vs. months or years) and could use much larger tolerances and longer time-steps. I guess this, if possible, would save computing time, which is the goal here. So, my question is: is there any way to set at least two different tolerances for different groups of equations in a system of ODEs? Or any other way to get around this problem? I use MATLAB R2017b.

Answer 1

[Jan]

No, there is no way to avoid the application of the small time steps to all components of the trajectory using the ODE integrators of Matlab. The only exception is, that the integrations can be separated, if the different components do not interact with each other.

This is a common problem at simulation large system, e.g. the global weather, a organism with resolution of single cells, the mass distribution of the universe, or a vehicle simulation considering high frequency vibrations of small elements. If you want to simulate the whole system, the part with the highest frequency determines the step size.

It is possible to create integrators, which identify high frequency parts dynamically. Then these parts are integrated with a smaller steps until the changes reach a limit, in which the other components are affected. The implementation is hard and needs some heuristics, because the identification of the low and high frequency parts require to simulate the system at first...

Using a vector of tolerances is equivalent to scaling the components, e.g. if a certain component is measured in meter or nanometer. This influences, how the different components are used to determine the step size, but again the steps are applied to all components.

[Francisco de Castro]

I'm confused now. Torsten's suggestion seemed promising and, in fact, you can assign a vector of values to AbsTol in ODESET: opt= odeset('AbsTol',[1e-4:1e-2:1e-1]) However, Jan is always right, so does this mean that the solver will ignore the extra values in AbsTol? Or what would it do with them? Maybe use only the smallest? F.

[Jan]

 However, Jan is always right,...

;-) Let's try this: "I'm left."

The vector provided as absolute or relative tolerance is considered for the step size control: E.g. ODE45 calculates each step with two methods (with the orders 4 and 5). The step size is chosen such, that the difference between the results is a little bit smaller than the given tolerance:

x4 = <result from 4th order>
x5 = <result from 5th order>
absDiff = abs(x4 - x5)
relDiff = absDiff ./ max(x4, x5)
if any(absDiff > absTol) || any(relDiff > relTol)
  ... reject the step, reduce the step size

Now these comparisons accept either and absTol and relTol as scalar, or as a vector of the same size as the trajectory. This means, that the step size control reacts with a different sensitivity to the different components, but you get still the same steps size for all components.

Answer 2

[CARLOS RIASCOS]

Hello friend, you could try, generate for example two systems of equations one with slow dynamics and another with fast dynamics, for the slow dynamics use a value of big h and for fast dynamics a value of small h, also consider in changing the limits of tspan, because if you are analyzing a dynamic that evolves in micro-seconds should not have limits like the code here example but rather tspan = 0: h_small: 0.01.

If you must make a model as a whole, you can see if the fast dynamics stabilize and if so you could even consider them as constants, but everything depends on the type of problem you have. I hope I've helped.

    h=1;
    tspan = 0:h:5;
    y0 = 0;
    [t,y] = ode45(@(t,y) 2*t, tspan, y0);
    plot(t,y)

[Francisco de Castro]

Thanks, but your suggestion would not work. All the variables must be in the same system since many of them interact with each other. I don't understand your comment about tspam, and in your example code there is nothing remotely related to the question.

Answer 3

[Francisco de Castro]

Thanks Torsten, your suggestion may be the solution. If you put it as an separate answer I'll mark it as accepted. F

Answer 4

[Torsten]

AbsTol can be a vector of the same size as the solution vector. Thus you can prescribe absolute tolerances for each solution component separately.

Best wishes

Torsten.