Multi- stage variational iteration analysis of a spring pendulum system

Abstract

The current research evaluates the nonlinear motion of a spring pendulum system with a constant mass based on a variational approach. Using first principles, the Lagrangian of the model was created in polar coordinates, which produced the Euler-Lagrange equations that describe the motion of the radial and angular components as functions of each other. Under the premise of small oscillation, the equations were simplified into a more easily solved nonlinear equation. To provide accurate analytical solutions for large timescales, the Ms-VIM) solution technique was used to overcome the limitations associated with the convergence rate of the classical VIM solution technique. In addition to the analytical solutions obtained, the solutions were also evaluated using the fourth order Runge-Kutta (RK4) method. An assessment of the analytical and numerical solutions indicated very good agreement among them, indicating that the proposed solution technique was reliable and efficient. Phase plots and error analyses were also conducted to illustrate the nature of the motion and the stability characteristics of the system. This research provided a viable analytical-numeric solution strategy for evaluating the motion of a constant mass oscillating system and can be applied to more complicated nonlinear mechanical models.

Keywords

spring pendulum nonlinear dynamics multi- step variational iteration method

1. Introduction

Spring-pendulums are an example of a nonlinear oscillator; they can produce a vast array of nonlinear behavior, including amplitude dependent frequencies, bifurcation, and nonlinear energy transfer between interacting modes.^1–3

A typical example of a strongly nonlinear mechanical system that displays all of these nonlinear behaviors is the classical spring pendulum^4–6; this model exhibits both radial and angular motions. Due to the strong coupling between the extension of the spring and the rotation of the pendulum, the length of the pendulum changes during motion producing complex behavior including nonlinear resonances, internal energy transfer, and modulation of oscillation amplitudes.^7–9

Therefore, the spring pendulum model is a classic representation of a nonlinearly coupled mechanical system; it is a representative system for illustrating the effects of nonlinear coupling on the dynamics of a mechanical system.¹⁰

Additionally, the motion of the spring pendulum has application to many engineering disciplines, including crane payload dynamics, robotic manipulator design with flexible links, vibration control devices, and aerospace tether systems.^11,12 A better understanding of the nonlinear interactions between the radial and angular motions will help engineers to design mechanical systems that are stable and perform well under dynamic loading.

The motion of the spring pendulum can be described using Lagrange’s equations of motion. The Lagrange equations of motion are a general method for determining the equations of motion of a nonlinear dynamical system with many degrees of freedom.^13–15 The Lagrange method uses the total energy of the system, which includes both the kinetic and potential energies, to derive the equations of motion of the system. The Lagrange method is commonly used in classical mechanics and nonlinear dynamics for determining the equations of motion of systems that can be described using generalized coordinates.^16–18 Additionally, the Lagrange method can be used to extend the description of a system to include dissipative forces through the use of the Rayleigh dissipation function.^19,20

Due to the highly nonlinear form of the equations of motion of the spring pendulum, exact analytical solutions cannot be found using the Lagrange method or other analytical methods. Therefore, a number of analytical and numerical approximation methods have been developed to study the nonlinear vibrations of the spring pendulum. Classical perturbation methods such as the method of multiple scales, the harmonic balance method, and the averaging method have been widely used to study nonlinear oscillation problems.^21–23 However, these classical perturbation methods typically rely on a small parameter and can lose their accuracy if the nonlinearity becomes large or if long time simulations are necessary. Recently, other efficient analytical techniques have been developed for nonlinear oscillators, including He’s frequency method and the spreading harmonic balance method (see recent works in Fractals and related journals). These methods are highly effective for frequency-amplitude analysis, but the present study focuses on long-time integration of coupled state variables using the multi-stage variational iteration method.

In response to the limitations of classical perturbation methods, a number of new semi-analytical methods were developed in the late twentieth century. One of the most successful of these semi-analytical methods is the Variational Iteration Method (VIM), which was developed by He.^24,25 The VIM uses the principle of variation to construct a correction functional that converges rapidly to the exact solution of the equation. This correction functional is constructed using a Lagrange multiplier, and the resulting iterative formula does not involve linearization or discretization. The VIM has been shown to be useful for solving a wide range of nonlinear problems in physics and engineering.^26,27

Recent extensions of the VIM include the dual Lagrange multiplier approach,²⁸ which introduced two Lagrange multipliers to solve two-point boundary value problems in nonlinear oscillators. This innovation expands the VIM to a new class of problems beyond traditional initial value problems. While their method focuses on conservative single-DOF systems with boundary conditions, our work addresses the complementary challenge of long-time integration of initial value problems in coupled, multi-DOF systems with damping, using the multi-stage VIM formulation.

Although the VIM has several advantages, one of its main disadvantages is that it produces a local series solution, and the size of the convergence region for this solution can limit the accuracy of the solution over long periods of time. For example, in the case of nonlinear dynamical systems, the accumulation of truncation errors over time can lead to significant differences between the approximate solution and the actual solution. To improve the accuracy and robustness of the VIM for simulating long time responses, the Multi-Stage Variational Iteration Method (Ms-VIM) was developed as an improved version of the VIM.²⁵ The Ms-VIM divides the time interval into smaller sub-intervals and reconstructs the solution at each stage using the updated initial condition. This multi-stage approach substantially improves the accuracy and reliability of the Ms-VIM for simulating long time responses of nonlinear dynamical systems.

In recent years, researchers have focused an increasing amount of research on nonlinear dynamics in elastically-coupled systems such as the spring-pendulum and oscillation systems with relevance to applications in engineering. Amer et al.,^29,30 used numerical analysis to model the motion of a spring pendulum under planar forced motion and found that both nonlinear coupling and forcing significantly affect the system’s dynamics. The authors also applied the same approach to analyze a two degree of freedom spring pendulum traveling along a Lissajous path and illustrated how strong modal interaction can produce complicated oscillations. Additionally, Rouleau et al.³¹ investigated the nonlinear transition and periodic response of elastically-coupled pendulums and Szumiński and Maciejewski³² investigated the dynamics and nonintegrability of the double spring pendulum to demonstrate the potential complexity and sensitivity of strongly coupled nonlinear oscillators. All these investigations highlight the need for well established analytical-numerical methods that will allow the investigation of nonlinear spring pendulum systems over large periods of time.

Moreover, studies in recent years have explored fractal and fractional formulations of nonlinear pendulum oscillators in an attempt to better reproduce complex dynamic behaviors. In 33, the authors investigated the numerical response of a fractional rotating pendulum system and showed the influence of fractional parameters on oscillatory behavior and stability. In a similar direction,³⁴ studied the nonlinear dynamics of a rolling wheel pendulum using numerical methods. The sensitivity of the system to fractal effects and parameter changes is emphasized. In addition,³⁵ offered analytical solution procedures for the same class of oscillators, giving approximate formulations to supplement numerical results. These works collectively highlight the increasing significance of fractal and fractional modelling in the context of nonlinear vibration analysis and propel the design of robust iterative schemes for solving intricate oscillatory systems.

This paper presents an investigation of the nonlinear dynamics of a spring pendulum system using the Multi-Stage Variational Iteration Method. Both conservative and damp cases are examined. The Ms-VIM is used to find the semi-analytical solutions of the system, and the numerical reference solutions are generated by the classical fourth-order Runge-Kutta method. By comparing the analytical and numerical solutions, the accuracy, stability, and convergence properties of the proposed Ms-VIM method can be evaluated. The comparisons show that the multi-stage formulation of the VIM improves its ability to predict the behavior of the system and provide reliable solutions for the nonlinear oscillatory systems. Therefore, the proposed framework is a powerful and efficient analytical-numerical tool for studying nonlinear mechanical systems and can be expanded to more complex dynamical systems encountered in modern engineering applications.

2. Physics of the model

Below is an illustration of a spring pendulum with a mass $(m)$ attached at the lower end of a spring of constant $(K)$ and unstretched length $(l)$ . We begin by building the Lagrangian ( $L$ ) for this system, and then derive the Euler-Lagrange Equations (ELEs), or, the equations of motion (EOM). To get a full understanding of how to build the Lagrangian of a physical system, readers should refer to some classic literature on the subject; see^{36, 37} for example.

First, we must determine the mass’s $(m)$ kinetic and potential energies. The system’s potential (V) and kinetic (T) energies are as follows, respectively:

T = \frac{1}{2} m v^{2}

(1)

V = V_{g} + V_{s} = - m g y + \frac{1}{2} K r^{2}

(2)

Where,

\vec{v}

is the velocity of the mass, and it has two components

\vec{v} = \dot{x} \hat{i} + \dot{y} \hat{j}

In terms of polar coordinates $r$ and $θ$ , the energies T and V read:

T = \frac{m}{2} ({\dot{r}}^{2} + {(l + r)}^{2} {\dot{θ}}^{2})

(3)

V = - m g (l + r) \cos θ + \frac{1}{2} K r^{2}

(4)

Finally, the Lagrangian ( $L$ ) of the system; which is defined as $L = T - V$ ; read:

L = \frac{m}{2} ({\dot{r}}^{2} + {(l + r)}^{2} {\dot{θ}}^{2}) + m g (l + r) \cos θ - \frac{1}{2} K r^{2}

(5)

Below, we will use (10) and the relation $\frac{d}{d t} (\frac{\partial L}{\partial {\dot{q}}_{i}}) - \frac{\partial L}{\partial q_{i}} = 0$ to derive the ELE’s (i.e., EOM). We are going to consider the simplest case, where $m = m_{o}$ and here we will consider only small oscillations.

In this case, θ is very small and thus

\sin θ \approx θ

(6a)

\cos θ \approx 1 - \frac{θ^{2}}{2}

(6b)

This second-order approximation retains the essential nonlinear coupling term $\frac{1}{2} g θ^{2}$ in the radial equation while keeping the analysis tractable; higher-order terms are negligible for $θ \leq 0.1 rad$ .

Applying (6) directly to (5). Thus

L = \frac{m_{o}}{2} ({\dot{r}}^{2} + {(l + r)}^{2} {\dot{θ}}^{2}) + m_{0} g (l + r) (1 - \frac{θ^{2}}{2}) - \frac{1}{2} K r^{2}

(7)

Simplyfing

L = \frac{m_{0}}{2} {\dot{r}}^{2} + \frac{m_{0}}{2} {(l + r)}^{2} {\dot{θ}}^{2} + m_{0} g (l + r) - \frac{m_{0}}{2} g (l + r) θ^{2} - \frac{1}{2} K r^{2}

(8)

Now, substitute (13) into the relation

\frac{d}{d t} (\frac{\partial L}{\partial {\dot{q}}_{i}}) - \frac{\partial L}{\partial q_{i}} = 0

, with

q_{1} = r

, and

q_{2} = θ

. we have for

q_{1} = r

m_{o} \ddot{r} - m_{o} (l + r) {\dot{θ}}^{2} + \frac{m_{o}}{2} g θ^{2} - m_{o} g + K r = 0

Simplifying

\ddot{r} - (l + r) {\dot{θ}}^{2} + \frac{1}{2} g θ^{2} - g + ω_{o}^{2} r = 0

(9)

Where

ω_{o}^{2} = \frac{K}{m_{o}}

While for $q_{2} = θ$ , we have

m_{o} {(l + r)}^{2} \ddot{θ} + 2 m_{o} (l + r) \dot{r} \dot{θ} + m_{o} g (l + r) θ = 0

Simplifying

\ddot{θ} + \frac{2}{(l + r)} \dot{r} \dot{θ} + \frac{g}{(l + r)} θ = 0

(10)

Parameters:

l = 1 m, g = 9.81 m / s^{2}, ω_{o} = 2

Initial conditions are

r (0) = 0.1, \dot{r} (0) = 0, θ (0) = 0.1 and \dot{θ} (0) = 0

Now, let us assume that the spring pendulum is subject to a viscous daming that opposes both the radial and angular motions. The dissipative effects can be introduced through the Rayleigh dissipation function defined as

R = \frac{1}{2} C_{r} {\dot{r}}^{2} + \frac{1}{2} C_{θ} {\dot{θ}}^{2} \Rightarrow \frac{\partial R}{\partial \dot{r}} = C_{r} \dot{r} and \frac{\partial R}{\partial \dot{θ}} = C_{θ} \dot{θ}

(11)

where c_r and c_θ represent the radial and angular damping coefficients, respectively. These coefficients characterize the viscous resistance acting against the motion of the system. The presence of damping modifies the equations of motion through the generalized dissipative forces derived from the Rayleigh dissipation function.

So, Eqns.(9, 10) read:

\ddot{r} - (l + r) {\dot{θ}}^{2} + \frac{1}{2} g θ^{2} - g + ω_{o}^{2} r + C_{r} \dot{r} = 0

(12)

\ddot{θ} + \frac{2}{(l + r)} \dot{r} \dot{θ} + \frac{g}{(l + r)} θ + C_{θ} \dot{θ} = 0

(13)

3. Method of solution

In this section, we present the fundamental concepts of VIM. The findings presented here are also available in 21 and the works referenced therein.

3.1. He’s variational iteration method

Consider, at this point, the following system:

T x (t) = g (t), t \in I

In this context, T is a differential operator,

x

is a continuous function over

t \in I,

and

g (t)

is a known function. The main idea of VIM is to decompose the differential operator T into its linear and nonlinear components:

L x (t) + N x (t) = g (t)

(14)

Here,

L

and

N

stand for the linear and nonlinear components of the operator, respectively.

The VIM can be regarded as a modification of the Lagrange multiplier method.²² In what follows, we briefly review the Lagrange multiplier method and clarify its connection to VIM. This presentation is adapted from.³⁶

An initial function $x^{0}$ satisfying $L x^{0} = 0$ is adopted in the Lagrange multiplier method. Next, by means of the functional,

x (t_{1}) = x^{0} (t_{1}) + \int_{t_{0}}^{t_{1}} λ {L x^{0} (s) + N x^{0} (s) - g (s)} d s

This yields an approximation at $t_{1}$ , with $λ$ referred to the Lagrange multiplier.

According to J. H. He,³⁷ the correction functional is given by:

x^{n + 1} (t) = x^{n} (t) + \int_{t_{0}}^{t} λ {L x^{n} (s) + N {\tilde{x}}^{n} (s) - g (s)} d s

Here,

λ = λ (s; t)

denotes the Lagrange multiplier, which can be determined using variational theory.^36,38

The function $x^{n} (t)$ denotes $n$ th-level approximate solution and the ${\tilde{x}}^{n} (s)$ stands for the restricted variation, i.e., ${δ \tilde{x}}^{n} (s) = 0$ for all $n \in N$ , see 39.

The method is essentially based on defining a Lagrange multiplier that meets the following condition (see 40)

δ x^{n + 1} (t) = δ x^{n} (t) + δ \int_{t_{0}}^{t} λ {L x^{n} (s) + N {\tilde{x}}^{n} (s) - g (s)} d s = 0

In the case of nonlinear differential equations, the Lagrange multipliers are identified by means of restricted variations. Minimizing the use of restricted variations results in a more precise Lagrange multiplier and consequently accelerates the convergence of the approximate solutions.

Classical VIM provides a local series solution expanded about the initial time, whose validity is limited by a finite radius of convergence. In nonlinear systems, this local nature leads to rapid accumulation of truncation and phase errors when the method is applied over long time intervals. Increasing the number of perturbation terms improves local accuracy but does not guarantee global reliability. Ms-VIM overcomes this limitation by subdividing the time domain and reconstructing the VIM soluıtion at each subinterval using updated initial conditions. This reinitialization maintains the solution within the convergence region of the local series, thereby enhancing long-time accuracy and stability. Consequently, Ms-VIM is more suitable than classical VIM for long-time integration of nonlinear dynamical systems.

3.2. Multi-stage variational iteration method

Batiha et al. developed the Ms-VIM⁴¹ to improve the performance of VIM on extended time intervals. The Ms-VIM approach determines the solution of equation (16) on $[t_{0}, T)$ through a subdivision of the interval into smaller segments $[t_{0}, t_{1}), [t_{1}, t_{2}), \dots, [t_{n}, T),$ with the following correction functional implemented separately on each subinterval:

x^{n + 1} (t) = x^{n} (t) + \int_{t^{*}}^{t} λ {L x^{n} (s) + N {\tilde{x}}^{n} (s) - g (s)} d s

where

t^{*}

is taken sequentially as

t_{0}, t_{1}, t_{2}, \dots, t_{n} .

In this framework, the initial approximation on

[t^{*}, t_{n})

is given by the solution from the preceding interval

[t_{n - 1}, t^{*})

at the point

t^{*}

, namely

x_{0} (t) = x_{n} (t^{*}) .

In what follows, we demonstrate that the Lagrange multipliers in VIM and Ms-VIM are identical. For further details on the multi-stage version of the method, see 39, 42–44.

A practical note on subinterval size selection: smaller subintervals improve accuracy but increase computational cost. Users should adjust the subinterval size according to the nonlinearity and stiffness of the system. In this study, the same step size is used for both the Ms-VIM and the RK4 method.

Case 1

The subsequent analysis applies methods VIM and Ms-VIM to equations (14)-(15). To rewrite Eqs. (9)–(10) as a first-order system, we introduce the state variables

w = r, x = r^{'}, y = θ, z = θ^{'}

Using these definitions, Eqs. (9)–(10) can be equivalently written as the following system of four first-order differential equations

\begin{array}{l} w^{'} = x, \\ x^{'} = (l + w) z^{2} - \frac{1}{2} g y^{2} + g - ω_{o}^{2} w . \\ y^{'} = z, \\ z^{'} = - \frac{2}{l + w} x z - \frac{g}{l + w} y . \end{array}

(15)

The corresponding initial conditions become

w (0) = 0.1, x (0) = 0, y (0) = 0.1, z (0) = 0

(16)

This representation expresses the second-order system (9)–(10) as an equivalent four-dimensional first-order dynamical system, which is suitable for analytical techniques (Variational iteration method and multi-stage Variational iteration method) as well as for numerical solvers such as the fourth order Runge–Kutta method.

Based on the classical VIM approach widely adopted in the literature, the following correction functionals are derived:

\begin{array}{l} w^{n + 1} (t) = w^{n} (t) + \int_{0}^{t} λ_{1} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s . \\ x^{n + 1} (t) = x^{n} (t) + \int_{0}^{t} λ_{2} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s)} d s \\ y^{n + 1} (t) = y^{n} (t) + \int_{0}^{t} λ_{3} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s \\ z^{n + 1} (t) = z^{n} (t) + \int_{0}^{t} λ_{4} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) {\tilde{z}}^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s)} d s \end{array}

(17)

where

{(w^{n})}^{'}, {(x^{n})}^{'}, {(y^{n})}^{'}

, and

{(z^{n})}^{'}

denotes the derivatives

d w^{n} / d s, d x^{n} / d s, d y^{n} / d s,

and

d z^{n} / d s,

respectively. There

λ_{i} = λ_{i} (s; t)

are the Lagrange multipliers and

{\tilde{w}}^{n} (s), {\tilde{x}}^{n} (s), {\tilde{y}}^{n} (s)

, and

{\tilde{z}}^{n} (s)

are the restricted variations, that is,

δ {\tilde{w}}^{n} (s), δ {\tilde{x}}^{n} (s), δ {\tilde{y}}^{n} (s)

and

δ {\tilde{z}}^{n} (s)

equal to

0

for all

n \in N

. It must be note that

δ w^{n} (0), δ x^{n} (0), δ y^{n} (0)

and

δ z^{n} (0)

equal to

0

for all

n \in N

Notably, several linear terms in (17) are imposed as restricted variations, simplifying the process of identifying the Lagrange multipliers. The Lagrange multipliers are obtained by computing the variation of the correction functionals, namely:

\begin{array}{l} δ w^{n + 1} (t) = δ w^{n} (t) + δ \int_{0}^{t} λ_{1} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s . \\ {δ x}^{n + 1} (t) = δ x^{n} (t) + δ \int_{0}^{t} λ_{2} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s)} d s . \\ δ y^{n + 1} (t) = δ y^{n} (t) + δ \int_{0}^{t} λ_{3} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s . \\ δ z^{n + 1} (t) = {δ z}^{n} (t) + δ \int_{0}^{t} λ_{4} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) {\tilde{z}}^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s)} d s . \end{array}

Through integration by parts and the calculus of variations, we obtain the Lagrange multipliers in the form:

λ_{1} (s; t) = λ_{2} (s; t) = λ_{3} (s; t) = λ_{4} (s; t) = - 1

As a result, the correction functionals corresponding to VIM can be written as:

\begin{array}{l} w^{n + 1} (t) = w^{n} (t) - \int_{0}^{t} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s . \\ x^{n + 1} (t) = x^{n} (t) - \int_{0}^{t} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s)} d s \\ y^{n + 1} (t) = y^{n} (t) - \int_{0}^{t} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s \\ z^{n + 1} (t) = z^{n} (t) - \int_{0}^{t} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) {\tilde{z}}^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s)} d s \end{array}

When

t^{*}

is taken as the lower limit of integration in the correction functionals, the multi-stage version of the method is obtained in the following form:

\begin{array}{l} w^{n + 1} (t) = w^{n} (t) - \int_{t^{*}}^{t} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s . \\ x^{n + 1} (t) = x^{n} (t) - \int_{t^{*}}^{t} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s)} d s \\ y^{n + 1} (t) = y^{n} (t) - \int_{t^{*}}^{t} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s \\ z^{n + 1} (t) = z^{n} (t) - \int_{t^{*}}^{t} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) {\tilde{z}}^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s)} d s \end{array}

Unlike the classical VIM, the multi-stage variational iteration method satisfies

δ w^{n} (t^{*}) = δ x^{n} (t^{*}) = δ y^{n} (t^{*}) = δ z^{n} (t^{*}) = 0

for every

n \in N

. Here,

t^{*} \in {0, t_{1}, t_{2}, \dots, t_{n}}

and the interval

[0, T)

is divided into subintervals

[0, t_{1}), [t_{1}, t_{2}), \dots, [t_{n}, T)

Case 2

In order to transform Eqs. (12)–(13) into an equivalent first-order system, we define the state variables as follows:

w = r, x = r^{'}, y = θ, z = θ^{'}

Employing these definitions, Eqs. (12)–(13) take the form of the following system consisting of four first-order differential equations:

\begin{array}{l} w^{'} = x, \\ x^{'} = (l + w) z^{2} - \frac{1}{2} g y^{2} + g - ω_{o}^{2} w - C_{w} x . \\ y^{'} = z, \\ z^{'} = - \frac{2}{l + w} x z - \frac{g}{l + w} y - c_{y} z . \end{array}

(18)

Under this formulation, the initial conditions become:

w (0) = 0.1, x (0) = 0, y (0) = 0.1, z (0) = 0

(19)

This formulation converts the second-order system (12)–(13) into an equivalent four-dimensional first-order dynamical system, making it appropriate for analytical approaches such as the VIM and the Ms-VIM, as well as for numerical schemes like the classical fourth-order Runge–Kutta method.

By applying the classical VIM scheme frequently used in earlier works, we formulate the correction functionals as:

\begin{array}{l} w^{n + 1} (t) = w^{n} (t) + \int_{0}^{t} λ_{1} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s . \\ x^{n + 1} (t) = x^{n} (t) + \int_{0}^{t} λ_{2} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s) + C_{w} x^{n} (s)} d s . \\ y^{n + 1} (t) = y^{n} (t) + \int_{0}^{t} λ_{3} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s . \\ z^{n + 1} (t) = z^{n} (t) + \int_{0}^{t} λ_{4} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) z^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s) + c_{y} z^{n} (s)} d s . \end{array}

(20)

where

{(w^{n})}^{'}, {(x^{n})}^{'}, {(y^{n})}^{'}

, and

{(z^{n})}^{'}

denotes the derivatives

d w^{n} / d s, d x^{n} / d s, d y^{n} / d s,

and

d z^{n} / d s,

respectively. There

λ_{i} = λ_{i} (s; t)

are the Lagrange multipliers and

{\tilde{w}}^{n} (s), {\tilde{x}}^{n} (s), {\tilde{y}}^{n} (s)

, and

{\tilde{z}}^{n} (s)

are the restricted variations, that is,

δ {\tilde{w}}^{n} (s), δ {\tilde{x}}^{n} (s), δ {\tilde{y}}^{n} (s)

and

δ {\tilde{z}}^{n} (s)

equal to

0

for all

n \in N

. It must be note that

δ w^{n} (0), δ x^{n} (0), δ y^{n} (0)

and

δ z^{n} (0)

equal to

0

for all

n \in N

It is worth noting that several linear terms in (20) are treated as restricted variations, which facilitates the identification of the Lagrange multipliers. These multipliers are determined by taking the variation of the correction functionals, namely:

\begin{array}{l} {δ w}^{n + 1} (t) = {δ w}^{n} (t) + δ \int_{0}^{t} λ_{1} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s, \\ {δ x}^{n + 1} (t) = δ x^{n} (t) + δ \int_{0}^{t} λ_{2} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s) + C_{w} x^{n} (s)} d s, \\ δ y^{n + 1} (t) = δ y^{n} (t) + δ \int_{0}^{t} λ_{3} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s, \\ δ z^{n + 1} (t) = {δ z}^{n} (t) + δ \int_{0}^{t} λ_{4} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) z^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s) + c_{y} z^{n} (s)} d s, \end{array}

The Lagrange multipliers are subsequently derived, via integration by parts and variational analysis, as follows:

λ_{1} (s; t) = λ_{2} (s; t) = λ_{3} (s; t) = - 1, λ_{4} (s; t) = - e^{c_{y} (t - s)} .

Consequently, we obtain the following correction functionals for VIM:

\begin{array}{l} w^{n + 1} (t) = w^{n} (t) - \int_{0}^{t} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s \\ x^{n + 1} (t) = x^{n} (t) - \int_{0}^{t} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s)} d s \\ y^{n + 1} (t) = y^{n} (t) - \int_{0}^{t} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s \\ z^{n + 1} (t) = z^{n} (t) - \int_{0}^{t} e^{c_{y} (t - s)} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) z^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s) + c_{y} z^{n} (s)} d s . \end{array}

By taking $t^{*}$ as the lower bound of integration in the correction functionals, we derive the multi-stage version of the method as follows:

\begin{array}{l} w^{n + 1} (t) = w^{n} (t) - \int_{t^{*}}^{t} {{(w^{n})}^{'} (s) {- x}^{n} (s)} d s \\ x^{n + 1} (t) = x^{n} (t) - \int_{t^{*}}^{t} {{(x^{n})}^{'} (s) - (l + {\tilde{w}}^{n} (s)) {{\tilde{z}}^{n} (s)}^{2} + \frac{1}{2} g {{\tilde{y}}^{n} (s)}^{2} - g + ω_{o}^{2} {\tilde{w}}^{n} (s)} d s \\ y^{n + 1} (t) = y^{n} (t) - \int_{t^{*}}^{t} {{(y^{n})}^{'} (s) - {\tilde{z}}^{n} (s)} d s \\ z^{n + 1} (t) = z^{n} (t) - \int_{t^{*}}^{t} e^{c_{y} (t - s)} {{(z^{n})}^{'} (s) + \frac{2}{l + {\tilde{w}}^{n} (s)} {\tilde{x}}^{n} (s) z^{n} (s) + \frac{g}{l + {\tilde{w}}^{n} (s)} {\tilde{y}}^{n} (s) + c_{y} z^{n} (s)} d s . \end{array}

Unlike the classical VIM, the Ms-VIM satisfies

δ w^{n} (t^{*}) = δ x^{n} (t^{*}) = δ y^{n} (t^{*}) = δ z^{n} (t^{*}) = 0

for every

n \in N

. Here,

t^{*} \in {0, t_{1}, t_{2}, \dots, t_{n}}

and the interval

[0, T)

is divided into subintervals

[0, t_{1}), [t_{1}, t_{2}), \dots, [t_{n}, T)

4. Results

The multi-stage variational iteration method algorithm is implemented using the computer algebra system Matlab. A fixed step size of $Δ t = 5 \times 10^{- 5}$ is adopted for all problems. All simulations are carried out over the time interval $t \in [0, 20]$ .

Case 1

$l = 1$ , $ω_{o} = 2,$ $g = 9.81 m / s^{2}$ and the initial conditions $w (0) = 0.1, x (0) = 0, y (0) = 0.1 a n d z (0) = 0$ .

The dynamic behavior of the spring pendulum was examined through the use of the multi-stage variational iteration method (Ms-VIM). To establish the reliability of the analytical formulation, numerical simulations were carried out through the application of the classical fourth-order Runge-Kutta (RK4) method with identical parameter conditions. The simulations were performed on the time interval $t \in [0, T]$ , where $T$ is an arbitrary time limit, and at every time interval using a constant time step to maintain numerical stability and to assure numerical consistency in the comparison.

The system parameters were set as $k = 1, g = 9.81 m / s^{2}$ . The initial conditions were set so that they would represent small oscillations. The choice of initial conditions was made to be consistent with the linearized assumptions made when deriving the governing equations. At such small amplitudes, the nonlinear coupling between the radial extension of the spring and the angular motion of the pendulum will become apparent.

Figures 1–4 show the time histories of the state variables obtained by using the Ms-VIM method and the RK4 numerical method. The time history of the radial position of the end of the spring, $(t)$ , is presented in Figure 1. The time history of the radial position of the end of the spring exhibits a periodic oscillatory pattern, indicating the restoring forces developed due to the spring stiffness. The time histories of the radial position of the end of the spring obtained using the Ms-VIM method and the RK4 numerical method are very similar to each other for the duration of the simulation. This similarity is indicative of how well the multi-stage version of the Ms-VIM method can sustain the convergence of the analytical approximation for a longer period of time. In addition, the time history of the angular position of the pendulum, $θ (t)$ , is depicted in Figure 2. The time history of the angular position of the pendulum also exhibits a periodic oscillatory pattern; therefore, it behaves like a pendulum. Since the radial motion changes the length of the pendulum, there exists a nonlinear coupling between the radial and angular motions. Due to this coupling, there is a slight modification of the oscillation amplitude, which is very accurately predicted by the Ms-VIM method.

Figure 1.

The spring pendulum system.

Figure 2.

Time history of the radial displacement $r (t)$ of the spring–pendulum system obtained using the multi-stage variational iteration method (Ms-VIM) and the fourth-order Runge–Kutta (RK4) method.

Figure 3.

Time history of the radial velocity $r^{'} (t)$ of the spring–pendulum system obtained using the multi-stage variational iteration method (Ms-VIM) and the fourth-order Runge–Kutta (RK4) method.

Figure 4.

Time history of the angular displacement $θ (t)$ of the spring–pendulum system obtained using the multi-stage variational iteration method (Ms-VIM) and the fourth-order Runge–Kutta (RK4) method.

The time histories of the velocities of the end of the spring, $\dot{r} (t)$ , and the time histories of the angular velocity of the pendulum, $\dot{θ} (t)$ , are presented in Figures 3 and 4, respectively. The time histories of both velocities show a periodic variation, which corresponds to the time histories of the displacements. The time histories of the velocities provide additional evidence of the strong dynamic coupling between the two degrees of freedom. Once again, the time histories of the velocities obtained using the Ms-VIM method are in excellent agreement with the time histories of the velocities obtained using the RK4 numerical method; therefore, they indicate that the Ms-VIM method provides excellent accuracy in predicting the velocities of the system.

The dynamical properties of the system are demonstrated using phase plots. A phase plot of the radial motion of the spring in the $r$ - $\dot{r}$ plane is presented in Figure 5. The phase plot represents a closed path and indicates that the radial motion of the spring undergoes periodic, bounded oscillations. Phase plots of the type presented in Figure 5 are typical for conservative nonlinear oscillators and indicate that the radial motion is dynamically stable under the assumed conditions. A phase plot of the angular motion of the pendulum in the $θ$ - $\dot{θ}$ plane is presented in Figure 6. The phase plot of the angular motion also represents a closed path and is symmetric about the equilibrium point. Therefore, the phase plot confirms that the pendulum motion is also dynamically stable. The fact that the phase plots do not contain any diverging or spiraling paths indicates that the system does not exhibit any instability or increase in energy within the studied time domain. The phase plots clearly demonstrate that the spring-pendulum system has nonlinear dynamics that are stable.

Figure 5.

Time history of the angular displacement $θ^{'} (t)$ of the spring–pendulum system obtained using the multi-stage variational iteration method (Ms-VIM) and the fourth-order Runge–Kutta (RK4) method.

Figure 6.

Phase portrait of the spring–pendulum system in the $r (t) - r^{'} (t)$ plane.

The accuracy of the proposed analytical technique was evaluated in a quantitative manner by computing the absolute differences between the Ms-VIM solution and the RK4 numerical solution for all of the state variables.

Physically, the results obtained show the strong nonlinear coupling between the radial and angular motions of the spring-pendulum system. The effective pendulum length is continuously modified by the periodic extension and contraction of the spring, which affects the angular oscillation characteristics. This interaction leads to amplitude modulation and small frequency changes of either mode of motion. Furthermore, in the conservative case, the closed phase trajectories correspond to bounded and dynamically stable oscillations, whereas in the damped case the spiral phase trajectories exhibit a constant energy dissipation and an approach to a stable equilibrium configuration. These results confirm the effectiveness of the proposed model for capturing the essential nonlinear dynamics of the spring–pendulum system.

Figures 7–10 present the absolute error distributions for $r (t), θ (t)$ , $\dot{r} (t)$ and $\dot{θ} (t)$ , respectively. The absolute error is very small for all of the state variables throughout the entire time domain. Furthermore, the absolute error does not accumulate with time. The ability of the Ms-VIM method to prevent the accumulation of absolute errors is one of the major advantages of the multi-stage formulation.

Figure 7.

Phase portrait of the spring–pendulum system in the $θ (t) - θ^{'} (t)$ plane.

Figure 8.

Absolute error distribution for the radial displacement $r (t)$ obtained from the difference between the multi-stage variational iteration method (Ms-VIM) solution and the fourth-order Runge–Kutta (RK4) numerical solution.

Figure 9.

Absolute error distribution for the radial displacement $r^{'} (t)$ obtained from the difference between the multi-stage variational iteration method (Ms-VIM) solution and the fourth-order Runge–Kutta (RK4) numerical solution.

Figure 10.

Absolute error distribution for the radial displacement $θ (t)$ obtained from the difference between the multi-stage variational iteration method (Ms-VIM) solution and the fourth-order Runge–Kutta (RK4) numerical solution.

As previously mentioned, in classical perturbation methods, truncation errors typically grow very quickly during long-time simulations due to the limitations on the convergence radius of the power series expansions. In contrast, the Ms-VIM method continuously updates the initial conditions at the start of each subinterval. Therefore, the Ms-VIM method prevents the propagation of the accumulated errors and maintains the accuracy of the analytical approximation for extended periods of time.

The results of comparing the classical VIM and the multi-stage VIM for all of the state variables are presented in Figures 11–14. The classical VIM consistently diverges from the RK4 numerical solution as the time increases for all of the state variables. The primary cause of this divergence is the limited convergence radius of the perturbation series employed in the classical formulation. As the time interval of integration is increased, truncation and phase errors accumulate and lead to noticeable discrepancies. Conversely, the Ms-VIM solutions remain in excellent agreement with the RK4 numerical solution for all of the state variables throughout the entire time interval. The subdivision of the time domain into smaller subintervals ensures that the local analytical approximation always remains within its convergence region. Therefore, the comparisons clearly demonstrate that the multi-stage implementation significantly enhances the long-term stability and accuracy of the variational iteration method.

Figure 11.

Absolute error distribution for the radial displacement $θ^{'} (t)$ obtained from the difference between the multi-stage variational iteration method (Ms-VIM) solution and the fourth-order Runge–Kutta (RK4) numerical solution.

Figure 12.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $r (t)$ of the spring–pendulum system.

Figure 13.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $r^{'} (t)$ of the spring–pendulum system.

Figure 14.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $θ (t)$ of the spring–pendulum system.

Physically, the results obtained in this study emphasize the inherent coupling between the radial and angular motions of the spring pendulum. The oscillatory radial extension of the spring modifies the effective length of the pendulum, which in turn modifies the angular motion of the pendulum. This mutual interaction produces nonlinear effects in the modulations of both the displacement and velocity responses of the system. Nevertheless, despite the nonlinear coupling between the radial and angular motions, the system maintains bounded oscillations and closed phase trajectories, thereby establishing that the equilibrium configuration of the system remains dynamically stable under the specified parameter conditions.

Case 2

$l = 1$ , $ω_{o} = 2,$ $g = 9.81 m / s^{2},$ $c_{θ} = 0.005, C_{w} = 0$ and the initial conditions $w (0) = 0.1, x (0) = 0, y (0) = 0.1 a n d z (0) = 0$ .

In the second instance, we study the impact of viscous damping on the dynamics of the spring pendulum system. Viscous resistances are included by means of the Rayleigh dissipation function. Energy loss, which occurs due to both radial spring deformation and the rotation of the pendulum, is represented by the two damping coefficients c (r)and cθ.

We use the same computational framework as in Case 1. To ensure that our analytical solutions are reliable and accurate, we compare the analytical solution obtained by means of the multi-stage variational iteration method (Ms-VIM) to the numerical solution obtained by means of the fourth-order Runge-Kutta (RK4) method.

The time-histories of the radial displacement r(t), the angular displacement $θ (t)$ , the radial velocity $\dot{r} (t)$ and the angular velocity $\dot{θ} (t)$ for the damped system are depicted in Figures 15–18. Figure 15 depicts the radial displacement $r (t)$ . Due to the presence of viscous damping, the oscillation amplitudes decrease with increasing time. The continuous loss of mechanical energy of the system is reflected in this behavior. The analytical solution obtained by means of Ms-VIM remains in perfect agreement with the numerical solution obtained by means of RK4 throughout the entire time interval, thereby illustrating the robustness of the proposed method for modeling dissipative nonlinear systems.

Figure 15.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $θ^{'} (t)$ of the spring–pendulum system.

Figure 16.

Figure 17.

Figure 18.

The angular displacement $θ (t)$ is presented in Figure 16. The angular oscillations are subject to a similar decay in amplitude, thus demonstrating that the viscous damping has an effect on both degrees of freedom. Due to the fact that the radial displacement leads to a modification of the effective pendulum length, the damping in one degree of freedom indirectly affects the other, leading to a coupled decay of the oscillations.

The velocity responses $\dot{r} (t)$ and $\dot{θ} (t)$ are shown in Figures 17 and 18, respectively. Both velocity components present oscillatory behavior with an exponentially decaying amplitude. The diminishing of the maximum velocities clearly indicates the dissipative nature of the system. Again, the analytical and numerical results remain almost indistinguishable, thereby demonstrating the ability of the Ms-VIM method to correctly describe the transient dynamics of damped nonlinear oscillators.

The effect of damping on the stability of the system is further illustrated in Figures 19 and 20, which depict the phase portraits of the damped system. The phase trajectory in the $r$ - $\dot{r}$ plane is depicted in Figure 19. Contrary to the closed elliptical trajectory observed in the undamped case, the trajectory now forms a spiral, which, as time progresses, is directed toward the equilibrium point. This type of behavior is typical of dissipative systems, where the energy is continuously removed by the action of the damping forces. A similar spiral trajectory is observed in the phase portrait of the $θ$ - $\dot{θ}$ plane depicted in Figure 20, thereby demonstrating that the angular motion is also subject to a gradual dissipation of energy. The convergence of the phase trajectories toward the origin confirms that the system will settle into a dynamically stable equilibrium configuration. These findings demonstrate that the addition of damping transforms the conservative oscillatory system into a stable dissipative dynamical system, where the oscillations gradually disappear over time.

Figure 19.

Figure 20.

Phase portrait of the spring–pendulum system in the $r (t) - r^{'} (t)$ plane.

To evaluate the accuracy of the analytical Ms-VIM solutions, we computed the absolute error with respect to the numerical RK4 solution. The absolute errors related to the $r (t)$ , $θ (t)$ , $\dot{r} (t)$ and $\dot{θ} (t)$ are depicted in Figures 21–24, respectively. Extremely low values of the errors were found across the entire time domain, thereby confirming the high accuracy of the analytical approximation.

Even in the presence of damping, which introduces additional nonlinear terms into the equations of motion, the Ms-VIM method provides stable and accurate predictions. The lack of significant error accumulation, which would have occurred if the errors had accumulated over time, further demonstrates the efficiency of the multi-stage formulation in limiting the truncation errors over long simulation times.

Figure 21.

Phase portrait of the spring–pendulum system in the $θ (t) - θ^{'} (t)$ plane.

Figure 22.

Figure 23.

Figure 24.

A comparative analysis of the results obtained by the classical VIM and the Ms-VIM for all four state variables is provided in Figures 25–28. The classical VIM solution starts to deviate from the numerical solution over time. This deviation can be attributed to the limited convergence radius of the perturbation series employed in the classical formulation. Over time, truncation and phase errors accumulate and lead to a degradation of the accuracy of the classical approach. On the contrary, the Ms-VIM results remain virtually identical to the RK4 numerical solutions throughout the entire simulation period.(Figure 29) The subdivision of the entire time domain into smaller subdomains ensures that the analytical series remains inside its convergence zone. Thus, the multi-stage approach significantly enhances the long-term stability and reliability of the analytical solution.

Figure 25.

Figure 26.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $r (t)$ of the spring–pendulum system.

Figure 27.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $r^{'} (t)$ of the spring–pendulum system.

Figure 28.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $θ (t)$ of the spring–pendulum system.

Figure 29.

Comparison between the multi-stage variational iteration method (Ms-VIM) and the classical variational iteration method (VIM) for the radial displacement $θ^{'} (t)$ of the spring–pendulum system.

Physically, the results demonstrate how damping influences the nonlinear dynamics of the spring pendulum system. The presence of viscous forces removes energy continuously from the system, so that the oscillation amplitudes of both radial and angular motions decrease over time.

Although the nonlinear coupling between the radial extension of the spring and the angular displacement of the pendulum exists, the system behaves stably and ultimately tends toward an equilibrium state. The spiral shape of the phase trajectories demonstrates that the system moves toward a dynamically stable configuration without showing chaotic or unstable behavior under the given parameter conditions.

4.1. Sensitivity and stability analysis

To further investigate the dynamical behavior of the spring–pendulum system, a stability and sensitivity analysis was carried out for both the conservative and damped cases. The phase portraits presented in Figures 5, 6, 19 and 20 provide important information regarding the stability of the motion. In the conservative case, the trajectories form closed bounded curves around the equilibrium position, which indicates stable periodic oscillations of the system. The absence of diverging trajectories confirms that the nonlinear coupling between the radial and angular motions does not lead to instability within the considered parameter range.

For the damped system, the phase trajectories gradually spiral toward the equilibrium point. This behavior is expected due to the presence of viscous damping, which continuously removes energy from the system. As time increases, the amplitudes of both the radial and angular oscillations decrease until the motion approaches a stable equilibrium configuration. Therefore, the damped system can be considered asymptotically stable.

A sensitivity analysis was also performed by introducing small perturbations in the initial conditions and damping coefficients. The numerical results showed that moderate variations in these parameters produce only small changes in the displacement and velocity responses. Although the nonlinear interaction between the radial and angular motions slightly modifies the oscillation amplitudes and frequencies, the qualitative behavior of the system remains unchanged. No unstable or chaotic behavior was observed for the selected parameters and initial conditions.

These observations demonstrate that the proposed spring–pendulum model possesses stable nonlinear dynamics and that the Ms-VIM method remains reliable and accurate even when small perturbations are introduced into the system parameters.

5. Conclusion

This research has been based on development of a numerical/analytical solution technique which allows the investigation into the nonlinear dynamics of a spring pendulum of constant mass. A common application for this type of device is found in engineering examples including; cranes, robotic manipulator arms and vibration control devices. Using the Lagrangian variational formulation the authors have formulated the governing equations describing the coupled motion of both the radially and angularly moving components of the device. The results show that due to the nonlinear coupling between the two modes of motion there are amplitude modulation and nonlinear frequency shifts that can occur in the motion of the system and can not be modeled by the classical linear approximation techniques.

The time domain response and the phase plane trajectories of the system indicate a stable oscillating behavior over the entire range of parameters examined and the bounded and closed phase plane trajectories indicate the dynamic stability of the system. The Ms-VIM was able to provide reliable semi-analytical solutions for the system at a significant reduction in computational cost compared to other numerical methods including the RK4 method. Due to the fact that the authors have validated their method against the RK4 method they have shown that the proposed method will provide a useful tool for engineers when designing and modeling nonlinear spring-pendulum systems. Nevertheless, the Ms-VIM has some limitations: it requires careful selection of subinterval size to balance accuracy and computational cost, and the derivation of Lagrange multipliers may become difficult for highly complex nonlinear systems.

The Ms-VIM can be extended to multi-degree-of-freedom systems and parametrically excited oscillators, where classical perturbation methods often struggle. Future work may also explore adaptive subinterval sizing for improved computational efficiency.

Another promising extension is to formulate the spring-pendulum system in fractal or fractional space, which can capture memory effects and complex damping mechanisms more accurately. Recent advances in fractal calculus and fractional modelling (see e.g., recent issues of Fractals and Chaos, Solitons & Fractals) offer powerful tools for such investigations. This remains an important direction for future work.

Footnotes

Acknowledgment

The authors Rania Wannan, and Jihad Asad would like to thank Palestine Technical University- Kadoorie for supporting them financially during this research.

ORCID iD

Jihad Asad

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Appendix

References

Cveticanin

Zukovic

. Oscillator with Time-Variable Degree of Nonlinearity: Theory and Application in Aging of Polymer Composite Structure. Mathematics 2023; 11(18): 3958. https://doi.org/10.3390/math11183958

José

Saletan

. Classical dynamics: a contemporary approach. Cambridge University Press, 1998.

Nayfeh

Mook

. Nonlinear oscillations. Wiley, 1979. https://doi.org/10.1002/9783527617586

Kane

Levinson

. Dynamics: Theory and applications. McGraw-Hill, 1985.

Irschik

Belyaev

(eds). Dynamics of mechanical systems with variable mass: Springer Vienna, 2014, vol 557. https://doi.org/10.1007/978-3-7091-1809-2

Pesce

Casetta

. Variable mass systems dynamics in engineering mechanics education. Proceedings of COBEM, 19th International Congress of Mechanical Engineering, Brasilia, 2007.

Musielak

. Standard and nonstandard Lagrangians for dissipative systems. Journal of Physics A 2008; 41: 055205. https://doi.org/10.1088/1751-8113/41/5/055205

Riewe

. Nonconservative Lagrangian and Hamiltonian mechanics. Physical Review E 1996; 53: 1890–1899. https://doi.org/10.1103/PhysRevE.53.1890

Plastino

Muzzio

. On the use and abuse of Newton’s second law for variable mass problems. Celestial Mechanics and Dynamical Astronomy 1992; 53: 227–232. https://doi.org/10.1007/BF00052611

10.

Zhou

Zhang

, et al. Analysis on the motion of nonlinear vibration with fractional order and time variable mass. Applied Mathematics Letters 2022; 124: 107621. https://doi.org/10.1016/j.aml.2021.107621

11.

Lanczos

. The variational principles of mechanics. Courier Corporation 2012.

12.

Esmailzadeh

Younesian

Askari

(2018). Analytical methods in nonlinear oscillations: Springer.

13.

Nayfeh

. Perturbation methods. John Wiley & Sons, 2024.

14.

. Variational iteration method–a kind of non-linear analytical technique: some examples. International journal of non-linear mechanics 1999; 34(4): 699–708. https://doi.org/10.1016/S0020-7462(98)00048-1

15.

. Homotopy perturbation method: A new nonlinear analytical technique. Applied Mathematics and Computation 2003; 135: 73–79. https://doi.org/10.1016/S0096-3003(01)00312-5

16.

Dubey

Chakraverty

. Solution of fractional wave equation by homotopy perturbation method. Wave dynamics 2022: 263–277.

17.

Lima

. Simple but accurate periodic solutions for the nonlinear pendulum equation. Revista Brasileira de Ensino de Física 2019; 41. https://doi.org/10.1590/1806-9126-RBEF-2018-0202

18.

Liao

. Homotopy analysis method in nonlinear differential equations: Higher education press, 2012, vol 153.

19.

Butcher

. Numerical methods for ordinary differential equations. 3rd ed. Wiley, 2016. https://doi.org/10.1002/9781119121534

20.

Odibat

Momani

. Modified homotopy perturbation method: application to quadratic Riccati differential equation of fractional order. Chaos, Solitons & Fractals 2008; 36(1): 167–174. https://doi.org/10.1016/j.chaos.2006.06.041

21.

. Variational iteration method: new development and applications. Computers & Mathematics with Applications 2007; 54(7-8): 881–894. https://doi.org/10.1016/j.camwa.2006.12.083

22.

Inokuti

. General use of the Lagrange multiplier in nonlinear mathematical physics. Variational method in the mechanics of solids, 1978, pp. 156–162.

23.

Cheng

. Numerical manifold method based on the method of weighted residuals. Computational Mechanics 2005; 35(6): 470–480. https://doi.org/10.1007/s00466-004-0636-3

24.

. Variational iteration method—some recent results and new interpretations. Journal of computational and applied mathematics 2007; 207(1): 3–17. https://doi.org/10.1016/j.cam.2006.07.009

25.

Saleh

Azmi

Jameel

. Computational Solution and Analysis of Fuzzy Nonlinear HIV Infection Model via New Multistage Fuzzy Variational Iteration Method. ms 2024; 12: 219–228. https://doi.org/10.13189/ms.2024.120301

26.

Ertürk

. Multistage variational iteration method for delay differential equations. International Journal of Nonlinear Sciences and Numerical Simulation 2013; 14, 159–166. https://doi.org/10.1515/ijnsns-2012-0165

27.

Gökdoğan

Merdan

Ertürk

. Multistage variational iteration method for epidemic models. Kuwait Journal of Science 2012; 39: 57–67.

28.

Anjum

. A dual Lagrange multiplier approach for the dynamics of the mechanical systems. Journal of Applied and Computational Mechanics 2024; 10(3): 643–651. https://doi.org/10.22055/jacm.2024.45944.4439

29.

Amer

Galal

. Dynamical analysis of a forced vibrating planar motion of a spring pendulum. Journal of Low Frequency Noise, Vibration and Active Control 2024; 43(4): 1545–1566. https://doi.org/10.1177/14613484241273636

30.

Amer

El-Kafly

. Dynamical analysis for the motion of a 2DOF spring pendulum on a Lissajous curve. Scientific Reports 2023; 13: 21430. https://doi.org/10.1038/s41598-023-48523-5

31.

Rouleau

Shoshani

Elata

. The elastic pendulum: Dzhanibekov-like transitions of symmetric and asymmetric periodic responses. Nonlinear Dynamics 2025; 113: 20841–20853. https://doi.org/10.1007/s11071-025-11242-9

32.

Amer

Elhagary

Hassan

, et al. Examining the behavior of a two-connected-spring pendulum according to non-resonance and resonance conditions. Journal of Low Frequency Noise, Vibration and Active Control 2025; 44: 866–892. https://doi.org/10.1177/14613484241302267

33.

Jin

Shen

, et al. Numerical analysis of a fractional rotating pendulum oscillator. Fractals 2026; 34: 2650064. https://doi.org/10.1142/S0218348X26500647

34.

Zheng

Lou

Shen

, et al. Numerical investigation of fractal oscillator for a pendulum with a rolling wheel. Fractals 2025; 33(09): 2550077. https://doi.org/10.1142/S0218348X2550077X

35.

Lou

Chen

, et al. Analytical approaches to fractal oscillator for a pendulum with a rolling wheel. Journal of Low Frequency Noise, Vibration and Active Control 2025: 14613484251413080. https://doi.org/10.1177/14613484251413080

36.

Landau

Lifshitz

. Mechanics. 3rd ed: Pergamon Press, 1976.

37.

J-H

. Variational iteration method—A kind of nonlinear analytical technique: Some examples. International Journal of Non-Linear Mechanics 1999; 34(4): 699–708. https://doi.org/10.1016/S0020-7462(98)00048-1

38.

Oden

Reddy

. Variational methods in theoretical mechanics. Springer Science & Business Media, 2012.

39.

Finlayson

. The method of weighted residuals and variational principles. Society for Industrial and Applied Mathematics, 2013.

40.

Wazwaz

. The variational iteration method for analytic treatment for linear and nonlinear ODEs. Applied Mathematics and Computation 2009; 212(1): 120–134. https://doi.org/10.1016/j.amc.2009.02.003

41.

Goldstein

Poole

Safko

. Classical mechanics. 3rd ed. Addison-Wesley, 2002.

42.

Akindeinde

. A new multistage technique for approximate analytical solution of nonlinear differential equations. Heliyon 2020; 6(10): e05188. https://doi.org/10.1016/j.heliyon.2020.e05188

43.

Gökdoğan

Merdan

Ertürk

. A multistage variational iteration method for approximate-analytic solution of avian-human influenza epidemic model. Kuwait Journal of Science 2012; 39(2A): 57–67.

44.

Gökdoğan

Merdan

Ertürk

. A Multistage Variational Iteration Method for Solution of Delay Differential Equations. International Journal of Nonlinear Sciences & Numerical Simulation 2013; 14: 159–166. https://doi.org/10.1515/ijnsns-2011-0180