Stability of fixed points in an approximate solution of the spring-mass running model

We consider a classical spring-mass model of human running which is built upon an inverted elastic pendulum. Based on our previous results concerning asymptotic solutions for large spring constant (or small angle of attack), we construct analytical approximations of solutions in the considered model. The model itself consists of two sets of differential equations - one set describes the motion of the centre of mass of a runner in contact with the ground (support phase), and the second set describes the phase of no contact with the ground (flight phase). By appropriately concatenating asymptotic solutions for the two phases we are able to reduce the dynamics to a one-dimensional apex to apex return map. We find sufficient conditions for this map to have a unique stable fixed point. By numerical continuation of fixed points with respect to energy, we find a transcritical bifurcation in our model system.


Introduction
Running is one of the most fundamental means of bipedal animal locomotion. Superficially, it may seem mundane. However, it is a result of complex interactions between neural and muscular systems [6,13]. Remarkable efficiency and stability of animal running over rough terrain attracts the attention of scientists throughout various disciplines spanning from biomechanics, medicine, applied mathematics, to robotics [3,5,16]. Running can also be considered from several other points of view that are not of biomechanical nature. For instance, in competitive sports, one is usually interested in optimization problems based on traversing a given distance in the shortest time. A classical model of Keller [17], which was based on early research of Hill [15], provides an interesting account of this optimization. Some modern approaches to that subject can be found in [28,21,1].
To understand the very process of running, one usually constructs a model that can capture its relevant properties. This can be done on a spectrum of various complexity levels. Here, we are concerned with the, so-called, spring-mass running model for bouncing gaits that has been introduced by Blickhan [4] and thoroughly investigated by MacMahon and Cheng [18]. This model can be considered conceptual in the sense that it is a simple realisation of the very motion of a hopping leg. It consists of a spring-loaded inverted pendulum (SLIP) that swings during the phase of contact with the ground (support phase) and is then ejected upwards into the aerial phase (flight phase). The exceptional accuracy of that model is based on the validity of the reductionist approach. It stems from dimensional analysis of fundamental physical quantities describing the gait of various, not necessarily bipedal, animals. Some experimental studies concerning the SLIP model have been conducted for example in [7,8,14]. This model, due to its robustness, has also been extensively used in robotics to design and construct various legged robots (see the seminal paper by Raibert [22] and a modern review [2]). Furthermore, in the literature one can find a number of important generalizations of the SLIP model such as varying attack angle [9], introducing control [24,27,26], damping [23], and multi-legged hoppers [10].
In spite of being a relatively simple model, SLIP does not posses a closed form analytical solution. Therefore, in order to study its properties one usually uses quantitative dynamical systems analysis or numerical methods. To gain more insight, an approximate analytical solutions can be obtained. This has been done in the work of Geyer, Seifarth, and Blickhan [11] which is one of the basis of our subsequent results. Moreover, some approaches to analysis of the current model based on approximate solutions have been conducted in [25,12,23] and in our earlier works [20,19]. There, we have conducted a rigorous Poicaré-Lindstedt asymptotic analysis based on the perturbation theory in the regime of large spring stiffness and small angle of attack. Moreover, we have also found some expansions in the limit of large horizontal velocity. A thorough review of mathematical studies of various locomotion models has been conducted in [16], to which we refer the interested reader.
In this paper, we obtain analytical solutions to the model by concatenating our approximate solutions in the support phase with the solutions of the flight phase, with the latter phase representing parabolic trajectories. Using these solutions we obtain a onedimensional Poincaré map that takes an apex in the flight into the next apex. We then investigate the existence of fixed points of the Poincaré map so constructed and find sufficient conditions for the stability of fixed points. This shows that this model can describe a stable running gait. Our approach is similar to the one conducted in [11] where the authors use a different approximation strategy stemming from the principles of mechanics, which, however, obscures the error one makes by necessarily approximating the solutions. In contrast to that, our approach is based on rigorous asymptotic analysis that clearly formulates the region of validity of the model [20,19]. Moreover, as our numerical studies indicate, the fixed point of the Poincaré map undergoes a transcritical bifurcation with increasing energy. A thorough investigation of this phenomenon, as well as the investigations of the existence of fixed points corresponding to asymmetric solutions, will be the subject of our future studies.
The rest of the paper is outlined as follows. In Sec. 2 we introduce the governing equations of the spring-mass model of running. Then, in Sec. 3, we introduce approximate solutions of our model equations obtained by means the Poincaré -Lindstedt method. These equations are then used in Sec. 4 to obtain a one-dimensional map which is also analyzed in this section. Finally, Sec. 5 concludes the paper. √ X 2 + Y 2 is the radial while θ is the angular position of the point of mass, and additionally ∆θ is angle swept during stance. Moreover, we assume that θ T D = −α, where α ∈ (0, π/2).

Spring-mass running model
In this section we will derive nonlinear equations, which describe running during the stance phase (see Fig. 1 for notation). Using nondimensional polar coordinates (L, θ), the Lagrange function of the stance phase is given by where the parameter K denotes the nondimensional spring stiffness. So the two Euler-Lagrange equations are as follows where P = L 2 θ denotes the nondimensional angular momentum.
Observe that the dimensionless form of the mechanical energy during the contact phase denoted by E s (see (1)) is given by On the other hand, the total energy E s given by equation (3), for small angles θ, takes an alternative form Consequently, the Lagrangian (1) is also modified. We should note that for small angles approximation, change of the angular momentum of the point mass we consider is zero, that is From equations (4) and (5), we now obtain a set of conditions that relate the system state at take-off to the state at touch-down. During take-off and touch-down, denoted by TO and TD, respectively, we have the radial symmetry with rest length L = 1 at each phase transition and so from (4) we have then Furthermore, after integration (5), the following equation θ T O/T D = cL −2 T O/T D = c holds, where c is some real constant. In our case: forward locomotion the constant c is positive.
We can now write Note that in the full system, in addition to (6), (8) and (9), we require (see also Fig. 1) where only the angle ∆θ swept in stance cannot simply be expressed by the state at the time of touch-down. Moreover, it should be observed that for small angles θ the conservation of energy and angular momentum do not enforce symmetry in the stance phase i.e. through the arbitrary value of θ parameter, asymmetric solutions are also allowed.
In the next section we will derive the formula for ∆θ.

Constructing approximate solutions
An asymptotic analysis of the main equations (2) with the use of the Poincaré -Lindstedt series was carried out in [20]. The solution of the problem is based on the perturbative expansion related to the significant spring stiffness (K → ∞). To simplify matters we set = 1/ √ K ( → 0 + ) and introduce In what follows we will operate on the order of O( 2 ) as → 0 + truncating higher order terms. All equations should then be understood as approximations in this asymptotic sense. For clarity we will retain using the equality sign instead of more rigorous approximate equality after performing the truncation.
• Radial motion during stance L(t) (see (23) in [20]): The L(t) approximation is as follows where the initial conditions have the form for α > 0 with X and Y are the horizontal and vertical velocities made dimensionless by the factor (acceleration of gravity × length) 1/2 , which are called Froude numbers in fluid mechanics. In addition, it is worth noting that from (3) where θ d refers to the angular momentum P at touch-down (T D).
The general solution (12) of the radial motion L(t) during stance describes a sinusoidal oscillation around L = 1 + A( ) with amplitude B( ) and frequency −1 ω( ), But, as shown in Fig. 2 on the left, the solution only holds for L(t) ≤ 1.
In formula (12), the parameter t ranges from 0 to the contact time t C given by where and → 0 + . Observe that L(0) = L(t C ) = 1 and L(t C /2) = 1 − ∆L M AX (see Fig. 2). The maximum leg deflection during the stance phase, denoted by ∆L M AX , is given by the difference of the amplitude B( ) of the radial motion L(t) and the shift A( ) of the touch-down position, i.e.
So • Angular motion during stance θ(t) (see (24) in [20]): The θ(t) approximation with the initial conditions (13) is as follows Observe that θ(0) = θ T D = −α and θ(t C ) = θ T O = ∆θ + θ T D = ∆θ − α (see Fig. 2 on the right). • Angle swept during stance: From formulas (20) and (17) we obtain the angle ∆θ swept during stance where C( ) is already defined in (18). Therefore, the angle swept ∆θ is uniquely defined by the system state at touch-down (L T D , θ T D , θ T D ) and the spring stiffness parameter K. Moreover, the landing angle θ T D = −α affects ∆θ also by determining the radial and angular landing velocity. Finally, let's observe, that the Taylor expansion of (21) to the second order, about = 0, is given by where → 0. Since ∆θ approximation given by formula (21) is of the order 2 , all terms of the order O( 3 ) in formula (22) are omitted.
An important case of stability analysis, considered also by Geyer, Seyfarth and Blickhan (cf. [11]), is In the special case we get the shift A( ) = 0 and the leg deflection parabola is equal to the amplitude i.e. ∆L M AX = B( ). Moreover, the angular velocity at touch-down is θ d = √ cos α and the contact time in formula (17) is equal to 2π /(2 − 2 θ 2 d ). Putting C( ) = 0 into the equation (21) or equivalently θ d = √ cos α into the equation (22), we get (see also (15)) Parameter combinations {α, E s (α), K(α, E s )} leading to a symmetrical stance phase, i.e. ∆θ = 2π , will be given, for the case of C( ) = 0 in Sec. 4.4. By setting the stance phase limits to symmetry, we derive in the next section approximations of the dimensionless leg stiffness K.

Symmetric solutions
By limiting ourselves to the second order expansion (22), the solution ∆θ = 2α comes down to the quadratic equation for √ K. A physically reasonable solution gives the following expression to approximate the spring stiffness where Using α → 0 + in (26), the second component tends to zero and we get 1 The same approximation of K as (27) was obtained in [20] and [29]. Moreover, in [19] it was shown that the leading order of the expansion of K as In the symmetric case we have a particular boundary value problem to solve. Let (θ(t, K), L(t, K)) be the solution of the system (2) with (13). Find K * and the smallest time t * > 0 satisfying θ(t * , K * ) = α, L(t * , K * ) = 1.
The above problem can be easily solved numerically using the shooting method, as was done for example in [18]. Graphical illustrations confirming the validity (27) as an approximation of the solution to the boundary problem (28) are presented in [20] and [29]. In Fig. 3, on the left, we can see what the approximations K given by (25) and K given by (27) look like. Obviously, if α → 0 + then K * → ∞ and both approximations of K * are very good. However, for moderate values of parameter α, K is slightly better than K. Moreover the right side of Fig. 3 illustrates the ratio of K * to both approximations K and K for varying α. Indeed, the approximation K behaves better than K in the presented range of the parameter α, i.e. 0.01 ≤ α ≤ 0.25. It turns out that the approximation K works very good for α < 0.05 while K is accurate only for α < 0.01. 1 Since the angular velocity at touch-down θ d (α) may depend on the angle α, to be asymptotically consistent we should write where θ d (0) = X . For simplicity we use the formula (27).  (25), and the second K approximation, calculated from (27), for varying α from 0.1 to 0.35. On the right: ratio of K * to K and K for varying α from 0.01 to 0.25. Here X = 1 and Y = 0.1.

Analytical apex return map
At the beginning of this section we will deduce the analytical solution to calculate dependencies between two subsequent apex heights. Identification of the system during the phase transitions: from flight to stance and from stance to flight will allow the derivation of the apex return map Y i+1 (Y i ). In Sec. 2 we showed that assuming small angles, the system state at take-off (T O) is related to the state at touch-down (TD) through a set of the equations (see (6), (8), (9) and (10)) Starting with the correct initial values, the mapping between the height Y i of the apex i and the touch-down state in dimensionless polar coordinates is expressed as where E s is the corresponding dimensionless energy of the system prior to touch-down, given from (3) by E s = (θ 2 d + L 2 d )/2 + cos α. Since the Froude numbers X and Y depend on Y i , then the radial velocity L d and the angular velocity θ d also depend on Y i (see (14) or (30)). Then, an additional mapping is required between the take-off state and the height Y i+1 of the apex i + 1, which is By using both mappings (30) and (31), the apex return map function Y i+1 (Y i ) can be constructed and takes the following form after some simplifications where ∆θ i = ∆θ(θ d , L d ) is given by (21) or (22). To mark the periodicity for this running model, we introduced the index i to the notation ∆θ i . So ∆θ i is the angle swept during stance between two subsequent flight phases, when the apex height Y i transforms into Y i+1 . Now you can see how ∆θ i = ∆θ(θ d , L d ) depends on Y i . From formulas (14) and (30) we have Moreover, observe that (see (15)) If the apex height Y i+1 is less then landing height (Y i+1 < cos α), the leg will get into the ground, and the stumble will occur. So the apex return map only exists otherwise. Regardless of the angle swept during stance, from equations (29) the kinetic energy 1 2 ((L ) 2 + L 2 (θ ) 2 ) is equal at touch-down and take-off. For asymmetric contact phases, On the other hand, during symmetric stance phases (Y T D = Y T O ) the corresponding shifts in energy at touch-down and take-off compensate each other, so the system energy E s is constant. To restore the conservative nature of the model, change resulting from asymmetry, needs to be corrected. Thus, the appropriate horizontal velocity in (31) must be given during take-off so that the kinetic energy (E s − Y i+1 ) reduces the lack of potential energy at the end of stance phase. For the new apex hight Y i+1 from (32) the adjusted E s is used.

Existence of fixed points
In the following section we need to check whether apex states Y i+1 restricted by symmetric contracts can be found. Solution (31) with condition (29) for θ T O = α leads to the corresponding apex height of fixed points where θ * d and L * d represent the values of θ d and L d given by (33) at the fixed point Y * . Since θ * d and L * d are linked by the relationship (34), the apex height of fixed points can be expressed in the following equivalent form Substituting (36) and (37) into equations (33) yields and also It is easy to see that θ * d cos α + L * d sin α > 0. So, it follows that equations (38) and (39) hold if From (15) and (40) we have Hence the system energy E s must satisfy the inequality E s ≥ E min For the system energy E s = E min s , the apex height of fixed point Y * is cos α, which is equal to the landing height. As the energy increases (E s > E min s ), the apex height increases, but it never exceeds the upper boundary, i.e. E s .
If the attack angle α is fixed, then the stiffness approximation given by (25) and (26) is the lowest for the minimum energy system given by (42), i.e. where Whence, we proved that there is a fixed point under a certain assumption: where α ∈ (0, π/2). Then there exist initial conditions for the stance phase Additionally we get where ∆ min = π 2 (cos α − 4α sin α) + 32α sin α.
The fixed points can only exist if a minimum energy E min s and a minimum stiffness K min are exceeded.
Considering θ * d = √ cos α case, if α → 0, the minimal system energy E min s is approaching the value of E s = 1.5. To get a dimensional overview, the system energy then has the value: mgl 0 E s ≈ 1125 J, where body mass m = 75 kg, gravitational acceleration g = 9.81 m/s 2 , and leg length l 0 = 1 m. If additionally the initial apex height is set to Y 0 = 1, then in this case it corresponds to the initial velocity x = X √ gl 0 = 2(E s − Y 0 )gl 0 , which is 3.13 m/s.

Stability of fixed points
Let us denote mapping (32) as Y i+1 = f (Y i ). Stable solutions of fixed point Y * fulfill the following condition |f (Y * )| < 1.
So to prove stability, at least one the parameter set (α, E s , θ * d ) leading to solutions of Y * satisfying (49) should be identified. Starting with (32), we get by using ∆θ i = 2α for symmetric contact phases. Later in the paper, we take the notation Since the expression in square brackets in (50) always remains positive, then (49) transforms into the following condition for the angle swept during stance which shows that higher or lower apex heights must be compensated appropriately larger or smaller amount of angular sweep ∆θ i . However, for stability reasons, the positive derivative d i ∆θ * cannot be greater than 2 sin α + 2 (Y * − cos α)(E s − Y * ) −1 . Then we substitute the equations fixed point Y * (36) and (37) satisfying the condition (45) and denoting Next, from (22) it follows that where K is given by (25).
and applying (34) we have and from (55) can be further deduced with D (α, E s , θ * d ) given by (52).
with d i ∆θ * given by the formula (56), then Y * given by (36) or (37) is a stable fixed point.

Remark 2.
When the angular velocity at touch-down for the fixed point Y * equals θ * d = √ cos α, then the condition (57) from Theorem 2 takes the form where from (52) If we additionally assume that K is large enough to be replaced by K = π 2 cos α 4α 2 (see (27)), then the formula for d i ∆θ * given by (56) takes the form

Stability regions for special case
On the one hand, the special case is a mathematical simplification compared to the general one. On the other hand, it reflects a typical running speed in humans. For small angles of attack α, the horizontal Froude number X relates to the angular velocity at touch-down with X ≈ θ d . Hence the case C( ) = 0 describes running for which the horizontal velocity X is approximately 1. For a leg length of one meter, this means human jogging speed about 3.13 m/s, which exactly corresponds to the initial velocity for a minimum energy of 1.5 and an initial apex height of 1. consisting of a combination of parameters leading to a stable fixed point. Due to Remark 2 we obtain the following E s (α) -region: where E − s and E + s are determined from (58) and (60), and E min s = 1 2 cos α + cos α is given by the condition (46), and also fixed α ∈ (π/36, π/6). It is easy to verify that E − s < E min s for each parameter α within the range under consideration. Thus, the lower energy limit of the system is E min s . Moreover, using (25) and (26) we can find the following K(α, E s )region for stable solutions where K min is given by (47) and (48), K + = K(α, E + s ), and also fixed α ∈ (π/36, π/6). Based on these results (see (61) and (62)), parameter combinations leading to stable fixed points are indicated in Fig. 4 as dark areas. For example, an angle of attack α = π/9, marked in Fig. 4 with a vertical line, necessitates a minimum system energy E min s = 1.4718 (shown in Fig. 4) and the lower stability constraint corresponds to a system energy E − s = 1.411 below this minimum (not shown in Fig. 4). On the other hand, the system energy related to the upper constraint E + s = 1.528 (shown in Fig. 4) still exceeds the critical level. Moreover, on the right side of Figure 4, we can check that for α = π/9: 18.1013 = K min < K < K + = 19.5960.
Taking into account the assumption of small values of the angle α, the range of the approximate solution is always bound to a spring stiffness exceeding physiologically reasonable values. Since angles of attack less than π/18 are unlikely to happen in a real human run with jogging speed, then the stiffness K determined from the symmetrical model (see (25) and (26)) cannot be a reliable physiological value. For example, taking α = π/36 and E s = 1.5, the dimensionless stiffness K(π/36, 1.5) is 322 (see Fig. 4). Compared to experiments, where typical stiffness values are in the range of 40 -80 (e.g. [29]), the obtained K = 322 is from 4 to 8 times greater.  Figure 5: Stability of spring-mass running. The return map function Y i+1 (Y i ) is given by (32) for the parameter set: α = π/9, E s = 1.48 and K = 18.3575 with marked the stable fixed point (black circle) and the unstable fixed point (white circle).
The return map function Y i+1 (Y i ) (see (32)) in Fig. 5 is shown for the parameter set α = π/9, E s = 1.48 and K = K(π/9, 1.48) = 18.3575, which belongs to the calculated regions of parameter combinations producing stable fixed point solutions. As predicted, the return map has a stable fixed point Y * (see (36) or (37)) attracting neighboring apex states within a few steps (arrow traces in the magnified region in Fig. 5). The stable fixed point Y * is calculated from the formula (36) or (37) and for the parameter set: α = π/9, E s = 1.48, θ * d = cos (π/9), and L * d = 2 × 1.48 − 3 × cos (π/9), it is 0.9399, exactly as shown in Fig. 5. Furthermore, the return map is characterized by an additional, unstable fixed point representing the upper limit of the basin of attraction of the stable one. Here, the basin of attraction contains all apex heights from the landing height Y i = cos α to the second, unstable fixed point (white circle). Both fixed points are merely observations from plotting equation (32). However, from Fig. 5 it is obtained that, if the system energy is not adequately selected (E s > E + s ), the fixed point given by (36) or (37) is unstable.

Stability regions for generalized symmetrical cases
In this section, we will deviate from the case of θ * d = √ cos α and present the regions of stability in a general approach. The upper and lower limits of the stability regions, denoted in the figures by a solid and a dashed line, respectively, are obtained from Theorem 2. Fig. 6 shows the stability regions of angular velocity θ * d depending on the changing energy (with a given α) and the changing angle of attack (with a given E s ). For a fixed angle of attack, the angular velocity θ * d increases with the increase in energy, while for a fixed energy, the angular velocity θ * d decreases with the increase in the angle of attack. For sprints, running technique changes significantly. This is due to higher θ * d values and smaller α values. While the angles of attack are smaller than in jogging, the take-off    Figure 7: On the left, the stability region of the energy E s (α), while on the right, the stability region of the dimensionless stiffness K(α, E s ), determined by (25) and (26), with θ * d = 1.5 for varying angles of attack α from 10 • (= π/18) to 30 • (= π/6). The upper and lower limits are respectively E + s and E − s on the left, and K + = K(α, E + s ) and K − = K(α, E − s ) on the right.
angles are significantly larger. This introduces an asymmetry for which the model is not suitable, so it does not make sense to analyse the stability for this running technique. In this section physiologically valid values of leg stiffness are considered (i.e. α greater than π/18). Fig. 7 illustrates the stability regions of energy E s and stiffness K(α, E s ), defined by (25) and (26), for the changing angle of attack in the case of a fixed angular velocity. As can be seen in Fig. 7, the energy E s increases with increasing angle of attack at constant angular velocity, while the stiffness K of course decreases.
If the angular velocity θ * d is constant, the energy increases with increasing angle of attack α. During the contact phase, for a larger angle of attack α, more energy is required for the mass point to stabilize at the same level. On the other hand, for a fixed α, more energy is needed to increase the angular velocity θ * d . The stiffness of the spring depends both on the angle α and the energy E s . The greater the energy, the greater the stiffness, while the stiffness decreases with increasing angle of attack. Let us note, that the effect of  α on leg stiffness is much stronger and a typical situation is a decrease in stiffness with an increase in both the angle of attack and energy. It is easy to conclude that with greater leg stiffness (small α), the spring deflects less. In Section 3, we observed that small deflection is also influenced significantly by large angular velocity θ * d (parameter A( ) is positive). The stability of the model enforces the symmetry of solutions, which for the small α, required for approximate solutions, gives unreasonable values of the running parameters. Thus, the compromise requires to analyze the stability for the attack angles in the range from π/18 to π/6.

Transcritical bifurcation
In the current section, we will present numerical evidence for the occurrence of a transcritical bifurcation in our mapping given by equation (32). In particular, by considering E s as a bifurcation parameter, we will show that fixed point solutions collide in a transcritical bifurcation and exchange stability. Let us focus in Fig. 8, which is a one-parameter bifurcation diagram, where we depict the existence and stability of fixed point solutions versus energy. The solid line in the figure indicates stable solutions, and the dashed line unstable ones.
Let us start from the left of the figure with the lowest energy level E s = 1.51. At this value of the bifurcation parameter, there is a pair of fixed point solutions, the one on the  Figure 9: Numerical depiction of transcritical bifurcation occurring under variation of bifurcation parameter E s in mapping (32). The angle of attack α is set to α = π/12, while the stiffness K obtained from (25)  lower branch is stable and the one on the upper branch is unstable. Increasing the value of energy parameter E s , these two solutions move closer together until they collide for E + s = 1.550 with Y * = 0.9754. Increasing E s further, past the value of E s = E + s we can see that the stable branch continues past the bifurcation point as an unstable branch, and the unstable one as the stable branch. We thus have a typical scenario of a transcritical bifurcation. To see this further, in Fig. 9, we numerically depict mapping (32) for three distinct values of the bifurcation parameter. Namely for E s < E + s , then for E s = E + s and, finally for E s > E + s . Let us consider first the graph shown by black dotted curve for E s = 1.515 < E + s . We can see in the graph two fixed points with the stable (black solid dot) one below the unstable one (black empty circle), which agrees with the bifurcation diagram presented earlier in Fig. 8. Increasing E s to E s = E + s = 1.550 the two fixed points collide -see the red solid curve in the figure. Notice that the graph is quadratically tangent to line Y i+1 = Y i at the bifurcation point Y * (red solid dot). Finally, increasing E s to E s = 1.575 > E + s allows us to see the exchange of stability of the fixed points -see the black dashed curve in Fig. 9 -the unstable fixed point (black empty circle), say Y u still lies above fixed point Y * , that is Y u > Y * , and the stable one (black solid dot), say Y s , lies below Y * , that is Y s < Y * . From these three graphs we can see that the defining and nondegeneracy conditions for a transcritical bifurcation are satisfied. Namely, we observe a quadratic tangency of the mapping at fixed point Y * at E s = E + s , and unfolding with respect to parameter E s , which shows the exchange of stability of the fixed points.

Conclusions
The analysis presented in the paper shows that we can construct a useful and relatively simple model as an approximation to the more complex spring-mass description of running. In [11], the authors used a similar strategy to reduce the spring-mass model to a one-dimensional map. However, their method of constructing a one-dimensional mapping as a reduced approximate model is different than ours. That is, we are led by rigorous application of scaling and perturbation analysis conducted in [20], which allows us to neglect several small terms and remain asymptotically consistent. This rigorous approximation leads to an integrable system that can be thoroughly analysed. In particular, we have proved that under some natural conditions the apex to apex mapping has a stable fixed point which concurs with realistic situations. We conduct stability analysis of this fixed point to determine the regions of its existence in parameter space. We also determine the condition for the existence of fixed points expressed as an inequality linking the system's energy with the angle of attack and angular velocity. This result is stated in Theorem 1.
We also perform a numerical continuation of fixed points with respect to energy. This has allowed us to discover that the stable fixed point undergoes a transcritical bifurcation. A thorough investigation of this behaviour as well as the analysis of asymmetric solutions, which are of special interest in the modelling of running, is the subject of our further studies.