Analytical and Numerical Approximations to Some Coupled Forced Damped Duffing Oscillators

Abstract

In this investigation, two different models for two coupled asymmetrical oscillators, known as, coupled forced damped Duffing oscillators (FDDOs) are reported. The first model of coupled FDDOs consists of a nonlinear forced damped Duffing oscillator (FDDO) with a linear oscillator, while the second model is composed of two nonlinear FDDOs. The Krylov–Bogoliubov–Mitropolsky (KBM) method, is carried out for analyzing the coupled FDDOs for any model. To do that, the coupled FDDOs are reduced to a decoupled system of two individual FDDOs using a suitable linear transformation. After that, the KBM method is implemented to find some approximations for both unforced and forced damped Duffing oscillators (DDOs). Furthermore, the KBM analytical approximations are compared with the fourth-order Runge–Kutta (RK4) numerical approximations to check the accuracy of all obtained approximations. Moreover, the RK4 numerical approximations to both coupling and decoupling systems of FDDOs are compared with each other.

1. Introduction

Nonlinear differential equations of the second- and third-order have attracted great attention from many researchers due to their many applications in various fields of science especially in engineering problems and in modeling robotics motion [1,2]. For instance, the Duffing equation (DE) is a second-order differential equation, which has spread widely due to its many applications in physics, chemistry, and in modeling many engineering problems [3,4]. It is typically regarded as a paradigm for nonlinear dynamics of dissipative systems. The Duffing equation (DE)/Duffing oscillator (DO) has been used successfully to describe and model a variety of physical and chemical processes, and several engineering problems, such as reinforcing springs, beam buckling, nonlinear electronic circuits, among others [3]. Given a physical system that has some form of potential energy associated with it, its dynamics in the neighborhood of a potential minimum can always be approximated by a linear oscillator. The oscillator restorative force coefficient incorporates the curvature of the potential well and is the first term of the Taylor series in which it can develop the position-dependent force acting on the already mentioned physical system. Taking more terms of the series improves the approximation. The next term of the Taylor series that brings about qualitative changes is the one that goes with the cube of the position (the quadratic term only implies a change in the equilibrium position). Narasimha [5] stated the general equations of the dynamics of a thin wire. The center point of a vibrating wire whose ends are embedded can be described as a one-dimensional effective oscillator. It was shown that for this oscillator the cubic term should not be neglected under resonance conditions. The micro clamped–clamped oscillator should share this characteristic with the vibrating wire as it is also a thin body tensed between anchors, therefore, its dynamics are governed by the following forced damped Duffing oscillator (FDDE)

$$m\ddot{x}+\gamma \dot{x}+kx+h{x}^3=fcos\omega t,$$

(1)

where m represents the pendulum mass, $\gamma$ gives the coefficient of the damping force (damping parameter), and k denotes the effective stiffness constant while f and $\omega$ indicate the amplitude and frequency of the driving force. The new term is $h{x}^3$, which represents the nonlinear relationship of the increase in the restoring force as the tension in the beams increases due to their elongation. Again, the equation is dimensionless. The unit of time is $1/{\omega }_0=\sqrt{m/k}$ and the position is $f/k$. There are now two free parameters, namely, $ϵ$ and $\eta$, defined by

$$\ddot{x}+ϵ\dot{x}+x+\eta x^3=cos\omega t,$$

(2)

with $ϵ=\gamma /\sqrt{mk}$ and $\eta =h{f}^2/k^3$, where $ϵ$ indicates the normalized friction coefficient and $\eta$ represents the coefficient of the nonlinear term. It is considered that these are the two parameters that define the system. We want to find the motion of the oscillator when it is excited with frequency $\omega$.

One of the desirable properties of a clock-like oscillator is that it has a stable frequency. The linear oscillator has a stable frequency near the resonance maximum, but the DO loses that property, particularly when fluctuations in amplitude affect the frequency of the system. The dependency of the frequency on the phase has a promising local maximum (zero slope), but according to the range in which frequency fluctuations can be tolerated, it may not be enough.

What can be interpreted from the measurements is that micro-oscillators behave as coupled Duffing oscillators (DOs) with another higher-frequency oscillator. In the current investigation, we propose a model consisting of a nonlinear oscillator that experiences a forced external coupling and a coupling due to the second oscillator. This second oscillator is linear and sustains its movement only by coupling with the main oscillator. The dynamics of this system can be modeled by the pair of equations

$$\left\{\begin{array}{c}{m}_1\ddot{x}+{\gamma }_1\dot{x}+k_{1}x+h{x}^3=fcos\omega t+J_{1}y,\\ m_2\ddot{y}+{\gamma }_2\dot{y}+k_{2}y=J_{2}x.\end{array}\right.$$

(3)

The two oscillators have different natural frequencies (given by the ratio $k_i/m_i={\omega }_i^2$). The width of their resonance curves is given by ${\gamma }_i$. The coupling, in principle, is asymmetric due to $J_1\ne J_2$. The system (3) may be rewritten in the following initial-value problem (i.v.p.) form, known as, coupled forced damped Duffing oscillators (FDDOs)

$$\left\{\begin{array}{c}{\mathcal{R}}_1\equiv \ddot{x}+2{\nu }_1\dot{x}+\alpha x+\beta x^3-\gamma cos\omega t-{ϵ}_{1}y=0,\hfill \\ {\mathcal{R}}_2\equiv \ddot{y}+2{\nu }_2\dot{y}+\mu y-{ϵ}_{2}x=0,\hfill \end{array}\right.$$

(4)

with the initial conditions (ICs)

$$\left\{\begin{array}{c}x\left(0\right)=x_0,y\left(0\right)=y_0,\\ \dot{x}\left(0\right)={\dot{x}}_0\mathrm{and}\dot{y}\left(0\right)={\dot{y}}_0,\end{array}\right.$$

(5)

where ${\nu }_1={\gamma }_1/\left(2m_1\right)$, $\alpha =k_1/m_1$, $\beta =h/m_1$, $\gamma =f/m_1$, ${ϵ}_1=J_1/m_1$, ${\nu }_2={\gamma }_2/\left(2m_2\right)$, $\mu =k_2/m_2$, and ${ϵ}_2=J_2/m_2$. We assume that the coupling is weak, that is, $\left|{ϵ}_1\right|\ll 1$ and $\left|{ϵ}_2\right|\ll 1$.

Furthermore, in our investigation, we try to find an explicit approximation to the following coupled FDDOs [6]

$$\left\{\begin{array}{c}{\mathbb{k}}_1\equiv \ddot{x}+2{\nu }_1\dot{x}+\alpha x+\beta x^3+k_1\left(y-x\right)-H_1\left(t\right)=0,\\ {\mathbb{k}}_2\equiv \ddot{y}+2{\nu }_2\dot{y}+\lambda y+\mu y^3+k_2\left(x-y\right)-H_2\left(t\right)=0,\end{array}\right.$$

(6)

with the ICs

$$\left\{\begin{array}{c}x\left(0\right)=x_0,y\left(0\right)=y_0,\\ \dot{x}\left(0\right)={\dot{x}}_0\mathrm{and}\dot{y}\left(0\right)={\dot{y}}_0,\end{array}\right.$$

(7)

where $\left(\epsilon ,\delta \right)$ indicate the damping coefficients of both first and second oscillators, $\left(H_1\left(t\right),H_2\left(t\right)\right)$ represent the driving forces, $\left(k_1,k_2\right)$ express the linear coupling stiffness parameters, $\left(\alpha ,\lambda \right)$ denote the linear stiffness parameters, and $\left(\beta ,\mu \right)$ represent the parameters of nonlinear terms. Furthermore, $H_{1,2}\left(t\right)$ represent any time-dependent functions.

Some studies have been conducted around coupled Duffing-type oscillators [7,8,9,10,11,12], but an explicit solution to the coupled Duffing-type oscillators is rarely found in the published papers. In our study, we use the Krylov–Bogoliubov–Mitropolsky (KBM) method [7,8,9] for analyzing the two models of coupled FDDO including the first model (4) and the second model (6). The proposed method is successful in analyzing many nonlinear oscillators such as a generalized Van der Pol oscillator, a pendulum with variable length, and a Rayleigh Equation [13]. Furthermore, this method as well as the ansatz method are applied to analyze the (un)forced pendulum–cart system oscillators [14] and many other oscillators [15,16,17,18,19,20,21]. All mentioned studies and many others prove the efficiency and accuracy of the KBM method in obtaining high-accuracy approximations to many strong nonlinear oscillators. Motivated by these investigations, we can use this method to find high-accurate approximations to many nonlinear coupled systems such as the two mentioned models of coupled oscillators (4) and (6).

2. Decoupling the Two Models of Coupled FDDOs

Before applying the KBM technique for analyzing the coupled systems (4) and (6), we first must decouple systems (4) and (6) to two individual FDDOs.

2.1. Decoupling the First Model (4) to Two Individual FDDOs

For decoupling system (4), the following linear solution to the system (4) is introduced

$$\left\{\begin{array}{c}x\left(t\right)=u+v,\hfill \\ y\left(t\right)=ru+sv,\hfill \end{array}\right.$$

(8)

where $\left(r,s\right)$ represent the transformation parameters and $u\equiv u\left(t\right)$ and $v\equiv v\left(t\right)$. We choose r and s so that the system (4) has the new form

$$\left\{\begin{array}{c}\ddot{u}=-\left(2\epsilon \dot{u}+Pu+R{u}^3\right)+{ϝ}_1\left(t\right),\\ \ddot{v}=-\left(2\delta \dot{v}+Qv+S{v}^3\right)+{ϝ}_2\left(t\right),\end{array}\right.$$

(9)

where ${ϝ}_1\left(t\right)=f_{1}cos\omega t$ and ${ϝ}_2\left(t\right)=f_{2}cos\omega t$ or any other time-dependent functions. Observe that this new system is linearly decoupled.

Inserting Equations (8) and (9) into the linearized form of system (4), i.e., in the following form

$$\left\{\begin{array}{c}\ddot{x}+2{\nu }_1\dot{x}=-\alpha x+{ϵ}_{1}y+H_1\left(t\right)\hfill \\ \ddot{y}+2{\nu }_2\dot{y}=-\mu y+{ϵ}_{2}x+H_2\left(t\right),\hfill \end{array}\right.$$

(10)

we have

$$\begin{array}{cc}\hfill & {\mathcal{R}}_1={ϝ}_1\left(t\right)+{ϝ}_2\left(t\right)-H_1\left(t\right)+\left(-P+\alpha +r{ϵ}_1\right)u\hfill \\ \hfill & +\left(-Q+\alpha +s{ϵ}_1\right)v-R{u}^3-S{v}^3+\beta {\left(u+v\right)}^3\hfill \\ \hfill & +2\left(-\epsilon +{\nu }_1\right)\dot{u}+2\left(-\delta +{\nu }_1\right)\dot{v}=0,\hfill \end{array}$$

(11)

and

$$\begin{array}{cc}{\mathcal{R}}_2=\hfill & r{ϝ}_1\left(t\right)+s{ϝ}_2\left(t\right)-H_2\left(t\right)+\left(-{ϵ}_2+\mu s-sQ\right)v\hfill \\ \hfill & +\left(-{ϵ}_2+\mu r-Pr\right)u-sS{v}^3-rR{u}^3\hfill \\ \hfill & +2r\left(-\epsilon +{\nu }_2\right)\dot{u}+2s\left({\nu }_2-\delta \right)\dot{v}=0,\hfill \end{array}$$

(12)

where $H_1\left(t\right)=\gamma cos\omega t$ and $H_2\left(t\right)=0$ as defined in Equation (4). After vanishing the coefficients of $u^0$, $v^0$, u, v, $u^3$, $v^3$, $\dot{u}$, and $\dot{v}$ in Equations (11) and (12), we have

$$\left\{\begin{array}{c}{ϝ}_1\left(t\right)+{ϝ}_2\left(t\right)-H_1\left(t\right)=0,\hfill \\ -H_2\left(t\right)+r{ϝ}_1\left(t\right)+s{ϝ}_2\left(t\right)=0,\hfill \\ \left(-P+\alpha +r{ϵ}_1\right)=0,\hfill \\ \left(-Q+\alpha +s{ϵ}_1\right)=0,\hfill \\ \left(-2\epsilon +2{\nu }_1\right)=0,\hfill \\ \left(2{\nu }_{2}s+-2\delta s\right)=0,\hfill \\ \left(-{ϵ}_2+\mu s-sQ\right)=0,\hfill \\ \left(-{ϵ}_2+\mu r-Pr\right)=0.\hfill \end{array}\right.$$

(13)

By solving system (13), we finally get

$$\left\{\begin{array}{c}\epsilon ={\nu }_1,\delta ={\nu }_2,\hfill \\ R=S=\beta ,\hfill \\ P=\alpha +r{ϵ}_1,\hfill \\ Q=\alpha +s{ϵ}_1,\hfill \\ {ϝ}_1\left(t\right)=\frac{-s{H}_1\left(t\right)+H_2\left(t\right)}{r-s},\hfill \\ {ϝ}_2\left(t\right)=\frac{r{H}_1\left(t\right)-H_2\left(t\right)}{r-s},\hfill \end{array}\right.$$

(14)

and

$$\left\{\begin{array}{c}r=\frac{-\alpha +\sqrt{\Delta }+\mu }{2{ϵ}_1},\\ s=\frac{-\alpha -\sqrt{\Delta }-\mu }{2{ϵ}_1}.\end{array}\right.$$

(15)

Note that the values of $\left(r,s\right)$ given in Equation (15) are valid only for $\Delta ={(\alpha -\mu )}^2-4{ϵ}_1{ϵ}_2>0.$

The original coupled system of FDDOs (4) is decoupled to two individual FDDOs as shown in Equation (9) with some residual errors and with the following new ICs

$$\left\{\begin{array}{c}u\left(0\right)=\frac{y_0-s{x}_0}{r-s},v\left(0\right)=\frac{r{x}_0-y_0}{r-s},\\ \dot{u}\left(0\right)=\frac{{\dot{y}}_0-s{\dot{x}}_0}{r-s}\mathrm{and}\dot{v}\left(0\right)=\frac{r{\dot{x}}_0-{\dot{y}}_0}{r-s}.\end{array}\right.$$

(16)

The residual errors for this case read

$$\left\{\begin{array}{c}{\mathcal{R}}_1=\beta \left(3u^{2}v+3v^{2}u\right)+2\left({\nu }_1-\delta \right)\dot{v},\\ {\mathcal{R}}_2=-sS{v}^3-rR{u}^3+2r\left({\nu }_2-\epsilon \right)\dot{u}.\end{array}\right.$$

Now, it is enough to solve any one of the FDDOs (9)

$$\left\{\begin{array}{c}\ddot{w}+2\kappa \dot{w}+{\mathsf{\Omega }}^{2}w+\sigma w^3=fcos\omega t,\\ w\left(0\right)=w_0\mathrm{and}\dot{w}\left(0\right)={\dot{w}}_0,\end{array}\right.$$

(17)

for arbitrary parameters to find approximations to the coupled system (4).

2.2. Decoupling the Second Model (6) to Two Individual FDDOs

Let us consider the following alternative form for the two coupled FDDOs

$$\left\{\begin{array}{c}{\mathbb{k}}_1\equiv \ddot{x}+2{\nu }_1\dot{x}+\alpha x+\beta x^3+k_1\left(y-x\right)-H_1\left(t\right)=0,\\ {\mathbb{k}}_2\equiv \ddot{y}+2{\nu }_2\dot{y}+\lambda y+\mu y^3+k_2\left(x-y\right)-H_2\left(t\right)=0.\end{array}\right.$$

(18)

For decoupling system (18), the following linear solutions are assumed

$$\left\{\begin{array}{c}x\left(t\right)=u+v\\ y\left(t\right)=ru+sv,\end{array}\right.$$

(19)

where $u\equiv u\left(t\right)$ and $v\equiv v\left(t\right)$ follow the FDDOs as

$$\left\{\begin{array}{c}\ddot{u}=-\left(2\epsilon \dot{u}+Pu+R{u}^3\right)+{ϝ}_1\left(t\right),\\ \ddot{v}=-\left(2\delta \dot{v}+Qv+S{v}^3\right)+{ϝ}_2\left(t\right).\end{array}\right.$$

(20)

Here, $\left(r,s\right)$ represent the transformation parameters and $\left(\epsilon ,\delta ,P,Q,R,S,{ϝ}_1\left(t\right),{ϝ}_2\left(t\right)\right)$ are undermined parameters.

Inserting Equation (19) into system (18) and using system (20), we finally obtain

$$\begin{array}{c}{\mathbb{k}}_1={ϝ}_1\left(t\right)+{ϝ}_2\left(t\right)-H_1\left(t\right)+\left(\alpha +k_{1}r-k_1-P\right)u\hfill \\ +\left(\alpha +k_{1}s-k_1-Q\right)v+\left(-S+\beta \right)v^3+\left(-R+\beta \right)u^3\hfill \\ +\beta \left(3v{u}^2+3u{v}^2\right)+2({\nu }_1-\delta )\dot{v}+2({\nu }_1-\epsilon )\dot{u},\hfill \end{array}$$

(21)

and

$$\begin{array}{c}{\mathbb{k}}_2=r{ϝ}_1\left(t\right)+s{ϝ}_2\left(t\right)-H_2\left(t\right)+\left(k_2-k_{2}r+\lambda r-rP\right)u\hfill \\ +\left(\lambda s-k_{2}s+k_2-sQ\right)v+r\left(r^2\mu -R\right)u^3+s\left(s^2\mu -S\right)v^3\hfill \\ +\mu \left(3r^2{u}^{2}sv+3s^2{v}^{2}ru\right)+2r({\nu }_2-\epsilon )\dot{u}+2s({\nu }_2-\delta )\dot{v}.\hfill \end{array}$$

(22)

Now, after vanishing the coefficients of $u^0$, $v^0$, u, v, $u^3$, $v^3$, $\dot{u}$, and $\dot{v}$ in Equations (21) and (22), we obtain

$$\left\{\begin{array}{c}\epsilon ={\nu }_1,\delta ={\nu }_2,\hfill \\ R=\beta ,S=\mu s^2,\hfill \\ P=-k_1+r{k}_1+\alpha ,\hfill \\ Q=-k_1+s{k}_1+\alpha ,\hfill \\ {ϝ}_1\left(t\right)=\frac{-s{H}_1\left(t\right)+H_2\left(t\right)}{r-s},\hfill \\ {ϝ}_2\left(t\right)=\frac{r{H}_1\left(t\right)-H_2\left(t\right)}{r-s},\hfill \end{array}\right.$$

(23)

and

$$\left\{\begin{array}{c}r=\frac{\lambda +k_1-k_2-\alpha +\sqrt{\Delta }}{2k_1},\hfill \\ s=\frac{\lambda +k_1-k_2-\alpha -\sqrt{\Delta }}{2k_1}.\hfill \end{array}\right.$$

(24)

Note that the values of $\left(r,s\right)$ given in Equation (24) are valid only for $\Delta =\left(\lambda +k_1-k_2-\alpha \right){}^2+4k_1{k}_2>0.$

Finally, the original coupled system of FDDOs (18) is decoupled to individual two FDDOs as shown in Equation (20) with some residual errors and with the following new ICs

$$\left\{\begin{array}{c}u\left(0\right)=\frac{y_0-s{x}_0}{r-s},v\left(0\right)=\frac{r{x}_0-y_0}{r-s},\\ \dot{u}\left(0\right)=\frac{{\dot{y}}_0-s{\dot{x}}_0}{r-s}\mathrm{and}\dot{v}\left(0\right)=\frac{r{\dot{x}}_0-{\dot{y}}_0}{r-s}.\end{array}\right.$$

(25)

Accordingly, the residual errors read

$$\left\{\begin{array}{c}{\mathbb{k}}_1=\left(-S+\beta \right)v^3+\beta \left(3v{u}^2+3u{v}^2\right)+2({\nu }_1-\delta )\dot{v},\\ {\mathbb{k}}_2=r\left(r^2\mu -R\right)u^3+\mu \left(3r^2{u}^{2}sv+3s^2{v}^{2}ru\right)+2r({\nu }_2-\epsilon )\dot{u},\end{array}\right.$$

(26)

Thus, it is enough to solve the FDDO (17) for arbitrary parameters to find some approximations to the coupled system (18).

3. KBM Technique for Analyzing the FDDO (17)

We first solve the following unforced damped DE:

$$\left\{\begin{array}{c}\ddot{w}+2\kappa \dot{w}+{\mathsf{\Omega }}^{2}w+\sigma w^3=0,\\ w\left(0\right)=w_0\mathrm{and}\dot{w}\left(0\right)={\dot{w}}_0.\end{array}\right.$$

(27)

According to the KBM method, the following perturbed problem is constructed:

$$\left\{\begin{array}{c}{\mathbb{N}}_1\equiv \ddot{w}+{\mathsf{\Omega }}^{2}w+p\left(2\kappa \dot{w}+\sigma w^3\right)=0,\\ w\left(0\right)=w_0\mathrm{and}\dot{w}\left(0\right)={\dot{w}}_0,\end{array}\right.$$

(28)

where p indicates the perturbed parameter $0<p<<1$.

Assume the solution to the i.v.p. (28) is given by the following ansatz form

$$w\left(t\right)=c_{1}acos\psi +p\left(d_1\left(a\right)cos3\psi +c_{2}asin\psi \right),$$

(29)

where the functions $a\equiv a\left(t\right)$ and $\psi \equiv \psi \left(t\right)$ are defined by

$$\left\{\begin{array}{c}\dot{a}=pA\left(a\right),\\ \dot{\psi }=\mathsf{\Omega }+p{\mathsf{\Phi }}_1\left(a\right).\end{array}\right.$$

(30)

Inserting solution (29) into the i.v.p. (28) yields

$$\begin{array}{c}{\mathbb{N}}_1=\frac{p}{4}[\left(a^3{c}_1^3-32{\mathsf{\Omega }}^2{d}_1\left(a\right)\right)cos\left(3\psi \right)\hfill \\ +a{c}_1\left(3a^2{c}_1^2\sigma -8\mathsf{\Omega }{\mathsf{\Phi }}_1\left(a\right)\right)cos\left(\psi \right)\hfill \\ -8c_1\mathsf{\Omega }sin\left(\psi \right)(A\left(a\right)+a\kappa )\left)\right].\hfill \end{array}$$

(31)

Equating to zero the coefficients of $sinj\psi$ and $cosj\psi$ ($j=0,1,3,...$), an algebraic system is obtained and by solving this system, we have

$$\begin{array}{cc}{d}_1\left(a\right)\hfill & =\frac{a^3{c}_1^3}{32{\mathsf{\Omega }}^2},\hfill \\ {\mathsf{\Phi }}_1\left(a\right)\hfill & =\frac{3a^2{c}_1^2\sigma }{8\mathsf{\Omega }},\hfill \\ A\left(a\right)\hfill & =-\kappa a.\hfill \end{array}$$

(32)

Inserting the values of $\left(a,\psi \right)$ given in Equation (32) into Equation (30) and solving the obtained equations using the ICs $a\left(0\right)=1$ and $\psi \left(0\right)=0$, we get

$$\left\{\begin{array}{c}a=e^{-p\kappa t},\\ \psi =\mathsf{\Omega }t+\frac{3c_1^2\sigma }{16p\kappa \mathsf{\Omega }}\left(1-e^{-2p\kappa t}\right).\end{array}\right.$$

(33)

For $p=1$, the following approximate solution (first-order approximation) to the unforced damped DE (27) is obtained

$$\begin{array}{cc}w\left(t\right)\hfill & =c_1{e}^{-\kappa t}cos\left[\mathsf{\Omega }t+\frac{3c_1^2\sigma }{16\kappa \mathsf{\Omega }}\left(1-e^{-2\kappa t}\right)\right]\hfill \\ \hfill & +c_2{e}^{-\kappa t}sin\left[\mathsf{\Omega }t+\frac{3c_1^2\sigma }{16\kappa \mathsf{\Omega }}\left(1-e^{-2\kappa t}\right)\right]\hfill \\ \hfill & +\frac{a^3{c}_1^3}{32{\mathsf{\Omega }}^2}cos\left[3\left(\mathsf{\Omega }t+\frac{3c_1^2\sigma }{16\kappa \mathsf{\Omega }}\left(1-e^{-2\kappa t}\right)\right)\right].\hfill \end{array}$$

(34)

Now, by applying the ICs given in the i.v.p. (27), the values of $c_{1,2}$ are obtained

$$\begin{array}{cc}{c}_2\hfill & =\frac{3c_1^3\sigma \kappa +32c_1\kappa {\mathsf{\Omega }}^2+32{\dot{w}}_0{\mathsf{\Omega }}^2}{12c_1^2\sigma \mathsf{\Omega }+32{\mathsf{\Omega }}^3},\hfill \\ \hfill & \frac{\sigma }{32{\mathsf{\Omega }}^2}{c}_1^3+c_1-w_0=0.\hfill \end{array}$$

(35)

Remark 1.

We have the following third-order approximation to the i.v.p. (28)

$$\begin{array}{c}w\left(t\right)=a{c}_{1}cos\left(\psi \right)+\frac{1}{3}(p+p^2+p^3)a{c}_{2}sin\left(\psi \right)\hfill \\ +p{W}_1+p^2{W}_2+p^3{W}_3,\hfill \end{array}$$

(36)

with

$$\begin{array}{cc}\dot{a}\hfill & =-pa\kappa +p^2{a}^3\frac{3c_1^2\sigma \kappa }{8{\mathsf{\Omega }}^2}-p^3{a}^5\frac{195c_1^4{\sigma }^2\kappa }{512{\mathsf{\Omega }}^4},\hfill \\ \dot{\psi }\hfill & =\mathsf{\Omega }+p{a}^2\frac{3c_1^2\sigma }{8\mathsf{\Omega }}-p^2{a}^4\frac{15c_1^4{\sigma }^2+128{\kappa }^2{\mathsf{\Omega }}^2}{256{\mathsf{\Omega }}^3}\hfill \\ \hfill & +p^3\frac{a^2\left(369a^4{c}_1^6{\sigma }^3-4608c_1^2\sigma {\kappa }^2{\mathsf{\Omega }}^2+1024c_2^2\sigma {\mathsf{\Omega }}^4\right)}{24576{\mathsf{\Omega }}^5},\hfill \end{array}$$

where the values of coefficients $W_i$$\left(i=1,2,3\right)$ are given in Appendix A.

It remains to solve the FDDO (17). For this purpose, we assume the solution is defined by the ansatz form

$$w\left(t\right)=z\left(t\right)+e_{1}cos\omega t+e_{2}sin\omega t+e_{3}cos3\omega t+e_{4}sin3\omega t,$$

(37)

where the function $z\equiv z\left(t\right)$ represents a solution to the unforced damped DE

$$\left\{\begin{array}{c}\ddot{z}+2\nu \dot{z}+{\mathsf{\Omega }}^{2}z+\sigma z^3=0,\\ z\left(0\right)=w_0-e_1\mathrm{and}\dot{z}\left(0\right)={\dot{w}}_0-\omega e_2,\end{array}\right.$$

(38)

where $e_i$$\left(i=1,2,3,4\right)$ are undetermined constants.

Inserting solution (37) into (17) yields

$$\begin{array}{c}\ddot{w}+2\kappa \dot{w}+{\mathsf{\Omega }}^{2}w+\sigma w^3-fcos\omega t=\\ S_{1}cos\left(\omega t\right)+S_{2}sin\left(\omega t\right)+S_{3}cos\left(3\omega t\right)+S_{4}sin\left(3\omega t\right)+\mathrm{h}.\mathrm{o}.\mathrm{t},\end{array}$$

(39)

where the values of the coefficients $S_i$$\left(i=1,2,3,4\right)$ are defined in Appendix B.

For $S_i=0$, we get

$$\left\{\begin{array}{c}8e_2\nu \omega +3e_1^3\sigma +3e_3{e}_1^2\sigma +3e_2^2{e}_1\sigma +6e_3^2{e}_1\sigma +6e_4^2{e}_1\sigma +6e_2{e}_4{e}_1\sigma -3e_2^2{e}_3\sigma -4e_1{\omega }^2+4e_1{\mathsf{\Omega }}^2-4f=0,\hfill \\ -2e_1\left(3e_2{e}_3\sigma +4\nu \omega \right)+e_2\left(3e_2^2\sigma +6e_3^2\sigma +6e_4^2\sigma -3e_2{e}_4\sigma -4{\omega }^2+4{\mathsf{\Omega }}^2\right)+3\left(e_2+e_4\right)e_1^2\sigma =0,\hfill \\ 24e_4\nu \omega +e_1^3\sigma +6e_3{e}_1^2\sigma -3e_2^2{e}_1\sigma +3e_3^3\sigma +3e_3{e}_4^2\sigma +6e_2^2{e}_3\sigma -36e_3{\omega }^2+4e_3{\mathsf{\Omega }}^2=0,\hfill \\ -24e_3\nu \omega -e_2^3\sigma +6e_4{e}_2^2\sigma +3e_4^3\sigma +3e_3^2{e}_4\sigma +3e_1^2\left(e_2+2e_4\right)\sigma -36e_4{\omega }^2+4e_4{\mathsf{\Omega }}^2=0.\hfill \end{array}\right.$$

(40)

By solving system (40), the values of coefficients $e_i$ can be obtained.

Example 1.

Let us now consider the following numerical example to the coupled system of FDDOs (42)

$$\left\{\begin{array}{c}\stackrel{..}{x}+0.2\dot{x}+3x+2x^3=0.1cos\left(t\right)+0.1y\hfill \\ \stackrel{..}{y}+0.2\dot{y}+2y=-0.1x,\hfill \\ x\left(0\right)=\dot{x}\left(0\right)=0and y\left(0\right)=\dot{y}\left(0\right)=0.\hfill \end{array}\right.$$

(41)

The decoupling system for system (41) reads

$$\left\{\begin{array}{c}x\left(t\right)=u+v,\hfill \\ y\left(t\right)=0.0990195u-10.099v,\hfill \end{array}\right.$$

(42)

with

$$\left\{\begin{array}{c}\ddot{u}+0.2\dot{u}+3.0099u+2u^3=0.099029cos\left(t\right),\\ \ddot{v}+0.2\dot{v}+1.9901v+2v^3=0.000970966cos\left(t\right),\\ u\left(0\right)=\dot{u}\left(0\right)=0\mathrm{and}v\left(0\right)=\dot{v}\left(0\right)=0.\end{array}\right.$$

The KBM explicit-form approximations to the decoupled system (41) read

$$\begin{array}{cc}x\hfill & =0.0050279sin\left(t\right)+1.90864\times {10}^{-6}sin\left(3t\right)\hfill \\ \hfill & +0.0496444cos\left(t\right)+9.55958\times {10}^{-6}cos\left(3t\right)\hfill \\ \hfill & -1.19996\times {10}^{-6}{e}^{-0.3t}cos\left({\chi }_3\right)\hfill \\ \hfill & -e^{-0.1t}\left(\begin{array}{c}0.000201709sin\left({\chi }_1\right)+0.00559604sin\left({\chi }_2\right)\\ 0.000942229cos\left({\chi }_1\right)+0.0487105cos\left({\chi }_2\right)\end{array}\right),\hfill \end{array}$$

(43)

and

$$\begin{array}{cc}y\hfill & =-0.00144313sin\left(t\right)+1.88989\times {10}^{-7}sin\left(3t\right)\hfill \\ \hfill & -0.00469313cos\left(t\right)+9.46579\times {10}^{-7}cos\left(3t\right)\hfill \\ \hfill & -1.18819\times {10}^{-7}{e}^{-0.3t}cos\left({\chi }_3\right)\hfill \\ \hfill & +e^{-0.1t}\left(\begin{array}{c}0.00203707sin\left({\chi }_1\right)-0.000554117sin\left({\chi }_2\right)\\ +0.00951559cos\left({\chi }_1\right)-0.00482329cos\left({\chi }_2\right)\end{array}\right),\hfill \end{array}$$

(44)

with

$$\begin{array}{cc}{\chi }_1\hfill & =1.41071t-2.35997\times {10}^{-6}\left(e^{-0.2t}-1\right),\hfill \\ {\chi }_2\hfill & =1.73491t-0.00512861\left(e^{-0.2t}-1\right),\hfill \\ {\chi }_3\hfill & =5.20472t-0.0153858e^{-0.2t}+0.0153858.\hfill \end{array}$$

Now, we can apply two methodologies for analyzing the coupled system of FDDOs (41). In the first methodology, we can use the KBM analytical approximation (34) for studying the properties of the decoupled system oscillators (42) and then make a comparison among the obtained results and the RK4 approximations to the coupling system (41) as illustrated in Figure 1. In the second methodology, we can directly analyze the decoupling system of oscillators (42) using the RK4 approximations (here we do not use approximation (34)) and compare the numerical results with the RK4 approximations to the coupling system (41) as shown in Figure 2. Furthermore, the maximum residual distance error for the two methodologies is estimated as follows

$$\left\{\begin{array}{c}{L}_x=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(41\right)}-\mathrm{KBM}\mathrm{approxi}.\right|=0.00196831,\\ L_y=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(41\right)}-\mathrm{KBM}\mathrm{approxi}.\right|=0.0002707,\\ L_x=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(41\right)}-{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(42\right)}\right|=0.00216587,\\ L_y=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(41\right)}-{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(42\right)}\right|=0.000361187.\end{array}\right.$$

(45)

It is noted from both Figure 1 and Figure 2 and the maximum error given in Equation (45) that both analytical and numerical approximations are compatible with each other, which confirms the high accuracy of the analytical solutions.

Example 2.

Let us now consider the following numerical example to the coupled system of FDDOs (6)

$$\left\{\begin{array}{c}\ddot{x}+0.2\dot{x}+3x+x^3+0.1\left(y-x\right)-0.1cos\left(0.2t\right)=0,\\ \ddot{y}+0.2\dot{y}+4y+1y^3+0.3\left(x-y\right)-0.2cos\left(0.2t\right)=0,\\ x\left(0\right)=\dot{x}\left(0\right)=0\mathrm{and}y\left(0\right)=\dot{y}\left(0\right)=0.\end{array}\right.$$

(46)

The decoupling system for system (46) reads

$$\left\{\begin{array}{c}x\left(t\right)=u+v,\hfill \\ y\left(t\right)=8.3589u-0.358899v,\hfill \end{array}\right.$$

(47)

with

$$\left\{\begin{array}{c}\ddot{u}+0.2\dot{u}+3.73589u+u^3=0.0270584cos\left(0.2t\right),\\ \ddot{v}+0.2\dot{v}+2.86411v+0.128808v^3=0.0729416cos\left(0.2t\right),\\ u\left(0\right)=\dot{u}\left(0\right)=0\mathrm{and}v\left(0\right)=\dot{v}\left(0\right)=0.\end{array}\right.$$

The KBM explicit form approximations to decoupled system (46) read

$$\begin{array}{cc}x\hfill & =0.000444959sin\left(0.2t\right)-2.19375\times {10}^{-8}sin\left(0.6t\right)\hfill \\ \hfill & +0.0331427cos\left(0.2t\right)-2.49289\times {10}^{-7}cos\left(0.6t\right)\hfill \\ \hfill & +e^{-0.3t}\left(-1.87858\times {10}^{-7}cos\left({\eta }_3\right)-3.28121\times {10}^{-9}cos\left({\eta }_4\right)\right)\hfill \\ \hfill & +e^{-0.1t}\left(\begin{array}{c}-0.00156899sin\left({\eta }_1\right)-0.000386924sin\left({\eta }_2\right)\\ -0.025822cos\left({\eta }_1\right)-0.00732025cos\left({\eta }_2\right)\end{array}\right),\hfill \end{array}$$

(48)

and

$$\begin{array}{cc}y\hfill & =0.000530975sin\left(0.2t\right)-9.32518\times {10}^{-9}sin\left(0.6t\right)\hfill \\ \hfill & +0.0519219cos\left(0.2t\right)-1.6307\times {10}^{-7}cos\left(0.6t\right)\hfill \\ \hfill & +e^{-0.3t}\left(6.7422\times {10}^{-8}cos\left({\eta }_3\right)-2.74273\times {10}^{-8}cos\left({\eta }_4\right)\right)\hfill \\ \hfill & +e^{-0.1t}\left(\begin{array}{c}0.000563108sin\left({\eta }_1\right)-0.00323426sin\left({\eta }_2\right)\\ +0.00926748cos\left({\eta }_1\right)-0.0611892cos\left({\eta }_2\right)\end{array}\right),\hfill \end{array}$$

(49)

with

$$\begin{array}{cc}{\eta }_1\hfill & =1.69237t-0.0000951546\left(e^{-0.2t}-1\right),\hfill \\ {\eta }_2\hfill & =1.93285t-0.0000519824\left(e^{-0.2t}-1\right),\hfill \\ {\eta }_3\hfill & =5.0771t-0.000285464\left(e^{-0.2t}-1\right),\hfill \\ {\eta }_4\hfill & =5.79854t-0.000155947\left(e^{-0.2t}-1\right).\hfill \end{array}$$

Both analytical approximations using the KBM method for the decoupling system (47) and the RK4 numerical approximation to coupled system (46) were compared with each other as shown in Figure 3. Moreover, the RK4 numerical approximations to the coupling system (46) and decoupling system (47) are introduced in Figure 4. Over and above, the maximum residual distance error for the two methodologies was estimated using the same parameters of Figure 3 as follows

$$\left\{\begin{array}{c}{L}_x=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(46\right)}-\mathrm{KBM}\mathrm{approxi}.\right|=0.000230824,\\ L_y=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(46\right)}-\mathrm{KBM}\mathrm{approxi}.\right|=0.0004614,\\ L_x=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(46\right)}-{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(47\right)}\right|=0.0000764837,\\ L_y=\underset{0\le t\le 40}{max}\left|{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(46\right)}-{\left.\mathrm{RK}4\right|}_{\mathrm{For}\mathrm{system}\left(47\right)}\right|=0.000301353.\end{array}\right.$$

(50)

It is evident that both analytical and numerical solutions exhibit a high degree of accuracy and stability for a long time.

4. Conclusions

The coupled forced damped Duffing oscillators (FDDOs) were analyzed analytically and numerically using highly accurate approximate methods. In our investigation, two different models for the coupled FDDOs were analyzed in detail. For the analytical approximation to the proposed models, the Krylov–Bogoliubov–Mitropolsky (KBM) method was used for solving the coupled FDDOs. However, to apply the mentioned method, the coupled FDDOs were reduced to a decoupled system of two individual standard FDDOs. After that, the proposed method was applied to find some approximations for both unforced and forced damped Duffing oscillators in terms of trigonometric functions. Some numerical examples using the KBM approximations were analyzed and discussed. Moreover, the KBM analytical approximations to the decoupled system and the RK4 numerical approximations to the coupled and decoupled systems were compared with each other to verify the high accuracy of all obtained approximations. Furthermore, the maximum residual error for all obtained approximations was estimated on the whole time domain. It was found that both analytical and numerical approximations were compatible with each other, which confirmed the high accuracy of the obtained approximations.