Solution and asymptotic analysis of a boundary value problem in the spring-mass model of running

TL;DR

This paper rigorously analyzes a boundary value problem in the spring-mass model of running, deriving asymptotic solutions and proving the existence and uniqueness of the solution, with implications for understanding stiffness dependence.

Contribution

It provides a rigorous asymptotic analysis of the boundary value problem in the spring-mass running model, including existence, uniqueness, and approximation of stiffness.

Findings

Existence and uniqueness of the solution to the boundary value problem.

Asymptotic estimates for the stiffness based on initial conditions.

Numerical validation of theoretical results.

Abstract

We consider the classic spring-mass model of running which is built upon an inverted elastic pendulum. In a natural way, there arises an interesting boundary value problem for the governing system of two nonlinear ordinary differential equations. It requires us to choose the stiffness to ascertain that after a complete step, the spring returns to its equilibrium position. Motivated by numerical calculations and real data we conduct a rigorous asymptotic analysis in terms of the Poicar\'e-Lindstedt series. The perturbation expansion is furnished by an interplay of two time scales what has an significant impact on the order of convergence. Further, we use these asymptotic estimates to prove that there exists a unique solution to the aforementioned boundary value problem and provide an approximation to the sought stiffness. Our results rigorously explain several observations made by other researchers concerning the dependence of stiffness on the initial angle of the stride and its velocity. The theory is illustrated with a number of numerical calculations.

Tables

Table 1. Table 1: Typical values of all appearing nondimensional physical parameters. Data based on [ 16 ] and calculated for various animals.

Symbol	Description	Typical value
$\alpha$	Angle of attack	$0.2-0.8$
$U$	Horizontal Froude number	$0.8-2.6$
$V$	Vertical Froude number	$0.05-0.5$

Equations

$$\left\{\begin{array}[]{l}m\dfrac{d^{2}x}{dt^{2}}=k\left(l_{0}-\sqrt{x^{2}+y^{2}}\right)\sin\theta,\vspace{6pt}\\ m\dfrac{d^{2}y}{dt^{2}}=k\left(l_{0}-\sqrt{x^{2}+y^{2}}\right)\cos\theta-mg,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\dfrac{d^{2}x}{dt^{2}}=\dfrac{k}{m}x\left(\dfrac{l_{0}}{\sqrt{x^{2}+y^{2}}}-1\right),\vspace{6pt}\\ \dfrac{d^{2}y}{dt^{2}}=\dfrac{k}{m}y\left(\dfrac{l_{0}}{\sqrt{x^{2}+y^{2}}}-1\right)-g.\end{array}\right.$$

$$x(0)=-l_0\sin \alpha ,\frac{dx}{dt}(0)=u,y(0)=l_0\cos \alpha ,\frac{dy}{dt}(0)=v,$$

$$\left\{\begin{array}[]{l}L\dfrac{d^{2}\theta}{dt^{2}}+2\dfrac{dL}{dt}\dfrac{d\theta}{dt}=\sin\theta,\vspace{6pt}\\ \dfrac{d^{2}L}{dt^{2}}-\left(\dfrac{d\theta}{dt}\right)^{2}L=K\left(1-L\right)-\cos\theta,\end{array}\right.$$

$$K=\frac{k{l}_0}{mg}.$$

$$\theta (0)=-\alpha ,\frac{d\theta }{dt}(0)={\theta }_d,L(0)=1,\frac{dL}{dt}(0)=-L_d,$$

$${\theta }_d=U\cos \alpha -V\sin \alpha ,L_d=U\sin \alpha +V\cos \alpha ,$$

$$U=\frac{u}{\sqrt{g{l}_0}},V=\frac{v}{\sqrt{g{l}_0}}.$$

$$\theta (t^{\ast },K^{\ast })=\alpha ,L(t^{\ast },K^{\ast })=1.$$

$$ϵ=\frac{1}{\sqrt{K}},$$

$$\tau =\frac{t}{ϵ},{\tau }^{+}=\omega (ϵ)\tau =\left(1+ϵ{\omega }_1+{ϵ}^2{\omega }_2+...\right)\frac{t}{ϵ}.$$

$$\begin{array}{r}\omega (ϵ{)}^2{L}^{\prime \prime }+\left(1-\omega (ϵ{)}^2({\theta }^{\prime }{)}^2\right)L=1-{ϵ}^2\cos \theta ,L{\theta }^{\prime \prime }+2L^{\prime }{\theta }^{\prime }=\omega (ϵ{)}^{-2}{ϵ}^2\sin \theta ,\end{array}$$

$$L=L_0+ϵL_1+{ϵ}^2{L}_2+...,\theta ={\theta }_0+ϵ{\theta }_1+{ϵ}^2{\theta }_2+...,$$

$$\begin{array}{rl}\omega (ϵ{)}^2\left(L_0^{\prime \prime }+ϵL_1^{\prime \prime }+{ϵ}^2{L}_2^{\prime \prime }+...\right)& +\left(1-\omega (ϵ{)}^2{\left({\theta }_0^{\prime }+ϵ{\theta }_1^{\prime }+{ϵ}^2{\theta }_2^{\prime }+...\right)}^2\right)\left(L_0+ϵL_1+{ϵ}^2{L}_2+...\right)\\ & =1-{ϵ}^2\cos \theta ,\end{array}$$

$$\begin{array}{rl}\left(L_0+ϵL_1+{ϵ}^2{L}_2+...\right)\left({\theta }_0^{\prime \prime }+ϵ{\theta }_1^{\prime \prime }+{ϵ}^2{\theta }_2^{\prime \prime }+...\right)+2& \left(L_0^{\prime }+ϵL_1^{\prime }+{ϵ}^2{L}_2^{\prime }+...\right)\left({\theta }_0^{\prime }+ϵ{\theta }_1^{\prime }+{ϵ}^2{\theta }_2^{\prime }+...\right)\\ & =\omega (ϵ{)}^{-2}{ϵ}^2\sin \left({\theta }_0+ϵ{\theta }_1+{ϵ}^2{\theta }_2+...\right).\end{array}$$

$$\begin{array}[]{lll}\epsilon^{0}:&L^{\prime\prime}_{0}+\left(1-(\theta_{0}^{\prime})^{2}\right)L_{0}=1,&L_{0}(0)=1,\quad L^{\prime}_{0}(0)=0,\\ \epsilon^{1}:&L^{\prime\prime}_{1}+\left(1-(\theta_{0}^{\prime})^{2}\right)L_{1}=-2\omega_{1}L_{0}^{\prime\prime}+2\left(\theta_{0}^{\prime}\theta_{1}^{\prime}+\omega_{1}(\theta_{0}^{\prime})^{2}\right)L_{0},&L_{1}(0)=0,\quad L^{\prime}_{1}(0)=-L_{d},\\ ...&&\end{array}$$

$$\begin{array}[]{lll}\epsilon^{0}:&L_{0}\theta_{0}^{\prime\prime}+2L_{0}^{\prime}\theta_{0}^{\prime}=0,&\theta_{0}(0)=-\alpha,\quad\theta^{\prime}_{0}(0)=0,\\ \epsilon^{1}:&L_{0}\theta^{\prime\prime}_{1}+2L_{0}^{\prime}\theta_{1}^{\prime}=-L_{1}\theta_{0}^{\prime\prime}-2L_{1}^{\prime}\theta_{0}^{\prime},&\theta_{1}(0)=0,\quad\theta^{\prime}_{1}(0)=\theta_{d},\\ ...&&\end{array}$$

$$L_0^2{\theta }_0^{\prime }=0,$$

$$L_1^{\prime \prime }+L_1=0,{\theta }_1^{\prime \prime }=0,$$

$$\begin{array}[]{lll}\epsilon^{2}:&L^{\prime\prime}_{2}+L_{2}=-2\omega_{1}L^{\prime\prime}_{1}+\theta_{d}^{2}-\cos\alpha,&L_{2}(0)=0,\quad L^{\prime}_{2}(0)=L_{d}\omega_{1},\\ \epsilon^{3}:&L^{\prime\prime}_{3}+L_{3}=&\\ &-\left(2\omega_{2}+\omega_{1}^{2}\right)L^{\prime\prime}_{1}-2\omega_{1}L_{2}^{\prime\prime}+\theta_{d}^{2}L_{1}+2\theta_{d}\left(\theta_{2}^{\prime}+2\omega_{1}\theta_{d}\right)-(\sin\alpha)\theta_{1},&L_{3}(0)=0,\quad L^{\prime}_{3}(0)=L_{d}(\omega_{2}-\omega_{1}^{2}),\\ ...&&\end{array}$$

$$\begin{array}[]{lll}\epsilon^{2}:&\theta_{2}^{\prime\prime}+2\theta_{d}L_{1}^{\prime}=-\sin\alpha,&\theta_{2}(0)=0,\quad\theta_{2}^{\prime}(0)=-\theta_{d}\omega_{1},\\ \epsilon^{3}:&\theta_{3}^{\prime\prime}+L_{1}\theta_{2}^{\prime\prime}=-2\left(L_{1}^{\prime}\theta_{2}^{\prime}+\theta_{d}L_{2}^{\prime}\right)+(\cos\alpha)\theta_{1}+\omega_{1}\sin\alpha,&\theta_{3}(0)=0,\quad\theta^{\prime}_{3}(0)=-\theta_{d}(\omega_{2}-\omega_{1}^{2}),\\ ...&&\end{array}$$

$${\omega }_2=-\frac{1}{2}{\theta }_d^2,$$

$$\begin{array}{r}L(t)\approx \tilde{L}(t):=1-ϵL_d\sin {\tau }^{+}(t)-{ϵ}^2\left(\cos \alpha -{\theta }_d^2\right)\left(1-\cos {\tau }^{+}(t)\right),\\ \text{where}{\tau }^{+}(t)=\tilde{\omega }(ϵ)\frac{t}{ϵ},\tilde{\omega }(ϵ)=1-\frac{{\theta }_d^2}{2}{ϵ}^2,ϵ\ll 1.\end{array}$$

$$\theta (t)\approx \tilde{\theta }(t):=-\alpha +ϵ{\theta }_d{\tau }^{+}(t)-\frac{1}{2}{ϵ}^2(\sin \alpha )({\tau }^{+}(t){)}^2+2{ϵ}^2{L}_d{\theta }_d\left(1-\cos {\tau }^{+}(t)\right),$$

$$L(\tau )=1-ϵ{\lambda }_1(\tau )-{ϵ}^2{\lambda }_2(\tau )-...,$$

$$\theta (t,\tau )={\theta }_0(t,\tau )+ϵ{\theta }_1(t,\tau )+{ϵ}^2{\theta }_2(t,\tau )+...$$

$$\begin{array}[]{llll}\theta_{0}(0,0)=-\alpha,&\dot{\theta}_{0}(0,0)=\theta_{d},&\theta^{\prime}_{0}(0,0)=0,&\\ \theta_{i}(0,0)=0,&\dot{\theta}_{i}(0,0)=0,&\theta^{\prime}_{i}(0,0)=0,&i>0,\end{array}$$

$$\begin{array}{rl}& \left(1-ϵ{\lambda }_1-...\right)\left({\ddot{\theta }}_0+{ϵ}^{-2}{\theta }_0^{\prime \prime }+2{ϵ}^{-1}{\dot{\theta }}_0^{\prime }+ϵ\left({\ddot{\theta }}_1+{ϵ}^{-2}{\theta }_1^{\prime \prime }+2{ϵ}^{-1}{\dot{\theta }}_1^{\prime }\right)+{ϵ}^2\left({\ddot{\theta }}_2+{ϵ}^{-2}{\theta }_2^{\prime \prime }+2{ϵ}^{-1}{\dot{\theta }}_2^{\prime }\right)\right)\\ & -2\left({\lambda }_1^{\prime }+ϵ{\lambda }_2^{\prime }+...\right)\left({\dot{\theta }}_0+{ϵ}^{-1}{\theta }_0^{\prime }+ϵ{\dot{\theta }}_1+{\theta }_1^{\prime }+...\right)=\sin \left({\theta }_0+ϵ{\theta }_1+{ϵ}^2{\theta }_2+...\right).\end{array}$$

$$\begin{array}[]{ll}\epsilon^{-2}:&\theta_{0}^{\prime\prime}=0,\\ \epsilon^{-1}:&2\dot{\theta}^{\prime}_{0}-\lambda_{1}\theta^{\prime\prime}_{0}+\theta^{\prime\prime}_{1}-2\lambda_{1}^{\prime}\theta_{0}^{\prime}=0,\\ \epsilon^{0}:&\ddot{\theta}_{0}+2\dot{\theta}^{\prime}_{1}-\lambda_{1}\left(2\dot{\theta}^{\prime}_{0}+\theta^{\prime\prime}_{1}\right)-\lambda_{2}\theta^{\prime\prime}_{0}-2\lambda^{\prime}_{1}\left(\dot{\theta}_{0}+\theta^{\prime}_{1}\right)-2\lambda^{\prime}_{2}\theta^{\prime}_{0}+\theta^{\prime\prime}_{2}=\sin\theta_{0},\\ ...\end{array}$$

$${\ddot{\theta }}_0+{\theta }_2^{\prime \prime }-2{\lambda }_1^{\prime }{\dot{\theta }}_0=\sin {\theta }_0.$$

Full text

**Keywords: singular perturbation theory, boundary value problem, Poincaré-Lindstedt series, elastic pendulum, running

**

AMS Classification: 34E10, 34B15

1 Introduction

Running is the fundamental way of rapid legged locomotion for terrestrial animals and due to its naturalness and everyday occurrence, it seems that there is nothing unusual in it. However, this way of movement requires a complex and accurate collaboration of neural, motor, and muscular systems with respect to the changing terrain [15]. The usual and common distinction between walking and running is that the latter contains an aerial phase during which the animal has no contact with the ground. This working definition is sufficient for us, however as research shows, it is too narrow to include certain animals or conditions of locomotion (see [10] for a broader classification based on the energetic concepts). Running is not just walking with a higher speed and there is a remarkable transition of one mean of motion to the other [2]. Legged locomotion of various animals has been investigated vigorously through the decades and it merges biology, engineering and mathematics into one successful endeavour. This topic was investigated by Aristotle [30] while the first biomechanical treatment was given by a seventeenth century Italian physiologist and mathematician Giovanni Borelli [9, 26]. The Reader can find interesting modern surveys concerning scientific accounts of locomotion in [7, 15, 36].

For humans, running gains another important meaning apart of being simple way of locomotion - namely, the sports. As a primal form of movement, running accompanies man from the very beginning. It is difficult not to agree with the fact that running is probably the simplest and the most natural sport that exist [3]. This is why it forms a basis for many other disciplines. In sport science it is quite common to analyse and to attempt to describe some movements that are specific to the considered discipline [17]. Many investigations lead to a better understanding of human performance during races of different lengths [14, 5, 38] which, in turn, provided better insights on improved training methods [13]. Mathematically, competitive running was described by Keller in his variational model [25] based on a physiological findings of Hill [21]. Some further generalizations and analysis are given in [32, 6, 37].

Mathematical modelling is an important part of the biomechanics. In this paper we are concerned with the so-called spring-mass model of running which is based on an inverted elastic pendulum (see the seminal papers [8, 27]). However, some earlier attempts accurately described walking with a similar construction utilizing inverted pendulum [29]. These first investigations lead to a proliferation of interesting concepts, models, and methods that helped to design walking robots [11]. Further generalizations of the spring-mass model are based on including additional legs, dampers and segments [19, 34, 35, 4, 33]. The two-legged version has an interesting bifurcation structure [28]. A very thorough review of models of legged locomotion is given in [24].

In this paper we revisit the conceptual spring-mass model. Our focus is to solve the naturally arising boundary-value problem for the spring stiffness via the asymptotic analysis. However, being essentially an elastic pendulum the mathematical description of the system consists of two nonlinear ordinary differential equations which possess a rich geometrical structure [18, 22] and chaotic behaviour [12, 1].

The paper is structured as follows. In Section 2 we state derive the model and state the main boundary value problem. By a numerical calculations we motivate that a perturbative expansion with respect to the large spring stiffness is relevant for the solution of the problem. Section 3 contains asymptotic analysis of main equations with the use of Poincaré-Lindstedt series. The material is divided into two parts: one gives heuristic derivation of the perturbation solutions while the other rigorously justifies them. In Section 4 we solve the initially stated boundary value problem and provide an approximation for its solution. We close the paper with several numerical calculations verifying our theory.

2 Model statement and motivation

The spring-mass model of running assumes that each leg can be described as an inverted elastic pendulum. For the grounded phase of the jump, we assume the situation depicted on Fig. 1.

Let $(x(t),y(t))$ denote the Cartesian coordinates of the point mass $m$. Balancing respective components of gravity and stiffness we can write

[TABLE]

where $l_{0}$ is the equilibrium length of the spring and $k$ is the stiffness. Plugging in the Cartesian formulas for the angle we obtain

[TABLE]

The initial conditions are the following

[TABLE]

where $u$ and $v$ are, respectively, horizontal and vertical velocities. It is very convenient to cast the governing equations (2) into nondimensional polar form. To this end, we naturally scale $x$ and $y$ with respect to $l_{0}$ and choose the pendulum time scale $\sqrt{g/l_{0}}$. Moreover, we introduce the nondimensional spring length $L=\sqrt{x^{2}+y^{2}}/l_{0}$ and the polar angle $\theta$. Eventually, the polar form of (2) is the following

[TABLE]

where the only nondimensional parameter (spring stiffness) present in the equations is given by

[TABLE]

The initial conditions for the polar coordinate system have the form

[TABLE]

with

[TABLE]

where we have defined the horizontal and vertical Froude numbers

[TABLE]

In Tab. 1 we have collected all nondimensional parameters appearing in the model. Notice that usually $U$ is of order of unity, while $V$ and $\alpha$ are small.

By the standard theory of ordinary differential equations the system (4) with (6) possesses a unique solution in the vicinity of $t=0$. On the other hand, from the point of view of applications a question of completely different nature is relevant. Since the inverted elastic pendulum models the forwardly hoping leg we are faced with a peculiar boundary value problem to solve.

Problem 1**.**

Let $(\theta(t;K),L(t,K))$ be the solution of the system (4) with (6). Find $K^{*}$ and the smallest time $t^{*}>0$ satisfying

[TABLE]

This means that the spring stiffness $K$ has to be determined to ensure that during the first cycle the spring will return to the equilibrium length precisely at the time for which the pendulum travels to the angle $\alpha$. This represents the grounded phase of the step, i.e. the leg completes the full cycle before jumping into the aerial phase (see Fig. 1).

The above problem can easily be solved numerically using the shooting method as was done for example in [27]. To illustrate this, we apply a numerical solver based on the forth order Runge-Kutta method for solving the initial value problem (4) and (6) with a given $K$. Then, the point $t^{*}$ is found such that $\theta(t^{*})=\alpha$. Next, $L(t^{*})$ is compared with $1$ and according to the difference the value of $K$ is corrected via the secant method for the next iteration. The loop continues until the required accuracy is attained.

Numerically found values of $K^{*}$ are depicted on Fig. 2 with respect to the initial angle of attack $\alpha$ for several values of $U$ and vice-versa. We see that in general $K^{*}$ is a moderately large parameter especially for small angles but also for realistic regime of constants (see Tab. 1). The dependence on $\alpha$ and $U$ is monotone and we can anticipate that for $U$ of orders of unity $K^{*}$ has a quadratic component of $U$. In the following sections we will prove these claims along with finding asymptotic expansions of $L$ and $\theta$ and proving existence and uniqueness for Problem 1.

3 Perturbation theory

In this section we provide an asymptotic analysis for solutions of the system (4) with (6) for large $K$. First, to get a hint how the possible asymptotic solution may look like we proceed with the formal singular perturbation theory. Then, we provide rigorous proofs concerning the order of approximation.

3.1 Formal expansions

The above numerical results (compare [27]) indicate that for a wide range of realistic initial conditions the appropriate stiffness being the solution of the Problem 1 is large. This suggests that we can gain a meaningful insight by expanding (4) for $K\rightarrow\infty$. To this end we will use the standard Poincaré-Lindstedt method (see [23]).

We will see that solutions live on different time scales. Since $K$ is the factor of $L$ in the second equation in (4) we expect that the dynamics takes place on the fast $\sqrt{K}t$ scale. On the other hand, the equation for the angle indicates that the main time scale for $\theta$ is the slow one $t$. Furthermore, due to the coefficient $\left(d\theta/dt\right)^{2}$ the period of oscillations is modulated by the slower time $t$. However, as we will shortly see, the solution of Problem 1 occurs on the fast scale and it will be more convenient to expand both variables with respect to that.

First, we will use the Poincaré-Lindstedt series to find the asymptotic behaviour of $L$ and $\theta$ as $K\rightarrow\infty$. To simplify matters we set

[TABLE]

hence we are looking for an expansion as $\epsilon\rightarrow 0$. Being lead by the above discussion we introduce the following fast time scales

[TABLE]

Now, substituting $\tau^{+}$ from (11) into (4) we obtain

[TABLE]

where now $L=L(\tau^{+})$, $\theta=\theta(\tau^{+})$, and prime denotes the derivative with respect to $\tau^{+}$. If we make the following formal asymptotic expansions

[TABLE]

the equations become

[TABLE]

and

[TABLE]

Next, collecting respective coefficients of $\epsilon^{n}$ in (12) we obtain the following chain of differential equations

[TABLE]

and

[TABLE]

To save space, we have written only the two first orders of $\epsilon$ since further equations complicate its form very quickly. When we obtain these initial approximation we will see that subsequent equations simplify considerably.

Starting from the $\epsilon^{0}$ equations we can multiply the one for $\theta_{0}$ by $L_{0}$, integrate, and obtain the conservation of angular momentum

[TABLE]

where we have used the fact that $\theta_{0}^{\prime}(0)=0$. This can only be satisfied if $\theta_{0}(\tau^{+})=-\alpha$. Therefore, the $L_{0}$ equation yields the solution $L_{0}(\tau^{+})=1$ and hence, the leading order solutions are constant. The $\epsilon^{1}$ equations are now simplified

[TABLE]

which quickly can be solved to obtain $L_{1}(\tau^{+})=-L_{d}\sin\tau^{+}$ and $\theta_{1}(\tau^{+})=\theta_{d}\tau^{+}$. These simple initial approximations simplify further asymptotic equations. The are the following

[TABLE]

and

[TABLE]

Notice that the $\epsilon^{2}$ order $L$-equation will not be forced by a resonant term only if we take $\omega_{1}=0$. Similarly, in the next order equation we can eliminate the secular terms when we take

[TABLE]

since $L^{\prime\prime}_{1}=-L_{1}$. Next, solving for $L_{2}$ yields the (formal) behaviour of the solution up to the second order

[TABLE]

On the other hand, now we can go back to the $\epsilon^{2}$ order equation for the angle and solve it immediately to obtain the approximation for $\theta$

[TABLE]

where $\tau^{+}(t)$ is defined in (23).

Having obtained the above approximation in the fast $\tau^{+}$ scale it is interesting to analyse the angle equation in (4) in the slow $t$ scale. Notice that since $L(t)\approx 1$ for $\epsilon\ll 1$ in the first approximation the fast time scale $\tau^{+}$ enters the equation only through the damping term. In that case we obtain a second order equation with a quickly varying coefficient which could be tackled by the homogenization theory (see for ex. [31]). However, in our case in order to find the leading order of $\theta$ it is more convenient to use the general multiple-scales method.

We start by assuming that $\theta=\theta(t,\tau)$. Since we are interested only in the leading order form of the asymptotic expansion we do not have to use the strained time scale $\tau^{+}$ which would introduce higher order terms. As an expansion for $L$ we use

[TABLE]

while the angle is expanded as follows

[TABLE]

The exact form of $\lambda_{1}$ can be inferred from (23). The initial conditions can be translated into the expansion as

[TABLE]

where, similarly as above, the dot indicates the derivative with respect to $t$ while prime is the derivative with respect to $\tau$. The equation can now be expanded to yield

[TABLE]

Comparing various powers of $\epsilon$ yields an array

[TABLE]

The first equation immediately gives us $\theta_{0}(t,\tau)=C_{0}(t)+D_{0}(t)\tau$. But from the initial conditions (27) we have $D_{0}(t)=0$ and hence $\theta^{\prime}_{0}=0$. As anticipated, the leading order term does not depend on the fast time scale $\tau$. Similarly, for the $\epsilon^{-1}$ equation initial conditions yield $\theta^{\prime}_{1}=0$ and hence $\theta_{1}(t,\tau)=C_{1}(t)$.

The first non-trivial behaviour comes from the $\epsilon^{0}$ equation which reduces to

[TABLE]

Now, since by assumption $t$ and $\tau$ are independent variables, we can average over the fast time scale to extract information only about the slow evolution. Integrating above yields

[TABLE]

where overline denotes the average over $\tau$ changing from [math] to $2\pi$. Since $\lambda_{1}$ is $2\pi$-periodic up to $O(\epsilon)$ due to (23), we have $\overline{\lambda^{\prime}_{1}}=0$. Moreover, since $\tau$ enters the equation only through the periodic terms we anticipate that $\theta_{2}$ is $\tau$ periodic and hence the average of its derivatives vanishes. This leaves

[TABLE]

which is the leading order equation for the evolution of $\theta_{0}$. This approximation suffices for our needs and we will not continue the multiple scales analysis. We therefore claim that

[TABLE]

where $\widetilde{\theta_{0}}$ is given by the solution of (32) (it can be solved analytically in terms of the Jacobi amplitude function $am$). Notice also that if we put $t=\epsilon\tau$ and expand the solution of (32) with respect to $\epsilon\rightarrow 0$ for fixed $\tau$ the first two terms correspond to the nonperiodic part of (24) with $\tau^{+}$ replaced by $\tau$. Moreover, if we make the reasonable small angle assumption $\sin\theta=\theta+O(\theta^{3})$ we will have

[TABLE]

Of course we can continue the multiple scales approach and obtain higher order terms. For this program to be successful, we should also include the $t$-scale expansion of the pendulum length $L$. We will not pursue this topic here since approximations soon become very complicated and are not needed in what follows.

The above analysis has been intended to be formal yet illustrative to clearly state the possible form of the singular asymptotic expansion and the interplay of multiple time scales. We now proceed to the rigorous proofs of the above results.

3.2 Rigorous proofs

First, we need the following result which is a generalization of Grönwall-Bellman’s lemma. Its generalized version has been proven in [20] but here, for completeness, we include a simplified proof of the case we need.

Lemma 1**.**

Let $\psi=\psi(t)$ and $f=f(t)$ be continuous and positive functions. Assume that the following inequality holds

[TABLE]

for $C>0$. Then

[TABLE]

In particular, where $f\equiv D=$const. we have

[TABLE]

Proof.

First, put $y(t)=\int_{0}^{t}(t-s)\psi(s)ds$. Then, it follows that $y^{\prime\prime}=y$ and hence from the assumption

[TABLE]

If we add and subtract $\sqrt{C}y^{\prime}$ and multiply both sides by $e^{\sqrt{C}t}$ we arrive at

[TABLE]

That is to say

[TABLE]

where we have used the fact that by the definition we have $y(0)=0$ and $y^{\prime}(0)=0$. Once again, multiplying by the integrating factor $e^{-\sqrt{C}t}$ we have

[TABLE]

The last integration gives

[TABLE]

and finally by evaluating the inner integral

[TABLE]

The assertion easily follows by the assumption and definition of $y$. The proof is complete. ∎

We are ready to prove the main result of this paper.

Theorem 1** (Fast time scale asymptotics).**

Let $(L(\tau^{+}),\theta(\tau^{+}))$ be the solution of (12). Then, the following asymptotic behaviour holds

[TABLE]

uniformly for $\tau^{+}\leq T<\infty$, where $\widetilde{L}$ and $\widetilde{\theta}$ are defined in (23) and (24).

Proof.

We begin by finding the asymptotic behaviour of $L$. Write

[TABLE]

and change the time scale into the strained fast time $\tau^{+}=\widetilde{\omega}(\epsilon)\epsilon^{-1}t$ (with a slight abuse of notation). According to (4) and (23) the equation for the remainder $\lambda$ has the form

[TABLE]

where, as before, the differentiation with respect to $\tau^{+}$ is denoted with a prime while the derivative with respect to $t$ is denoted with a dot, and

[TABLE]

Moreover, the initial conditions for the remainder are zero: $\lambda(0)=\dot{\lambda}(0)=0$. Notice that we have deliberately retained the $t$-derivative of $\theta$ evaluated at a point $\epsilon\widetilde{\omega}^{-1}\tau^{+}$.

The equation for $\lambda$ can be easily transformed into an integral equation with the use of the Green’s function

[TABLE]

Out claim is that $f(\tau^{+},\epsilon)=O(\epsilon^{3})$ uniformly with respect to $\tau^{+}\leq T$. To prove it, first notice that by (23) we have

[TABLE]

Since, by uniform continuity on compact intervals $\dot{\theta}(\epsilon\widetilde{\omega}^{-1}\tau^{+})^{2}=\theta_{d}^{2}+O(\epsilon)$, the first term in the definition of $f$, namely (47), is $O(\epsilon^{4})$ as $\epsilon\rightarrow 0$. Now, the terms in the bracket in (47) can be written as

[TABLE]

Note that $\cos\alpha-\cos\theta(\epsilon\widetilde{\omega}^{-1}\tau^{+})=O(\epsilon)$. Again, we use uniform continuity of $\dot{\theta}$ and conclude that the terms (50) are $O(\epsilon)$. Combining the previous two estimates we conclude that

[TABLE]

for some constant $D>0$. Now, from the integral equation (48) we can infer, for $C>0$, that

[TABLE]

where we have used the fact that $\dot{\theta}$ is bounded and $|\sin(u-v)|\leq|u-v|$. By Lemma 1 it follows that

[TABLE]

And the first part of the assertion is proved.

The second part proceeds analogously by writing

[TABLE]

for $\widetilde{\theta}$ defined in (24). In the $\tau^{+}$ scale the governing equation (12) yields the evolution of the remainder

[TABLE]

where

[TABLE]

We can transform (55) into its equivalent integral form by multiplying it by $L$ and integrating twice keeping in mind vanishing initial conditions. Hence,

[TABLE]

where the kernel is defined by

[TABLE]

First, we will extract the meaningful information from $g$. To this end use the first part of the theorem to write $L=1-\epsilon L_{d}\sin\tau^{+}+O(\epsilon^{2})$ uniformly with respect to $\tau^{+}\leq T$ as $\epsilon\rightarrow 0$ to obtain

[TABLE]

Now, the terms with $\cos\tau^{+}$ cancel leaving

[TABLE]

since all other terms are uniformly bounded. Whence, the integral equation (57) becomes

[TABLE]

Next, by a simple estimate we have

[TABLE]

which, by the fact that $\widetilde{\theta}+\alpha=O(\epsilon)$, implies

[TABLE]

Moreover, due to asymptotic expansion of $L$ the kernel can be written as

[TABLE]

uniformly for $\tau^{+}\leq T$ as $\epsilon\rightarrow 0$. Therefore, combining (63) and (64) with (61) we arrive at

[TABLE]

for some constants $E,F>0$ and we have used the fact that all the $O(\epsilon^{3})$ terms are uniform with respect to $\tau^{+}\leq T$. Invoking Lemma 1 finally yields

[TABLE]

This ends the proof. ∎

From the above proof we can immediately spot a place when the usual trade-off between the order of approximation and interval length can be made. Notice that in both final estimates (53) and (66) we could even allow for $\tau^{+}\leq\epsilon^{-1}T$ and still uniformly bound the exponential term. Therefore, the expansions should be valid on longer $\epsilon$-expanding intervals. As we will see this is only partially true and we loose one order of convergence.

Corollary 1**.**

Let $(L(t),\theta(t))$ be the solution of (12). Then, the following asymptotic behaviour holds

[TABLE]

on an expanding interval $\tau^{+}\leq\epsilon^{-1}T<\infty$, where $\widetilde{L}$ and $\widetilde{\theta}$ are defined in (23) and (24).

Proof.

It suffices to reanalyse the proof of Theorem 1. When estimating the size of (47) we used the assumption that $\epsilon\widetilde{\omega}\tau^{+}\rightarrow 0$ when $\epsilon\rightarrow 0$ uniformly with respect to $\tau^{+}$ on compact intervals. Now, since $\epsilon\tau^{+}$ stays bounded in that limit we can only conclude that $|f(\tau^{+},\epsilon)|\leq C\epsilon^{2}$. The final estimate (53) now reads

[TABLE]

The proof for $\theta$ has to be changed exactly in the same way. We obtain $g(\tau^{+},\epsilon)=\epsilon^{2}\sin\alpha+O(\epsilon^{2})$ and continue the reasoning accordingly. ∎

We thus see the interplay between the two time scales. The pendulum length $L$ oscillates on the $\tau^{+}$ scale with the period modulated by the evolution of $\theta$ on the slow scale. Since $\theta(\widetilde{\omega}^{-1}\epsilon\tau^{+})$ is uniformly continuous only on compact subsets of $\tau^{+}$, the asymptotic expansion looses one order to account for that on $\epsilon$-expanding $\tau^{+}$ intervals. It is also a very well known fact that without coupling between $L$ and $\theta$ the Poincaré-Lindstedt series would approximate the solutions with full order on $\epsilon$-expanding intervals.

Having in mind the above discussion we can prove the leading-order asymptotic expansion of $\theta$ on the slow $t$ scale.

Theorem 2** (Slow time scale asymptotics).**

Let $\theta=\theta(t)$ be solution of (4). Then,

[TABLE]

uniformly for $t\leq T$, where $\widetilde{\theta}_{0}$ is defined in (33).

Proof.

The integral equation for $\theta$ can be obtained by multiplying the first equation in (4) by $L$ which brings up the angular momentum

[TABLE]

The above can be integrated twice and manipulated to yield

[TABLE]

where the kernel is given by (58). Now, write

[TABLE]

where $\widetilde{\theta}_{0}(t)$ is a solution of (32) with the original initial conditions (6) and hence $\psi(0)=\dot{\psi}(0)=0$. Then, by the above argument we have

[TABLE]

Due to Corollary 1 we have $L(t)=1-\epsilon L_{d}\sin\left(\widetilde{\omega}\epsilon^{-1}t\right)+\lambda(t)$, where $\lambda=O(\epsilon^{2})$ as $\epsilon\rightarrow 0$ uniformly for $t\leq T$ since a compact interval for $t$ is $\epsilon$-expanding for $\tau^{+}$. Hence

[TABLE]

where we used the fact that $\ddot{\theta}_{0}=\sin\theta_{0}$.

The $\psi$ term is associated with $g_{1}$ and we can estimate it as follows

[TABLE]

where $E_{1}>0$. Further, for some constants $E_{2},F_{2}>0$, the next term is simply

[TABLE]

by the fact that $t\leq T$. Lastly, we have $\dot{\lambda}=O(\epsilon)$ and by integration by parts we obtain

[TABLE]

The term in the brackets vanishes because $G(t,t,\epsilon)=0$ which leaves us with

[TABLE]

where $E_{3}>0$. Combining our results for $g_{i}$, $i=1,2,3$, we now have

[TABLE]

for some constants $E,F>0$. Invoking Lemma 1 yields

[TABLE]

The proof is complete. ∎

We have thus found the exact asymptotic expansion of $L$ and $\theta$ on two time scales. Now, we will proceed to applying these results to solving the boundary value problems stated at the beginning of this paper.

4 Boundary value problem

Armed with above results we will proceed to reanalyse Problem 1. First, we will prove that we can always find its unique solution at least for sufficiently small initial angles $\alpha$.

Theorem 3**.**

There exists a number $\alpha_{0}>0$ such that Problem 1 has a unique solution for $|\alpha|<\alpha_{0}$. Moreover,

[TABLE]

and

[TABLE]

Proof.

We will work on the fast $\tau^{+}$ scale. Since $L=\widetilde{L}+O(\epsilon^{3})$ the time $\tau^{*}$ of the first return of $L$ to its initial condition satisfies

[TABLE]

as $\epsilon\rightarrow 0$. We can see that to the leading order $\tau^{*}\approx\pi$. This observation can be made more accurate. Let $\mu$ be the solution of the following equation

[TABLE]

then, by classical theory we will have $\tau^{*}=\mu+O(\epsilon^{2})$ as $\epsilon\rightarrow 0$. Since $\mu>0$ we have $\sin\mu=\sqrt{1-\cos^{2}\mu}$ and the above is a quadratic equation in $\cos\mu$. The solution is

[TABLE]

Hence, the skipped terms in (83) are of the same order as the difference $\tau^{*}-\mu$. Therefore,

[TABLE]

and (81) follows.

Now, since we know $\tau^{*}$ we require that $\theta(\tau^{*})=\alpha$. Multiplying $\theta$ equation (12) by $L$, integrating twice, and using the kernel (58) we obtain

[TABLE]

We have to show that there exists a number $\epsilon(\alpha)$ for which the above has a solution. To this end define

[TABLE]

Observe that $F(0,0)=0$ and

[TABLE]

When we put $\epsilon=0$ only the first term above survives and hence

[TABLE]

Therefore, by the Implicit Function Theorem it follows that there exists a number $\alpha_{0}$ and a function $\epsilon:(-\alpha_{0},\alpha_{0})\rightarrow\mathbb{R}$, such that $F(\alpha_{0},\epsilon(\alpha_{0}))=0$. The boundary value problem has thus a unique solution.

The last part of the proof is to find an approximation to the solution. Since $\tau^{*}\approx\pi+O(\epsilon^{2})$ we can use $\widetilde{\theta}$ to determine $\epsilon(\alpha)$ when we truncate $O(\epsilon^{2})$ terms. We have

[TABLE]

Solving and remembering that $\epsilon=K^{-1/2}$ yields (82) and the proof is complete. ∎

After proving our results it is required to find the conditions for $\widetilde{L}$, $\widetilde{\theta}$ and $\widetilde{K}^{*}$ to be good approximations of the corresponding quantities. First of all, the perturbation expansion of $L$ and $\theta$ have been obtained under the assumption that $K\rightarrow\infty$ (i.e. $\epsilon\rightarrow 0$) with all other parameters fixed. We can a posteriori verify the assumptions on these. This can be done by noting that the subsequent terms in the expansions (23) and (24) have to be of higher order. That is to say, we require that $\epsilon^{2}\theta_{d}^{2}\ll\epsilon L_{d}$ and $\epsilon^{2}L_{d}\theta_{d}\ll\epsilon\theta_{d}$. Expressing this in terms of $K$ yields the consistency condition

[TABLE]

For example, the above is satisfied if $L_{d},\theta_{d}=O(1)$. Moreover, from (82) we read that the requirement of $K\gg 1$ simultaneously keeping the conditions (92) satisfied forces $\alpha\ll 1$. Hence, the angle should be small what is also consistent with the data. Further, if we go back to (7) and express $L_{d}$ and $\theta_{d}$ in terms of the Froude numbers $U$ and $V$, we see that $L_{d},\theta_{d}\approx U$ since $V\leq U$ for fixed $\alpha$. Whence, (92) reduce to

[TABLE]

Moreover, from the above we can infer about the validity of the approximation (82) which can be plugged into above to conclude that

[TABLE]

which is consistent with the small angle assumption. Referring to Tab. 1 we see that our asymptotic approximations are valid in the realistic regime of parameters.

In [27] several approximations of $K^{*}$ have been proposed. Authors claimed that for slower velocities, the dependence of $K^{*}$ on $U$ should be quadratic while for faster velocities, linear. In both of these cases Authors gave very complex fitted empirical formulas which closely reproduced the numerical results. Our approximation (82) gives a systematic explanation of the leading order behaviour of $K^{*}$ for small angles. It can also be treated as an approximation of the quadratic part of dependence on $U$. Note however, that some other components of the velocity might be missing and finding them is a subject of our future work. Furthermore, in the cited work Authors introduce the so-called effective vertical stiffness $K_{vert}$ which reduce to $K$ when the subject is required to jump vertically (i.e. $\alpha=0$ and $V=0$). Authors heuristically derive that

[TABLE]

where $Q$ is an unknown constant dependent on the contact time. Notice the similarity to our systematically devised result (82). This suggest that for small velocities and angles, the two stiffness parameters behave in a similar fashion.

We will illustrate our results with several numerical simulations. On Fig. 3 we can see an exemplary verification of Theorem 1. Error is calculated on the fast $\tau^{+}$ scale by choosing a compact interval $[0,\pi]$ and comparing solutions of (12) with their approximations (23) and (24) at $\tau^{+}=\pi$ for the worst case. Observe that on the log-log scale the plots become parallel to the superimposed $K^{-3/2}$ line indicating the correct order of convergence. Note also that we have used the original variables $K$, $U$ and $V$.

The next example concerns the validity of (82) as the approximation to the solution of Problem 1. Since, as we noted above, the natural assumption for its accuracy is $\alpha\ll 1$, we compare $K^{*}$ with $\widetilde{K}^{*}$ for different values of the angle. Results are given on Fig. 4. We can see that both values are close to each other and their ration converge to $1$ when $\alpha\rightarrow 0^{+}$. Note, however, that although decently accurate, $\widetilde{K}^{*}$ is only the leading order approximation for $K^{*}$. Finding the subsequent corrections is the aim of our future work.

5 Conclusion

We have solved a nonlinear boundary value problem that is very natural to modelling legged locomotion. It appeared that the equations live on multiple time scales, however, the problem had its solution on the fast scale $\tau^{+}$. Having proved the validity of asymptotic expansions we had used them in applying the Implicit Function Theorem to grant the existence and uniqueness of solution to Problem 1. It is also worth to mention that the approximation of the stiffness (82) is consistent with all of the previously speculated features of the numerical solution. We have thus justified several claims about its behaviour for a realistic regime of parameters.

The further work is will be based on considering expansions for larger velocities and finding the subsequent terms in the expansion for $K^{*}$ with respect to small $\alpha$. As was noted in [27] the stiffness starts to depend linearly for large values of $U$. We also plan to justify this claim analytically.