When solving Ordinary Differential Equations in time using, say R, is common to find ourselves in a part of the parameter space that does not allow for correct solutions. For simple models it is easy to derive parameter conditions in a way that we only end up with appropriate solutions.
Consider a complex ecological system with many species, that is, multiple ODEs. It can be difficult or impossible to arrive analytically at the intervals for parameters that lead to proper states (i.e. all species abundances should be positive for example). We can try to choose random parameter combinations until we get results that make sense, but often times we get hyper-exponential blow-ups in growth.
Is there any numerical or analytical trick that can help in avoiding these solutions?
Hello Alfonso,
do you mean that your ODEs system will diverge for some choices of the parameters, while it should lead to a fixed point (your final state) when the parameters belong to certain areas? Or is it just a numerical divergence? In the latter case use implicit methods with higher order to solve the equations.
If you're looking for the area of stability, you can use a systematic grid algorithm by stepping each parameter in a given range. However it takes long! You stop any simulation when its solution becomes greater than a predefined bound.
Gabriel
For me, too, it is a bit unclear if your question is about divergence in numerical time integration, or about determination of stability boundaries.
Adding to what Gabriel Thomas wrote, if the question is about stability boundaries, then as an alternative to a brute force grid algorithm you could also consider numerical continuation and bifurcation analysis. There exist a number of freely available software implementations. One of these is MatCont, co-developed by Yuri Kuznetsov from Utrecht. (If you need pointers, feel free to let me know.)
Afonso Dimas Martins,
Can you offer an example of such a system of ODEs, perhaps in the literature?
As Gabriel Thomas and Sebastiaan Janssens say, there may be no efficient general methods for determining the area of stability, therefore it may be necessary to find a method that uses some properties of your equations. And I agree that one should first be certain that the divergence is inherent to the system of ODEs and not a problem of the solver.
I could add that if you know some parameters for which the system converges, start from them to explore the parameters space.
Take care.
see for example
Article On Blow-Up Solutions of Differential Equations with Poincaré...
I suggest you see the article
On the Computation of Blow-up by A.M. Stuart and M.S. Floater.
Assuming that the blow up is due to the actual nature of the system and not to the numerical integration (if numerical approximation is performed by standard ODE packages, e.g. MATLAB solvers, time step adaptation and/or use of stiff solvers generally take care of most numerical instabilities) one basic possibility is always to linearize your system close to the typical initial state and study its eigenvalues. If you have eigenvalues with large positive real parts for certain parameter ranges, it is clear that blow up will ensue.
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics
- Statistical Software
- R