How to find simple nonlocal stability and resilience measures

Stability of dynamical systems is a central topic with applications in widespread areas such as economy, biology, physics and mechanical engineering. The dynamics of nonlinear systems may completely change due to perturbations forcing the solution to jump from a safe state into another, possibly dangerous, attractor. Such phenomena cannot be traced by the widespread local stability and resilience measures, based on linearizations, accounting only for arbitrary small perturbations. Using numerical estimates of the size and shape of the basin of attraction, as well as the systems returntime to the attractor after given a perturbation, we construct simple nonlocal stability and resilience measures that record a systems ability to tackle both large and small perturbations. We demonstrate our approach on the Solow–Swan model of economic growth, an electro-mechanical system, a stage-structured population model as well as on a high-dimensional system, and conclude that the suggested measures detect dynamic behavior, crucial for a systems stability and resilience, which can be completely missed by local measures. The presented measures are also easy to implement on a standard laptop computer. We believe that our approach will constitute an important step toward filling a current gap in the literature by putting forward and explaining simple ideas and methods, and by delivering explicit constructions of several promising nonlocal stability and resilience measures.


Introduction
Understanding stability of dynamical systems is important for widespread areas of research such as control theory, economy, biology, physics, hydrodynamics and mechanical engineering. As computers and computational tools have become more and more efficient, simulations and numerical investigations of realistic, high-dimensional, advanced mathematical models have become easier and also increasingly popular. As such advanced models account for numerous factors and their interplay, the dynamics are often nontrivial, nonlinear and coexisting attractors (bistability) often exist or may be hard to rule out. Such systems algebraic structure are usually also complicated and classical mathematical analysis is therefore often difficult. This calls for further research in mathematical analysis, as well as in numerical methods, for finding efficient ways of quantifying the stability of dynamical systems.
The analysis of the stability of attractors (e.g., equilibrium points, limit cycles, quasiperiodic or strange attractors) in dynamical systems naturally split into local analysis and nonlocal analysis. The local stability approach is usually based on linearizations and yields information in a small neighborhood of the attractor, saying little or nothing about the systems behavior a bit away from the attractor. Measuring stability in nonlinear dynamical systems using local methods (such as eigenvalues of the Jacobian matrix at an equilibrium) delivers information of how the system reacts only on arbitrary small perturbations. A perturbation of a given size may push a locally stable state into another attractor having a completely different, possibly dangerous, behavior. For example, researchers believe that the crash of the aircraft YF-22 Boeing in April 1992 was caused by a sudden switch to an unsafe attractor [1,2]. Figure 1 shows three systems having completely different abilities to withstand perturbations, but exactly the same local stability and local resilience.
The local stability analysis should therefore be complemented by, or replaced by, a nonlocal approach that considers properties of the basin of attraction (henceforth basin) to the attractor under investigation. The size of basin (its volume) constitutes a natural candidate for a nonlocal stability measure as a large basin indicates that the system comes back to the attractor with a high probability, given a random perturbation. It is also easy to numerically estimate the size of the basin. However, if the distance from the attractor to the boundary of the basin is short in some direction, then a small perturbation in this direction can push the system to another, possibly dangerous, attractor, even though the basin is large, see Fig. 1c. Therefore, the shape of basin is another natural candidate to be included in a (a) (b) (c) Fig. 1 Three identical balls rest at equilibrium on three different stands. All stands have the same shape in a neighborhood of the equilibrium and hence a local approach would rank the systems in a-c as equally stable and equally resilient. The basin of attraction is largest in (a), followed by (c) while it is smallest in (b). Therefore, measuring only the size of the basin of attraction would rank b as least stable even though it is easier to push the ball off the stand in c as the equilibrium in c is closer to the boundary of the basin. Hence, both the size and the shape of the basin of attraction should be included when measuring stability and resilience of a nonlinear dynamical system nonlocal measure construction. It is not trivial how to construct measures based on the shape of the basin, but we may conclude that a large and convex basin, having the attractor in the middle, should indicate high stability of the system. This is because in such case a large perturbation can be given to the system, in any direction, without pushing the system to another attractor.
In the field of mathematics, it is classical to attempt to find Lyapunov functions in order to prove nonlocal stability of attractors in dynamical systems. Such method gives analytical estimates of both the size and the shape of the basin. However, it is often a difficult task to find suitable Lyapunov functions, making this approach insufficient, especially when dealing with high-dimensional complicated systems of equations. If a Lyapunov function is not available, another investigation of the basin has to be considered in order to understand the nonlocal stability of the attractor. The aim of this paper is to deliver simple methods for doing this by considering several nonlocal stability and resilience measures built upon numerical estimates of the basin and on the returntime of trajectories corresponding to perturbations from the attractor.
Keeping in mind that works involving nonlocal stability approaches can be hard to find as papers generally put focus on the results of a stability analysis rather than the stability measure itself, a literature survey indicates that nonlocal stability measures are rarely used. Most works connecting to applied dynamical systems, e.g., in theoretical ecology and mechanical engineering, seem to rely solely on a local investigation of stability. This is surprising because even though a mathematical approach is difficult, today's computers can often calculate estimates of the basin within a couple of minutes. Lundström [3] and Lundström and Aidanpää [4] used two nonlocal stability measures when investigating stability of an electric generator. The authors refer to the systems ability to absorb perturbations as robustness of the system. Their first stability measure estimates the basin's size, reflecting the probability that the generator system returns to the "safe" attractor given a perturbation. A similar approach was used by Menck et al. [5] to prove properties of neutral networks and power grids; they refer to the nonlocal stability measure as basin stability and give also several examples of systems in which local stability measures fail to detect crucial destabilizing effects. Lundström and Aidanpää [4] also used a stability measure based upon the shape of the basin through approximation of the smallest pertur-bation needed to push the system from the safe attractor into another "unsafe" attractor. A similar idea was presented by Klinshov et al. [6] as the stability threshold, and algorithms for calculating it in the setting of general attractors can be found in [6] and Kerswell et al. [7].
A first aim of this paper is to expand on the above ideas by constructing and evaluating simple nonlocal measures of stability using fundamental properties of the basin such as its size and shape. To do so, we consider two nonlocal measures accounting for both large and small perturbations. The first measure, denoted by P, estimates the size of the basin and thereby answers the natural question; what is the probability that the system returns to the attractor given a perturbation? The second measure, denoted by D, builds on the shape of the basin and concerns the question; what is the smallest perturbation needed to push the solution to another attractor? The measure D calls for suitable ways to compare distances in the state-space, and, focusing on mechanical systems, we construct the version D energy reflecting the least amount of energy needed to push the systems solution into another attractor. We also suggest the version D rel , considering relative distances, and show how to use a version of it when evaluating harvesting strategies in stage-structured populations.
A second aim of this paper is to include another natural candidate for nonlocal stability into the measures construction, i.e., the returntime of a trajectory to the attractor given a perturbation. Doing so lead us into the resilience of a system and construction of resilience measures, an important concept increasingly used in ecology (see, e.g., [8][9][10][11][12][13]). There are at least two different definitions of resilience of a system in the literature. The first definition, and the more traditional, concentrates on stability near an equilibrium, where resistance to disturbance and speed of return to the equilibrium, following a perturbation, are used to measure the property [14][15][16]. This view of resilience provides one of the foundations for economic theory as well and may be termed engineering resilience [17]. The second definition emphasizes conditions far from any equilibrium, where instabilities can flip a system into another regime of behavior, i.e., to another stability domain [18]. In this case, the measurement of resilience is the magnitude of disturbance that can be absorbed before the system changes its structure by changing the variables and processes that control behavior. This second view may be termed ecological resilience [19] and clearly involves properties of the basin of attraction.
Even though neither of the above definitions of resilience restrict to small perturbations, widespread technics for calculating resilience are based solely on local analysis by considering, e.g., the eigenvalues of the Jacobian matrix at equilibrium. See e.g. Neubert and Caswell [8], Arnoldi et al. [11] and the references therein for more on standard local resilience measures, alternative resilience measures and their relations as well as discussions concerning shortcomings of the local approach in general. For further discussion on the use of local resilience in ecology and the fact that it can be difficult to assess from an empirical point of view, see Haegeman et al. [12].
Mitra et al. [20] suggest an integrative measure of resilience based on both local and nonlocal features, and Lundström et al. [21] use a simple nonlocal resilience measure for evaluating harvesting strategies of age-and stage-structured populations. However, our impression is that, in analogue with the literature on pure stability measures discussed above, the nonlocal approach is under-used also in the setting of resilience. As a consequence, valuable information of a systems ability to sustain perturbations may be unrecorded: Recall that all three systems in Fig. 1 have identical resilience according to any measure based on a local linearization (such as the largest eigenvalue) near the equilibrium.
This motivates our second aim: to find simple nonlocal resilience measures based on the attractors basin as well as on the returntime for trajectories starting within the basin (and thus returns to the attractor under investigation). We will present four nonlocal resilience measures, accounting for large as well as small perturbations, of which the first two steams from seeing resilience as the reciprocal of the returntime, i.e., as a rate of return given a perturbation. Indeed, our simple measure R will answer the question; what is the expected rate of return of the system given a random perturbation? We also keep track of a "worst-case resilience," or slowest recovery from a set of perturbations, through the measure R worst answering; what is the slowest rate of return of the system given a random perturbation?
We proceed by introducing the concept basin-time as the subset of the basin of attraction from which all trajectories return to the attractor within a certain time limit. Based upon the basin-time, we build two more resilience measures. The first is P τ which concerns the question what is the probability that the system returns from a random perturbation in time τ ? The second is D τ , what is the least perturbation from which the system will not return in time τ ? The above questions should be relevant for wide ranges of applications, e.g., in biology, fisheries management, economy, and management of financial systems. Our resilience measures are applicable in the setting of simple linear systems as well as in high-dimensional nonlinear systems. All suggested stability and resilience measures are easy to implement on a standard laptop computer, and together, they yield an explaining picture of the systems dynamics in a "large" neighborhood of the attractor.
In Sect. 2, we present the main ideas and definitions of the nonlocal stability and resilience measures, while in Sect. 3 we demonstrate our approach on four independent mathematical models; the Solow-Swan model of economic growth, a simple electromechanical model, a stage-structured population model and a high-dimensional model related to evolution of a multidimensional phenotype. In the applications, we discuss what can be learned by each measure, and why. We also consider relations between the measures as well as their relation to linear measures. In Sect. 4, we discuss the achieved results, how to apply the simple measures efficiently to more general situations, some topics for future research, and finally we give some concluding words.

Methods
We focus on models governed by a dynamical system on the standard form Here, f (x, t) is a given function defined on R d ×R taking values in R d , d is the dimension, andx ∈ R d gives the initial condition. However, our ideas are general and the following measure constructions have potential for generalization to systems described by a wider range of mathematical tools as well, such as systems of partial differential equations, systems combining partial and ordinary differential equations, stochastic differential equations and cellular automaton.
To set up the measures, the fist step is to choose a suitable set of perturbations to test the system for. This choice may preferably reflect what the system may be exposed to in reality. Perturbations can be taken deterministically or randomly from a predefined probability distribution; the measures presented below can handle any reasonable choice of perturbations. A perturbation is inferred through an initial condition of (1), and thus the set of perturbations is a set of initial conditions.
For each initial condition in the set of chosen perturbations, the corresponding trajectory is numerically integrated in order to see if it returns to the attractor (of which stability and/or resilience are to be investigated) or not. All initial conditions that returned are saved as "safe" in a set I safe and those that did not return as "unsafe" in a set I unsafe . The time needed for the trajectory to recover the attractor is also saved as the returntime T (x i ) for each safe initial condition x i ∈ I safe . Let N safe , N unsafe and N tot denote the number of the safe, the unsafe and the total tested initial conditions, respectively. By "return to the attractor," we mean that the trajectory has entered a small predefined neighborhood of the attractor. The choice of this (very) small neighborhood should be based on what is considered as normal, or perfect, operation of the system.

Nonlocal stability measures
We choose to measure nonlocal stability through the size of basin using the simple measure P estimated bŷ The measure P reflects the probability that the system returns to the attractor given a perturbation from the pre-specified set of initial conditions. We have 0 ≤ P ≤ 1; if P = 0, then no trajectories return and all perturbations are unsafe, and in case all trajectories return then P = 1. The similar measure was considered by Lundström and Aidanpää [4] and by Menck et al. [5]. We measure stability through the shape of basin of an attractor by considering the shortest distance from the attractor to the boundary of the basin. We denote this measure by D and estimate it aŝ where A is the attractor of which stability is to be investigated and dist(x, A) denotes the distance from x to A defined in a suitable way. The measure D reveals the smallest perturbation needed to push the system to another attractor. A similar idea was considered by Lundström and Aidanpää [4], Klinshov et al. [6] and Kerswell et al. [7]. In Sect. 3, we will demonstrate our approach to estimate D in the setting of equilibria, and in Sect. 4 we further discuss how one can estimate our suggested measures for some general attractors in a simple way.
It is important and often nontrivial to find a suitable distance function, or norm, dist(x, A), for the specific application of the measure D. Standard Euclidean distance can be used due to its simplicity, but is seldom motivated as different dimensions in the statespace may have completely different meanings in terms of perturbations, e.g., in mechanical systems one often ends up with comparing velocity-dimensions to displacement-dimensions. To deal with such case, we consider an energy-norm to construct the version D energy , which finds the least amount of energy needed to push the system from its current state into another attractor. Indeed, a perturbation can be given through a displacement, a velocity impulse, or a combination of these. The energy-norm is defined as the potential energy needed to perform the displacement perturbation plus the kinetic energy given through the velocity impulse. We further explain and demonstrate this idea in Sect. 3.2. When dealing with biomass in population dynamics, we found it relevant to use a relative distance function and introduce the version D rel , considering the size of the perturbation in relation to the biomass at the original state. We explain and demonstrate this idea in Sect. 3.3. As the measure D searches for a particular perturbation, "the most dangerous one," one may consider a distribution of perturbations that puts more weight near the boundary of the basin (if such information is available) in order to estimate it efficiently.

Nonlocal resilience measures
By building resilience on the recovery rate (the reciprocal of the returntime), we consider the simple resilience measure R through the estimatê The measure R yields the expected rate of return of the system given a random perturbation. It accounts for the basin of attraction by giving zero resilience to initial conditions in I unsafe of which trajectories not return to the attractor. The parameter t > 0 is added to avoid division by zero and explosion for very short returntimes and should be given a small value relative to the systems timescale (we will use t = 1 if nothing else is mentioned). We consider also the worst-case resilience, R worst , as the reciprocal of the slowest recovery (longest returntime) through the estimatê The measure R worst reflects the slowest rate of return of the system given a random perturbation. It is often relevant to restrict perturbations further in order to specify this measure. This can be done by replacing the set of all safe perturbations, I safe , by a refined set L ⊂ I safe of perturbations. We give discussions and examples of this in Sects. 3.1 and 3.3.
To proceed, we introduce the concept basin-time as the subset of the basin of attraction from which all trajectories return within a time limit, t ≤ τ . Based upon the basin-time, we build two more resilience measures, P τ and D τ , by simply replacing the basin of attraction by the basin-time in the constructions of the nonlocal stability measures P and D. To calculate these measures, one has to, for a given τ , save all initial conditions that returned in time τ in a set I τ safe and those that not returned in time τ in a set I τ unsafe . Let N τ safe , N τ unsafe and N τ tot denote the number of initial conditions in each set, respectively. We estimate the resilience measures P τ and D τ througĥ In analogue with the stability measures P and D, we note that the measure P τ is based upon the size of the basin-time, while D τ is based upon the shape of the basin-time, and that properties and extensions of the stability measures trivially generalize to these resilience-type measures. In contrary to the stability measures, the resilience measures P τ and D τ will be applicable also in the setting of simple linear systems having only one global attractor (in which always P = 1 and D = constant) as they record also if the dynamics becomes faster/slower by accounting for the returntime. The measure P τ reflects the probability that the system returns from a random perturbation in time τ , and D τ yields the smallest perturbation from which the system will not return in time τ . The presented approach involves some degrees of freedom in terms of measure-parameters, e.g., the choice of a suitable set of perturbations, the choice of a distance function for D and the choice of time limit τ for the measures building on the basin-time. Even though these choices may alter the results, they should not cause much of a problem since one is always interested in how a measure reacts on a change of a system, and thereby a relative change of the value of the measure for fixed values of the measureparameters.
When working with randomly inferred perturbations, an efficient sampling method may be useful for selecting the initial conditions from the prescribed probability distribution. One may thereby obtain sets of initial conditions that are well spread over the desired part of the state-space, which improves the "covering" of the state-space and reduces the variance of the estimates of the measures. In such way, the number of initial conditions that need to be tested in order to obtain good enough estimates may be reduced. In Sect. 3.2, we explain a simple way to achieve this by using the local pivotal method by Grafström et al. [22] when sampling the perturbations.
In the following sections, we will frequently compare our results to the standard local resilience measure given by the real part of the eigenvalue of the Jacobian matrix at equilibrium with largest real part. We denote this measure by −λ max .

Results
In this section, we compare and discuss properties of the suggested nonlocal stability and resilience measures on four different models: the Solow-Swan model of economic growth, a simple electro-mechanical model, a stage-structured population model and a highdimensional model related to evolution. Through the examples, we also discuss relations to a local approach for measuring stability and resilience. Each model is considered in an independent section, and the reader may start reading the application closest to her/his interest.

The Solow-Swan model of economic growth
In this section, we consider the simple one-dimensional differential equation is the production function and s is a constant fraction of the income saved per worker, and C x is the investment required to maintain capital per worker, then the differential equation (2) is the well known Solow-Swan model of economic growth [23,24]. With this model, Solow and Swan showed that, under suitable assumptions on F(x), the capital per worker function x(t) reaches an equilibrium after sufficient time. Indeed, the function F(x) is often assumed to grow faster than maintenance C x as x is small, and slower than C x as x is large, see Fig. 2. Under such assumption, the Solow-Swan model (2) has a globally stable equilibrium (E) when saving equals maintenance, and an unstable equilibrium at the origin: to the left of E, F(x) > C x and therefore x is growing until it reaches E, and to the right of E, F(x) < C x and therefore x is decreasing until it reaches E.
To test the measures of stability and resilience presented in Sect. 2, we will now consider different variations in the function F = F(x) that make the "safe" equilibrium E less stable. Indeed, we consider four variants of perturbing the function F, denoted by F a , F b , F c and F d , see Fig. 2a-d, and explain how each measure captures the corresponding consequences of the stability and the resilience in each case.
To describe the four cases of stressing the dynamics, we first note that if F changes to F a or F b , as shown in Fig. 2a, b, then the Solow-Swan model still has the same equilibria and the global attractor is still equilibrium E. The basin of attraction for E is unchanged in both cases, but the returntime to E will be longer (for all perturbations in case F a and large perturbations in case F b ) since the functions F a and F b are closer to the line C x than F is. In case of F a , the local curvature near E has changed, but the curvatures of F, F b , F c and F d are identical in a small neighborhood of E.
In the cases of F c and F d , the Solow-Swan model has, besides equilibrium E, a locally stable equilibrium Fig. 2 The saved output per worker is given by the function F (green, dashed) of which we consider four variants F a , F b , F c and F d (red, dotted). The maintenance Cx is linear and illustrated by the black solid line. In a, b the economy converges to the globally stable equilibrium E given any perturbation. In c, d trajectories starting to the left of the unstable equilibrium E 1 will reach the origin E 0 , while trajectories starting to the right of E 1 will reach E.
at the origin (E 0 ) and an unstable equilibrium (E 1 ), see Fig. 2c, d. Observe that now E is only locally stable, and the economy can therefore converge to two attractors, E 0 and E. The new equilibrium E 1 produces a border of the basins for E 0 and E; all trajectories starting to the left of E 1 will reach the origin, while trajectories starting to the right of E 1 will reach E. We will focus our discussion on all measures considered in Sect. 2, as well as their relation to local measures. We do not specify the distance function for D, but the reader may think of it as the Euclidean distance for simplicity. The exact choice of distribution for perturbations is not important for the below results, but we assume that perturbations are chosen so that tested initial conditions are distributed with positive probability over the whole interval displayed in Fig. 2. Smaller intervals centered at E can also be considered without altering the main results in this section.

Local measures fail to detect changes to F b , F c and F d
Since the curvatures of F, F b , F c and F d coincide in a small neighborhood of E, any local measure, such as −λ max , based on the eigenvalues of the Jacobian matrix, will fail to detect these variations of F. Indeed, the same is true for any function having similar curvature as F near E. This underlines the non-sufficiency of considering only local stability measures; they will miss all changes of the system except those arbitrary close to E. Moreover, a local analysis of stability and resilience will not reveal the crucial fact that the modeled economy can converge to two completely different stable states in cases of F c and F d ; the origin, corresponding to no capital per worker, and E. If F changes to F a , then the local curvature also changes and a local measure will successfully record a decrease in resilience.

The measures P and D detect changes to F c and F d but fail in case of F a and F b
The size and the shape of the basin of E will be the same in the three cases of F, F a and F b , and therefore the measures P and D do not work. Indeed, D is the distance from E to the origin and P = 1 in all three cases, even though trajectories will return faster to E in the case of F than with F a and F b .
If F changes toward F d through F c , meaning that the unstable equilibrium E 1 moves from the origin toward E, then the size and shape of the basin of E changes and therefore both measures P and D (D is simply the distance between E and E 1 ) will record a decrease in stability, recall Fig. 2c, d. In the case of F d , when equilibrium E 1 is "very" close to E, then D gives the best warning signal that the economy may jump to a completely different stable state (the origin) due to a small perturbation. The size of basin measure P will not give such a clear warning signal since the basin is still large to the right of equilibrium E. Indeed, if we consider in addition that E 1 → E, then D → 0 and P → c for some constant c > 0.
The resilience measures R, R worst , P τ and D τ detect all considered variations of F In contrary to a local resilience measure that can only work in case F a , all suggested nonlocal resilience measures record the crucial changes of the dynamics for all considered variations of F. When F changes to F a and F b , the returntime to E will be longer since the functions F a and F b are closer to the line C x than F is. Therefore, the measures R and R worst , built directly on the returntime, record the destabilization of E. Moreover, in contrary to the basin of equilibrium E, the basin-time of E will, for any suitable choice of τ , decrease in size as F changes to F a or F b . The reason is that fewer initial conditions will reach E in time τ when the returntime becomes longer. Hence, the resilience measures P τ and D τ also record the destabilization of E as F changes to F a or F b .
When F changes toward F d through F c , then, in addition to the slower convergence of the trajectories toward E, the size of the basin of E also decreases as the new equilibrium E 1 comes closer to E. Because D τ is bounded from above by the distance from E 1 to E, and since R worst → 0 immediately if E 1 comes into the range of pre-specified perturbations, these two measures give the clearest warning signals as E 1 moves toward E. Measures R and P τ also records a decrease in resilience, but as they account for that E attracts all trajectories from the right in a similar way for all three cases F, F c and F d , they will not approach zero as E 1 → E. For R worst to work properly in this situation the set of predefined perturbations, L, should not include perturbations in the whole range from E to the origin, as that would force the measure down to zero too early as soon as E 1 is born.
We have seen that all four nonlocal resilience measures R, R worst , P τ and D τ are suitable for analyz- Fig. 3 The electro-mechanical machine consists of a wagon with mass m connected to a linear spring with stiffness k and a linear damper with damping factor c. The spring always forces the wagon to x = 0 where the spring force is zero. The wagon is also affected by a nonlinear magnetic force due to a magnet placed at x = a ing the economy described by the Solow-Swan model. Which of these four measures to use depends upon the question to be answered. The questions related to these measures, recall Sects. 1 and 2, should all be relevant for applications in economy.

A simple electro-mechanical system
In this section, we consider an electro-mechanical machine consisting of a wagon with mass m, free to move without friction horizontally in the x-direction. The wagon is connected to a linear spring with stiffness k and a linear damper with damping factor c. The spring always forces the wagon to x = 0 where the spring force is zero. The wagon is also affected by a nonlinear magnetic force from a magnet placed at x = a, see Fig. 3.
We assume that the magnetic force is given by k m (x − a) −2 where the parameter k m is a magnetic stiffness. The governing equation of motion is gives the location of the wagon. By introducing another dimension for the velocity, y = y(t) = dx dt , we obtain a two-dimensional nonlinear dynamical system on the form (1): for x < a and y ∈ R n . System (3) has (if k is large enough) two equilibria, of which one is unstable and the other is locally stable (E) and constitutes the "safe" attractor that we will investigate. There is also always an "unsafe" attractor at x = a where the wagon slams into the magnet and the machine damages. More advanced electro-mechanical systems having similar dynamic properties as system (3) arise, e.g., when modeling electric generators and electric motors. One is then typically concerned with at least four-dimensional systems of equations governing the motion of the rotor center (see the Jeffcott-rotor model) and with nonlinear magnetic forces acting between the rotor and the stator.
We will investigate stability and resilience properties of the locally stable equilibrium E as the spring performance decays, using measures suggested in Sect. 2. In a first case, we consider only a decrease in the spring stiffness k. In a second case, we assume also that the spring breaks (k is set to zero) if the spring is deformed to quickly, i.e., if the wagon moves faster than at a certain speed, y(t) ≥ y limit for some t > 0. In such case, the wagon slams into the magnet at x = a with probability 1.

Setting up relevant measures
In order to specify the set of perturbation to which the machine is to be tested for, one should consider, e.g., behavior of connected machines or dynamical systems as well as other relevant information that can be achieved through the environment of the system. One may consider only displacements, or only velocity impulses, or combinations of them. A velocity impulse may be due to a "smash" on the wagon, while a displacement may be due to a "smash" on the fundament of the machine. The less information you have, the more general set of perturbations may be necessary to include. As we have little information in this example, we chose to simply test the machine for perturbations normally distributed in both displacement and velocity and sample our initial conditions from a twodimensional normal distribution centered at E with standard deviation 5 in both the x-and y-direction. To numerically integrate trajectories from their initial conditions to the attractor, we used MATLAB's ode-solver ODE45 with standard tolerance settings. We integrated each trajectory until it either reached the small neighborhood given by a ball of radius 0.01 centered at E, or reached a state where x(t) > a − 0.01.
As the dynamics of the electro-mechanical system (3) is typically fast and the work of engineers should naturally focus on if the machine damages or not, rather than the time of recover given a perturbation, we invoke only our simplest resilience measure R and put focus on the pure stability measures. Indeed, as system (3) has more than one attractor, both the size and the shape of basin are crucial to investigate in order to understand stability and we therefore include the measures P and D in our investigation. We consider two versions of the latter, one using simply Euclidean distance (D Euc ) to compare displacement to velocity (x to y), and one using an energy-norm to find the distance (D energy ). Using Euclidean distance may seem naive since it compares velocity to displacement in an unmotivated way, but we chose to include D Euc anyway due to its simplicity.
More natural is the measure D energy that will find the least amount of energy needed to push the wagon from the locally stable equilibrium E into the "unsafe" attractor at x = a. This is certainly relevant information concerning stability of the machine. Giving the wagon an initial velocity y 0 implies an input of the kinetic energy of m y 2 0 /2. To find the potential energy input given by an initial displacement x 0 , we calculate the work needed for moving the wagon from E to the initial state x 0 . We assume that the wagon is moved slowly and that damping is neglected. We also assume that no energy is transferred back in cases when the wagon moves freely, due to the magnetic force, further away from E = (x eq , 0). The potential energy input given by an initial displacement from x eq to x 0 is then given by and so every initial condition (x 0 , y 0 ) is associated with the energy We can now define the energy-normed distance and hence our measure D energy as While D Euc is estimated by the radius of the smallest possible circle, centered at E, not lying entirely in  Figure 4 shows some initial conditions, trajectories as well as the circle estimating D Euc and the level curve estimating D energy . Besides giving a stability measure, the construction of D Euc and D energy reveals the most critical direction of a perturbation (in the sense of the corresponding measure) given by the red and blue arrows in Fig. 4.
The measure D energy gives the best warning before the machine damages Figure 5 shows the measures P, D Euc , D energy , R and −λ max for decreasing spring stiffness k. At k ≈ 0.054 the locally stable equilibrium E disappears in a fold bifurcation, leaving x = a as the only attractor and the machine damages as the wagon slams into the magnet for any initial condition. The results for no speed limit are given by the black solid curves, while dotted red curves correspond to the case with speed limit y limit = 2. Figure 5d gives also the local measure −λ max (dashed blue), which, due to its local nature, is unable to record effects of the speed limit and is therefore given by the same curve in both the case with and without speed limit.
The basin of attraction to the safe equilibrium E is relatively large until just before spring stiffness k reaches the critical value of the bifurcation. Thus, even though E is very close to the boundary of the basin, so that only a small perturbation to the wagon will damage the machine, the measure P will not give a clear warning, see Fig. 5a. This case is similar to Fig. 1c in Introduction; considering only the size of basin is simply not enough in such situation. However, one should note that the system is considerably damped, and by decreasing damping the measure P would tell more. This is because less damping allows for more oscillating motions which would force the basin to decrease smoother with decreasing spring stiffness k.
As the measure of the shape of basin, D Euc , records the distance to the boundary of the basin, it gives a slightly better warning signal, at least in the case without speed limit, c.f. Fig. 5a, b. However, in the case with speed limit the measure D Euc is, similar to P, nearly constant until just before the crash. This is because the Fig. 5 The considered measures as functions of the spring stiffness k in the case with speed limit (red, dotted) and without speed limit (black, solid). a P, b D Euc , c D energy and d R and −λ max × 0.12 (blue, dashed). A total number of 55,000 initial conditions were tested (1000 for each value of k) and each simulation (with and without speed limit) elapsed in 20 min. Parameters are m = c = k m = 1, a = 5 and y limit = 2.  Fig. 4b is bounded by the speed limit until just before the spring stiffness k reaches its critical value. Figure 5c shows the more natural energy-based measure of the shape of basin, D energy , that records the least amount of energy that can push the wagon into the magnet at x = a. This measure approaches zero in a linear way and delivers an early warning signal, in both the cases with and without speed limit, before the crash. One reason for this is that D energy puts focus on the "most dangerous" direction of perturbation toward the magnet at x = a, in contrary to D Euc that uses Euclidean distances; in terms of energy it is easiest to perturb the wagon in the positive x-direction, see the level curves in Fig. 4. To stress generality of the measure D energy , we point out that even though analytical expressions as (4) for the energy level curves may be difficult to find in general, numerical calculations of the energy of a displacement should be possible also for complicated systems.
The resilience measures R and −λ max shown in Fig. 5d record both a warning before the spring stiffness k reaches its critical value and the machine damages. Observe, however, that as R drops significantly as the speed limit is imposed, the local measure −λ max is unable to take into account the nonlinear effect of the speed limit. This is because the local behavior of the system is identical for all possible values of y limit . This highlights again the unsufficiency of measuring resilience only through local measures; clearly the system becomes less resilient as y limit → 0, and a suitable resilience measure should record this fact.
We conclude that, as in the case of the Solow-Swan model considered in Sect. 3.1, our nonlocal approach is able to record the considered destabilization of the system. Indeed, we believe that the measure D energy is a promising candidate for estimating stability of mechanical systems, and we further discuss its generalizations to more general situations of "small attractors" in Sect. 4.

Variance reduction through the local pivotal method
When perturbations are imposed randomly by drawing initial conditions from a given probability distribution, an efficient sampling technique may be used in order to sample the perturbations well spread over the state-space. Indeed, we have to sample a relatively large number of initial conditions in order to get a "good  6 The measures P (a) and R (b) as functions of the spring stiffness k with variance reduction through the local pivotal method (black, solid) and without variance reduction (red, dotted). 100 initial conditions, sampled from a two-dimensional normal distribution centered at E with standard deviation 5, are examined for each value of spring stiffness k. A total number of 5500 initial conditions were tested, and the simulation elapsed in 3 min. Parameters are set to m = c = k m = 1, a = 5 and y limit = ∞. (Color figure online) enough covering" of the desired part of the state-space, as well as small enough variances of our estimates, and each sampled trajectory is expensive in terms of computational time. In this section, we will show how the local pivotal method by Grafström et al. [22] easily can be applied in order to sample well-spread sets of initial conditions. This results in efficient variance reduction of our estimates, and thus less initial conditions need to be tested, resulting in shorter computational time.
The local pivotal method generates well-spread samples from discrete populations. To obtain a discrete population, we randomly draw a large number, say N , of independent points from the assumed probability distribution. We then sample the n desired initial conditions, out of the N points, n << N , using the local pivotal method. The local pivotal method is very fast, it can be run in any space-dimension, and an open source R-package named Balanced sampling, for taking the sample, is available [25]. Figure 6 shows the stability measure P and the resilience measure R when spring stiffness k decreases with and without using the local pivotal method for sampling the perturbations. In these simulations, we have used N = 10,000 and n = 100. A significant improvement can be observed. Compare also to Fig. 5a, d in which 1000 initial conditions were tested for each value of the spring stiffness k. The simulations presented in Fig. 6 elapsed in 3 min, of which the sampling through the local pivotal method elapsed in few seconds. (The corresponding simulation in Fig. 5 elapsed in 20 min.) By estimating the variance for spring stiffness k = 0.5 when 100 initial conditions were tested, we conclude that the variance was reduced by a factor of approximately 7, for both measures, when using the local pivotal method.

A stage-structured population model
In this section, we consider a stage-structured consumer-resource population model introduced by de Roos et al. [26]. The model is based upon the assumption that individuals are composed into two stages: juveniles and adults, depending only on their size. Juveniles are born with size s born and grow until they reach the size s max at which they cease to grow, mature, and become adults. Juveniles use all available energy for growth and maturation, while adults do not grow and instead invest all their energy in reproduction. The total juvenile biomass is denoted by J = J (t), while the total adult biomass is denoted by A = A(t). Both juveniles and adults forage on a shared resource R = R(t). The net biomass production per unit biomass for juveniles and adults equals the balance between ingestion and maintenance rate T according to where q describes the difference in ingestion rates between juveniles and adults, σ represents the efficiency of resource ingestion, H is half-saturation of the consumers, and the maximum juvenile and adult ingestion rates per unit biomass equal I max and q I max , respectively. The model incorporates stage-selective harvesting by allowing for separate harvesting rates on juveniles (h J ) and on adults (h A ), in addition to the background mortality rates d J and d A . Resource turnover rate is denoted by r , and R max is the maximum resource density. The model consists of the following three-dimensional dynamical system: where The function v(x) describes the maturation rate by determining how fast juvenile biomass is transferred into adult biomass. We adopt the parameter values H = T = r = 1, I max = 10, d J = d A = 0.1, q = 0.85, σ = 0.5, R max = 2, s born /s max = 0.01 [26,27].
Concerning the dynamics of the model, extensive numerical investigations indicate that there exists always a unique global attractor, and, depending on the parameter values, this attractor is either a positive equilibrium (J, A, R) = (J eq , A eq , R eq ) or an extinction equilibrium given by (J, A, R) = (0, 0, R max ) [27]. This means that the basin of attraction is always the whole R 3 + , and therefore the size and shape of basin will not deliver any information. In fact, we always have P = 1 and D = constant. However, as we will see in the following subsections, we obtain the sought after information by recording the returntime through the resilience measures given in Sect. 2.
Our first aim here is to show how to use the resilience measures presented in Sect. 2 to investigate the consequences of the populations ability to recover the globally stable equilibrium, given a perturbation, as harvesting pressure increases. Our second aim for this section is to show how to apply the same measures of resilience in order to compare the efficiency of different harvesting strategies when accounting for both yield and resilience.

Setting up relevant measures
For this application, we consider perturbations of sizes relative to the actual size of the population biomass at equilibrium. This choice is mainly based on our believes that it is reasonably that a fraction, rather than an absolute amount, of the population is eliminated through, e.g., illegal harvesting, tough weather conditions, sudden diseases or due to some other reason. Further motivations for considering perturbations of relative size is given by the fact that the populations equilibrium biomass can vary drastically as function of the harvesting rates and other parameters, making absolute measures less relevant.
We will consider two approaches of imposing perturbations on system (6). In a first case including resource, we test the system for perturbations in all three dimensions; i.e., juveniles, adults and the resource. This approach is relevant when studying how the system reacts when both the population and the resource are perturbed, e.g., by sudden environmental changes. In a second case excluding resource, we leave the resource at equilibrium and test the system for perturbations only in the two dimensions of juveniles and adults. This is relevant when studying how the system reacts when the population, but not (explicitly) the resource, is perturbed, e.g., by sudden diseases, by illegal harvesting or by errors in an implementation of a harvesting strategy. We proceed by specifying the four resilience measures R, R worst , P τ and D τ , given in Sect. 2, for the approach including resource as well as for the approach excluding resource.
Including resource We impose perturbations randomly by drawing initial conditions from a threedimensional normal distribution centered at the equilibrium (J eq , A eq , R eq ) with standard deviation given by 1 2 J eq , A eq , R eq . Initial conditions including nonpositive values are excluded. For the measures P τ and D τ based on the size and the shape of the basin-time, we have to specify τ . In general, to find a suitable value of τ one should recall the sought after information and that P τ estimates the probability that the population returns to equilibrium before time τ given a random perturbation, while D τ estimates the smallest perturbation that can force the population to not return in time τ . The value of τ should also be in the range of typical returntimes of the system, and for our application we will use τ = 5 for both measures. To find a relevant version of the "shape of basin-time measure", D τ , we follow the above reasoning of relative distances to the population size and define For the resilience measure R worst , estimating "worst resilience using the slowest recovery," we use a smaller set of perturbations (L) than for the other measures. In particular, we restrict perturbations to a smaller neigh-borhood of the attractor by taking L as the subset of the above defined initial conditions x = (x 1 , x 2 , x 3 ) satisfying Doing so, we ensure that a perturbation for R worst never removes more than half the juvenile biomass, or half the adult biomass, from the population; indeed, if J eq = A eq = R eq = 1, then (8) is a sphere centered at (1, 1, 1) with radius 1 2 . Excluding resource In this case, initial conditions are drawn from a two-dimensional normal distribution centered at (J eq , A eq ) with standard deviation 1 2 J eq , A eq . The measure D τ is defined through (7) but without the third term under the square rot. Similarly, the set L for R worst is defined as those initial conditions that end up in the ellipse given by (8) but without the resource term. Figure 7 shows the basin-time, some trajectories, the ellipse (8) defining the set L (blue), as well as the ellipse giving the estimate for the measure D τ rel (red) in the case excluding resource.

(a)
To numerically integrate trajectories from their initial conditions to the attractor, we used MATLAB's ode-solver ODE45 with standard tolerance settings. In both the case including resource and excluding resource, we integrated each trajectory until it reached the small neighborhood given by the ellipse (8) but with right-hand side 0.01.
All resilience measures R, P τ , D τ rel and R worst give similar results Figure 8 shows the suggested measures P τ , D τ rel , R and R worst together with −λ max for increasing harvesting pressure. We assume here equal harvesting rates on juveniles and adults, i.e., h J = h A . The case with resource is given by black solid curves, while dotted red curves correspond to the case without resource. In Fig. 8a, the dashed blue curve gives the local resilience measure −λ max .
All five measures agree on a similar result; the resilience first increases with harvesting pressure, up to a maximum at approximately h J = h A = 0.8, after which it decreases to its lowest level before the population goes extinct. All five measures of resilience behave as similar "smooth" functions of the harvesting rates, which indicates that system (6) has a "weak" nonlinearity; the linearized system probably agrees well with the original one in a relatively large neighborhood of the investigated equilibrium.
When comparing the approaches of including and excluding the resource in the measures constructions, we see that the curves corresponding to the approach including resource are lower than the curves for resource excluded. This is true mainly for relatively high harvesting rates and for the measure R, see Fig. 8a. Thus, excluding resource will in this situation be in favor of a slightly higher harvesting pressure.
Comparing harvesting strategies using the measures R, R worst , P τ and D τ rel In this section, we use the proposed measures of resilience to compare the strategy of harvesting only adult biomass (adult harvesting), h J = 0, h A > 0, to harvesting both juvenile and adult biomass at the same rates (equal harvesting), h J = h A > 0. For a given harvesting strategy, we can find the yield (catch) as Fig. 9 The yield as functions of the five measures of resilience for increasing harvesting pressure. Equal harvesting (black, solid) and adult harvesting (green, dotted). a R and −λ max × 1 3 (blue, dashed), b R worst , c P τ and d D τ rel . The dots correspond to the strategy h J = h A = 0.8 which seems to be one of the best possible strategies. In this simulation, 2000 initial conditions are examined for each harvesting rate combination. We consider here the case including the resource in the measure constructions.

(c)
To find out which strategy that is most efficient when accounting for both yield and resilience, we plot the yield as functions of the five resilience measures P τ , D τ rel , R, R worst and −λ max in Fig. 9. It is clear, in terms of all resilience measures, that equal harvesting performs better than pure adult harvesting. This is true since for any given yield (close to the maximum yield that can be obtained) equal harvesting gives much higher resilience, in terms of all measures, than adult harvesting. In fisheries management, harvesting of juvenile fish is a debated topic, and the present result strengthens previous results by Lundström et al. [21] arguing for equal harvesting. In particular, in that paper the authors evaluate harvesting strategies with respect to both yield and conservation using two population models, of which one is model (6). They use four measures to quantify conservation, and their resilience measure is similar to the presented measure R.
We have seen that all nonlocal resilience measures R, R worst , P τ and D τ , as well as the local measure −λ max , give similar results in the sense of recommending the same harvesting strategy. This does not mean that the presented nonlocal approach is superfluous. Indeed, it tells that the system behaves "smooth and linearly" near the attractor and that we may not see sudden jumps or other unexpected behavior in the dynamics. This is valuable information that cannot be obtained through a local measure based on a linearized system. Which of the suggested nonlocal measures to include in a study depends upon the question to be answered. All questions related to these measures (see Sects. 1 and 2) should be relevant for a variety of applications in biology.

Applications to a high-dimensional system
We end the results section by demonstrating that our approach works fine also in the setting of highdimensional models in which chaotic attractors, as well as several coexisting equilibria and periodic solutions, may arise. Indeed, we consider the dynamical system where d is the dimension and b i j and a i jk are constants with values drawn from a normal distribution with mean zero and variance 1. Systems similar to (9) arise in biology when describing the evolution of a multidimensional phenotype [28]. The system is dissipative, which means that there is a bounded region in the state-space to which all trajectories enter after sufficient time. We will consider the two cases of dimension d = 20 and dimension d = 50. If −λ is large enough, then system (9) has a locally stable equilibrium at the origin. We will implement the stability measure P as well as the resilience measure R for evaluating changes in the stability and the resilience of system (9) when λ varies. When d = 20, system (9) has at least three attractors; the equilibrium at the origin, another equilibrium (E 1 ) as well as a periodic orbit, see Fig. 10a. In this case, we investigate how the stability of equilibrium E 1 is lost when λ decreases and more and more trajectories attracts to the origin (Fig. 10c). For each tested parameter value, we randomly draw 200 independent initial conditions from a 20-dimensional normal distribution centered at E 1 with standard deviation 10, and integrate until the trajectory reaches a distance 0.01 from E 1 or until an end-time 10 is reached. We consider 40 different values of the parameter λ, and thus we test a total of 8000 initial conditions. The elapsed time of this simulation was 3 h.
When d = 50, system (9) has, in addition to the equilibrium at the origin, a chaotic attractor, see Fig. 10b. In this case, we investigate how the stability of the origin is lost when λ increases (Fig. 10d). For each tested parameter value, we randomly draw 200 initial conditions, from a 50-dimensional normal distribution centered at the origin with standard deviation 2, and integrate until the trajectory reaches a distance 0.01 from the equilibrium or until end-time 1. We consider 20 parameter values, and thus we test a total of 4000 initial conditions. On the chaotic attractor, the ode-solver takes small time-steps resulting in long computational time. The elapsed time of this simulation, performed with MATLAB's ode-solver ODE45 on standard tolerance settings, was 8 h. To speed up the simulations, one may produce the code in, e.g., JAVA or C++, incorporating further optimizations of the ODE-solver, see, e.g., Mulansky [29] and references therein.
In this section, we have chosen to present implementations of the measures P and R. However, all other measures suggested in Sect. 2 can be implemented successfully as well, with similar elapsed time, on system (9).

Discussion
We have considered two nonlocal stability measures and four nonlocal resilience measures based on properties of the systems basins, on the introduced concept basin-time, as well as on the returntime to the attractor. We showed how to quantify stability and resilience using these measures on four different dynamical systems modeling economy, electro-mechanics, stagestructured populations and the evolution of a multidimensional phenotype. We compared and discussed our results in relation to local measures through the paper, and explained which measures that give what information in each case. We concluded that, in contrary to a local approach, our measures give a good understanding of the stability properties of the attractor in all examples considered. Due to this fact and due to its simplicity, we believe that the presented methods have a large potential to increase the understanding of stability in a wide range of research fields involving applied dynamical systems. We underline that an implementation of our approach reveals, with a high probability, if there is another unexpected attractor in the range of the considered perturbations. Indeed, as our method implies testing the system for a large number of initial conditions even hidden attractors (see, e.g., [1]) can automatically be found (when producing Figs. 5 and 8 more than 50,000 initial conditions where tested in each system).
When deciding which measures to include in a stability analysis, we first note that measuring the size and the shape of the basin should always be done in some way when dealing with nonlinear systems. Thus, if no analytical estimates of the basin is available, measures P and a suitable version of D should always be relevant. In case the returntime is important, a resilience measure should be included. Which of the resilience measures to use depends on the question to be answered. We point out that one may get estimates of all suggested measures from the same samples of perturbations. It is therefore almost no more expensive to estimate all measures than just one single measure.
It seems difficult to construct an overall "best measure" in the general setting of a wide range of applications, but it is of course possible to build integrated measures that records the relevant information in even more general settings. However, that may ruin the simplicity in the measures construction, leading to results that are harder to interpret and implementations that are more complicated. Therefore, we choose to stay with simple measure constructions and recommend instead to use several simple measures simultaneously.

"Small" attractors can be handled simply as an equilibrium
We have demonstrated our approach on four models in which the stability and the resilience were investigated on the simplest attractor, an equilibrium. The presented ideas are, however, applicable also in cases of general attractors, not only in the setting of an equilibrium. In such general case, it may be more difficult to implement the ideas as it is usually not easy, e.g., to calculate the distance from the attractor to the boundary of the basin, and to determine when a trajectory has returned to the attractor. However, all the suggested measures are simple to calculate for some complicated attractors. Indeed, if the attractor is "small" in the sense that it can be contained in a ball B(y, r ), with radius r and center y, where r is small in comparison with the perturbations one wish to test the system for, then one can use a similar procedure as if the attractor was an equilibrium at y. This is independent of the type of the attractor; it can be periodic, quasiperiodic or even chaotic. To explain this, recall that when the attractor is an equilibrium at y, then we may use a ball B(y, r ), r very small, to define the neighborhood determining when a trajectory has returned to the equilibrium and the simulation can be stopped. In the setting of an arbitrary attractor contained in a ball B(y, r ), we can proceed similarly; when the trajectory has entered B(y, r ) and stayed there a "sufficient" time, then we can assume that the attractor has been reached. The returntime is thereafter determined through the moment when the trajectory entered the ball B(y, r ) for the last time.
This is an important point because these situations occur frequently in applications of dynamical systems. Indeed, equilibria often bifurcate into "small" but complicated attractors due to variations in model param-eter values. To explain a typical such situation, consider a "perfect" electric generator in which the rotor is balanced and no other external forcing on the rotor is present. In such situation, the rotor rotates at a stable equilibrium. However, in reality there are always imperfections such as mass imbalance on the rotor, shape deviations of the rotor and the stator as well as external forces transferred from nearby machines such as turbines. When adding the effects from such phenomena to the perfect rotor model, the stable equilibrium will most probably bifurcate into more complicated dynamic behavior, e.g., mass unbalance implies simple periodic motion. Engineers have to keep track of these effects when designing the machine. Typically, large amplitudes are unwanted and the generator should be constructed so that the attractor is "small" in the above sense, i.e., so that the rotor stays close to the previously existing equilibrium (the attractor in case of the perfect machine) under operation. As our approach can handle these situations in a simple way, we believe that the suggested measures will constitute useful tools for engineers when designing new generators and motors as well as when performing restore work and understanding behavior of existing machines. Simple versions of P and D have already been used in this manner, see Lundström [3] and Lundström and Aidanpää [4].

Alternative measures and topics for future research
We have considered two nonlocal stability measures and four nonlocal resilience measures. As we believe in a value of straightforward interpretations of the measures in terms of application, our choice of suggested measures relies heavily on the simplicity in the construction. Anyway, we choose to define the resilience measures R and R worst through the reciprocal of the returntime, see Sect.
2. An even more straightforward way to record the systems ability to return fast would be to consider the returntime explicitly. For bistable systems however, one then has to integrate the effect of the "unsafe" perturbations, of which trajectories not return to the attractor, in a suitable way. By taking the reciprocal of the returntime and consider resilience as rate of return, we solved this by giving zero resilience to unsafe perturbations. Another reason for considering the resilience, in place of the returntime explicitly, is the more direct connection to the literature and established resilience measures such as the real part of the largest eigenvalue. In connection to this discussion, we mention that the function 1/(T (x) + t ), used in the definitions of R and R worst , may be replaced by another decreasing function as well, such as e.g. e −T (x) , decreasing faster to zero and therefore punish more for long returntimes, or (T (x) + t ) −a for some exponent a > 0. If the attractor is globally stable, one may also take the mean of returntimes before taking the reciprocal [21]. There are many suggestions on alternatives to standard resilience measures in the literature. Neubert and Caswell [8] consider several alternative measures (e.g., reactivity and maximum possible amplification) to resilience together with local methods for calculating them. They demonstrate their measures on several examples and discuss also disadvantages of local methods for calculating stability and resilience measures. Extensions of the work of Neubert and Caswell have been done by, e.g., Arnoldi et al. [11,13]. These papers consider small perturbations and put focus on local measures. Mitra et al. [20] consider an integrative quantifier of resilience based on both local and nonlocal theory. Their measure builds upon the three aspects of resilience by Walker et al. [30]: latitude, resistance and precariousness. Hellmann et al. [31] consider the nonlocal measure survivability, reflecting the likelihood that the trajectory of a random initial condition leaves a certain desired (or allowed) subset of the systems state-space. A related paper by van Kan et al. [32] considers constrained basin stability which discusses measures related to P and P τ . Interesting topics for future research are to consider nonlocal versions of existing local measures and also to compare and discuss relations between existing measures and those considered here.
Even though our approach is already relatively fast (the simulations in Fig. 5 elapsed in 20 min, Fig. 8 in 1 h, and that of the high-dimensional system in Fig. 10 in three respective 8 h), it is a topic for future research to make the approach faster. All simulations presented in this paper were performed in MATLAB, using the ode-solver ODE45, on a standard laptop computer. To speed them up, one may produce an optimized code in, e.g., JAVA or C++, using a fast ODE-solver, see e.g. Mulansky [29] and the references therein. To increase efficiency through the way perturbations are inferred, we have considered well spread samples of initial conditions through the local pivotal method by Grafström et al. [22] in Sect. 3.2. Further improvements in this direction may be obtained by optimization of the quotient n/N , i.e., the quotient of the number of sampled initial conditions and the number of points in the discrete population to which the local pivotal method is applied. Importance sampling for choosing the set of initial conditions may also be useful and should be possible to implement through the local pivotal method as well. In general, improvements in the estimates by further refine the selection of the initial conditions, randomly or systematically, are interesting topics for future research.
Another interesting topic for future research would be to create an efficient approach for finding estimates of the neighborhood to attractors in nonlinear systems in which the linearized system approximates the original system well. One could, for example, classify dynamical systems based on the nature of their nonlinearity so that "weak" means that the linearization holds in a large neighborhood, and "strong" means that it holds only in a small neighborhood. Such results would complement the linear local stability approach by giving knowledge about the magnitude of perturbations that can be handled by the already well established local stability and resilience measures. A related field of research concerns how changes in a systems recovery rate, toward an attractor, can be used in order to predict changes in the size and shape of the basin of attraction. The recovery rate usually decreases as the basin shrinks to nil and a regime shift is approached, a phenomenon known in physics as "critical slowing down." This can be observed in both the electro-mechanical system considered in Sect. 3.2 and in the population model considered in Sect. 3.3. Indeed, from Figs. 4d and 8a-d, we see that the resilience decreases, according to all tested resilience measures, before stability is lost in a bifurcation. Van Nes and Scheffer [33] showed that, in six ecological models, recovery rates from small perturbations decrease as a regime shift was approached. For more on numerical and experimental approaches for detecting bifurcations before they occur by exploiting the critical slowing down phenomenon, we refer to Lim and Epureanu [34], Ghadami and Epureanu [35,36] and Van de Leemput et al. [37].

Concluding words
We have presented two nonlocal stability measures and four nonlocal resilience measures in order to comple-ment the widespread local stability-measures approach based on linearizations. The measures are simple and easy to implement on a standard laptop computer. All measures give answers to questions which should be relevant for a variety of applications.
We have considered the two nonlocal stability measures P and D, based on the size and shape of the basin of attraction, respectively. These measures account for both large and small perturbations; P reflects the probability that the system returns to the attractor given a perturbation, while D estimates the smallest perturbation needed to push the solution to another attractor. The measure D calls for suitable ways to compare distances in the state-space and, for mechanical systems, we constructed the version D energy reflecting the least amount of energy needed to push the systems solution into another attractor. We also introduced the version D rel considering relative distances when dealing with population dynamics.
We proceeded by considering the systems returntime to the attractor given a perturbation and defined the simple nonlocal resilience measures R and R worst . In contrary to local resilience measures, based on linearizations and eigenvalues of the Jacobian matrix, these measures account for large as well as small perturbations; R gives the expected rate of return given a random perturbation, while R worst reflects the slowest rate of return given a random perturbation. Moreover, we introduced the concept basin-time as the subset of the basin of attraction from which all trajectories return to the attractor in time τ . By replacing the basin of attraction with the smaller set basin-time in the constructions of the stability measures P and D, we derived two more resilience measures, P τ and D τ . The measure P τ reflects the probability that the system returns from a random perturbation in time τ , while D τ estimates the least perturbation from which the system will not return in time τ .
We demonstrated our approach on four mathematical models (the Solow-Swan model of economic growth, a system modeling an electro-mechanical machine, a stage-structured population model and a high-dimensional model) by "stressing" the dynamics in each model through variations in parameter values. During these stressed scenarios, we discussed which dynamic behavior that can be recorded, and why, by each measure. For example, all nonlocal resilience measures detect all considered variations in the Solow-Swan model, while local measures failed to record the destabilization in three out of four stressed cases. The measure D energy gave the clearest warning signal before the electro-mechanical machine damaged due to a reduction of a spring stiffness, and all resilience measures agreed on similar results as harvest pressure increased in the stage-structured population model; they all gave arguments in favor of equal harvesting of juveniles and adults, compared to pure adult harvesting. In our examples we have seen, and in general we believe, that the three concepts size of basin, shape of basin and returntime constitute an important fundament when building stability and resilience measures. We can also conclude that no measure is "the best" measure through all applications and thus complementary measures, such as P measuring the size of basin and D measuring the shape of basin, should be considered simultaneously.
Although the presented approach is more timeconsuming than classical linear methods, it was no problem to run it on a 50-dimensional dynamical system showing complicated dynamics. We have also demonstrated how to improve the estimates of our measures, and thereby reduce the computational time, by taking well-spread samples of initial conditions through the local pivotal method.
Even though valuable information of a systems stability can easily (using the presented approach and today's computers) be recorded from properties of the basin of attraction and the returntime, the use of nonlocal measures of stability and resilience is today rare in the literature. We believe that the simple ideas and measure constructions presented here will constitute an important step in order to quantify stability and resilience through nonlocal methods and thereby help to fill this gap in the literature.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestric ted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.