instability ODE45/ODE23 Reltol

I have found an article with a mathematical model for adsorption on spherical particles in a batch reactor. It is quite straightforward, so I wanted to implement it in matlab. The authors of the article used Runge-Kutta 4th order to solve the system of 2 coupled ODE's. That's why I tried ode45 and ode23.

All parameters and data were published in the article, so it should be easy to solve the system. This, however, turned out to be not that simple, at first I wanted the RelTol to be quite low, say 1e-4. Matlab gave me the warning that the integration tolerances could not be met, so I tried less tight tolerances and ended up with 0.1 for ode45 and 0.282 for ode23, which are unacceptable high tolerances.

For most other tolerances the solutions became unstable, my questions to you are:

• What could be the reason for the instability? I checked the ODE's for a million times now, they can't be wrong anymore...
• Is there a standard way for determining the RelTol in such a way that it gives stable solutions?
• I would like to give you insight in the script I wrote, can I upload it somewhere

Thanks in advance for your answer!

[Arnaud Miege]

What have you set AbsTol to? Typically an order of magnitude smaller than RelTol is a good rule of thumb. You might also want to try *tightening* the tolerances rather than relaxing them. It's counter-intuitive but it can help.

[Arnaud Miege]

Also, if the system is stiff, you may want to try ode15s or ode23t.

Answer 1

[Jan]

AbsTol means, that local discretization error is controlled by an absolute limit:

abs(y1 - y2) < AbsTol

with y1 and y2 are two solution obtained by methods, e.g. orders. If y1 is smaller near to 0, the absolute error looses its meaning. Imagine the trajectory y lives on the interval 10E-64 to 10e-70. In such cases the relative tolerance is more meaningful:

norm(y1 - y2) / min(abs(y1), abs(y2)) < RelTol

Anyhow, even this fails for y1 = 0 and y2 = 1e-70 due to a division by zero... Therefore an addition threshold is needed to dertmine the local errors. For ODE45 this threshold is AbsTol/RelTol - a reasonable, but arbitrary choice. The chpoise of the NORM can be importamt also, e.g. if the components of y have different magnitudes: 1. if y(t) is [1e100, 1e-100], 2. if y(t1) is [1, 1] and y(t2) is [1e100, 1].

Conclusion: The re-caculation of an integration explained in a paper is not trivial! Ask the author, which integrator has been used.

[Adriaan]

Thanks for the quick answers, I already tried to tighten the tolerances, because I want an accurate solution in the first place. This however caused my computer to display the annoying word 'busy' for several hours...

I tried the stiff solvers as well, but they happen to be also very sensitive for the set tolerance/parameters combination. As I understand from the story above of Jan Simon, adding an AbsTol is only relevant when you already have an idea which local discretization error the solver could handle and if the absolute values are not hard to handle for matlab (like values near zero).

At last, I have sent an e-mail to the author of the article, but as it is a publication of 2003 he is not very eager to mail back. I appreciate your help, after the weekend I will dig deeper into the Norm thing, maybe that helps...