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

Abstract

1 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 (see Dickinson et al. (2000) and Gordon et al. (2017)). 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 (cf. Biewener & Patek (2018), Collins et al. (2005) and Holmes et al. (2006)). 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 (1973), which was based on early research of Hill (1925), provides an interesting account of this optimization. Some modern approaches to that subject can be found in Aftalion & Martinon (2019), Pritchard (1993) and Woodside (1991).

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 (1989) and thoroughly investigated by McMahon & Cheng (1990). This model can be considered conceptual in the sense that it is a simple realization 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 e.g. in Farley et al. (1991, 1993) and He et al. (1991). 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 (1986) and a modern review Aguilar et al. (2016)). Furthermore, in the literature one can find a number of important generalizations of the SLIP model such as varying attack angle (cf. Gan et al. (2018)), introducing control (cf. Sato & Buehler (2004); Takahashi et al. (2017) and Shahbazi et al. (2016)), damping (cf. Saranli et al. (2010)) and multi-legged hoppers (cf. Gan & Remy (2014)).

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, approximate analytical solutions can be obtained. This has been done in the work of Geyer et al. (2005) 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 Schwind & Koditschek (2000), Ghigliazza et al. (2005), Saranli et al. (2010) and in our earlier works (see Płociniczak & Wróblewska (2020) and Okrasińska-Płociniczak & Płociniczak (2020)). There, we have conducted a rigorous Poincaré–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 Holmes et al. (2006), 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 one-dimensional 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 Geyer et al. (2005), 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 order of the error made when using approximations (see Płociniczak & Wróblewska (2020) and Okrasińska-Płociniczak & Płociniczak (2020)). 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. To sum up, the main new results of this paper are: the construction of a reduced mapping for fixed point analysis by means of asymptotic solutions with precise error estimates on physical parameters such as stiffness and attack angle, the stability analysis of the fixed point solutions and the discovery and analysis of a transcritical bifurcation in the reduced mapping. It is important to note that this model has practical applications in the interpretation of real human running. The results presented in this article have been utilized for analysing the data from the experiment described in Wróblewska et al. (2023). In particular, in Wróblewska et al. (2023), we apply the model to experimental data in view of determining the efficiency of running by constructing model-based measure for energy levels linked with the existence of stable running gaits in humans.

The rest of the paper is outlined as follows. In Section 2 we introduce the governing equations of the spring-mass model of running. Then, in Section 3, we introduce approximate solutions of our model equations obtained by means of the Poincaré–Lindstedt method. These equations are then used in Section 4 to obtain a one-dimensional map which is also analysed in this section. Finally, Section 5 concludes the paper.

2 Spring-mass running model

In this section we will derive non-linear equations, which describe running during the stance phase (see Fig. 1 for notation). We naturally scale spring length with respect to |$1/l_0$| and choose the pendulum time scale |$\sqrt{g/l_0}$|⁠, where |$l_0$| is the initial leg length of the leg (i.e. the distance of the human centre of mass from the ground in a standing position). Using non-dimensional polar coordinates |$(L,\theta )$|⁠, the Lagrange function of the stance phase is given by

Fig. 1.

Schematic of the spring-mass model for running. Dimensionless parameters: |$(X,Y)$|—Cartesian coordinates of the point mass, |$(L,\theta )$|—polar coordinates, where |$L=\sqrt{X^2+Y^2}$| is the radial, while |$\theta $| is the angular position of the point of mass, and additionally |$\varDelta \theta $| is angle swept during stance. Moreover, we assume that during touch-down and take-off, denoted by TD and TO, respectively, the value of |$\theta $| is |$\theta _{{\scriptstyle{\text{TD}}}}=-\alpha $| and |$\theta _{{\scriptstyle{\text{TO}}}}$|⁠, where |$\alpha \in (0,\pi /2)$|⁠.

Open in new tabDownload slide

where the parameter |$K$| denotes the non-dimensional spring stiffness. So the two Euler–Lagrange equations are as follows:

where |$P=L^2\theta ^{\prime}$| denotes the non-dimensional angular momentum. Moreover, the prime sign (⁠|${}^{\prime} $|⁠) represents the total derivative with respect to the dimensionless time of the pendulum.

Observe that the dimensionless form of the mechanical energy during the contact phase denoted by |$E_s$| (see (2.1)) is given by

In what follows, we consider small angle approximation in the equations of motion (2.2), which translates on the total energy (2.4). That is, the total energy |$E_s$| given by equation (2.3), for small angles |$\theta $|⁠, takes an alternative form

We should add here that the small angle approximation in the equations of motion (2.2) is not a linearization, but it corresponds to the vector field in which terms that are small are ignored (where the meaning of small should be understood as it is used in asymptotics). The leading order term of |$\sin \theta $|⁠, for small |$\theta $|⁠, is |$0$|⁠, but the first non-zero term is |$\theta $|⁠, which is small, and hence it is neglected in the vector field equations. Strictly speaking, as it is shown in Płociniczak & Wróblewska (2020), the solutions which we use here are expansions for |$\theta (t)$| and |$L(t)$| in terms of the inverse of the square root of |$K$|⁠, which will be defined as the small parameter, we replace |$\cos \theta $| with |$\cos \alpha $| and |$\sin \theta $| with |$\sin \alpha $|⁠, for some fixed |$\alpha $|⁠, in the vector field equations to get the accuracy up to and including |$O(\epsilon ^2)$| in our solutions, where |$\epsilon =1/\sqrt{K}$|⁠. Since we are interested in small |$\alpha $| these two approximations may be understood as considering the system vector field where |$\cos \theta = 1$| and |$\sin \theta = 0$|⁠. However, if one considers |$\alpha $| as not being small, but such that |$K$| remains the large parameter, the asymptotic expansions still hold up to and including the order |$O(\epsilon ^2)$|⁠. In this latter case, there is a constant change of the angular momentum, but the constraints on the boundary values for symmetric solutions, namely equations (2.6), (2.8) and (2.9) still hold.

Consequently, the Lagrangian (2.1) is also modified, and hence for small angle approximation, the change of the angular momentum of the point mass we consider is zero, i.e.

From equations (2.4) and (2.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 (see Fig. 1), we have the radial symmetry with rest length |$L=1$| at each phase transition

and so from (2.4) we have then

The subscript TO/TD means that the equations are valid in both TD and TO time positions in stance. Furthermore, after integrating (2.5), we get |$\; \theta ^{\prime}_{{\scriptstyle{\text{TO/TD}}}}=cL^{-2}_{{\scriptstyle{\text{TO/TD}}}}=c$|⁠, where |$c$| is some real constant. In our case: forward locomotion, the constant |$c$| is positive. Thus

and |$\;2\left (\widetilde{E}_s-1\right ) - c^2= \left (L^{\prime}\right )^2_{{\scriptstyle{\text{TO/TD}}}},\;$| and hence |$|L^{\prime}_{{\scriptstyle{\text{TO}}}}|=|L^{\prime}_{{\scriptstyle{\text{TD}}}}|$|⁠. We can now write

In addition to (2.6), (2.8) and (2.9), we require (see also Fig. 1)

where only the angle |$\varDelta \theta $| swept in stance cannot simply be expressed by the state at the time of touch-down. It should be observed that for small angles |$\theta $|⁠, the conservation of energy and angular momentum do not enforce symmetry in the stance phase, i.e. through the arbitrary value of |$\theta $| state at take off, asymmetric solutions are also allowed. In other words, solutions for which |$|\theta _{TD}|\not = |\theta _{TO}|$|⁠.

In the next section we will derive the formula for |$\varDelta \theta $|⁠.

3 Constructing approximate solutions

An asymptotic analysis of the main equations (2.2) with the use of the Poincaré–Lindstedt series was carried out in Płociniczak & Wróblewska (2020). This solution (see Płociniczak & Wróblewska (2020)) is based on the perturbative expansion related to the significant spring stiffness (⁠|$K \rightarrow \infty $|⁠). To simplify matters we set |$\epsilon =1/\sqrt{K}$|⁠, (⁠|$\epsilon \rightarrow 0^{+}$|⁠) and introduce

In what follows we will operate of the order of |$O(\epsilon ^2)$| as |$\epsilon \rightarrow 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.

Fig. 2.

|$L$| and |$\theta $| approximations during the stance phase. Here |$\alpha =0.1$|⁠, |$X^{\prime}=1$|⁠, |$Y^{\prime}=0.1$| and |$K=15$|⁠.

Open in new tabDownload slide

• Radial motion during stance |$L(t)$| (see (23) in Płociniczak & Wróblewska (2020)): The |$L(t)$| approximation is as follows:

  $$ \begin{align}& L(t)= 1-\epsilon L_d\sin{\left(\epsilon^{-1} \widetilde\omega(\epsilon)t\right)} + \epsilon^2\left(\theta_d^2-\cos{\alpha}\right)\left(1-\cos{\left(\epsilon^{-1} \widetilde\omega(\epsilon)t\right)}\right), \end{align} $$

  (3.2)

  where the initial conditions have the form for |$\alpha>0$|  

  $$ \begin{align}& \begin{array}{ll} \theta(0)=\theta_{{\scriptstyle{\text{TD}}}}=-\alpha, & L(0)=L_{{\scriptstyle{\text{TD}}}}=1, \\ \theta^{\prime}(0)=\theta^{\prime}_{{\scriptstyle{\text{TD}}}}=\theta_d, & L^{\prime}(0)=L^{\prime}_{{\scriptstyle{\text{TD}}}}=-L_d, \end{array} \end{align} $$

  (3.3)

  with

  $$ \begin{align}& \theta_d = X^{\prime} \cos\alpha + Y^{\prime} \sin\alpha, \quad L_d = X^{\prime} \sin\alpha - Y^{\prime} \cos\alpha. \end{align} $$

  (3.4)

   |$X^{\prime}$| and |$Y^{\prime}$| are the horizontal and vertical velocities made dimensionless by the factor (acceleration of |$\text{gravity} \times \text{length})^{1/2}$|⁠, which are called Froude numbers in fluid mechanics. In addition, it is worth noting that from (2.3)

  $$ \begin{align}& L_d=|L^{\prime}_{{\scriptstyle{\text{TD}}}}|=\sqrt{2E_s-\theta_d^2-2\cos{\alpha}}, \end{align} $$

  (3.5)

  where |$\theta _d$| refers to the angular momentum |$P$| at touch-down (TD).
  The general solution (3.2) of the radial motion |$L(t)$| during stance describes a sinusoidal oscillation around |$L=1+A(\epsilon )$| with amplitude |$B(\epsilon )$| and frequency |$\epsilon ^{-1} \widetilde \omega (\epsilon )$|⁠, where

  $$ \begin{align}& A(\epsilon)=\epsilon^2\left(\theta_d^2-\cos{\alpha}\right) \quad \text{and} \quad B(\epsilon)=\epsilon\sqrt{\epsilon^{2}(\theta_d^2-\cos{\alpha})^2 + L_d^2}. \end{align} $$

  (3.6)

  Thus, using the notation (3.6), the approximation (3.2) can be represented as follows:

  $$ \begin{align}& L(t)= 1+A(\epsilon) - \epsilon L_d\sin{\left(\epsilon^{-1} \widetilde\omega(\epsilon)t\right)}-A(\epsilon)\cos{\left(\epsilon^{-1} \widetilde\omega(\epsilon)t\right)}. \end{align} $$

  (3.7)

  But, as shown in Fig. 2 on the left, the solution only holds for |$L(t)\leq 1$|⁠.
  In formula (3.2), the parameter |$t$| ranges from 0 to the contact time |$t_C$| given by

  $$ \begin{align}& t_{C}=\frac{2\epsilon}{\widetilde\omega(\epsilon)}\arccos{\frac{A(\epsilon)}{B(\epsilon)}}=\frac{4\epsilon}{2-\epsilon^2\theta_{d}^2}\arccos{C(\epsilon)}, \end{align} $$

  (3.8)

  where

  $$ \begin{align}& C(\epsilon)=\frac{A(\epsilon)}{B(\epsilon)}=\frac{\epsilon\left(\theta_d^2-\cos{\alpha}\right)}{\sqrt{\epsilon^{2}(\theta_d^2-\cos{\alpha})^2 + L_d^2}}, \end{align} $$

  (3.9)

  and |$\epsilon \rightarrow 0^{+}$|⁠. Observe that |$L(0)=L(t_C)=1$| and |$L(t_C/2)=1-\varDelta L_{MAX}$| (see Fig. 2). The maximum leg deflection during the stance phase, denoted by |$\varDelta L_{MAX}$|⁠, is given by the difference of the amplitude |$B(\epsilon )$| of the radial motion |$L(t)$| and the shift |$A(\epsilon )$| of the touch-down position, i.e.

  $$ \begin{align}& \varDelta L_{MAX}= B(\epsilon) - A(\epsilon). \end{align} $$

  (3.10)

  So restriction to small values of |$L-1$| is adequately formulated by |$B(\epsilon )-A(\epsilon ) = O(\epsilon ) \ll 1$| as |$\epsilon \rightarrow 0^{+}$|⁠. For illustrations, we refer to Fig. 2. In the case of the data from Fig. 2, we get |$A(\epsilon ) = -0.0016$| and |$B(\epsilon ) = 0.0515$|⁠, where |$\epsilon =1/\sqrt{15}$|⁠. Hence |$\varDelta L_{MAX}=0.0531$|⁠, and we can check, that |$\varDelta L_{MAX}=1-L(t_C/2)$|⁠, where |$t_C=0.8558$|⁠. In addition, note that the shift |$A(\epsilon )$| may have also a positive result, when |$\theta _d>\sqrt{\cos{\alpha }}$|⁠, then it will move |$1$| below |$1+A(\epsilon )$| and reduce the maximum spring compression |$\varDelta L_{MAX}$|⁠.
• Angular motion during stance |$\theta (t)$| (see (24) in Płociniczak & Wróblewska (2020)): The |$\theta (t)$| approximation with the initial conditions (3.3) is as follows:

  $$ \begin{align}& \theta(t)= - \alpha +\theta_d\widetilde\omega(\epsilon)t-\frac{1}{2}\left(\widetilde\omega(\epsilon)t\right)^2\sin{\alpha} +2\epsilon^2L_d\theta_d\left(1-\cos{\left(\epsilon^{-1} \widetilde\omega(\epsilon) t\right)}\right). \end{align} $$

  (3.11)

  Observe that |$\theta (0)=\theta _{{\scriptstyle{\text{TD}}}}=- \alpha $| and |$\theta (t_C)=\theta _{{\scriptstyle{\text{TO}}}}=\varDelta \theta +\theta _{{\scriptstyle{\text{TD}}}} = \varDelta \theta - \alpha $| (see Fig. 2 on the right).
• Angle swept during stance: From formulas (3.11) and (3.8) we obtain the angle |$\varDelta \theta $| swept during stance

  $$ \begin{align}& \varDelta \theta= \theta(t_C) +\alpha =2\epsilon\theta_d \times\left[\arccos{C(\epsilon)} + 2\epsilon L_d\left(1-C^2(\epsilon)\right)\right] - 2\epsilon^2\sin{\alpha}\left(\arccos{C(\epsilon)}\right)^2, \end{align} $$

  (3.12)

  where |$C(\epsilon )$| is already defined in (3.9). Therefore, the angle swept |$\varDelta \theta $| is uniquely defined by the system state at touch-down |$(L^{\prime}_{{\scriptstyle{\text{TD}}}}, \theta _{{\scriptstyle{\text{TD}}}}, \theta ^{\prime}_{{\scriptstyle{\text{TD}}}})$| and the spring stiffness parameter |$K$|⁠. Moreover, the landing angle |$\theta _{{\scriptstyle{\text{TD}}}}=- \alpha $| affects |$\varDelta \theta $| also by determining the radial and angular landing velocity.For the data in Fig. 2 we get |$\theta (t_C)= \theta _{{\scriptstyle{\text{TO}}}} = 0.7334$|⁠. Hence |$\varDelta \theta =\theta _{{\scriptstyle{\text{TO}}}} - \theta _{{\scriptstyle{\text{TD}}}}=0.8334$|⁠, which agrees with the value calculated in formula (3.12).Finally, let us observe that the Taylor expansion of (3.12) to the second order, about |$\epsilon =0$|⁠, is given by

  $$ \begin{align}& \varDelta \theta= \pi\theta_d \epsilon + \left[4\theta_d L_d - 2\frac{\theta_d}{L_d}\left(\theta_d^2-\cos{\alpha}\right) - \frac{1}{2}\pi^2\sin{\alpha}\right]\epsilon^2, \end{align} $$

  (3.13)

  where |$\epsilon \rightarrow 0$|⁠. Since |$\varDelta \theta $| approximation given by formula (3.12) is of the order |$\epsilon ^2$|⁠, all terms of the order |$O(\epsilon ^3)$| in formula (3.13) are omitted.

An important case of stability analysis, considered also by Geyer et al. (2005), is

In the special case we get the shift |$A(\epsilon )=0$| and the leg deflection parabola is equal to the amplitude, i.e. |$\varDelta L_{MAX} =B(\epsilon )$|⁠. Moreover, the angular velocity at touch-down, i.e. |$\theta _d=\sqrt{\cos{\alpha }}$|⁠, is identical to the frequency of the pendulum, and the contact time from the formula (3.8) is equal to |$2\pi \epsilon /(2-\epsilon ^2\theta _{d}^2)$|⁠. Putting |$C(\epsilon ) = 0$| into the equation (3.12) or equivalently |$\theta _d=\sqrt{\cos{\alpha }}$| into the equation (3.13), where |$\epsilon =1/\sqrt{K}$|⁠, we get (refer to equation (3.5) as well)

Parameter combinations |$\{\alpha , E_s(\alpha ), K(\alpha ,E_s)\}$| leading to a symmetrical stance phase, i.e. |$\varDelta \theta =2\alpha \;$|⁠, will be given, for the case of |$C(\epsilon ) = 0$| in Section 4.4.

All the above approximations work well for large values of the stiffness |$K$| and small values of the attack angle |$\alpha $|⁠. In Table 1 we present typical values of the human run parameters. The table was prepared based on the experiment data from Wróblewska (2020) as well as data from previous papers: McMahon & Cheng (1990), Arampatzis et al. (1999) and Geyer et al. (2005). Notice that usually |$X^{\prime}$| is of the order of unity, while |$Y^{\prime}$| and |$\alpha $| are small. Proper running is combined with touching down under the centre of mass, which guarantees small values of |$\alpha $| (cf. Romanov (2014)). On the other hand, very small |$\alpha $| angles imply an asymmetric run. Additionally, typical values of dimensionless leg stiffness |$K$| range from around 12 to around 70 (cf. Arampatzis et al. (1999)). Thus, |$\epsilon $| in a real run does not exceed 0.289. The last parameter appearing in Table 1 is the maximum leg deflection (for comparison see (3.10)). As can be seen, the parameter |$\varDelta L_{max}$| reaches the smallest of the values presented in Table 1. Finally, if we calculate the mechanical energy from the formula |$E_s=(\theta _d^2 + L_d^2)/2 + \cos \theta $| (see (2.3)) for the values of |$\alpha $|⁠, |$X^{\prime}$| and |$Y^{\prime}$| from Table 1, it is from 1.2 to 3.6.

By setting the stance phase limits to symmetry, we derive in the next section approximations of the dimensionless leg stiffness |$K$|⁠.

3.1 Symmetric solutions

By limiting ourselves to the second-order expansion (3.13), the solution |$\varDelta \theta = 2\alpha $| comes down to the quadratic equation for |$\sqrt{K}$|⁠. A physically reasonable solution gives the following expression to approximate the spring stiffness:

where

Using |$\alpha \rightarrow 0^+$| in (3.17), the second component tends to zero and we get ¹

The same approximation of |$K$| as (3.18) was obtained in Płociniczak & Wróblewska (2020) and Wróble–wska [2020]. Moreover, in Okrasińska-Płociniczak & Płociniczak (2020) it was shown that the leading order of the expansion of |$K$| as |$\alpha \rightarrow 0^+$| is indeed |$(\pi X^{\prime})^2/(2\alpha )^2$|⁠.

We highlight here that the stiffness |$K$|⁠, given by equation (3.16), is a function of the state at touch down. That is, |$K = K(\alpha ,\theta _d,L_d)$|⁠. In Section 4.4, we consider a special case of fixed point solutions when |$\theta _d = \sqrt{\cos \alpha }$|⁠. In this case, if the asymptotic approximation of |$K$|⁠, i.e. |$\tilde{K}$|⁠, is being taken into consideration, we get

Note that in this case, |$\tilde{K}$| is a function of touch down angle |$\alpha $| and may be regarded as a constant since |$\alpha $| is a fixed parameter.

In the symmetric case we have a particular boundary value problem to solve. Let |$(\theta (t,K), L(t,K))$| be the solution of the system (2.2) with (3.3). Find |$K^*$| and the smallest time |$t^*>0$| satisfying

The above problem can be easily solved numerically using the shooting method, as was done e.g. in McMahon & Cheng (1990). Graphical illustrations confirming the validity (3.18) as an approximation of the solution to the boundary problem (3.19) are presented in Płociniczak & Wróblewska (2020) and Wróblewska (2020).

In Fig. 3, on the left, we can see what the approximations |$K$| given by (3.16) and |$\widetilde{K}$| given by (3.18) look like. Obviously, if |$\alpha \rightarrow 0^+$| then |$K^*\rightarrow \infty $| and both approximations of |$K^*$| are very good. However, for moderate values of parameter |$\alpha $|⁠, |$K$| is slightly better than |$\widetilde{K}$|⁠.

Fig. 3.

On the left: comparison between the numerically calculated value of |$K^*$| and its |$K$| approximation, calculated from (3.16), and the second |$\widetilde{K}$| approximation, calculated from (3.18), for varying |$\alpha $| from |$6^{\circ }(=\pi /30)$| to |$20^{\circ }(=\pi /9)$|⁠. On the right: ratio of |$K^*$| to |$K$| and |$\widetilde{K}$| for varying |$\alpha $| from |$1^{\circ }(=\pi /180)$| to |$14^{\circ }(=7\pi /90)$|⁠. Here |$X^{\prime}=1$| and |$Y^{\prime}=0.1$|⁠.

Open in new tabDownload slide

Moreover the right side of Fig. 3 illustrates the ratio of |$K^*$| to both approximations |$K$| and |$\widetilde{K}$| for varying |$\alpha $|⁠. Indeed, the approximation |$K$| behaves better than |$\widetilde{K}$| in the presented range of the parameter |$\alpha $|⁠, i.e. |$1^{\circ } \leq \alpha \leq 14^{\circ }$| degree. It turns out that the approximation |$K$| works very good for |$\alpha <3^{\circ }$|⁠, while |$\widetilde{K}$| is accurate only for |$\alpha <1^{\circ }$|⁠. In a real run, the typical angle of attack |$\alpha $| is about |$14^{\circ }$| (see Table 1), and the relative error for the approximation |$K$| (⁠|$(K^*-K)/K\cdot 100\%$|⁠) is 6%, while for |$\widetilde{K}$| (⁠|$(K^*-\widetilde{K})/\widetilde{K}\cdot 100\%$|⁠) it is 12% (see the right side of Fig. 3). Therefore, approximations could be applied for most of the values appearing from the real run, taking into account the size of the errors 6% and 12%. Figure 3 also shows that the relative error decreases rapidly as |$\alpha $| gets closer to |$0^+ $|⁠.

4 Stability of spring-mass running

4.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 Section 2 we showed that assuming small angles, the system state at take-off (TO) is related to the state at touch-down (TD) through a set of the equations (see (2.6), (2.8), (2.9) and (2.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 (2.3) by |$E_s = (\theta ^2_{d} + L^2_{d})/2 +\cos{\alpha }$|⁠. Since the Froude numbers |$X^{\prime}$| and |$Y^{\prime}$| depend on |$Y_i$|⁠, then the radial velocity |$L_d$| and the angular velocity |$\theta _d$| also depend on |$Y_i$| (see (3.4) or (4.2)). Then, an additional mapping is required between the take-off state and the height |$Y_{i + 1}$| of the apex |$i + 1$|⁠, which is

Using both mappings (4.2) and (4.3), the apex return map function |$Y_{i+1}(Y_i)$| can be constructed and takes the following form after some simplifications:

where |$\varDelta \theta _i = \varDelta \theta (\theta _d,L_d)$| is given by (3.12) or (3.13). To mark the periodicity for this running model, we introduced the index |$i$| to the notation |$\varDelta \theta _i$|⁠. So |$\varDelta \theta _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 |$\varDelta \theta _i= \varDelta \theta (\theta _d,L_d)$| depends on |$Y_i$|⁠. From formulas (3.4) and (4.2) we have

Moreover, observe that (see (3.5))

If the apex height |$Y_{i+1}$| is less then landing height (⁠|$Y_{i+1} < \cos{\alpha }$|⁠), 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 (4.1) the kinetic energy |$\frac{1}{2} ( (L^{\prime})^2 + L^2(\theta ^{\prime})^2 )$| is equal at touch-down and take-off. For asymmetric contact phases, since |$Y_{{\scriptstyle{\text{TD}}}} \neq Y_{{\scriptstyle{\text{TO}}}}$|⁠, there is a net change in system energy |$\varDelta E=Y_{{\scriptstyle{\text{TO}}}}-Y_{{\scriptstyle{\text{TD}}}}$|⁠. On the other hand, during symmetric stance phases (⁠|$Y_{{\scriptstyle{\text{TD}}}}=Y_{{\scriptstyle{\text{TO}}}}$|⁠) 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 (4.3)

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 height |$Y_{i+1}$| from (4.4) the adjusted |$E_s$| is used. We note here that the energy compensation expressed in equation (4.7) is non-unique. That is, a different manner of such a compensation may be used.

4.2 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 (4.3) with condition (4.1) for |$\theta _ {TO} = \alpha $| leads to the corresponding apex height of fixed points

where |$\theta ^*_d$| and |$L^*_d$| represent the values of |$\theta _d$| and |$L_d$| given by (4.5) at the fixed point |$Y^*$|⁠. Since |$\theta ^*_d$| and |$L^*_d$| are linked by the relationship (4.6), the apex height of fixed points can be expressed in the following equivalent form:

Substituting (4.8) and (4.9) into equations (4.5) yields

and also

It is easy to see that

So, it follows that equations (4.10) and (4.11) hold if

From (3.5) and (4.12) we have

Hence the system energy |$E_s$| must satisfy the inequality |$E_s \ge E_s^{min}$| with

For the system energy |$E_s = E_s^{min}$|⁠, the apex height of fixed point |$Y^*$| is |$\cos{\alpha }$|⁠, which is equal to the landing height. As the energy increases (⁠|$E_s> E_s^{min}$|⁠), the apex height increases, but it never exceeds the upper boundary, i.e. |$E_s$|⁠.

If the attack angle |$\alpha $| is fixed, then the stiffness approximation given by (3.16) and (3.17) is the lowest for the minimum energy system given by (4.14), i.e.

where

Whence, we proved that there is a fixed point under a certain assumption:

Considering |$\theta _d^*=\sqrt{\cos \alpha }$| case, if |$\alpha \rightarrow 0$|⁠, the minimal system energy |$E_s^{min}$| 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 \approx 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^{\prime}=X^{\prime}\sqrt{gl_0}=\sqrt{2(E_s-Y_0)gl_0}$|⁠, which is |$3.13 \: m/s$|⁠.

4.3 Stability of fixed points

Let us denote mapping (4.4) as |$Y_{i+1}=f(Y_i)$|⁠. Stable solutions of fixed point |$Y^*$| fulfil the following condition:

So to prove stability, at least one the parameter set |$(\alpha , E_s, \theta ^*_d)$| leading to solutions of |$Y^*$| satisfying (4.21) should be identified. Starting with (4.4), we get

using |$\varDelta \theta _{i} = 2\alpha $| for symmetric contact phases. Later in the paper, we take the notation |$\frac{d}{d\:Y_i}\varDelta \theta _i\Big \rvert _{Y^*}=d_i\varDelta \theta ^*$|⁠. Since the expression in square brackets in (4.22) always remains positive, then (4.21) 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 |$\varDelta \theta _i$|⁠. However, for stability reasons, the positive derivative |$d_i\:\varDelta \theta ^*$| cannot be greater than |$2\left [\sin{\alpha } + 2\sqrt{(Y^* - \cos{\alpha })(E_s - Y^*)}\right ]^{-1}$|⁠. Then we substitute the equations fixed point |$Y^*$| (4.8) and (4.9) satisfying the condition (4.17) and denoting |$2\sqrt{(Y^* - \cos{\alpha })(E_s - Y^*)}$| by |$D\left (\alpha ,E_s,\theta _d^*\right )$| we get

Next, from (3.13) it follows that

where |$d_i K^*=\frac{d}{d\:Y_i}K\Big \rvert _{Y^*}$| and |$K$| is given by (3.16). Moreover we get

and also

Resolving |$d_i \theta _d^*$| (see (4.5))

and applying (4.6) we have

and from (4.29) can be further deduced

and finally from (4.28)

with |$D\left (\alpha ,E_s,\theta _d^*\right )$| given by (4.24).

4.4 Stability regions for special case |$\theta _d^*=\sqrt{\cos \alpha }$|

On the one hand, the special case is a mathematical simplification compared with the general one. On the other hand, it reflects a typical running speed in humans. For small angles of attack |$\alpha $|⁠, the horizontal Froude number |$X^{\prime}$| relates to the angular velocity at touch-down with |$X^{\prime} \approx \theta _d$|⁠. Hence the case |$C(\epsilon ) = 0$| describes running for which the horizontal velocity |$X^{\prime}$| is |$\sim $|1. For a leg length of 1 m, 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.

In situations where the leg stiffness can be approximated by |$\widetilde{K}$| given in (3.18) in our special case of |$\theta _d^* = \sqrt{\cos{\alpha }}$| we essentially have |$\widetilde{K}=\widetilde{K}(\alpha )$|⁠. Therefore, the stiffness is constant with respect to |$Y_i$|⁠. This has an impact on the condition on the stability region. Hence, Theorem 4.2 has the following form:

Two curves: lower, when |$f^{\prime}(Y^*)=-1\;$| i.e. |$\;d_i\varDelta \theta ^* =2[\sin{\alpha } + D(\alpha ,E_s)]^{-1}$|⁠, and upper, when |$f^{\prime}(Y^*)=1\;$| i.e. |$\;d_i\varDelta \theta ^* =0$| (see 4.21) are the region constraints: |$E_s(\alpha )$| and |$K(\alpha ,E_s)$|⁠, consisting of a combination of parameters leading to a stable fixed point. When we approximate the stiffness with |$\widetilde{K}$|⁠, it is possible to find the |$E_s(\alpha )$| region analytically. That is, due to Remark 4.2 we have

where |$E_s^-$| and |$E_s^+$| are determined from (4.33) and (4.36), and |$E_s^{min}=\frac{1}{2\cos{\alpha }}+\cos{\alpha }$| is given by the condition (4.18), and also fixed |$\alpha \in \left (\pi /36,\pi /6\right )$|⁠. It is easy to verify that |$E_s^-<E_s^{min}$| for each parameter |$\alpha $| within the range under consideration. Thus, the lower energy limit of the system is |$E_s^{min}$|⁠. Moreover, using (3.16) and (3.17) we can find the following |$K(\alpha , E_s)$|—region for stable solutions:

where |$K^{min}$| is given by (4.19) and (4.20), |$K^+=K(\alpha ,E_s^+)$|⁠, and also fixed |$\alpha \in \left (\pi /36,\pi /6\right )$|⁠.

Based on these results (see (4.37) and (4.38)), parameter combinations leading to stable fixed points are indicated in Fig. 4 as dark areas. For example, an angle of attack |$\alpha =\pi /9$|⁠, marked in Fig. 4 with a vertical line, necessitates a minimum system energy |$E_s^{min}=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 Fig. 4, we can check that for |$\alpha =\pi /9$|⁠:

Fig. 4.

Parameter interdependence for fixed point solutions with |$\theta _d^*=\sqrt{\cos \alpha }$|⁠. The stability regions: |$E_s(\alpha )$| (on the left) and |$K(\alpha )$| (on the right) for spring-mass running predicted by the analytical approximations (4.37) and (4.38) for varying angles of attack |$\alpha $| from |$5^{\circ }(=\pi /36)$| to |$30^{\circ }(=\pi /6)$|⁠. The region of stable fixed points, with eigenvalues |$-1<f^{\prime}(Y^*) < 1$|⁠, is marked in grey. The region of unstable fixed points with eigenvalues |$f^{\prime}(Y^*)>1$| is marked in blue. |$f^{\prime}$| is the eigenvalue notation used in the legend.

Open in new tabDownload slide

Taking into account the assumption of small values of the angle |$\alpha $|⁠, the range of the approximate solution is always bound to a spring stiffness exceeding physiologically reasonable values. Since angles of attack less than |$\pi /18$| are unlikely to happen in a real human run with jogging speed, then the stiffness |$K$| determined from the symmetrical model (see (3.16) and (3.17)) cannot be a reliable physiological value. For example, taking |$\alpha =\pi /36$| and |$E_s=1.5$|⁠, the dimensionless stiffness |$K(\pi /36,1.5)$| is 322 (see Fig. 4). Compared with experiments, where typical stiffness values are in the range of 30–70 (see Wróblewska (2020) and Table 1), the obtained |$K = 322$| is from 4.5 to 10.5 times greater.

The return map function |$Y_{i+1}(Y_i)$| (see (4.4)) in Fig. 5 is shown for the parameter set |$\alpha =\pi /9$|⁠, |$E_s=1.48$| and |$K=K(\pi /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 (4.8) or (4.9)) attracting neighbouring 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 (4.8) or (4.9) and for the parameter set: |$\alpha =\pi /9$|⁠, |$E_s=1.48$|⁠, |$\theta _d^*=\sqrt{\cos{(\pi /9)}}$| and |$L_d^*=\sqrt{2\times 1.48-3\times \cos{(\pi /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{\alpha }$| to the second, unstable fixed point (white circle). Both fixed points are merely observations from plotting equation (4.4). 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 (4.8) or (4.9) is unstable.

Fig. 5.

Stability of spring-mass running. The return map function |$Y_{i+1}(Y_i)$| is given by (4.4) for the parameter set: |$\alpha =\pi /9$|⁠, |$E_s=1.48$| and |$K=18.3575$| with marked the stable fixed point (black circle) and the unstable fixed point (white circle).

Open in new tabDownload slide

4.5 Stability regions for generalized symmetrical cases

In this section, we will deviate from the case of |$\theta _d^*=\sqrt{\cos \alpha }$| 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 4.2. Figure 6 shows the stability regions of angular velocity |$\theta _d^*$| depending on the changing energy (with a given |$\alpha $|⁠) and the changing angle of attack (with a given |$E_s$|⁠). For a fixed angle of attack, the angular velocity |$\theta _d^*$| increases with the increase in energy, while for fixed energy the angular velocity |$\theta _d^*$| decreases with the increase in the angle of attack.

Fig. 6.

On the left, the stability region of the angular velocity |$\theta _d^*(E_s)$| with |$\alpha =\pi /9$| for varying energy |$E_s$| from |$1.45$| to |$3$|⁠. The horizontal line shows the case of |$\theta _d^*=\sqrt{\cos \alpha }=0.9694$|⁠. On the right, the stability region of the angular velocity |$\theta _d^*(\alpha )$| with |$E_s=1.8$| for varying angles of attack |$\alpha $| from |$10^{\circ }(=\pi /18)$| to |$30^{\circ }(=\pi /6)$|⁠. The upper and lower limits are |$\left (\theta _d^*\right )^+$| and |$\left (\theta _d^*\right )^-$|⁠, respectively. The region of stable fixed points, with eigenvalues |$-1<f^{\prime}(Y^*) < 1$|⁠, is marked in grey. The region of unstable fixed points with eigenvalues |$f ^{\prime}(Y^*)>1$| is marked in blue. |$f^{\prime}$| is the eigenvalue notation used in the legend.

Open in new tabDownload slide

For sprints, running technique changes significantly (see Bushnell (2004); Cunningham et al. (2013)). This is due to higher |$\theta _d^*$| values and smaller |$\alpha $| values. While the angles of attack are smaller than in jogging, the take-off angles are significantly larger (see Chan-Roper et al. (2012); Balandin et al. (2022)). 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. Moreover, in this section physiologically valid values of leg stiffness are considered (i.e. |$\alpha $| greater than |$\pi /18$|⁠).

Figure 7 illustrates the stability regions of energy |$E_s$| and stiffness |$K(\alpha ,E_s)$|⁠, defined by (3.16) and (3.17), 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 (see also Table 2).

Fig. 7.

On the left, the stability region of the energy |$E_s(\alpha )$|⁠, while on the right the stability region of the dimensionless stiffness |$K(\alpha , E_s)$|⁠, determined by (3.16) and (3.17), with |$\theta _d^*=1.5$| for varying angles of attack |$\alpha $| from |$10^{\circ }(=\pi /18)$| to |$30^{\circ }(=\pi /6)$|⁠. The upper and lower limits are, respectively, |$E_s^+$| and |$E_s^-$| on the left, and |$K^+=K(\alpha , E_s^+)$| and |$K^-=K(\alpha , E_s^-)$| on the right. The region of stable fixed points, with eigenvalues |$-1<f^{\prime}(Y^*) < 1$|⁠, is marked in grey. The region of unstable fixed points with eigenvalues |$f ^{\prime}(Y^*)>1$| is marked in blue. |$f^{\prime}$| is the eigenvalue notation used in the legend.

Open in new tabDownload slide

If the angular velocity |$\theta _d^*$| is constant, the energy increases with increasing angle of attack |$\alpha $|⁠. During the contact phase, for a larger angle of attack |$\alpha $|⁠, more energy is required for the mass point to stabilize at the same level. On the other hand, for a fixed |$\alpha $|⁠, more energy is needed to increase the angular velocity |$\theta _d^*$|⁠. The stiffness of the spring depends both on the angle |$\alpha $| 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 |$\alpha $| 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 |$\alpha $|⁠), the spring deflects less (see Table 2). In Section 3, we observed that small deflection is also influenced significantly by large angular velocity |$\theta _d^*$| (parameter |$A(\epsilon )$| is positive).

The stability of the model enforces the symmetry of solutions, which for the small |$\alpha $|⁠, required for approximate solutions, gives unreasonable values of the running parameters. Thus, the compromise requires to analyse the stability for the attack angles in the range from |$\pi /18$| to |$\pi /6$|⁠.

4.6 Transcritical bifurcation

In this section, we will present numerical evidence for the occurrence of a transcritical bifurcation in our mapping given by equation (4.4). It is also known that bifurcations occur in other approximate solutions (see Geyer et al. (2005)), but their occurrence in the full model has not yet been verified and is a subject for further research. 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 on 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.

Fig. 8.

One-parameter bifurcation diagram of (4.4). Note transcritical bifurcation of fixed point |$Y^* = 0.9754$| at |$E_s^+=1.550$|⁠. The angle of attack |$\alpha $| is set to |$\alpha =\pi /12$|⁠, while the stiffness |$K$| obtained from (3.16) increases from |$34.77$| for |$E_s = 1.51$| to |$38.73$| for |$E_s = 1.58$|⁠.

Open in new tabDownload slide

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 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 (4.4) 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 blue 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 green dashed curve in Fig. 9—the unstable fixed point (black empty circle), say |$Y_u$| still lies above fixed point |$Y^*$|⁠, i.e. |$Y_u> Y^*$|⁠, and the stable one (black solid dot), say |$Y_s$|⁠, lies below |$Y^*$|⁠, i.e. |$Y_s < Y^*$|⁠. From these three graphs we can see that the defining and non-degeneracy 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.

Fig. 9.

Numerical depiction of transcritical bifurcation occurring under variation of bifurcation parameter |$E_s$| in mapping (4.4). The angle of attack |$\alpha $| is set to |$\alpha =\pi /12$|⁠, while the stiffness |$K$| obtained from (3.16) is |$35.07$| for |$E_s = 1.515$|⁠, |$36.75$| for |$E_s=E_s^+ = 1.550$| and |$37.94$| for |$E_s=1.575$|⁠.

Open in new tabDownload slide

The physical consequence of the transcritical bifurcation in running activity is to lower the height of the stable point after exceeding the critical energy |$E_s^+$|⁠. There is more energy in the system than is needed for a given speed and must be lost on hitting the ground by changing velocity vectors. Such running becomes less economical in terms of energy consumption and is associated with a suboptimal running technique. In the considered example, the transcritical bifurcation occurs for the following set of parameters: |$\alpha =\pi /12$|⁠, |$E_s=1.550$|⁠, |$K=36.75$|⁠, |$\theta _d^*=1$|⁠, |$L_d^*=0.410$|⁠, so the horizontal and vertical Froude number is |$X^{\prime}=1.072$| and |$Y^{\prime}= 0.137$|⁠, respectively. These are reasonable values for a run (see Table 1) that is only slightly faster than |$\theta _d^*=\sqrt{\cos \alpha }$|⁠.

5 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 Geyer et al. (2005), 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 Płociniczak & Wróblewska (2020), 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 should note here that the full system is a system of two Hamiltonian equations with a holonomic constraint. The fact that the dynamics of such a system may be reduced to a system exhibiting asymptotic stability, as we show in the current paper, is in itself not obvious, and it is an interesting theoretical issue.

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 4.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. For the considered Poincaré map, we verified that at the bifurcation point, the analytical conditions for a transcritical bifurcation are satisfied (see the Apendix). Moreover, we checked that a transcritical bifurcation also occurs in the approximate solutions by Geyer et al. (2005) (not reported in Geyer et al. (2005)). However, the full SLIP model, by which we mean a switched system consisting of the two Hamiltonian systems (describing the support phase by equation (2.2) and the free flight phase) and a constraint on the touch-down angle, has not been investigated. It is not at all obvious if the transcritical bifurcation will be also observed there. We should add here that the small angle approximation for the vector field in the support phase has a significant importance for the construction of the reduced mapping in that the energy of the system has to be appropriately accounted for due to switchings from the support to the stance phase and vice versa. Putting it yet another way, due to the switched nature of the full model system, the implication of the small angle approximation on the dynamics of the full model system is not at all obvious. How possible bifurcations are perturbed by the reduction shown in this paper is an interesting issue which deserves separate investigations. The energy issue raised here is further explained in our recent work in Kowalczyk et al. (2023).

From the standpoint of mathematical modelling running an introduction of a modified model may be highly useful, which would better capture the fact that during switchings to and from the flight phase there is certainly present the dissipation of energy as well as the generation of energy by the tendon-muscle complex. Hence in the future work, model modifications will be focused on this point.

Conflict of interest

The authors declare that they have no conflict of interest.

Data availability

The datasets generated during and/or analysed during this study are available from the corresponding author on reasonable request.

Footnotes

References

A Appendix

Before we verify that the analytical conditions for the transcritical bifurcation are satisfied, we highlight here that parameter |$K$| must be considered as depending on both the state and the parameter variables, namely on the apex height |$Y$| and the energy system |$E$|⁠. We know that, in general, |$K = K (\alpha , \theta _d, L_d)$|⁠. Now, from fixed point equations (4.8) and (4.9), we can identify the relationship between |$(\alpha , \theta _d, L_d)$| and |$(Y,E)$|⁠. Hence, for fixed point solutions, we may equivalently express |$K (\alpha , \theta _d, L_d) $| as |$K(Y,E)$|⁠.

We are now ready to verify the conditions for a transcritical bifurcation at point |$(Y^*, E_s^+)$|⁠, discussed in Section 4.6, in terms of the derivatives of the Poincaré mapping |$f(Y,E)$| (see (4.4)), where the mapping is given by

We assume that at |$(Y,E)=(Y^*,E_s^+)$|⁠, (A.1) has a non-hyperbolic fixed point, i.e.

and

Now, in Fig. 8 we have two curves of fixed points passing trough |$(Y^*, E_s^+)$|⁠. So, by the implicit function theorem, it is necessary to have

Moreover, we need to get

Theorem 4.1 corresponds to the condition (A.2). The derivatives used under conditions (A.3), (A.4), (A.5) and (A.6) are shown below in their explicit form (for a comparison of the first derivative see (4.22)):

Hence, conditions (A.3), (A.4), (A.5) and (A.6) can be expressed as derivatives of the angle swept |$\varDelta \theta $| as follows:

Note that from the first line above it follows that at the bifurcation, the angle swept during the stance has a stationary point. To check whether this is indeed the case notice that (see (4.31))

So

and it can be shown that

where (see (4.5))

and

Moreover, the second partial derivatives (see second line of (A.8)) are different from zero at the bifurcation point |$(Y^*, E_s^+)=(0.9754,1.550)$| of the example in Section 4.6.