A Biologically Motivated Finite Difference Approach for Simulating Singularly Perturbed Vertical Motion in Human Gait

TL;DR

This paper introduces a novel finite difference numerical method, inspired by biology, for accurately simulating the vertical motion in human gait models involving singularly perturbed differential equations, addressing boundary layer challenges.

Contribution

The study develops a domain decomposition and time-rescaling based finite difference scheme that effectively solves singularly perturbed gait equations with second-order accuracy.

Findings

The method achieves second-order spatial accuracy.

Numerical experiments confirm stability and efficiency.

Framework applicable to complex biomechanical models.

Abstract

In this study, we present a simulation-based numerical method for solving a class of singularly perturbed second-order differential equations that come from a simplified biologically motivated model of human gait. Important physical factors such as gravity, damping, and leg stiffness are included in the model, which also depicts the vertical motion of the center of mass of the body during walking or running. Most of the time, standard numerical methods are ineffective in resolving boundary layer behavior that occurs due to the small perturbation parameter in the governing equation. We use a domain decomposition technique to divide the problem domain into inner and outer regions to tackle this difficulty. The boundary layer resolves the steep gradients. We applied a time-rescaling transformation to the inner region. Each subdomain is discretized, and the resulting tridiagonal systems are efficiently solved using the Thomas algorithm within the mixed finite difference framework. A detailed convergence analysis demonstrates second-order accuracy in space. The numerical results validate the proposed scheme's accuracy, stability, and efficiency through experiments based on modified human gait models. The framework serves as a fundamental tool for biomechanical simulation. The modeling is a foundation for future research, incorporating nonlinearities, time delays, and real-world scenarios data on how people walk.

Tables

Table 1. Table 1 : Data Table

	0.0009	0.009	0.001	0.01
0.0000	4.0000	4.0000	4.0000	4.0000
0.0200	4.7329	4.8068	4.7185	4.7185
0.0392	4.7104	4.7164	4.7079	4.7079
0.0584	4.6264	4.6268	4.6261	4.6261
0.0776	4.5389	4.5388	4.5389	4.5389
0.0968	4.4526	4.4525	4.4526	4.4526
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
0.3432	3.4801	3.4802	3.4801	3.4802
0.3530	3.4461	3.4462	3.4461	3.4463
0.3628	3.4125	3.4126	3.4125	3.4126
0.3726	3.3793	3.3793	3.3793	3.3794
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
0.9704	1.8587	1.8590	1.8587	1.8592
0.9802	1.8405	1.8433	1.8405	1.8443
0.9900	1.8229	1.8455	1.8229	1.8496
1.0000	2.0000	2.0000	2.0000	2.0000

Table 2. Table 2 : Data Table

	0.0009	0.009	0.001	0.01
0	1.0000	1.0000	1.0000	1.0000
0.0100	0.9806	0.9492	0.9805	0.9432
0.0198	0.4956	0.5065	0.4956	0.5086
0.0296	0.3314	0.3338	0.3314	0.3344
0.0394	0.2490	0.2497	0.2490	0.2499
0.0492	0.1994	0.1997	0.1994	0.1997
0.0590	0.1663	0.1664	0.1663	0.1664
0.0688	0.1426	0.1427	0.1426	0.1427
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
0.2844	0.0345	0.0345	0.0345	0.0345
0.2942	0.0333	0.0333	0.0333	0.0333
0.3040	0.0323	0.0323	0.0323	0.0323
0.3138	0.0313	0.0313	0.0313	0.0313
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
0.9606	0.0102	0.0102	0.0102	0.0102
0.9704	0.0101	0.0101	0.0101	0.0101
0.9802	0.0100	0.0100	0.0100	0.0100
0.9900	0.0099	0.0099	0.0099	0.0099
1.0000	0.0100	0.0100	0.0100	0.0100

Table 3. Table 3 : Data Table

	0.0009	0.009	0.001	0.01
0	9.6000	9.6000	9.6000	9.6000
0.0100	9.6421	6.6457	9.6259	6.3281
0.0198	9.6171	8.5673	9.6168	8.3379
0.0296	9.5239	9.1661	9.5239	9.0424
0.0394	9.4310	9.3100	9.4310	9.2510
0.0492	9.3390	9.2986	9.3390	9.2724
0.0590	9.2480	9.2348	9.2480	9.2237
0.0688	9.1578	9.1538	9.1578	9.1493
⋮	⋮	⋮	⋮	⋮
0.6078	5.3420	5.3422	5.3420	5.3423
0.6176	5.2899	5.2901	5.2899	5.2902
0.6274	5.2383	5.2385	5.2383	5.2386
0.6372	5.1872	5.1875	5.1872	5.1875
⋮	⋮	⋮	⋮	⋮
0.9312	3.8659	3.8660	3.8659	3.8661
0.9410	3.8282	3.8283	3.8282	3.8283
0.9508	3.7909	3.7908	3.7909	3.7907
0.9606	3.7539	3.7531	3.7539	3.7525
0.9704	3.7173	3.7128	3.7173	3.7105
0.9802	3.6811	3.6578	3.6811	3.6504
0.9900	3.6433	3.5263	3.6429	3.5088
1.0000	3.0000	3.0000	3.0000	3.0000

Equations

$$z_{tt}+\epsilon u(t)z_t+v(t)z=-g,$$

$$z(0)=f_1,z(t_f)=f_2.$$

$$T=\epsilon t\Rightarrow t=\frac{T}{\epsilon }.$$

$$z_t=\epsilon Z_T,z_{tt}={\epsilon }^2{Z}_{TT}.$$

$${\epsilon }^2{Z}_{TT}(T)+{\epsilon }^{2}U(T)Z_T(T)+V(T)Z(T)=-g,$$

$$Z(0)=f_1,Z(\epsilon t_f)=f_2.$$

$$Z_{TT}(T_i)$$

$$Z_T(T_i)$$

$${\epsilon }^2\left(\frac{Z_{i+1}-2Z_i+Z_{i-1}}{k^2}\right)+{\epsilon }^2{U}_i\left(\frac{Z_{i+1}-Z_{i-1}}{2k}\right)+V_i{Z}_i=-g,$$

$$C_i{Z}_{i-1}+A_i{Z}_i+B_i{Z}_{i+1}=-g,$$

$$C_i$$

$$A_i$$

$$B_i$$

$$z_{tt}(t_i)$$

$$z_t(t_i)$$

$$\frac{z_{i+1}-2z_i+z_{i-1}}{k^2}+\epsilon u_i\cdot \frac{z_{i+1}-z_{i-1}}{2k}+v_i{z}_i=-g,$$

$${\tilde{C}}_i{z}_{i-1}+{\tilde{A}}_i{z}_i+{\tilde{B}}_i{z}_{i+1}=-g,$$

$${\tilde{C}}_i$$

$${\tilde{A}}_i$$

$${\tilde{B}}_i$$

$$a_i{x}_{i-1}+b_i{x}_i+c_i{x}_{i+1}=d_i,i=1,2,\dots ,N-1,$$

$${\epsilon }^2{Z}_{TT}(T_i)+{\epsilon }^2{U}_i{Z}_T(T_i)+V_i{Z}_i=-g,$$

$$Z_{TT}(T_i)$$

$$Z_T(T_i)$$

$${\tau }_i(k)={\epsilon }^2\cdot \frac{k^2}{12}{Z}^{(4)}(T_i)+{\epsilon }^2{U}_i\cdot \frac{k^2}{6}{Z}^{(3)}(T_i)+O(k^4).$$

$${\tau }_i(k)=O(k^2).$$

$$PZ=R,$$

$$P=\left[\begin{array}{ccccc}{A}_1& B_1& 0& \cdots & 0\\ C_2& A_2& B_2& \cdots & 0\\ 0& C_3& A_3& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & B_{N-2}\\ 0& \cdots & 0& C_{N-1}& A_{N-1}\end{array}\right],$$

$$C_i$$

$$A_i$$

$$B_i$$

$$|A_i|>|C_i|+|B_i|.$$

$$\left|V_i-\frac{2{\epsilon }^2}{k^2}\right|>{\epsilon }^2\left(\frac{1}{k^2}+\frac{U_i}{2k}\right)+{\epsilon }^2\left(\frac{1}{k^2}-\frac{U_i}{2k}\right)=\frac{2{\epsilon }^2}{k^2}.$$

$$V_i>\frac{4{\epsilon }^2}{k^2}.$$

Full text

A Biologically Motivated Finite Difference Approach for Simulating Singularly Perturbed Vertical Motion in Human Gait

Shubhangini Gupta

email id: shubhangini_p240187ma@nitc.ac.in

Sourav Banerjee

email id: souravsbbanerjee@gmail.com

Tamal Pramanick

email id: tamal@nitc.ac.in, ∗Corresponding author

Abstract. In this study, we present a simulation-based numerical method for solving a class of singularly perturbed second-order differential equations that come from a simplified biologically motivated model of human gait. Important physical factors such as gravity, damping, and leg stiffness are included in the model, which also depicts the vertical motion of the center of mass of the body during walking or running. Most of the time, standard numerical methods are ineffective in resolving boundary layer behavior that occurs due to the small perturbation parameter in the governing equation. We use a domain decomposition technique to divide the problem domain into inner and outer regions to tackle this difficulty. The boundary layer resolves the steep gradients. We applied a time-rescaling transformation to the inner region. Each subdomain is discretized, and the resulting tridiagonal systems are efficiently solved using the Thomas algorithm within the mixed finite difference framework. A detailed convergence analysis demonstrates second-order accuracy in space. The numerical results validate the proposed scheme’s accuracy, stability, and efficiency through experiments based on modified human gait models. The framework serves as a fundamental tool for biomechanical simulation. The modeling is a foundation for future research, incorporating nonlinearities, time delays, and real-world scenarios data on how people walk.

Key words. Human gait model, Mixed finite difference method, Convergence, Singularly perturbation, Vertical motion

**Mathematics subject classifications. 65L11, 65L70, 65L80 **

1 Introduction

Mathematical modeling of human gait has emerged as a critical area of research due to its wide-ranging applications in clinical gait analysis, rehabilitation engineering, prosthetic limb development, humanoid robotics, and sports science. The most important aspect of this modeling is the vertical displacement of the body’s center of mass (COM), which is crucial for the energy efficiency and mechanical stability of locomotion. Simple mechanical models, like spring-mass systems, inverted pendulum models, and variations on these, have been used to successfully model key aspects of gait mechanics, such as impact forces, energy transfer, leg compliance, and ground response forces, see e.g., [8, 6].

Often using minor physical quantities as stiffness ratios, mass inertia, or damping coefficients, biomechanical models produce governing equations (cf. [2, 4]). These small parameters cause singular perturbations that occurs due to narrow areas near domain boundaries where the solution experiences fast variation to develop, which is named as boundary layers. Standard numerical approaches, which frequently fail to resolve the steep gradients precisely and may suffer from instability or loss of precision, find difficulty in these layers. When one models the motion of the center of mass (COM) while walking or running, where the dynamics require a complicated interaction between position, velocity and acceleration. One can readily model such biological processes as periodic dynamic systems. Often found near the ends of the domain, the resulting singularly perturbed differential or differential-difference equations are notable for their abrupt boundary layer behavior. Precisely capturing these characteristics calls for specific numerical techniques since ordinary discretization approaches often fail to preserve stability or provide sufficient accuracy in the presence of these layers [14, 5, 7].

In this study, we introduce a novel mixed finite difference approach to address the solution of a singularly perturbed differential-difference model that characterizes the intricacies of human gait. Our formulation is inspired by reduced-order nonlinear oscillatory models in biomechanics (see [10, 15]) and integrates structural insights from research on singularly perturbed delay differential equations (cf. [16, 12, 13, 9]). Our method diverges from shift operator-based or solely fitted schemes by using a terminal boundary point to partition the domain into separate inner and outer subregion. This enables us to convert the initial problem into an asymptotically comparable singular perturbation formulation. The use of time rescaling in the boundary layer area is a key new feature of our framework. This higher time resolution is especially useful for apps like tracking athletic motion, because it lets us get a clear picture of how quickly the solution changes without having to refine the mesh too much. In both subdomains, central difference approximations are used, and the Thomas algorithm is used to quickly solve the discrete systems that are left over, see e.g. [18, 19, 8].

The vertical motion of the human body’s center of mass during locomotion can be described by a second-order ordinary differential equation that incorporates key biomechanical properties such as leg stiffness and damping. Taking into account the small perturbation parameter $\varepsilon$, the governing equation is given by:

[TABLE]

where $u(t)=\varepsilon\mu(t)$ represents a damping function and $v(t)$ denotes the stiffness. The boundary conditions are $z(0)=f_{1}$, $z(t_{f})=f_{2}$. In the presence of the small parameter $\varepsilon$, the solution exhibits a boundary layer near the domain endpoints.

We use the suggested technique over several terminal point selections and validate its performance on several model instances, including those with perturbation- delay and advancement parameters (cf. [5]). This phenomena guarantees the computational robustness [17]. We show that the method is well-suited for stiff and layer-dominated regimes generally seen in biomechanical systems since it achieves second-order precision and uniform convergence with respect to the perturbation parameter $\varepsilon$.

To the best of authors knowledge, this is the first work to introduce a hybrid finite difference framework that combines domain decomposition, boundary layer scaling, and structural adaptation for solving singularly perturbed models in human gait. Numerical experiments are presented to validate the theoretical convergence behavior and demonstrate the method’s versatility and effectiveness. The work is arranged as follows. The mathematical model for the human gait is presented in the next section, namely Section 2, together with a time-rescaling technique and domain decomposition analysis. The discretization using the mixed finite difference approach and the algorithmic implementation is covered in Section 3, and Section 4 offers a theoretical convergence study. Numerical experiments verifying the accuracy and homogeneity of convergence of the technique are found in Section 5, and section 6 ends the work with a summary and recommendations for next work including extensions to nonlinear models with changing coefficients and control inputs.

2 Mathematical Modeling and Problem Formulation

The vertical motion of the human body’s center of mass during locomotion can be modeled using a second-order singularly perturbed differential equation incorporating biomechanical effects such as damping and leg stiffness. The governing equation is given by:

[TABLE]

where $\varepsilon$ is a small positive perturbation parameter ($0<\varepsilon\ll 1$), $u(t)=\varepsilon\mu(t)$ is the damping function, $v(t)$ represents leg stiffness, $g$ is the gravitational constant, $z(t)$ is the vertical displacement of the COM. The boundary conditions are:

[TABLE]

2.1 Time Rescaling in the Boundary Layer

To resolve the steep gradients near the boundary at $t=0$, we introduce a stretched variable:

[TABLE]

Let $Z(T)=z\left(\frac{T}{\varepsilon}\right)$, then the derivatives transform as:

[TABLE]

Substituting into equation (2.1), we obtain the transformed equation in the inner region:

[TABLE]

where $U(T)=\mu(T/\varepsilon)$ and $V(T)=v(T/\varepsilon)$, with modified boundary conditions:

[TABLE]

To accurately solve the singularly perturbed second-order differential equation representing human gait dynamics, we employ a mixed finite difference method based on domain decomposition and time rescaling. This method divides the domain into two regions, inner and outer regions. On the inner region a steep gradients (boundary layers) are expected to occur and a stretched variable is used to resolve this. On the other hand, the outer region, where the solution varies smoothly, standard discretization is sufficient.

3 Domain Decomposition using Mixed Finite Difference Scheme

We divide the computational domain $[0,t_{f}]$ into two subdomains to address boundary layer behavior:

• •

Inner region: $T\in[0,T_{p}]$, where $T_{p}=\varepsilon t_{p}$ captures the boundary layer, and

• •

Outer region: $t\in[t_{p},t_{f}]$, where the solution varies smoothly.

3.1 Discretization in the Inner Region

Let $T_{i}=ik$ be a uniform discretization of the interval $[0,T_{p}]$, where $k$ is the mesh step size. Denoting $Z(T_{i})=Z_{i}$, the central difference approximations are:

[TABLE]

Substituting these into equation (2.2), we obtain the finite difference scheme:

[TABLE]

which can be rearranged as:

[TABLE]

where:

[TABLE]

3.2 Discretization in the Outer Region

Similarly, for the outer region $t\in[t_{p},t_{f}]$ with uniform step size $k$ and grid points $t_{i}=t_{p}+ik$, denoting $z(t_{i})=z_{i}$, we use:

[TABLE]

Substituting into the original equation (2.1) gives:

[TABLE]

which we write as:

[TABLE]

where:

[TABLE]

These two systems form the backbone of the mixed finite difference scheme and will be solved using the Thomas algorithm in the subsequent section.

3.3 Thomas Algorithm for Tridiagonal Systems

The resulting linear systems in both the inner and outer regions are tridiagonal and can be efficiently solved using the Thomas algorithm (a simplified form of Gaussian elimination for tridiagonal matrices), cf. [1, 3]. This allows for a fast and memory-efficient solution with computational complexity $O(N)$.

Let the tridiagonal system be of the form:

[TABLE]

with boundary values $x_{0}$ and $x_{N}$ specified.

The Thomas algorithm proceeds in two steps:

Forward elimination: modify the coefficients to eliminate lower diagonal elements. 2. 2.

Back substitution: solve for $x_{i}$ starting from $i=N-1$ down to $1$.

This method is applied independently to the inner and outer systems and yields the approximate solution across the entire domain, ensuring high accuracy near the boundary layer and smooth regions.

4 Convergence and Stability Analysis

To validate the accuracy and robustness of the proposed mixed finite difference scheme, we analyze both the local truncation error and the global error behavior. This section provides a mathematical justification of the second-order convergence and stability of the method.

4.1 Local Truncation Error

We consider the finite difference scheme for the inner region:

[TABLE]

where central difference approximations are used:

[TABLE]

Substituting into the discrete scheme, we obtain the local truncation error $\tau_{i}(k)$:

[TABLE]

Assuming the derivatives $Z^{(3)}$ and $Z^{(4)}$ are bounded, we conclude:

[TABLE]

4.2 Matrix Form and Stability

The discretized system can be expressed in matrix-vector form:

[TABLE]

where $P\in\mathbb{R}^{(N-1)\times(N-1)}$ is tridiagonal, $Z$ is the vector of unknowns, and $R=(-g,-g,\dots,-g)^{T}$. The matrix $P$ is given by:

[TABLE]

where

[TABLE]

To ensure numerical stability, we require that the matrix $P$ is strictly diagonally dominant, i.e.,

[TABLE]

Substituting the expressions, we get:

[TABLE]

This is satisfied if:

[TABLE]

This condition is physically reasonable in biomechanical models where leg stiffness $V_{i}$ is always positive and bounded below. Hence, the matrix $P$ is invertible and the scheme is stable.

4.3 Norm Bound For The Global Error Estimate

Let $\bar{Z}$ be the exact solution vector at grid points, and $Z$ be the numerical solution. Define the global error vector:

[TABLE]

From the exact and numerical systems:

[TABLE]

where $\tau(k)$ is the local truncation error vector. Taking norms:

[TABLE]

Since $P$ is strictly diagonally dominant (as shown in equation (4.2)), the inverse $\|P^{-1}\|$ is bounded independently of $k$, and from (4.1), $\|\tau(k)\|=O(k^{2})$. Thus:

[TABLE]

This confirms that the proposed numerical method achieves second-order convergence with respect to the step size.

4.4 Final Estimate on the Error

Given the local truncation error is of the form

[TABLE]

we conclude

[TABLE]

Hence, the global error satisfies

[TABLE]

where $C=\|P^{-1}\|$ is bounded and independent of $k$.

The proposed finite difference method is second-order accurate in time. The matrix formulation ensures numerical stability due to diagonal dominance and monotonicity. This analysis supports the numerical observations and validates the convergence behavior shown in the examples.

5 Numerical Illustrations

In this section, we demonstrate the accuracy and robustness of the proposed mixed finite difference scheme by solving a series of singularly perturbed test problems that simulate vertical motion in human gait. Each problem contains an inner boundary layer, and we observe the method’s performance as the perturbation parameter $\varepsilon$ varies. The computations were performed for several values of $\varepsilon$, and the maximum absolute errors were recorded over uniform grids of increasing resolution. All systems were solved using the Thomas algorithm for tridiagonal systems [21, 22]. The results validate the theoretical second-order convergence established earlier.

Example 5.1** (Exponential Stiffness and Damping).**

Consider the singularly perturbed equation:

[TABLE]

This model includes both exponentially growing damping and exponentially decaying stiffness. The solution of (5.1) exhibits steep variation near $T=0$. The following table (cf. Table 1) describes the solution for various values of $\varepsilon$ with respect to different each time level and the solution for various values of $\varepsilon$ is depicted in Figure 1.

Example 5.2** (Linear Damping and Quadratic Stiffness).**

We now solve the following model:

[TABLE]

This problem exhibits strong stiffness for large $T$ values and is ideal for testing the robustness of the method under rapidly changing coefficients. The detailed numerical solutions are described in Table 2 and Figure 2.

Example 5.3** (Pure Stiffness-Driven System).**

The final test problem is:

[TABLE]

Here, damping is absent, and only the stiffness $e^{T}$ governs the system dynamics. The data table and the figure is given in the following (see Table 3 and Figure 3, respectively).

5.1 Observation and Discussion

Across all three examples we observe the following:

• •

The method captures boundary layer behavior efficiently without oscillations.

• •

Uniform second-order convergence is observed as $\varepsilon\to 0$.

• •

Results validate both the theoretical convergence rate and the scheme’s robustness.

6 Applications and Discussion

The presented numerical framework provides a powerful tool for analyzing singularly perturbed models of human gait. Its robustness in resolving boundary layers, coupled with second-order convergence and computational efficiency, makes it suitable for a wide range of practical and interdisciplinary applications.

6.1 Biomechanical Relevance of Stiffness Modeling

Leg stiffness is a key biomechanical property that governs energy efficiency, stability, and shock absorption during walking and running. It is regulated by the joint mechanics and muscle-tendon complex, and adapts dynamically to different locomotion tasks.

Our model accounts for time-varying stiffness and damping effects, both of which are biologically realistic. In particular:

• •

Positive stiffness is common in passive tissues and muscles under stretching or isometric contraction.

• •

Negative stiffness can occur transiently during rapid concentric contractions (e.g., hopping or sprinting), due to the nonlinear force–length relationship of active muscle fibers.

Despite local occurrences of negative stiffness at the contractile component level, the total leg stiffness during typical gait remains globally positive, ensuring mechanical stability. This complexity underscores the need for numerical methods that can handle rapid variations– captured here using domain decomposition and time rescaling.

6.2 Integration with Biomechanics and Robotics

The proposed method can serve as a core component for simulation platforms in the following areas:

• •

Clinical biomechanics: Analysis of gait disorders, planning of rehabilitation therapies, and prosthetic design.

• •

Sports science: Performance monitoring, fatigue analysis, and injury prevention in athletes.

• •

Humanoid robotics: Development of energy-efficient biped locomotion using biologically inspired stiffness control.

• •

Musculoskeletal modeling: Integration with finite element models for estimating muscle forces, joint torques, and strain profiles.

6.3 Computational Advantages

In addition to its biomechanical relevance, the method offers several computational benefits:

• •

Efficient solution of tridiagonal systems using the Thomas algorithm.

• •

Uniform accuracy across stiff regimes, avoiding oscillations common in standard finite difference schemes.

• •

Scalability through domain decomposition, which can be optimized to reduce computation time by approximately 40%.

6.4 Toward Intelligent Gait Analysis

Recent advancements in data-driven modeling and wearable sensing technology open the door for real-time gait monitoring. The current method lays the foundation for such developments by offering:

• •

A framework to integrate real-world sensor data (e.g., IMUs, pressure insoles),

• •

Compatibility with deep learning-based feature extraction (e.g., for gait recognition),

• •

Extensibility to incorporate nonlinear feedback and control in smart prosthetics and assistive devices.

6.5 Discussion Summary

The proposed scheme not only advances the numerical treatment of singularly perturbed equations but also establishes a connection between rigorous mathematics and applied biomechanics. Its relevance spans across theoretical, clinical, and engineering domains. With further development, it holds potential for enhancing intelligent gait systems and biomechanical modeling pipelines.

6.6 Conclusion and Future Work

In this study, we developed a robust and efficient numerical framework for simulating vertical motion in human gait governed by a singularly perturbed second-order differential equation. The method integrates domain decomposition with a time-rescaling transformation to resolve steep boundary layers effectively.

By discretizing the inner and outer regions using central difference schemes and solving the resulting tridiagonal systems via the Thomas algorithm, we achieved second-order accuracy with uniform convergence. Theoretical error analysis confirms this convergence rate, while numerical experiments validate the method’s effectiveness across a range of perturbation parameters and boundary-layer-dominated regimes.

The proposed method stands out due to:

• •

Its ability to handle rapid solution changes near boundaries without mesh refinement,

• •

Its modularity, allowing adaptation to a variety of biomechanical models,

• •

Its computational efficiency, suitable for near real-time simulations.

Future Work

Several promising directions can be pursued to extend the scope and impact of this work:

• •

Nonlinear modeling: Incorporate time-dependent or nonlinear stiffness and damping functions to reflect more complex gait dynamics.

• •

Gait phase transitions: Extend the model to simulate transitions between gait modes (e.g., walk-to-run, stance-to-swing).

• •

Data-driven integration: Fuse the numerical model with real-world gait datasets and wearable sensor data for subject-specific simulations.

• •

Adaptive meshing: Develop an adaptive mesh refinement strategy to optimize grid resolution near layers.

• •

Control and feedback: Explore integration with control-theoretic models for applications in smart prosthetics and robotics.

• •

Three-dimensional extension: Generalize the current 1D framework to simulate multi-DOF or 3D gait models using vector-valued formulations.

The hybrid structure of the proposed method provides a solid foundation for future interdisciplinary work in computational biomechanics, wearable systems, and intelligent gait analysis technologies.