Asymptotic Solution of a Boundary Value Problem for a Spring–Mass Model of Legged Locomotion

Running is the basic mode of fast locomotion for legged animals. One of the most successful mathematical descriptions of this gait is the so-called spring–mass model constructed upon an inverted elastic pendulum. In the description of the grounded phase of the step, an interesting boundary value problem arises where one has to determine the leg stiffness. In this paper, we find asymptotic expansions of the stiffness. These are conducted perturbatively: once with respect to small angles of attack, and once for large velocities. Our findings are in agreement with previous results and numerical simulations. In particular, we show that the leg stiffness is inversely proportional to the square of the attack angle for its small values, and proportional to the velocity for large speeds. We give exact asymptotic formulas to several orders and conclude the paper with a numerical verification.

human or dog is itself a complex interplay of motor, neural, and muscular systems (Dickinson et al. 2000;Gordon et al. 2017). Therefore, the research on movement is intense and spans many disciplines such as biomechanics, sports medicine, applied mathematics, and robotics (Biewener and Patek 2018;Daniels et al. 1978;Tibshirani 2005;Collins et al. 2005). Building a complete model of running can be prohibitively difficult and, thus, to proceed one may just focus on some particular features or simplify matters with a use of low-dimensional conceptual models (Dickinson et al. 2000). There are many approaches to that problem. For example, in sports one usually is concerned about finishing the race in a shortest time. Hence, a particular model can focus only on several relevant parameters such as runner's power, air drag, or race strategy. This was a basis for a conceptual model of Keller proposed in the second half of the twentieth century (Keller 1973) (which, in turn, stemmed from fundamental ideas of Hill forwarded about fifty years earlier Hill 1925). An interested reader can find some more modern accounts of competitive running models in Woodside (1991), Pritchard (1993), Aftalion (2017), Aftalion and Martinon (2019).
In this paper, we are mainly concerned not in race performance but rather the very mechanism of running. In particular, we analyse the classical spring-mass model that has been highly successful in describing the motion and dynamics of legged locomotion. The model has been firstly analysed in seminal papers (Blickhan 1989;McMahon and Cheng 1990) where authors compared theoretical predictions with real data. This model is based on an inverted elastic pendulum which acts in two phases: grounded (stance) and aerial. As was noted in Blickhan and Full (1993), despite of tremendous diversity of various locomotion modes, anatomical differences, and sizes of particular animals, the spring-mass model describes the fundamental mechanical principles of bouncing gaits very accurately. Therefore, this simple model is very robust and has been investigated in many studies from different vantage points: experimental (Farley et al. 1991(Farley et al. , 1993Dickinson et al. 2000) and biomechanical modelling (He et al. 1991;Geyer et al. 2005). Further, the spring-mass model was a starting point of many projects in robotics (where it is usually called the Spring Loaded Inverted-Pendulum (SLIP) model). From the first SLIP-based hopper (Raibert 1986) to more advanced modern robots (Rummel and Seyfarth 2008;Akinfiev et al. 2003;Armour et al. 2007), this seemingly simple dynamical system inspired the design of many fast legged robots. A very readable survey about locomotion in robotics can be found in Aguilar et al. (2016). Apart from that, the basic spring-mass model has been generalized and extended in many different directions. For example: multiple leg hoppers (Gan and Remy 2014;Geyer et al. 2006), dynamic swing leg motion when the attack angle is not constant (Gan et al. 2018), damping (Saranlı et al. 2010), and various control problems (Ghigliazza et al. 2005;Sato and Buehler 2004;Takahashi et al. 2017;Shahbazi et al. 2016). On the other hand, from the mathematical point of view, the spring-mass model possesses interesting dynamical properties such as complex bifurcation diagrams (Merker et al. 2015) or periodic gait stability issues (Geyer et al. 2005;Seipel and Holmes 2007;Hamzaçebi and Morgül 2017). A detailed numerical analysis of the dynamics has recently been conducted in Zaytsev et al. (2019). Moreover, in certain situations one aims to find approximate solutions of the considered model since they can be used in studying global dynamics under certain assumptions. This programme has been conducted for example in Geyer et al. (2005), Schwind and Koditschek (2000), Ghigliazza et al. (2005), Saranlı et al. (2010) and also in the prequel of our present investigations Płociniczak and Wróblewska (2020) which we conduct in a similar spirit. An interested reader will find a very thorough account of running models formulated as dynamical systems in Holmes et al. (2006).
The spring stiffness present in the model has to be adjusted accordingly in order to smoothly move from the grounded phase into the aerial. This gives rise to an interesting nonlinear boundary value problem. In our previous work (Płociniczak and Wróblewska 2020), we have proved that it has a unique solution for sufficiently small angles of attack (realistic assumption). We have also provided several asymptotic expansion of the solution for small angles. As a by-product, we have obtained a useful approximation of the important parameter being the leg stiffness. This work expands and completes several ideas born therein. While in Płociniczak and Wróblewska (2020), we were concerned mostly about various properties of solutions, here we are mainly concerned about the leg stiffness, its asymptotic form, and dependence on several other physical parameters. Since it is impossible to find an exact form of the stiffness, we have to use other methods to infer about its behaviour. In particular, we prove a rigorous result that gives a two-term asymptotic expansion of the stiffness for small angles. In the leading order, it coincides with our previous heuristics. Moreover, we consider the mathematically interesting singular case of large velocity in which we have to use a two-parameter perturbation. These results give very simple expressions that reveal some biomechanical properties of the running or hoping leg. Furthermore, numerical analysis shows that they are accurate and, hence, may be used in practical work. To sum up, this paper extends, formalizes, and improves the estimate on leg stiffness found in Płociniczak and Wróblewska (2020) along with additional material concerning large velocity asymptotics.
In the next section, we introduce the model and devise the main boundary value problem. The subsequent section deals with perturbation methods and asymptotic solution. There, we state our main results. Throughout the paper, we include some numerical examples that illustrate the theory.

Model and Problem Statement
The schematic of the model setting is depicted in Fig. 1. One leg is modelled as an inverted elastic pendulum with an axis on the ground. By equating forces in the Cartesian coordinates, we can write (for details, see Blickhan 1989;McMahon and Cheng 1990; Płociniczak and Wróblewska 2020) where l 0 is the equilibrium length of the spring and k is its stiffness. We supplement the above system with initial conditions where u and v are, respectively, horizontal and vertical velocities, and α is the angle of attack. As we noted in Płociniczak and Wróblewska (2020), it is much more convenient to analyse (1) in polar coordinates. Transforming and rescaling leads us to where L = x 2 + y 2 /l 0 , and the nondimensional stiffness is given by The initial conditions now have the form with the horizontal and vertical Froude numbers Data from McMahon and Cheng (1990), Farley et al. (1993) indicate that for a wide variety of animals we have therefore α and V can be considered small, while U is of order of unity.
In the above biomechanical setting, a certain boundary value problem naturally arises. It is to find the stiffness K * and the smallest time t * for which we have That is to say, we have to determine the stiffness for which the leg has travelled an angle 2α precisely at the time of returning to its original length (see Fig. 1). In Płociniczak and Wróblewska (2020), we have shown that this problem has a unique solution at least for small α and provided some approximations to K * and t * . Specifically, we have found that for U of order of unity we formally have In what follows, we show that the above approximations are indeed the leading order expansions as α → 0, more precisely we prove that Moreover, we also consider an interesting case of large U for which In this way, we prove that K * depends linearly on U for large velocities what has been observed before with numerical means only. Notice that above formulas are asymptotic also for small α. However, led by our analytical results we also show that K * can be very accurately approximated with the following quasi-empirical formula for all meaningful angles The above extrapolates the small-angle leading order asymptotic behaviour of K * into the larger interval α ∈ [0, 1].

Asymptotic Solution
We begin by noticing that the main equation (3) for L can be written as It follows that the nature of the solution depends on the relative sizes of K and dθ/dt. The latter, in turn, has a magnitude of order of U . Loosely speaking, the solutions are either of trigonometric or hyperbolic nature. This makes a profound difference for the solution of the problem (8). Below, we explore these two regimes perturbatively.
Because of a large number of different quantities, the scope of below two subsections is local. This means that all introduced variable names are valid only within the respective subsection. Specifically, x and y are different in each of them.

Small Angle
The solution of the boundary value problem (8) requires that K = K (α). We start our analysis with an observation that K (α)α 2 is bounded when α → 0 + . This makes the subsequent perturbation theory possible. Moreover, this result has important implications since previously the angular dependence of K * has been known only numerically.
Lemma 1 Let (K * (α), t * (α)) be the solution of the boundary value problem (8). Then, Proof Start by using a new time scale Then, the θ -equation in (3) can be transformed into the integral form by multiplying it by L and integrating twice. We have where G is the kernel whose explicit form can be found. Setting θ( τ ) = α, we obtain an implicit equation Now, in Płociniczak and Wróblewska (2020) we have shown that there exists a monotone relation K * = K * (α) for sufficiently small α. Moreover, K * (α) → ∞ and τ * (α) → π as α → 0 + . If we take the same limit in the equation above the second term on the right vanishes since, as can be checked with an elementary calculation, G is uniformly bounded for each α. Our previous results (Płociniczak and Wróblewska 2020) also indicate that in this limit L → 1 and hence, for the both sides to be equal to each other we should have K (α) ∼ π 2 U 2 α −2 /4 as α → 0 + . This finishes the proof.
Since we are interested in solving the boundary value problem (8), we will rescale the relevant variables in a natural way to simplify matters. First, set which are initial conditions for L and θ as in (5). Then, we put for new dependent variables x and y. As we can see x measures the spring departure from initial length while y is normalized angle. The magnitude of 1 − L is chosen according to the asymptotic analysis done in our previous work (Płociniczak and Wróblewska 2020) (but see also Geyer et al. 2005). As for the time scale, we choose the one balancing L with L in (3) where τ is the strained fast time variable introduced because of the frequency modulation. Notice that ω is terminated after the second-order term. We will not conduct the higher-order expansion. In these new variables, Eq.
(3) has the form with initial conditions Here, prime denotes the derivative with respect to the fast time τ . We can see that due to Lemma 1 the derivative y (0) is bounded as α → 0 + . Now, we make the following perturbation expansions for the solutions of ODEs and the unknown stiffness where we notice that it is more feasible to expand the square root of K . Plugging these into (21), collecting respective powers of α, and using initial conditions (22) yields the following systems The α 0 equation can immediately be solved Now, since x 0 (τ ) = − sin τ we can remove the secular terms from the α 1 equation only by taking ω 1 = 0. The corresponding solution will then have the form Further, the second-order correction to the frequency ω 2 can be determined from the α 2 equation Therefore, by removing the secular terms we have The process can be continued in a standard way however, as can be seen, equations become very cluttered. Nevertheless, we have all the expansions needed to determine k −1 and k 0 . Our systematic expansions can be compared with more heuristic approach in Geyer et al. (2006) where Authors find some accurate approximations, couple them with the aerial phase, and investigate the stability of such model. We require that K (α) ∼ K * (α) as α → 0 + that is to say, coefficients k i have to be chosen to asymptotically solve the boundary value problem (8). To this end, we make the further expansion for the time t * τ * (α) = τ 0 + τ 1 α + τ 2 α 2 + · · · .
Now, t i and k i can be found by perturbatively solving Using (26), (27), and collecting respective powers of α, we obtain a series of linear systems Solving the above leads us to It is now only a matter of returning to the original time scale t and squaring √ K (α). Therefore, we have proved the following result.
The numerical example of the above expansion is depicted in Fig. 2. We have plotted the relative error of the asymptotic approximation K * i of the numerical solution K * of (8). Here, i ∈ {−1, 0, 1} denotes the largest order of α taken in the expansion. Clearly, all errors converge to 0 for small angles and what could have been anticipated, K * 1 looses accuracy for larger values of α (say, α > 0.3). We can see that K * 0 provides a reasonable approximation even for larger angles.

Large Velocity
Similarly as in the previous section, we begin by finding the leading order of the expansion K = K (U ) when U → ∞. Here however, we determine only how fast K (U ) approaches infinity as U increases.

Lemma 2 Let (K * (U ), t * (U )) be the solution of the boundary value problem (8). Then K (U )U −1 is bounded when U → ∞.
Proof Similarly as in the proof of Lemma 1, we start by introducing a faster time scale Now, we have K = K (U ) and the initial derivatives have the form Since, by assumption (8) has a solution for U → ∞ the derivatives has to stay bounded in the limit. We have then two possibilities: with initial conditions Notice that the equation for θ can be rewritten as a conservation of angular momentum which, using the initial conditions, yields θ( τ ) = −α for all τ which is a contradiction since it is impossible to have θ = α. Therefore, we are left with U /K (U ) > 0 for U → ∞ what concludes the proof.
Thanks to the above result we know that K (U ) ∝ U for large U . This observation has been made before in McMahon and Cheng (1990) but only on the basis of numerical calculations. In what follows we will find the proportionality constant in this dependence. It appears that this problem is more difficult than the one considered before since the first-order perturbation equation does not posses an analytical solution. However, we will see that further expanding in the small α can help to obtain useful results.
We start by introducing appropriate scaling. Let We have thus rescaled the angle according to its range in the considered problem. The time scale has been chosen based on the linear approximation to θ for which θ(t) ≈ −α +(U cos α −V sin α)t. Then, θ(t) ≈ α when t ≈ 2α(U cos α −V sin α) −1 which gives the above chosen time scale for large U and small α. (1 − 2y))) , where we have denoted the small parameter and the prime denotes a derivative with respect to τ . The initial conditions (5) have the form (42) Therefore, we can see that both dependent variables and their derivatives are now of order of unity.
We are ready to make the perturbation expansion for which (40) becomes an array of equations 0 : (44) Notice that the above equations are nonlinear in all orders. However, the 0 system can be exactly integrated in a closed form. To see this notice that from the second equation we have the conservation of angular momentum From which it follows that Plugging it to the first equation, we obtain a single nonlinear ODE This equation has an analytical solution which can be obtained by reduction in order After returning to y 0 and integrating, we have Notice that this equation is degenerate in the sense that it does not contain the K 1 coefficient but, nevertheless, automatically solves (8). To wit, note that the leading order time of the solution is for which we have both x 0 (τ * ) = 1 and y 0 (τ * ) = 0. In order to gain some insight how does the solution of (40) behave for small positive we have to analyze the 1 equation in (44). Unfortunately, we cannot obtain its closedform solution. However, since α is small, we can at least try to find the subsequent perturbation approximation. To this end, put x 1 (τ ) = x 10 (τ )+x 12 (τ )α 2 +O(α 4 ), y 1 (τ ) = y 10 (τ )+y 12 (τ )α 2 +O(α 4 ) as α → 0.
This system is easy to be solved explicitly yielding Having these expansions in hand, we can proceed to the solution (8). To this end, we set and then require that as → 0 + and α → 0 + . The asymptotic solution of the above algebraic system is τ * 0 given by (50) and which coincides with (11) after returning to the original time scale. Therefore, we have proved the following result.
Notice that according to (56), we have K 1 ∝ V α −2 for small angles. Since our approximation is only the asymptotic leading order for angles and not the exact solution, it is interesting to extend this result numerically and find the empirical proportionality constant valid for a larger set of α. Let The ratio of K 1 to the numerical value of K 1 found by solving the boundary value problem (8) corresponding to Eq. (40) is plotted in Fig. 3. As we can see, the analysed quantity is almost independent of V . Whence, it is motivated to stipulate that for some G with G(0) = 1. Numerical calculations plotting K 1 / K 1 are shown in Fig. 4. We can see that the sought form of G(α) is approximately quadratic and the least squares best fit is given by with the determination coefficient R 2 = 0.999. The accuracy is thus very good for α ∈ [0, 1]. On Fig. 5, we plot the relative error of the approximation (11) to the solution of the problem (8) for the original Eq. (3). Notice that the approximations achieve an accuracy of several percent for U = O(10 2 ). Also, the corrected coefficient K p as in (58) performs very well even for large angles while the leading order term K 1 does not converge to 0 for such.

Conclusion
The perturbation analysis conducted above helped us to determine the asymptotic form of the solution of (8). As was shown, the physically important case of small angles, in principle could have been resolved for any required order. However, the complexity of expressions forced us to stop at two meaningful terms. Nevertheless, the rest can be obtained in a standard way with a use of some symbolic manipulation environment such as Wolfram Mathematica as was in our case. Moreover, found analytic formulas proved to be very accurate in their region of validity and explained how does the leg stiffness depends on the attack angle and velocities.
The case of large (horizontal) velocity also proved to be attackable with the perturbation theory. However, the leading order equations were nonlinear forcing higher orders not to be resolvable analytically. Despite this difficulty, we have been able to use the additional expansion in small angles and find leading order behaviour in two parameters. This let us to determine the form of dependence of the leg stiffness to other physical parameters and, thanks to that knowledge, to empirically extend the approximation for all meaningful angles.
Asymptotic analysis was only possible due to the fact that we firstly have found the leading order behaviours of K * (α) and K * (U ) in the respective limits. Moreover, the key step was to perturbatively look for expansions of the stiffness K * and the time t * . Of course, this programme can be carried over to higher orders with a symbolic manipulation software. However, our findings are sufficient to understand the fundamental relations between stiffness and various physical parameters. We hope that this research will contribute to the understanding of various gaits in locomotion and robotics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Tibshirani, R.: Who is the fastest man in the world? In: Anthology of Statistics in Sports, pp. 311-316. SIAM (2005) Woodside, W.: The optimal strategy for running a race (a mathematical model for world records from 50 m to 275 km). Math. Comput. Modell. 15(10), 1-12 (1991) Zaytsev, P., Cnops, T., David Remy, C.: A detailed look at the SLIP model dynamics: bifurcations, chaotic behavior, and fractal basins of attraction. J. Comput. Nonlinear Dyn. 14(8), 081002 (2019) Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.