Mathematical Model of the Movement of the Human Limbs

1. Introduction

The plane double gravitational pendulum serves as a valuable model for understanding the complexities arising from nonlinearity in physical systems. The system’s behavior is described by a set of coupled, nonlinear, second-order differential equations that defy analytical solution [1]. Numerical simulation provides a powerful tool for exploring its dynamics, but requires careful attention to numerical errors and energy conservation. The double pendulum’s sensitivity to initial conditions and its chaotic nature underscore the limitations of long-term predictability in nonlinear systems. Visualization tools, such as animations and phase space plots, are crucial for analyzing and interpreting the simulation results [2]. In addition to numerically solving the equations we have also used these methods. The numerical exploration of the double pendulum’s dynamics reveals its most compelling characteristic: its sensitivity to initial conditions. This sensitivity, a hallmark of chaotic systems, means that infinitesimally small differences in the initial angles or velocities of the pendulums can lead to essentially different trajectories over time [3,4]. The chaotic nature of the double pendulum can be most appealingly visualized through phase space plots. These plots, which map the system’s state over time, reveal intricate and often fractal-like structures. Unlike simpler systems with predictable periodic motion that result in closed, repeating loops in phase space, the double pendulum exhibits trajectories that wander seemingly randomly through the phase space, never repeating themselves exactly. This behavior highlights the system’s inherent unpredictability and its departure from the predictable dynamics of linear systems.

The current article demonstrates that it can be helpful in matters of biomechanics, which is related to the problems of human walking. We show that within the typical range and parameters of the human motion the double pendulum behaves relatively regularly and does not have essential chaotic behavior.

The pendulum systems are classical models of nonlinear dynamics [5]. We would like to stress that, generally, the study of the double pendulum is not merely an academic exercise; it holds significance in various scientific and engineering fields [6,7,8,9,10,11,12,13]. Understanding its complexities aids in comprehending the dynamics of multi-body systems, informing the design of robotic arms, satellite stabilization systems, financial modeling, and even contributing to our understanding of planetary motion. Furthermore, it turns helpful for apprehending questions related, as stated, to physics of human motion, and also climate modeling [14,15] (for which Edward Lorenz’s coined the famous term “butterfly effect”, in article titled “Deterministic Nonperiodic Flow” published in his 1963 paper in the Journal of the Atmospheric Sciences), fundamental physics, and even entertaining sport activities, like golf, namely for modeling the golf swing with double pendulums [16]. The understanding of coupled oscillators is crucial in analyzing the behavior of electrical circuits, mechanical vibrations, and even the rhythmic firing of neurons in the brain.

The planar double gravitational pendulum, a seemingly simple system composed of two pendulums connected in series and swinging under the influence of gravity, belies its deceptive simplicity. It embodies a profound intersection of classical mechanics, nonlinear dynamics, and chaos theory, serving as a potent example of how seemingly deterministic systems can exhibit highly unpredictable behavior. A large number of investigations, both theoretical and experimental, have been dedicated to the problem of determining the chaotic behavior of double pendulum [17,18,19].

2. The Model

In this section we describe a model of the limbs of the human body. In its simplest version the model represents a type of gravitational pendulum. They can be reduced to plane double or triple gravitational pendulum, i.e., we consider the corresponding limb of the human body to embrace two segments. Any of the branches of the pendulum are supposed to be axis symmetric and characterized by the corresponding densities ${\rho }_i\left(s\right)$, $i=1,2$, where s is defined so that it unicly determines the position of the corresponding cross-section of the given branch of the combined.

Of course, we need to determine the Lagrangian of the system. For that we need to clarify what is the kinetic and potential energies of the system in terms of the properly defined variables - see Figure 1. For the total kinetic energy K we obtain consequently

$$K=K_1+K_2$$

(1)

where $K_i$, $i=1,2$ are the kinetic energies of the segments. Thus, for $K_1$ one has

$$K_1={\displaystyle \frac{1}{2}}{\int }_{\left\{S_1\right\}}{\rho }_1\left(s_1\right)A\left(s_1\right)v_1^2(s_1,t)d{s}_1,$$

(2)

where $\left\{S_1\right\}$ is the set of allowed values of $s_1\in [0,L_1]$ and $A\left(s_1\right)$ is the area of the cross-section of segment 1 at position $s_1$ measured from the beginning of the segment. In what follows we will take the densities of a given segment to be uniform, i.e., ${\rho }_i\left(s\right)={\rho }_i$. The last actually reflects the experimental situation - data for the mass densities of a given segment differ from segment to segment, but within the specific segment they are assumed homogeneous. Then, taking into account that

$$v_1^2(s_1,t)={\dot{x}}_1^2(s_1,t)+{\dot{y}}_1^2(s_1,t),$$

(3)

and that

$$x_1(s_1,t)=s_{1}sin\left({\theta }_1\left(t\right)\right),\mathrm{and}{y}_1(s_1,t)=-s_{1}cos\left({\theta }_1\left(t\right)\right),$$

(4)

we immediately obtain

$${\dot{x}}_1(s_1,t)=s_{1}cos\left[{\theta }_1\left(t\right)\right]{\dot{\theta }}_1\left(t\right),\mathrm{and}{\dot{y}}_1(s_1,t)=s_{1}sin\left[{\theta }_1\left(t\right)\right]{\dot{\theta }}_1\left(t\right),$$

(5)

and, therefore

$$v_1^2(s_1,t)={\dot{x}}_1^2(s_1,t)+{\dot{y}}_1^2(s_1,t)=s_1^2{\dot{\theta }}_1^2\left(t\right),$$

(6)

which leads to

$$\begin{array}{c}\hfill K_1={\displaystyle \frac{1}{2}}{\rho }_1{\int }_{\left\{S_1\right\}}A\left(s_1\right)v_1^2(s_1,t)d{s}_1={\displaystyle \frac{1}{2}}{\rho }_1{\dot{\theta }}_1^2\left(t\right){\int }_{\left\{S_1\right\}}A\left(s_1\right)s_1^{2}d{s}_1={\displaystyle \frac{1}{2}}{\dot{\theta }}_1^2\left(t\right)I_1.\end{array}$$

(7)

Here

$$I_1={\rho }_1{\int }_0^{L_1}A\left(s_1\right)s_1^{2}d{s}_1$$

(8)

is a moment of inertia of segment “1” about an axes perpendicular to it and centered at its top. We have assumed that the length of the segment is $L_1$.

Now we shall deal with the kinetic energy $K_2$ of the second segment in the limb. Let $s_2\in [0,L_2]$ is the coordinate along the second segment. Then, looking at the geometry one derives

$$\begin{array}{ccc}\hfill x_2(s_2,t)& =& L_{1}sin{\theta }_1\left(t\right)+s_{2}sin{\theta }_2\left(t\right),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill y_2(s_2,t)& =& -L_{1}cos{\theta }_1\left(t\right)-s_{2}cos{\theta }_2\left(t\right),\hfill \end{array}$$

(10)

which leads to

$$\begin{array}{ccc}\hfill {\dot{x}}_2(s_2,t)& =& L_{1}cos{\theta }_1\left(t\right){\dot{\theta }}_1\left(t\right)+s_{2}cos{\theta }_2\left(t\right)\dot{{\theta }_2}\left(t\right),\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\dot{y}}_2(s_2,t)& =& -L_{1}sin{\theta }_1\left(t\right){\dot{\theta }}_1\left(t\right)+s_{2}sin{\theta }_2\left(t\right)\dot{{\theta }_2}\left(t\right),\hfill \end{array}$$

(12)

and, therefore, to

$$v_2^2(s_2,t)=L_1^2{\dot{\theta }}_1^2\left(t\right)+s_2^2{\dot{\theta }}_2^2\left(t\right)+2L_1{s}_{2}cos[{\theta }_1\left(t\right)-{\theta }_2\left(t\right)]{\dot{\theta }}_1\left(t\right){\dot{\theta }}_2\left(t\right).$$

(13)

Using the above expressions, for the kinetic energy $K_2$ we derive, consequently:

$$\begin{array}{ccc}\hfill K_2& =& {\displaystyle \frac{1}{2}}{\int }_{S_2}{\rho }_2\left(s_2\right)A\left(s_2\right)v_2^2(s_2,t)d{s}_2\hfill \\ & =& {\displaystyle \frac{1}{2}}{\rho }_2{\int }_{\left\{S_2\right\}}\left[L_1^2{\dot{\theta }}_1^2\left(t\right)+s_2^2{\dot{\theta }}_2^2\left(t\right)+2L_1{s}_{2}cos[{\theta }_1\left(t\right)-{\theta }_2\left(t\right)]{\dot{\theta }}_1\left(t\right){\dot{\theta }}_2\left(t\right)\right]\hfill \\ & =& {\displaystyle \frac{1}{2}}{L}_1^2{\dot{\theta }}_1^2\left(t\right)\left[{\int }_{\left\{S_2\right\}}{\rho }_{2}A\left(s_2\right)d{s}_2\right]+{\displaystyle \frac{1}{2}}{\dot{\theta }}_2^2\left(t\right)\left[{\int }_{\left\{S_2\right\}}{\rho }_{2}A\left(s_2\right)s_2^{2}d{s}_2\right]+L_{1}cos[{\theta }_1\left(t\right)-{\theta }_2\left(t\right)]{\dot{\theta }}_1\left(t\right){\dot{\theta }}_2\left(t\right)\left[{\int }_{\left\{S_2\right\}}{\rho }_{2}A\left(s_2\right)s_{2}d{s}_2\right]\hfill \\ & =& {\displaystyle \frac{1}{2}}{m}_2{L}_1^2{\dot{\theta }}_1^2\left(t\right)+{\displaystyle \frac{1}{2}}{I}_2{\dot{\theta }}_2^2\left(t\right)+L_1{m}_2{C}_{2}cos[{\theta }_1\left(t\right)-{\theta }_2\left(t\right)]{\dot{\theta }}_1\left(t\right){\dot{\theta }}_2\left(t\right),\hfill \end{array}$$

(14)

where $\left\{S_2\right\}\equiv [0,L_2]$,

$$\begin{array}{ccc}\hfill m_2& =& {\rho }_2{\int }_{\left\{S_2\right\}}A\left(s_2\right)d{s}_2,\hfill \\ \hfill I_2& =& {\rho }_2{\int }_{\left\{S_2\right\}}{s}_2^{2}A\left(s_2\right)d{s}_2,\hfill \\ \hfill C_2& =& {\displaystyle \frac{{\int }_{\left\{S_2\right\}}{s}_{2}A\left(s_2\right)d{s}_2}{{\int }_{\left\{S_2\right\}}A\left(s_2\right)d{s}_2}}.\hfill \end{array}$$

(15)

Obviously, $C_2$ is the coordinate of the center of mass of the second segment measured from its beginning.

Summarizing the above, we conclude that, if the limb consists of just two segments, its total kinetic energy $K=K_1+K_2$ equals to

$$K={\displaystyle \frac{1}{2}}\left[I_1+m_2{L}_1^2\right]{\dot{\theta }}_1^2\left(t\right)+{\displaystyle \frac{1}{2}}{I}_2{\dot{\theta }}_2^2\left(t\right)+L_1{m}_2{C}_{2}cos[{\theta }_1\left(t\right)-{\theta }_2\left(t\right)]{\dot{\theta }}_1\left(t\right){\dot{\theta }}_2\left(t\right).$$

(16)

Now we turn to the potential energy of the system V supposing the bodies to be acted upon by the gravity field of the Earth. Obviously

$$V=V_1+V_2,$$

(17)

where $V_1$ and $V_2$ are the potential energies of the segments. From the geometry depicted in Figure 1 one has

$$\begin{array}{ccc}\hfill V_1& =& g{\int }_0^{L_1}{\rho }_1\left(s_1\right)A\left(s_1\right)y\left(s_1\right)d{s}_1\hfill \\ & =& -g{\rho }_1{\int }_0^{L_1}A\left(s_1\right)s_{1}cos\left({\theta }_1\right)d{s}_1=-g{\rho }_{1}cos\left({\theta }_1\right){\int }_0^{L_1}A\left(s_1\right)s_{1}d{s}_1\hfill \\ & =& -gcos\left({\theta }_1\right)m_1{C}_1,\hfill \end{array}$$

(18)

where

$$C_1={\displaystyle \frac{{\int }_0^{L_1}{\rho }_1{s}_{1}A\left(s_1\right)d{s}_1}{{\int }_0^{L_1}{\rho }_{1}A\left(s_1\right)d{s}_1}}={\displaystyle \frac{{\int }_0^{L_1}{s}_{1}A\left(s_1\right)d{s}_1}{{\int }_0^{L_1}A\left(s_1\right)d{s}_1}}.$$

(19)

In a similar way, for $V_2$ one derives

$$\begin{array}{ccc}\hfill V_2& =& g{\int }_0^{L_2}{\rho }_2\left(s_2\right)A\left(s_2\right)y\left(s_2\right)d{s}_2\hfill \\ & =& g{\int }_0^{L_2}{\rho }_2\left(s_2\right)A\left(s_2\right)[-L_{1}cos\left({\theta }_1\right)-s_{2}cos\left({\theta }_2\right)]d{s}_2\hfill \\ & =& -g{L}_{1}cos\left({\theta }_1\right){\int }_0^{L_2}{\rho }_2\left(s_2\right)A\left(s_2\right)d{s}_2-gcos\left({\theta }_2\right){\int }_0^{L_2}{\rho }_2\left(s_2\right)s_{2}A\left(s_2\right)d{s}_2\hfill \\ & =& -g{m}_2\left[L_{1}cos\left({\theta }_1\right)+cos\left({\theta }_2\right)C_2\right].\hfill \end{array}$$

(20)

From Equations (18) and (20) we arrive at

$$V=V_1+V_2=-gcos{\theta }_1(C_1{m}_1+L_1{m}_2)-gcos{\theta }_2{C}_2{m}_2.$$

(21)

Once we know the kinetic K and the potential V energies, we derive the Lagrangian $\mathcal{L}$ of the system

$$\begin{array}{ccc}\hfill \mathcal{L}=K-V& =& {\displaystyle \frac{1}{2}}\left[I_1+m_2{L}_1^2\right]{\dot{\theta }}_1^2+{\displaystyle \frac{1}{2}}{I}_2{\dot{\theta }}_2^2+L_1{m}_2{C}_{2}cos[{\theta }_1-{\theta }_2]{\dot{\theta }}_1{\dot{\theta }}_2\hfill \\ & & +g(C_1{m}_1+L_1{m}_2)cos{\theta }_1+g{C}_2{m}_{2}cos{\theta }_2.\hfill \end{array}$$

(22)

One can also determine the Hamiltonian $\mathcal{H}$ of the system

$$\begin{array}{ccc}\hfill \mathcal{H}=K+V& =& {\displaystyle \frac{1}{2}}\left[I_1+m_2{L}_1^2\right]{\dot{\theta }}_1^2+{\displaystyle \frac{1}{2}}{I}_2{\dot{\theta }}_2^2+L_1{m}_2{C}_{2}cos[{\theta }_1-{\theta }_2]{\dot{\theta }}_1{\dot{\theta }}_2\hfill \\ & & -g(C_1{m}_1+L_1{m}_2)cos{\theta }_1-g{C}_2{m}_{2}cos{\theta }_2.\hfill \end{array}$$

(23)

Since we did not introduce explicit dependence on the time t, the Hamiltonian of the system is conserved and equal to the stationary energy of the system. According to the standard Hamiltonian formulation, the state of a system is characterized not only by its positions (generalized coordinates) ${\theta }_1$ and ${\theta }_2$, but also by its momenta $p_1$ and $p_2$. Hence, the Hamiltonian of the system depicted in Figure 1 is:

$$\begin{array}{c}H=K+V={\displaystyle \frac{I_2}{2{\Delta }_0}}{p}_1^2+{\displaystyle \frac{I_1+m_2{L}_1^2}{2{\Delta }_0}}{p}_2^2-{\displaystyle \frac{m_2{L}_1{C}_2\cos \left({\theta }_1-{\theta }_2\right)}{{\Delta }_0}}{p}_1{p}_2\\ -g\left[\left(m_1{C}_1+m_2{L}_1\right)\cos {\theta }_1+m_2{C}_2\cos {\theta }_2\right],\end{array}$$

(24)

where

$$p_1={\displaystyle \frac{\partial K}{\partial {\dot{\theta }}_1}},p_2={\displaystyle \frac{\partial K}{\partial {\dot{\theta }}_2}},\mathrm{and}{\Delta }_0=I_1{I}_2+m_2{I}_2{L}_1^2-m_2^2{L}_1^2{C}_2^2{\cos }^2\left({\theta }_1-{\theta }_2\right).$$

(25)

3. The Equations of Motion

From Equation (22) one obtains the corresponding Euler-Lagrange equations, governing the behavior of the two segments of the considered limb. According to the general theory the equation are to be obtained from

$${\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial \mathcal{L}}{\partial \dot{{\theta }_i}}}-{\displaystyle \frac{\partial \mathcal{L}}{\partial {\theta }_i}}=0,i=1,2.$$

(26)

In this way we obtain a set of two coupled, second-order, nonlinear differential equations:

$$\begin{array}{c}gsin\left[{\theta }_1\right](C_1{m}_1+L_1{m}_2)+\left(I_1+L_1^2{m}_2\right){\ddot{\theta }}_1+C_2{L}_1{m}_2\left\{{\dot{\theta }}_1{\dot{\theta }}_{2}sin[{\theta }_1-{\theta }_2]\right.\hfill \\ \left.+{\ddot{\theta }}_{2}cos[{\theta }_1-{\theta }_2]-{\dot{\theta }}_2\left[{\dot{\theta }}_1-{\dot{\theta }}_2\right]sin[{\theta }_1-{\theta }_2]\right\}=0.\hfill \end{array}$$

(27)

and

$$\begin{array}{c}g{C}_2{m}_{2}sin\left[{\theta }_2\right]+I_2{\ddot{\theta }}_2\hfill \\ +C_2{L}_1{m}_2\left\{\ddot{{\theta }_1}cos[{\theta }_1-{\theta }_2]-{\dot{\theta }}_1^{2}sin[{\theta }_1-{\theta }_2]\right\}=0.\hfill \end{array}$$

(28)

Equation (27) can be rewritten in the simpler form

$$\begin{array}{c}gsin\left[{\theta }_1\right](C_1{m}_1+L_1{m}_2)+\left(I_1+L_1^2{m}_2\right){\ddot{\theta }}_1\hfill \\ +C_2{L}_1{m}_2\left\{{\ddot{\theta }}_{2}cos[{\theta }_1-{\theta }_2]+{\dot{\theta }}_2^{2}sin[{\theta }_1-{\theta }_2]\right\}=0.\hfill \end{array}$$

(29)

Equations (28) and (29) form the system of two coupled nonlinear differential equations. When modeling the process of walking one shall keep in mind that ${\theta }_2\le {\theta }_1$, when ${\theta }_1>0$.

In the form of a dynamical model this system is

$$\begin{array}{c}{\ddot{\theta }}_1=a_1-a_2{\dot{\theta }}_1^2-a_3{\dot{\theta }}_2^2,\hfill \\ {\ddot{\theta }}_2=-a_4+a_5{\dot{\theta }}_1^2+a_6{\dot{\theta }}_2^2,\hfill \end{array}$$

(30)

where

$$\begin{array}{c}{a}_1={\displaystyle \frac{g}{{\Delta }_0}}\left[m_2^2{C}_2^2{L}_{1}cos\left({\theta }_1-{\theta }_2\right)sin{\theta }_2-I_2\left(C_1{m}_1+L_1{m}_2\right)sin\left({\theta }_1\right)\right],\hfill \\ a_2={\displaystyle \frac{m_2^2{C}_2^2{L}_1^{2}cos\left({\theta }_1-{\theta }_2\right)sin\left({\theta }_1-{\theta }_2\right)}{{\Delta }_0}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{m_2^2{C}_2^2{L}_1^{2}sin\left[2\left({\theta }_1-{\theta }_2\right)\right]}{{\Delta }_0}},\hfill \\ a_3={\displaystyle \frac{m_2{C}_2{I}_2{L}_{1}sin\left({\theta }_1-{\theta }_2\right)}{{\Delta }_0}},\hfill \\ a_4=\begin{array}{c}\hfill -\end{array}{\displaystyle \frac{g{C}_2{m}_2}{{\Delta }_0}}\left\{L_1\left[C_1{m}_{1}sin\left({\theta }_1\right)cos({\theta }_1-{\theta }_2)+L_1{m}_{2}cos\left({\theta }_1\right)sin({\theta }_1-{\theta }_2)\right]-I_{1}sin\left({\theta }_2\right)\right\},\hfill \\ a_5={\displaystyle \frac{m_2{C}_2{L}_{1}sin\left({\theta }_1-{\theta }_2\right)\left(L_1^2{m}_2+I_1\right)}{{\Delta }_0}},\hfill \\ a_6={\displaystyle \frac{1}{2}}{\displaystyle \frac{m_2^2{C}_2^2{L}_1^{2}sin\left[2({\theta }_1-{\theta }_2\right)]}{{\Delta }_0}},\hfill \\ {\Delta }_0=I_1{I}_2+m_2{L}_1^2\left[I_2-m_2{C}_2^2{cos}^2\left({\theta }_1-{\theta }_2\right)\right].\hfill \end{array}$$

(31)

The data needed for the evaluation of the above quantities in case of human body are presented in Table 1 on the example of the average Bulgarian males. The data are taken from Ref. [20]. There, actually, the relative location of the center of mass of the body segments are given, but the transformation to physical lengths is trivial. In Ref. [20] one provides data for principal moments of inertia $I_{xx},I_{yy}$ and $I_{zz}$, where the coordinate system is centered at the center of mass of the corresponding segment; the z axes is along the segment, while x and y axes for a plane orthogonal to it. Due to the supposed there symmetry, $I_{xx}=I_{yy}$. One way of acting, in order to determine $I_1$ and $I_2$ is that they are equal to

$$I_i=(I_{xx}+I_{yy}-I_{zz})/2+m_i{c}_1^2,$$

(32)

where we are using the Steiner’s theorem. Another, direct approach is to recall that the thigh and shank are modeled there as frustum of cones and to derive that for the segment $"i"$ one has

$$I_i={\displaystyle \frac{1}{30}}\pi L_i^3{\rho }_i\left(R_{i,1}^2+3R_{i,1}{R}_{i,2}+6R_{i,2}^2\right),i=1,2.$$

(33)

Let us perform dimensional analysis for the quantities shown in Equation (31). One has that $\left[{\Delta }_0\right]={(kg\times c{m}^2)}^2$. Then $\left[a_2\right]=\mathcal{O}\left(1\right)$, $\left[a_3\right]=\mathcal{O}\left(1\right)$, $\left[a_5\right]=\mathcal{O}\left(1\right)$, and $\left[a_6\right]=\mathcal{O}\left(1\right)$, i.e., these constants are dimensionless. Finally, $\left[a_1\right]=s^{-2}$ and $\left[a_4\right]=s^{-2}$. We have also to take into account that $g=9.8m/s^2=980cm/s^2$.

4. Biomechanical Data

In order to apply Equation (33), we have to take into account that, according to Ref. [21], for the thigh $R_{1,1}=9.5cm,R_{1,2}=5.5cm,{\rho }_1=1062kg/m^3$ while for the shank $R_{2,1}=5.5cm,R_{2,2}=4.0cm,$${\rho }_2=1088kg/m^3$.

Using the values shown in Table 1, one obtains (in standard SI units)

$${\Delta }_0=0.168-0.078{cos}^2({\theta }_1-{\theta }_2).$$

(34)

and

$$\begin{array}{c}{a}_1={\displaystyle \frac{56.522sin\left({\theta }_1\right)-19.216cos({\theta }_1-{\theta }_2)sin\left({\theta }_2\right)}{{cos}^2({\theta }_1-{\theta }_2)-2.147}},\hfill \\ a_2=-{\displaystyle \frac{sin\left[2({\theta }_1-{\theta }_2)\right]}{2\left[{cos}^2({\theta }_1-{\theta }_2)-2.147\right]}},\hfill \\ a_3=-{\displaystyle \frac{0.402sin({\theta }_1-{\theta }_2)}{{cos}^2({\theta }_1-{\theta }_2)-2.147}},\hfill \\ a_4={\displaystyle \frac{81.416sin\left({\theta }_1\right)cos({\theta }_1-{\theta }_2)+59.036cos\left({\theta }_1\right)sin({\theta }_1-{\theta }_2)-43.479sin\left({\theta }_2\right)}{{cos}^2({\theta }_1-{\theta }_2)-2.1473}},\hfill \\ a_5=-{\displaystyle \frac{5.335sin({\theta }_1-{\theta }_2)}{{cos}^2({\theta }_1-{\theta }_2)-2.147}},\hfill \\ a_6=-{\displaystyle \frac{sin\left[2({\theta }_1-{\theta }_2)\right]}{2\left[{cos}^2({\theta }_1-{\theta }_2)-2.147\right]}}.\hfill \end{array}$$

(35)

Note that $a_2=a_6$ for any ${\theta }_1$ and ${\theta }_2$.

The human gait is usually considered as a cycle decomposed into eight phases, see Figure 2. The phases can be categorized in various ways. Here we do that by using the ankles for the hip, knee, and ankle joint, for each of the eight phases of the human gait cycle. The angles are given, i.e. their typical values, in Table 2.

5. Numerical Results

The small free oscillations of a conservative linear system in the phase portrait have the form of closed trajectories in the neighborhood of a stable steady state. It is well-known that amplitudes of oscillations for linear systems (for nonlinear conservative ones also) depend on the initial conditions. On the other hand, its period is permanent and independent of initial conditions and the energy level of the system. The periodic oscillations in conservative systems are always orbitally stable, i.e. for very small perturbations in the initial conditions a periodical motion passes into another one as these two motions are very close. Solving numerically Equation (30) we obtain:

1)

For small initial angles and speeds we obtain the following solutions for ${\theta }_1$ and ${\theta }_2$:

2)

For intermediate initial angles and speeds we obtain the following solutions for ${\theta }_1$ and ${\theta }_2$:

Figure 3. (a) The figure shows the behavior of ${\theta }_1$ and ${\theta }_2$ for $t\in [0,20]$ and a specific choice of the initial conditions: ${\theta }_1\left(0\right)=5\pi /100$, ${\theta }^{\prime }\left[0\right]=0.1$, ${\theta }_2\left[0\right]=5\pi /100$, ${\theta }_2^{\prime }\left[0\right]=0.1$.(b) Phase diagram of the system under the same boundary conditions.

Figure 3. (a) The figure shows the behavior of ${\theta }_1$ and ${\theta }_2$ for $t\in [0,20]$ and a specific choice of the initial conditions: ${\theta }_1\left(0\right)=5\pi /100$, ${\theta }^{\prime }\left[0\right]=0.1$, ${\theta }_2\left[0\right]=5\pi /100$, ${\theta }_2^{\prime }\left[0\right]=0.1$.(b) Phase diagram of the system under the same boundary conditions.

Figure 4. (a) The figure shows the behavior of ${\theta }_1$ and ${\theta }_2$ for $t\in [0,20]$ and a specific choice of the initial conditions: ${\theta }_1\left(0\right)=5\pi /36$, ${\dot{\theta }}_1\left[0\right]=0.1$, ${\theta }_2\left[0\right]=5\pi /36$, ${\dot{\theta }}_2\left[0\right]=0.1$. (b) Phase diagram of the system under the same boundary conditions. It demonstrates that the movement is not purely periodic but there is also an element of chaos in it. We observe that the spread of the phase diagram is about three times larger than for the small angles. i.e., larger the initial angles larger is the spread of the phase diagram.

Figure 4. (a) The figure shows the behavior of ${\theta }_1$ and ${\theta }_2$ for $t\in [0,20]$ and a specific choice of the initial conditions: ${\theta }_1\left(0\right)=5\pi /36$, ${\dot{\theta }}_1\left[0\right]=0.1$, ${\theta }_2\left[0\right]=5\pi /36$, ${\dot{\theta }}_2\left[0\right]=0.1$. (b) Phase diagram of the system under the same boundary conditions. It demonstrates that the movement is not purely periodic but there is also an element of chaos in it. We observe that the spread of the phase diagram is about three times larger than for the small angles. i.e., larger the initial angles larger is the spread of the phase diagram.

6. Small Angles Approximation

In the current section we solve both analytically and numerically Equation (28) and Equation (29) form the system of two coupled nonlinear differential equations in the small angles approximation. Performing small angles expansion in these equation and retaining only linear terms, we obtain

$$\ddot{{\theta }_1}\left(t\right)=g{\displaystyle \frac{I_2\left[C_1{m}_1+L_1{m}_2\right]{\theta }_1\left(t\right)-C_2^2{L}_1{m}_2{\theta }_2\left(t\right)}{L_1^2{m}_2\left[C_2^2{m}_2-I_2\right]-I_1{I}_2}}$$

(36)

and

$$\ddot{{\theta }_2}\left(t\right)=g{m}_2{C}_2{\displaystyle \frac{{\theta }_2\left(t\right)\left[I_1+L_1^2{m}_2\right]-L_1{\theta }_1\left(t\right)\left[C_1{m}_1+L_1{m}_2\right]}{L_1^2{m}_2\left[C_2^2{m}_2-I_2\right]-I_1{I}_2}}.$$

(37)

Obviously, these two equations describe two linearly coupled oscillators. Looking for a solution in the form

$${\theta }_1\left(t\right)=A_{1}exp\left[i\omega t\right],and{\theta }_2\left(t\right)=A_{2}exp\left[i\omega t\right]=A_{1}rexp\left[i\omega t\right],withr=A_2/A_1,$$

(38)

we obtain

$$\omega =\pm \sqrt{g}{\displaystyle \frac{\sqrt{C_1{m}_1+L_1{m}_2}}{\sqrt{L_1{m}_2(C_2+L_1)+I_1}}},$$

(39)

and

$$r={\displaystyle \frac{(C_1{m}_1+L_1{m}_2)(C_2{L}_1{m}_2+I_2)}{C_2{m}_2\left[L_1{m}_2(C_2+L_1)+I_1\right]}}={\displaystyle \frac{{\omega }^2}{g}}{\displaystyle \frac{C_2{L}_1{m}_2+I_2}{C_2{m}_2}}.$$

(40)

Plugging in Equation (39) and Equation (40) the values of all quantities involved, we obtain

$$\omega =\pm 4.709andr=1.618.$$

(41)

The behavior of ${\theta }_1\left(t\right)$ and ${\theta }_2\left(t\right)$ is shown in Figure 5. One observes for very small angles periodic behavior. The changes in the behaviors of ${\theta }_1\left(t\right)$ and ${\theta }_2\left(t\right)$ are in phase with each other.

7. Comparison of Our Model Results to Real Walking

In this section, we will briefly discuss how our results compare to real walking. Walking is a fundamental human activity, essential for daily life, work, and social interaction. Its significance extends to well-being, enabling shopping, commuting, and maintaining connections with others. The study of human motion, with a particular focus on gait, has a long history, dating back to ancient times [23]. This research encompasses diverse areas, including neuromuscular disorders, joint replacement, sports injuries, and assistive devices [24,25]. Furthermore, the increasing interest in humanoid and wearable robots [26] underscores the continuing importance of understanding human locomotion. The results of our model of human walking indicate that for small initial angles (periodic movement), the period $T=2\pi /\omega \simeq 1.33$ seconds. One forward swing of the leg corresponds to a half cycle, which takes $T_{1/2}=0.67$ s. This is quite close to the number reported in Ref. [27,28], according to which at slow walk $T_{1/2}\simeq 0.60$ s. $T_{1/2}=0.67$ s corresponds to about 2400 step/hour. If we assume the average step length to be approximately $0.6$ m, we obtain a speed of approximately $1.44$ km/h. For the movements shown in Figure 3, $T_{1/2}=0.81$ s, while for Figure 4, $T_{1/2}=0.44$ s. These data correspond to $2.67$ km/h and $5.9$ km/h, respectively. The corresponding step lengths are, then, $2Lsin\left(\theta \right)$. Having that for $T_{1/2}=0.81$ s one has $\theta =2\times 5\pi /100\simeq {18}^{\circ }$ we obtain for the step length $(L_1+L_2)sin\left({18}^{\circ }\right)=0.27$ m and $(L_1+L_2)sin\left({50}^{\circ }\right)=0.68$ m for $T_{1/2}=0.44$ s. According to Ref. [28] [p. 127], see also Refs. [29,30], the speed for a slow walk is about $0.5m/s=1.8$ km/h. There one takes a step length of $0.3$ m, the time for a step $T_{1/2}\simeq 0.60$ s. For a fast walk, the speed is estimated about $2.0m/s=7.2$ km/h, and the step length of about, $0.7$ m, so the step time is about $0.35$ s. As we see all these numbers are very close to the numbers we obtained within our model.

8. Conclusions

The mathematical model described in the current article is based on a double pendulum representation of a limb (e.g., thigh + shank) with realistic mass–inertia parameters. The model uses measured anthropometric data (those for average Bulgarian males in the paper). Naturally, it can be tailored to specific populations or individuals, making it a bridge between abstract physics and personalized bio-mechanical solutions. Thus, although it’s a theoretical construct, it has clear real-world applications across several domains: i) gait analysis and rehabilitation; The model can simulate normal and pathological walking patterns, helping clinicians design targeted physiotherapy or post-surgery recovery plans. ii) prosthetics and orthotics design; By adjusting parameters (mass, length, inertia) to match a patient’s anatomy, engineers can optimize artificial limbs for more natural movement. iii) injury prevention; Understanding joint loads and motion dynamics can guide ergonomic recommendations for athletes and workers. Many of the biomechanical problems associated with the evolution of total knee and total hip replacements have been evaluated in terms of walking mechanics. Mathematical modeling of human walking is an approach to estimate the human joint kinematics and dynamics during walking. More specifically, mathematical modeling and analysis help us distinguish between healthy patterns and the presence of any pathological alterations. Summarizing, we conclude that in the current paper, a realistic mathematical model of human limb movement is proposed using a double pendulum that moves in a plane under the influence of gravity. The model simulates the complete limbs of the human body and can easily be scaled to study different groups (children, women and et al.) of human populations. The model is dynamic, and it has a realistic representation, such as a slow/fast healthy walk. The model can be, obviously, improved by introducing torques acting on the segments of the pendulum. This will lead to source terms on the right-hand side of Equation (28) and Equation (29). The problem will become even more interesting and complicated if these torques depend on the time [31]. Understanding torques enables inferences about the muscles that most likely contract at certain moments to cause the leg to swing in tandem with gravity [24]. In addition, one can also consider a double inverted pendulum [32,33,34,35]. As already stated above, the key characteristic of the model is its non-linearity and chaos. Small changes in initial conditions can lead to different outcomes. This makes it a powerful model for studying the complex control required for human movement [9,31,36]. The double pendulum is a cornerstone model for understanding how the human body moves, as it perfectly illustrates the mechanics of multi-segment chains. Much of our gait (walking/running) can be explained by passive pendulum-like motion, which is incredibly energy-efficient. The model helps us understand this: indeed, during the swing phase of walking or running, the leg behaves remarkably like a double pendulum.