Analytical Solutions of Some Strong Nonlinear Oscillators

Abstract

Oscillators are omnipresent; most of them are inherently nonlinear. Though a nonlinear equation mostly does not yield an exact analytic solution for itself, plethora of elementary yet practical techniques exist for extracting important information about the solution of equation. The purpose of this chapter is to introduce some new techniques for the readers which are carefully illustrated using mainly the examples of Duffing’s oscillator. Using the exact analytical solution to cubic Duffing and cubic-quinbic Duffing oscillators, we describe the way other conservative and some non conservative damped nonlinear oscillators may be studied using analytical techniques described here. We do not make use of perturbation techniques. However, some comparison with such methods are performed. We consider oscillators having the form x¨+fx=0 as well as x¨+2εẋ+fx=Ft, where x=xt and f=fx and Ft are continuous functions. In the present chapter, sometimes we will use f−x=−fx and take the approximation fx≈∑j=1Npjxj, where j=1,3,5,⋯N only odd integer values and x∈−AA. Moreover, we will take the approximation fx≈∑j=0Npjxj, where j=1,2,3,⋯N, and x∈−AA. Arbitrary initial conditions are considered. The main idea is to approximate the function f=fx by means of some suitable cubic or quintic polynomial. The analytical solutions are expressed in terms of the Jacobian and Weierstrass elliptic functions. Applications to plasma physics, electronic circuits, soliton theory, and engineering are provided.

1. Introduction

Both the ordinary and partial differential equations have an important role in explaining many phenomena that occur in nature or in medical engineering, biotechnology, economic, ocean, plasma physics, etc. [1, 2]. Duffing equation is considered one of the most important differential equations due to its ability for demonstrating the scenario and mechanism of various nonlinear phenomena that occur in nonlinear dynamic systems [3, 4, 5, 6, 7, 8, 9, 10, 11]. It is one of the most common models for analyzing and modeling many nonlinear phenomena in various fields of science such as the mechanical engineering [12], electrical engineering [13], plasma physics [14, 15], etc. Mathematically, the Duffing oscillator is a second-order ordinary differential equation with a nonlinear restoring force of odd power

$$\left\{\begin{array}{c}\ddot{x}+f\left(x\right)=0,\\ f\left(x\right)={\sum }_{i=1}^{\infty }{K}_i{x}^{2i-1},\end{array}\right.$$E1

where $f\left(-x\right)=-f\left(x\right)$ is a continuous function on some interval $\left[-A,A\right]$ with $f\left(0\right)=0,$$K_i$ is a physical coefficient related to the physical problem under study, and $i=\mathrm{1,2,3},\cdots \infty$. It is clear from Eq. (1) that there is no any friction/dissipation (this force arises either as a result of taking viscosity into account or the collisions between the oscillator and any other particle, etc.), and this only occurs in standardized systems such as superfluid (fluid with zero viscosity which it flows without losing any part from its kinetic energy sometimes like Bose–Einstein condensation) or the systems isolated from all the external force that resist the motion of the oscillator. The undamped Duffing equation [9] is considered one of the effective and good models for explaining many nonlinear phenomena that are created and propagated in optical fiber, Ocean, water tank, the laboratory and space collisionless and warm plasma (we will demonstrate this point below). As well known in fluid mechanics and in the fluid theory of plasma physics; the basic fluid equations of any plasma model can be reduced to a diverse series of evolution equations that can describe all phenomena that create and propagate in these physical models. For example, we can mention some of the most famous evolution equations that have been used to explain several phenomena in plasma physics and other fields of sciences; the family of one dimensional ($1-D$) korteweg–de Vries equation (KdV) and it is higher-orders, including the KdV, KdV-Burgers (KdVB), modified KdV (mKdV), mKdV-Burgers (mKdVB), Gardner equation or called Extended KdV (EKdV), EKdV-Burgers (EKdVB), KdV-type equation with higher-order nonlinearity. All the above mentioned equations are partial differential equations and by using an appropriate transformation, we can convert them into ordinary differential equations of the second orders. If the frictional force is neglected, some of these equations can be converted into the undamped Duffing equation with $f\left(x\right)\approx P\left(x\right)=K_{1}x+K_2{x}^3$ like the mKdV equation, the KdV equation can be transformed to the undamped Helmholtz equation with $f\left(x\right)\approx P\left(x\right)=K_{1}x+K_2{x}^2$ [16], the Gardner equation can be converted into the undamped H-D equation for $f\left(x\right)\approx P\left(x\right)=K_{1}x+K_2{x}^2+K_3{x}^3$ [17, 18], and so on the other mentioned equations.

However, these undamped models (without friction/dissipation) do not exist much in reality except under harsh conditions. In order to describe and simulate the natural phenomena that arise in many realistic physical models and dynamic systems, the friction/dissipation forces must be taken into account, as is the case in many plasma models and electronic systems. Accordingly, the following damped (non-conservative) Duffing equation will be devoted for this purpose

$$\ddot{x}+2\epsilon \dot{x}+f\left(x\right)=0.$$E2

If the frictional force does not neglect, so that all PDEs that have “Burgers $\equiv {\partial }_x^2\left(\cdot \right)$” term like KdVB-, mKdVB-, EKdVB-, KPB-, mKPB-, EKPB-, ZKB-, mZKB, EZKB-Eq. [1, 2], etc. can be transformed to damped Duffing equation ($\ddot{x}+2\epsilon \dot{x}+K_{1}x+K_2{x}^3=0$), damped Helmholtz equation ($\ddot{x}+2\epsilon \dot{x}+K_{1}x+K_2{x}^2=0$), and damped Duffing-Helmholtz equation ($\ddot{x}+2\epsilon \dot{x}+K_{1}x+K_2{x}^2+K_2{x}^3=0$). Eq. (2) without [19] and with [7, 20, 21] including damping term ($2\epsilon \dot{x}$) for $f\left(x\right)\approx P\left(x\right)=K_{1}x+K_2{x}^3$ has been investigated and solved analytically and numerically by many authors using different approaches in order to understand its physical characters [22, 23, 24, 25, 26, 27, 28].

Many authors investigated the (un)damped Duffing equation, (un)damped Helmholtz Eq. [16, 29, 30, 31], and undamped H-D equation. On the contrary, there is a few numbers of published papers about damped Duffing-Helmholtz equation [32, 33]. For example, Zúñiga [32] derived a semi-analytical solution to the damped Duffing-Helmholtz equation in the form of Jacobian elliptic functions, but he putted some restrictions on the coefficient of the linear term, and then obtained a solution that gives good results compared to numerical solutions. Also, it is noticed that Zúñiga solution [32] is very sensitive to the initial conditions. Gusso and Pimentel [33] obtained obtain improved approximate analytical solution to the forced and damped Duffing-Helmholtz in the form of a truncated Fourier series utilizing the harmonic balance method.

In this chapter, we display some novel semi-analytical (approximate analytical) solutions to the strong higher-order nonlinear damped oscillators of the following initial value problem (i.v.p)

$$\left\{\begin{array}{c}\ddot{x}+2\epsilon \dot{x}+\mathit{px}+{\mathit{qx}}^3+{\mathit{rx}}^5=\mathcal{F}\left(t\right),\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0,\end{array}\right.$$E3

and its family ($\epsilon =0$ or $r=0$ or $\epsilon =r=0$).

Our new semi-analytical solution to Eq. (3) is derived in terms of Weierstrass and Jacobian elliptic functions. Also, we will solve Eq. (3) numerically using Runge–Kutta 4th (RK4) and make a comparison between both the semi-analytical and numerical solutions. Moreover, as some realistic physical application to the problem (3) and its family will be investigated.

2. Duffing equation

Let us consider the standard (undamping) Duffing equation in the absence both friction ($2\epsilon \dot{x}$) and excitation ($\mathcal{F}\left(t\right)$) forces [34, 35]

$$\ddot{x}+\mathit{px}+{\mathit{qx}}^3=0,x=x\left(t\right),$$E4

which is subjected to the following initial conditions

$$x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0.$$E5

The general solution of Eq. (4) maybe written in terms of any of the twelve Jacobian elliptic functions.

For example, let us assume

$$x\left(t\right)=c_1\mathrm{cn}\left(\sqrt{\omega }t+c_2,m\right).$$E6

By inserting solution (6) in Eq. (4), we get

$$\ddot{x}+\mathit{px}+{\mathit{qx}}^3=\left(c_1^{3}q-2c_1\mathit{m\omega }\right){\mathrm{cn}}^3+\left(2c_1\mathit{m\omega }+c_{1}p-c_1\omega \right)\mathrm{cn},$$E7

where cn $=$ cn $\left(\sqrt{\omega }t+c_2,m\right).$

Equating to zero the coefficients of cn ${}^j$ gives an algebraic system whose solution gives

$$\omega =\sqrt{p+{\mathit{qc}}_1^2}\text{and}m=\frac{{\mathit{qc}}_1^2}{2\left(p+{\mathit{qc}}_1^2\right)}.$$E8

Thus, the general solution of Eq. (4) reads

$$x\left(t\right)=\mathrm{cn}\left(\sqrt{p+{\mathit{qc}}_1^2}t+c_2,\frac{{\mathit{qc}}_1^2}{2\left(p+{\mathit{qc}}_1^2\right)}\right).$$E9

The values of the constants $c_1$ and $c_2$ could be determined from the initial conditions given in Eq. (5).

Definition 1. The number $\mathrm{\Delta }={\left(p+{\mathit{qx}}_0^2\right)}^2+2q{\dot{x}}_0^2$ is called the discriminant of the i.v.p. (4) -(5). Below three cases will be discussed depending on the sign of the discriminant $\mathrm{\Delta }$.

2.1 First case: $\mathrm{\Delta }>0$

For $\mathrm{\Delta }>0$, the solution of the i.v.p. (4)–(5) is given by

$$x\left(t\right)=\sqrt{\frac{\sqrt{\mathrm{\Delta }}-p}{q}}\mathrm{cn}\left(\sqrt[4]{\mathrm{\Delta }}t-sign\left({\dot{x}}_0\right){\mathrm{cn}}^{-1}\left(\sqrt{\frac{q}{\sqrt{\mathrm{\Delta }}-p}}{x}_0,\frac{1}{2}-\frac{p}{2\sqrt{\mathrm{\Delta }}}\right),\frac{1}{2}-\frac{p}{2\sqrt{\mathrm{\Delta }}}\right).$$E10

Making use of the additional formula

$$\mathrm{cn}\left(x+y,m\right)=\frac{\mathrm{cn}\left(x,m\right)\mathrm{cn}\left(y,m\right)+\mathrm{sn}\left(x,m\right)\mathrm{dn}\left(x,m\right)\mathrm{sn}\left(y,m\right)\mathrm{dn}\left(y,m\right)}{1-m\mathrm{sn}\left(x,m\right)\mathrm{sn}\left(y,m\right)},$$E11

the solution (10) could be expressed as

$$x\left(t\right)=\frac{x_0\mathrm{cn}\left(\left.\sqrt{\omega }t\right|m\right)+\frac{{\dot{x}}_0}{\sqrt{\omega }}\mathrm{dn}\left(\left.\sqrt{\omega }t\right|m\right)\mathrm{sn}\left(\left.\sqrt{\omega }t\right|m\right)}{1+\frac{p+{\mathit{qx}}_0^2-\omega }{2\sqrt{\mathrm{\Delta }}}\mathrm{sn}{\left(\left.\sqrt[4]{\mathrm{\Delta }}t\right|m\right)}^2},$$E12

where

$$m=\frac{1}{2}\left(1+\frac{p}{\sqrt{\mathrm{\Delta }}}\right)\text{and}\omega =\sqrt[4]{\mathrm{\Delta }}.$$E13

Solution (12) is a periodic solution with period

$$T=4\left|\frac{K\left(m\right)}{\sqrt{\omega }}\right|.$$E14

Example 1.

Let us consider the i.v.p.

$$\left\{\begin{array}{c}{x}^{″}\left(t\right)+x\left(t\right)+x^3\left(t\right)=0,\\ x\left(0\right)=1\&x^{\prime }\left(0\right)=-1.\end{array}\right.$$E15

Using formula (10), the exact solution of the i.v.p. (15) reads

$$x\left(t\right)=-\sqrt{\sqrt{6}-1}\mathrm{cn}\left(\sqrt[4]{6}t+{\mathrm{cn}}^{-1}\left(-\frac{1}{\sqrt{\sqrt{6}-1}},\frac{1}{2}-\frac{1}{2\sqrt{6}}\right),\frac{1}{12}\left(6-\sqrt{6}\right)\right).$$E16

According to the relation (12)–(13), the exact solution of the i.v.p. (15) is also written as

$$x\left(t\right)=\frac{2\sqrt[4]{6}\mathrm{dn}\left(\sqrt[4]{6}t|\frac{1}{12}\left(6-\sqrt{6}\right)\right)\mathrm{sn}\left(\sqrt[4]{6}t|\frac{1}{12}\left(6-\sqrt{6}\right)\right)-2\sqrt{6}\mathrm{cn}\left(\sqrt[4]{6}t|\frac{1}{12}\left(6-\sqrt{6}\right)\right)}{\left(\sqrt{6}-2\right)\mathrm{sn}{\left(\sqrt[4]{6}t|\frac{1}{12}\left(6-\sqrt{6}\right)\right)}^2-2\sqrt{6}},$$E17

and its periodicity is given by

$$T=\frac{4K\left(\frac{1}{12}\left(6-\sqrt{6}\right)\right)}{\sqrt[4]{6}}\approx 3.27458.$$

In Figure 1, the comparison between the exact analytical solution (17) and the approximate numerical RK4 solution is presented. Full compatibility between the two analytical and numerical solutions is observed.

Figure 1.

A comparison between the analytical solution (17) and the approximate numerical RK4.

2.2 Second case: $\mathrm{\Delta }<0$

For $\mathrm{\Delta }<0,$ in this case $q<0$ and then, $\delta =\frac{p^2-\mathrm{\Delta }}{-q}>0$, $\delta \stackrel{\mathrm{def}}{=}\left(2p+{\mathit{qx}}_0^2\right)x_0^2+2{\dot{x}}_0^2>0$. Let us introduce the solution in the following form

$$x\left(t\right)=A-\frac{2A}{1+y\left(t\right)},$$E18

where $y=y\left(t\right)$ is a solution of some Duffing equation

$$y^{″}\left(t\right)+\mathit{my}\left(t\right)+{\mathit{ny}}^3\left(t\right),$$E19

with initial conditions

$$\begin{array}{l}y\left(0\right)=y_0=\frac{2A{\dot{x}}_0}{{\left(A-x_0\right)}^2}\\ y^{\prime }\left(0\right)={\dot{y}}_0=\frac{A+x_0}{A-x_0}.\end{array}$$E20

Inserting ansatz (18) into Eq. (4) and taking the below relation into account

$$y^{\prime }{\left(t\right)}^2={\dot{y}}_0+{\mathit{my}}_0^2+\frac{n}{2}{y}_0^4-{\mathit{my}}^2\left(t\right)-\frac{n}{2}{y}^4\left(t\right),$$E21

we get

$$\begin{array}{l}-A\left(A^{2}q+4{\mathit{my}}_0^2+2{\mathit{ny}}_0^4+p+4{\dot{y}}_0^2\right)\\ +A\left(3A^{2}q-2m-p\right)y\left(t\right)\\ -A\left(3A^{2}q-2m-p\right)y{\left(t\right)}^2\\ +A\left(A^{2}q-2n+p\right)y{\left(t\right)}^3=0.\end{array}$$E22

Equating the coefficients of $y^j\left(t\right)$ to zero, gives an algebraic system. A solution to this system gives

$$\begin{array}{l}m=\frac{1}{2}\left(-p+3A^{2}q\right),n=\frac{1}{2}\left(p+A^{2}q\right),\\ A=\sqrt[4]{\frac{\left(2p+{\mathit{qx}}_0^2\right)x_0^2+2{\dot{x}}_0^2}{-q}}=\sqrt[4]{\frac{\delta }{-q}}.\end{array}$$

Note that the i.v.p. (19)–(20) has a positive discriminant and it is given by

$${\left(m+{\mathit{ny}}_0^2\right)}^2+2n{\dot{y}}_0^2=\frac{\delta {\left(A-x_0\right)}^4\left(2A^4{\dot{x}}_0^2+\delta {\left(A^2+x_0^2\right)}^2\right)}{4A^8{x}_0^2}.$$

Then the problem reduces to the first case. Accordingly, the solution of the i.v.p. (4)–(5) maybe written in the form,

$$x\left(t\right)=A-\frac{2A}{1+B\frac{b_0\mathrm{cn}\left(\left.\sqrt{\omega }t\right|m\right)+b_1\mathrm{sn}\left(\left.\sqrt{\omega }t\right|m\right)\mathrm{dn}\left(\left.\sqrt{\omega }t\right|m\right)}{1+b_2{\mathrm{sn}}^2\left(\left.\sqrt{\omega }t\right|m\right)}},$$E23

where

$$\begin{array}{l}m=\frac{B^2\left(A^{2}q+p\right)}{2A^2\left(B^2+3\right)q+2\left(B^2-1\right)p},\\ \omega =\frac{1}{2}\left(A^2\left(B^2+3\right)q+\left(B^2-1\right)p\right),\\ b_0=\frac{A+{\mathrm{x}}_0}{\mathit{AB}-B{\mathrm{x}}_0},b_1=\frac{2A{\dot{\mathrm{x}}}_0}{B\sqrt{\omega }{\left(A-{\mathrm{x}}_0\right)}^2},\\ b_2=-\frac{2A{\mathrm{x}}_0\left({\mathrm{x}}_0-A\right)\left(p+q{{\mathrm{x}}_0}^2\right)+\omega \left(A+{\mathrm{x}}_0\right){\left(A-{\mathrm{x}}_0\right)}^2+4A{{\dot{\mathrm{x}}}_0}^2}{2\omega {\left(A-{\mathrm{x}}_0\right)}^2\left(A+{\mathrm{x}}_0\right)},\\ A=\sqrt[4]{\frac{\left(2p+{\mathit{qx}}_0^2\right)x_0^2+2{\dot{x}}_0^2}{-q}}=\sqrt[4]{\frac{\delta }{-q}},\\ B=\sqrt{\frac{A\left(2\sqrt{2}\sqrt{q\left(A^{2}q-p\right)}-3\mathit{Aq}\right)+p}{A^{2}q+p}}.\end{array}$$E24

The solution (23) is unbounded and its periodicity is given by

$$T=\left|\frac{4K\left(m\right)}{\sqrt{\omega }}\right|=\left|\frac{4K\left(1-m\right)}{m\sqrt{\omega }}\right|.$$E25

Example 2.

Let us assume the following i.v.p.

$$\left\{\begin{array}{c}{x}^{″}\left(t\right)+x\left(t\right)-x^3\left(t\right)=0,\\ x\left(0\right)=-1\&x^{\prime }\left(0\right)=-1.\end{array}\right.$$E26

The solution of the i.v.p. (26) according to the relation (23) reads

$$x\left(t\right)=1.31607-\frac{2.63215}{1+\frac{0.13647\mathrm{cn}\left(1.75396t|1.00353\right)-0.27976\mathrm{dn}\left(1.75396t|1.00353\right)\mathrm{sn}\left(1.75396t|1.00353\right)}{1.-1.00463\mathrm{sn}{\left(1.75396t|1.00353\right)}^2}},$$E27

and the periodicity of this solution is given by

$$T=9.57783.$$E28

Solution (27) is displayed in Figure 2.

Figure 2.

The profile of solution (27) is plotted against $t$.

2.3 Third case: $\mathrm{\Delta }=0$

If the discriminant vanishes $\left(\mathrm{\Delta }=0\right)$, then $q<0$ and the only solution of problem (4) with

$$x^{\text{'}}{\left(0\right)}^2={\dot{x}}_0^2=\frac{{\left(p+{\mathit{qy}}_0^2\right)}^2}{-2q},$$E29

reads

$$x\left(t\right)=\sqrt{-\frac{p}{q}}tanh\left[\sqrt{\frac{p}{2}}t\pm {tanh}^{-1}\left(x_0\sqrt{-\frac{q}{p}}\right)\right].$$E30

which may be verified by direct computation.

Remark 1. The solution of the i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+\mathit{px}+{\mathit{qx}}^3=0,\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)=0,\end{array}\right.$$E31

is given by

$$x\left(t\right)=x_0\mathrm{cn}\left(\sqrt{p+{\mathit{qx}}_0^2},\frac{{\mathit{qx}}_0^2}{2\left(p+{\mathit{qx}}_0^2\right)}\right).$$E32

Remark 2. For $p+\sqrt{p^2+2q}$ x_0 ${}^2>0$, then the solution of the i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+\mathit{px}+{\mathit{qx}}^3=0,\\ x\left(0\right)=0\&x^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E33

is given by

$$x\left(t\right)=\frac{\sqrt{2}{\dot{\mathrm{x}}}_0}{\sqrt{\sqrt{p^2+2q{{\dot{\mathrm{x}}}_0}^2}+p}}\mathrm{sn}\left(\sqrt{\frac{p+\sqrt{p^2+2q{{\dot{\mathrm{x}}}_0}^2}}{2}}t,-\frac{p^2+q{{\dot{\mathrm{x}}}_0}^2-\sqrt{p^2+2q{{\dot{\mathrm{x}}}_0}^2}p}{q{{\dot{\mathrm{x}}}_0}^2}\right).$$E34

Remark 3. According to the following identity

$$\mathrm{cn}\left(\sqrt{\omega }t,m\right)=1-\frac{S_0}{1+S_1\wp \left(t;g_2,g_3\right)},$$E35

with

$$\begin{array}{l}{S}_0=\frac{6}{\left(4m+1\right)},S_1=\frac{12}{\left(4m+1\right)\omega },\\ g_2=\frac{1}{12}\left(16m^2-16m+1\right){\omega }^2,\\ g_3=\frac{1}{216}\left(2m-1\right)\left(32m^2-32m-1\right){\omega }^3,\end{array}$$

the solution of the i.v.p. (4)–(5) could be written in terms of the Weierstrass elliptic function $\wp \equiv \wp \left(t;g_2,g_3\right)$. More precisely, if $\mathrm{\Delta }>0$ then

$$x\left(t\right)=A-\frac{A\left(\frac{4p}{3A^{2}q+p}+2\right)}{1+\frac{12}{3A^{2}q+p}\wp \left(t+t_0;g_2,g_3\right)},$$E36

with

$$\begin{array}{l}{t}_0={\wp }^{-1}\left(\frac{3A^{3}q+3A^2{\mathit{qx}}_0+5\mathit{Ap}+{\mathit{px}}_0}{12\left(A-x_0\right)};g_2,g_3\right),\\ g_2=\frac{1}{12}\left(-3A^4{q}^2-6A^2\mathit{pq}+p^2\right),\\ g_3=\frac{p}{216}\left(9A^4{q}^2+18A^2\mathit{pq}+p^2\right),\end{array}$$E37

and

$$A=\sqrt{\frac{-p\pm \sqrt{{\left(p+{\mathit{qx}}_0^2\right)}^2+2q{\dot{x}}_0^2}}{q}}=\pm \sqrt{\frac{-p\pm \sqrt{\mathrm{\Delta }}}{q}}.$$E38

The solution (36) is periodic with period

$$T=2{\int }_{\rho }^{+\infty }\frac{\mathit{dx}}{\sqrt{4x^3-g_{2}x-g_3}},$$E39

where $\rho$ is the greatest real root of the cubic $4x^3-g_{2}x-g_3=0$.

Remark 4. An approximate analytic solution of the i.v.p. (31) is given by

$$x\left(t\right)=\frac{x_0\sqrt{1+\lambda }cos\left(\mathit{wt}\right)}{\sqrt{1+\lambda {cos}^2\left(\mathit{wt}\right)}},$$E40

where

$$w=\frac{1}{2}\sqrt{\frac{5\left(p+{\mathit{qx}}_0^2\right){\lambda }^2+\left(12p+11{\mathit{qx}}_0^2\right)\lambda +8p+6{\mathit{qx}}_0^2}{3\lambda +2}}$$E41

and $\lambda$ is a root of the cubic

$$25\left(p+{\mathit{qx}}_0^2\right){\lambda }^3+\left(58p+59{\mathit{qx}}_0^2\right){\lambda }^2+2\left(16p+21{\mathit{qx}}_0^2\right)\lambda +8{\mathit{qx}}_0^2=0.$$E42

Example 3.

Let us consider the i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+x+10x^3=0,\\ x\left(0\right)=4,x^{\prime }\left(0\right)=0.\end{array}\right.$$E43

The approximate solution in trigonometric form is given by

$$x_{\mathrm{app}}\left(t\right)=\frac{3.34603cos\left(10.7542t\right)}{\sqrt{1-0.300255{cos}^2\left(10.7542t\right)}}.$$E44

The exact solution reads

$$x\left(t\right)=4\mathrm{cn}\left(\sqrt{161}t|\frac{80}{161}\right),$$E45

with period

$$T=\frac{4K\left(\frac{80}{161}\right)}{\sqrt{161}}.$$

The error on the interval $0\le t\le T$ equals $0.025.$

The comparison between the approximate analytic solution (44) and the exact analytic solution (45) is illustrated in Figure 3.

Figure 3.

A comparison between the approximate solution (44) and exact solution (45).

Remark 5. An approximate analytical solution of the i.v.p. (33) is given by

$$x_{\mathrm{app}}\left(t\right)=\frac{{\dot{x}}_{0}sin\left(\sqrt{\omega }t\right)}{\sqrt{\omega }\sqrt{1+\lambda {sin}^2\left(\sqrt{\omega }t\right)}},$$E46

where

$$\omega =-\frac{\sqrt{{\lambda }^2\left(64p^2-160q{\dot{x}}_0^2\right)+25p^2{\lambda }^4+80p^2{\lambda }^3-128q{\dot{x}}_0^2\lambda }+\mathit{p\lambda }\left(5\lambda +8\right)}{16\lambda }$$E47

and $\lambda$ is a solution of the quintic

$$\begin{array}{l}125p^2{\lambda }^5+10\left(79p^2+125q{\dot{x}}_0^2\right){\lambda }^4+40\left(43p^2+85q{\dot{x}}_0^2\right){\lambda }^3+\\ 8\left(196p^2+389q{\dot{x}}_0^2\right){\lambda }^2+64\left(8p^2+17q{\dot{x}}_0^2\right)\lambda +128q{\dot{x}}_0^2=0\end{array}$$E48

Example 4.

The approximate trigonometric solution of

$$\left\{\begin{array}{c}\ddot{x}+3x+5x^3=0,\\ x\left(0\right)=0,x^{\prime }\left(0\right)=1.\end{array}\right.$$E49

reads

$$x_{\mathrm{app}}\left(t\right)=\frac{0.499502sin\left(2.00199t\right)}{\sqrt{1-0.0817025{sin}^2\left(2.00199t\right)}}.$$E50

The exact solution is

$$x\left(t\right)=\sqrt{\frac{2}{3+\sqrt{19}}}\mathrm{sn}\left(\sqrt{\frac{1}{2}\left(3+\sqrt{19}\right)}t,\frac{1}{5}\left(-14+3\sqrt{19}\right)\right),$$E51

with period

$$T=4\sqrt{\frac{1}{5}\left(\sqrt{19}-3\right)}K\left(\frac{1}{5}\left(-14+3\sqrt{19}\right)\right)=3.1383.$$E52

The error on the interval $0\le t\le T$ equals $0.00018291.$

Figure 4 demonstrates the comparison between the approximate analytic solution (50) and the exact analytic solution.

Figure 4.

A comparison between the approximate solution (50) and exact solution (51).

3. An analytical solution of the undamped Duffing-Helmholtz Equation

The undamped Duffing-Helmholtz equation reads

$$\left\{\begin{array}{c}\ddot{x}+\mathit{px}+{\mathit{qx}}^2+{\mathit{rx}}^3=0,\\ x\left(0\right)=x_0\text{and}{x}^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E53

We will give a solution to the i.v.p. (53) in terms of Weierstrass elliptic functions. For solving this problem the following ansatz is considered

$$x\left(t\right)=A+\frac{B}{1+C\wp \left(t+t_0;g_2,g_3\right)},$$E54

where $\mathit{BC}\ne 0.$

Substituting the ansatz (54) into the ordinary differential equation (ode) $\mathbb{R}\equiv \ddot{x}+\mathit{px}+{\mathit{qx}}^2+{\mathit{rx}}^3=0$, gives

$$\frac{1}{2\left(1+C{\wp }^3\right)}\sum _{j=0}^3{K}_j{\wp }^j=0,$$E55

with

$$\begin{array}{l}{K}_3=2C^2\left(A^3\mathit{Cr}+A^2\mathit{Cq}+\mathit{ACp}+2B\right),\\ K_2=2C\left(3A^3\mathit{Cr}+3A^2\mathit{BCr}+3A^2\mathit{Cq}+2\mathit{\text{ABCq}}+3\mathit{ACp}+\mathit{BCp}-6B\right),\\ K_1=C\left(6A^{3}r+12A^2\mathit{Br}+6A^{2}q+6{\mathit{AB}}^{2}r+8\mathit{ABq}+6\mathit{Ap}+2B^{2}q-3{\mathit{BCg}}_2+4\mathit{Bp}\right),\\ K_0=A^{3}r+6A^2\mathit{Br}+2A^{2}q+6{\mathit{AB}}^{2}r+4\mathit{ABq}+2\mathit{Ap}+2B^{3}r+2B^{2}q-4{\mathit{BC}}^2{g}_3+{\mathit{BCg}}_2+2\mathit{Bp}.\end{array}$$

Equating the coefficients $K_j$ to zero will give us an algebraic system. Solving this system, we finally get

$$\begin{array}{l}B=-\frac{6A\left(A^{2}r+\mathit{Aq}+p\right)}{3A^{2}r+2\mathit{Aq}+p},C=\frac{12}{3A^{2}r+2\mathit{Aq}+p},\\ g_2=-\frac{1}{12}\left(3r^2{A}^4+4{\mathit{qrA}}^3+6{\mathit{prA}}^2-p^2\right),\\ g_3=\frac{1}{216}\left[\left(9{\mathit{pr}}^2-3q^{2}r\right)A^4+\left(12\mathit{pqr}-4q^3\right)A^3+\left(18p^{2}r-6{\mathit{pq}}^2\right)A^2+p^3\right],\end{array}$$E56

The values of $t_0$ and $A$ could be determined from the initial conditions $x\left(0\right)=x_0$ and $x^{\prime }\left(0\right)={\dot{x}}_0$ and

$$\ddot{x}\left(0\right)+\mathit{px}\left(0\right)+{\mathit{qx}}^2\left(0\right)+{\mathit{rx}}^3\left(0\right)=0.$$E57

We have

$$t_0=\pm {\wp }^{-1}\left(\frac{x_0-A-B}{C\left(A-x_0\right)};g_2,g_3\right).$$E58

The number $A$ is a solution to the quartic

$$3{\mathit{rA}}^4+4{\mathit{qA}}^2+6\mathit{pA}-\left(3{\mathit{rx}}_0^4+4{\mathit{qx}}_0^3+6{\mathit{px}}_0^2+6{\dot{x}}_0^2\right)=0.$$E59

Example 5.

The solution of the i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+x+2x^2+3x^3=0,\\ x\left(0\right)=1\text{and}{x}^{\prime }\left(0\right)=1,\end{array}\right.$$E60

according to the relation (54) is given by

$$x\left(t\right)=1.07627-\frac{2.72078}{1+0.762858\wp \left(t-0.148317;-\mathrm{7.16667,0.675926}\right)}.$$E61

In Figure 5, the comparison with the approximate analytic solution (61) and the approximate numerical solution using RK4 is investigated.

Figure 5.

A comparison between solution (61) and the approximate numerical solution using RK4.

The periodicity of solution (61) is given by

$$T=2{\int }_{0.0938538}^{\infty }\frac{1}{\sqrt{4x^3+7.16667x-0.675926}}\mathit{dx}=3.12129.$$

4. The solution of the forced undamped Duffing-Helmholtz equation

Suppose that the physical system to be studied is under the influence of some constant external/excitation force, so the standard Duffing-Helmholtz equation can be reformulated to the following constant forced Duffing-Helmholtz i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+\mathit{px}+{\mathit{qx}}^2+{\mathit{rx}}^3=F,\\ x\left(0\right)=x_0\text{and}{x}^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E62

For solving the i.v.p. (62), the following assumption is introduced

$$x\left(t\right)=y\left(t\right)+\zeta ,$$E63

where $\zeta$ is a solution to the cubic algebraic equation

$$r{\zeta }^3+q{\zeta }^2+\mathit{p\zeta }-F=0.$$E64

Substituting Eq. (63) into the i.v.p. (62), we have

$$y^{\text{'}\text{'}}\left(t\right)+\left(p+2\mathit{q\zeta }+3r{\zeta }^2\right)y\left(t\right)+\left(q+3\mathit{r\zeta }\right)y{\left(t\right)}^2+\mathit{ry}{\left(t\right)}^3=0.$$E65

Note that the constant forced Duffing-Helmholtz Eq. (62) has been reduced to the standard Duffing-Helmholtz Eq. (65) with the following new initial conditions

$$y\left(0\right)=x_0-\zeta \&y^{\prime }\left(0\right)={\dot{x}}_0.$$E66

Example 6.

Suppose that we have the following i.v.p. and we want to solve it

$$\left\{\begin{array}{c}\ddot{v}+2v-12v^2+v^3=4,\\ v\left(0\right)=1\&v^{\prime }\left(0\right)=1.\end{array}\right.$$E67

It is clear that the i.v.p. (67) is a constant forced Duffing-Helmholtz equation. The solution of this problem is given by

$$v\left(t\right)=15.8046-\frac{15.7714}{1+0.0322539\wp \left(0.761045-t;-\mathrm{9.41667,47.287}\right)}.$$E68

The comparison between the solution (68) and the RK4 solution is introduced in Figure 6.

Figure 6.

A comparison between the solution (68) and the RK4 solution.

The periodicity of solution (68) is given by

$$T=2{\int }_{1.93657}^{\infty }\frac{1}{\sqrt{4x^3+9.41667x-47.287}}\mathit{dx}=1.68202.$$

5. An approximate analytic solution of the forced damped Duffing-Helmholtz equation

Let us define the following i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+2\epsilon \dot{x}+\mathit{px}+{\mathit{qx}}^2+{\mathit{rx}}^3=F,\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E69

Suppose that

$$\underset{x\to +\infty }{lim}x\left(t\right)=d,\epsilon >0,$$E70

then the first equation in system (69) can be written as

$$\mathit{pd}+{\mathit{qd}}^2+{\mathit{rd}}^3=F.$$E71

For solving the i.v.p. (69), the following ansatz is assumed

$$x\left(t\right)=exp\left(-\mathit{\rho t}\right)y\left(f\left(t\right)\right),$$E72

with

$$f\left(t\right)=\frac{1-exp\left(-2\left(\epsilon -\rho \right)t\right)}{2\left(\epsilon -\rho \right)},$$E73

where the function $y\equiv y\left(t\right)$ represents the exact solution to the following i.v.p.

$$\left\{\begin{array}{c}{y}^{″}\left(t\right)+\left(3d^{2}r+2\mathit{dq}-2\mathit{\epsilon \rho }+p+{\rho }^2\right)y\left(t\right)+\left(3\mathit{dr}+q\right)y{\left(t\right)}^2+\mathit{ry}{\left(t\right)}^3=0,\\ y\left(0\right)=x_0-d\&y^{\prime }\left(0\right)={\dot{x}}_0+\rho \left(x_0-d\right).\end{array}\right.$$E74

Let us define the following residual

$$R\left(t\right)\equiv \ddot{x}\left(t\right)+2\epsilon \dot{x}\left(t\right)+\mathit{px}\left(t\right)+{\mathit{qx}}^2\left(t\right)+{\mathit{rx}}^3\left(t\right)-F,$$E75

and by applying the condition $R^{\prime }\left(0\right)=0$, we obtain

$$\begin{array}{l}4{\rho }^3-12{\mathit{\epsilon \rho }}^2+\left(3d^{2}r+3\mathit{dq}+3{\mathit{drx}}_0+8{\epsilon }^2+4p+5{\mathit{qx}}_0+6{\mathit{rx}}_0^2\right)\rho \\ -4\epsilon \left(d^{2}r+\mathit{dq}+{\mathit{drx}}_0+p+{\mathit{qx}}_0+{\mathit{rx}}_0^2\right)=0.\end{array}$$E76

By solving this equation we can get the value of $\rho$.

Example 7.

Let

$$\left\{\begin{array}{c}\ddot{x}+0.02\dot{x}+5x+2x^2+x^3=1/2,\\ x\left(0\right)=0.1\&x^{\prime }\left(0\right)=0.1.\end{array}\right.$$E77

The approximate analytic solution of the i.v.p. (77) reads

$$\begin{array}{cc}{x}_{\mathrm{app}}\left(t\right)& =0.0961263\\ & +e^{-0.0099959t}\left(\begin{array}{c}0.0429135\\ -\frac{0.252808}{1+2.13749\wp \left(0.631364-121944.\left(1-e^{-8.2\times {10}^{-6}t}\right);\mathrm{2.43588,0.737005}\right)}\end{array}\right).\end{array}$$E78

The distance error as compared to the RK4 numerical solution is given by

$$\underset{0\le t\le 5}{max}\left|x_{\mathrm{app}}\left(t\right)-x_{\mathrm{RK}4}\left(t\right)\right|=0.000944148.$$E79

Also, the comparison between solution (78) and RK4 solution is presented in Figure 7.

Figure 7.

A comparison between solution (78) and RK4 solution.

Remark 5. For the damped and constant forced Helmholtz equation

$$\left\{\begin{array}{c}\ddot{x}+2\epsilon \dot{x}+\mathit{px}+{\mathit{qx}}^2=F,\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E80

The value of $d$ can be determined from: $\mathit{pd}+{\mathit{qd}}^2=F.$ However, if this equation has no real solutions we can choose $d=0$.

Remark 6. Letting $q=0$, we obtain the damped and constant forced Duffing equation

$$\left\{\begin{array}{c}\ddot{x}+2\epsilon \dot{x}+\mathit{px}+{\mathit{rx}}^3=F,\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E81

In this case, the number $d$ must be a root to the cubic $\mathit{pd}+{\mathit{rd}}^3=F.$

6. Approximate analytic solution of the damped and trigonometric forced Duffing-Helmholtz equation

Let us define the following new i.v.p.

$$\left\{\begin{array}{c}\ddot{x}+2\epsilon \dot{x}+\mathit{px}+{\mathit{qx}}^2+{\mathit{rx}}^3=Fcos\left(\mathit{\omega t}\right),\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E82

We suppose that $q^2-4\mathit{pr}<0,$ and the following residual is defined

$$R\left(t\right)\equiv \ddot{x}\left(t\right)+2\epsilon \dot{x}\left(t\right)+\mathit{px}\left(t\right)+{\mathit{qx}}^2\left(t\right)+{\mathit{rx}}^3\left(t\right)-Fcos\left(\mathit{\omega t}\right).$$E83

Let us define the solution of i.v.p. (82) as follows

$$x\left(t\right)=exp\left(-\mathit{\rho t}\right)y\left(t\right)+c_{1}cos\left(\mathit{\omega t}\right)+c_{2}sin\left(\mathit{\omega t}\right),$$E84

where

$$\begin{array}{l}9F^2{r}^2{c}_1^3+96{\epsilon }^2\mathit{Fr}{\omega }^2{c}_1^2+\\ 4\left(64{\epsilon }^4{\omega }^4+16{\epsilon }^2{\omega }^6+3F^2\mathit{pr}-3F^{2}r{\omega }^2+16{\epsilon }^2{p}^2{\omega }^2-32{\epsilon }^{2}p{\omega }^4\right)c_1\\ -4F\left(-16{\epsilon }^2{\omega }^4+3F^{2}r+16{\epsilon }^{2}p{\omega }^2\right)=0.\end{array}$$E85

$$\begin{array}{l}-6144{\epsilon }^3\mathit{Fr}{\omega }^3+432F^2{r}^3{c}_2^3+2304\epsilon {\mathit{Fr}}^2\omega \left(p-{\omega }^2\right)c_2^2\\ +3072{\epsilon }^{2}r{\omega }^2\left(4{\epsilon }^2{\omega }^2+p^2-2p{\omega }^2+{\omega }^4\right)c_2=0.\end{array}$$E86

The function $y\equiv y\left(t\right)$ is a solution to the i.v.p.

$$\left\{\begin{array}{c}{y}^{″}\left(t\right)+2\epsilon y^{\prime }\left(t\right)+\tilde{p}y\left(t\right)+\mathit{qy}{\left(t\right)}^2+\mathit{ry}{\left(t\right)}^3=0,\\ y\left(0\right)=\left(x_0-c_1\right)\&y^{\prime }\left(0\right)=\left({\dot{x}}_0-\omega c_2\right).\end{array}\right.$$E87

where $\tilde{p}=\frac{1}{2}\left(2p+3{\mathit{rc}}_1^2+3{\mathit{rc}}_2^2-4\mathit{\epsilon \rho }+2{\rho }^2\right)$.

The value of $\rho$ can be determined from the following equation

$$\begin{array}{c}4{\rho }^3-12{\mathit{\epsilon \rho }}^2+\left(4p+8{\epsilon }^2-5{\mathit{qc}}_1+12{\mathit{rc}}_1^2+6{\mathit{rc}}_2^2+5{\mathit{qx}}_0-12{\mathit{rc}}_1{x}_0+6{\mathit{rx}}_0^2\right)\rho \\ -2\epsilon \left(2p-2{\mathit{qc}}_1+5{\mathit{rc}}_1^2+3{\mathit{rc}}_2^2+2{\mathit{qx}}_0-4{\mathit{rc}}_1{x}_0+2{\mathit{rx}}_0^2\right)=0.\end{array}$$E88

Example 8.

Let

$$\left\{\begin{array}{c}\ddot{x}+0.2\dot{x}+13x+x^2+x^3=0.25cos\left(0.5t\right),\\ x\left(0\right)=0\&x^{\prime }\left(0\right)=-0.2.\end{array}\right.$$E89

The approximate analytic solution of the i.v.p. (89) is given by

$$\begin{array}{ll}{x}_{\mathrm{app}}\left(t\right)& =e^{-0.100036t}\left(\begin{array}{c}0.0587764-\\ \frac{0.350893}{1+0.914739\wp \left(0.529144-13766.9\left(1-e^{0.0000726379t}\right);\mathrm{14.0404,10.1967}\right)}\end{array}\right)\\ & +0.000153771sin\left(0.5t\right)+0.0196062cos\left(0.5t\right)\end{array}$$E90

The distance error according to the RK4 numerical solution is calculated as

$$\underset{0\le t\le 60}{max}\left|x_{\mathrm{app}}\left(t\right)-x_{\mathrm{RK}4}\left(t\right)\right|=0.000671928.$$E91

Moreover, solution (90) is compared with RK4 solution as shown in Figure 8.

Figure 8.

A comparison between solution (90) and RK4 solution.

7. An analytic solution of cubic-quintic Duffing equation

Let us consider the following ordinary differential equation [36]

$$\ddot{x}+\mathit{\alpha x}+\beta x^3+\gamma x^5=0,x=x\left(t\right),$$E92

which is subjected to the following initial conditions

$$x\left(0\right)=x_0\text{and}{x}^{\prime }\left(0\right)={\dot{x}}_0.$$E93

Theorem 1.

a. Suppose that $x_0\ne 0$, then the solution of the i.v.p. (92)–(93) is given by

$$x\left(t\right)=x_0\frac{\sqrt{1+\lambda }v\left(t\right)}{\sqrt{1+\lambda v^2\left(t\right)}},$$E94

where the function $v\equiv v\left(t\right)$ is the solution to the following Duffing equation

$$\left\{\begin{array}{l}\ddot{v}+\mathit{pv}+{\mathit{qv}}^3=0,\\ v\left(0\right)≔v_0=1,\\ v^{\prime }\left(0\right)≔{\dot{v}}_0=\frac{{\dot{x}}_0}{x_0}\left(1+\lambda \right).\end{array}\right.$$E95

The values of the coefficients $p$ and $q$ are given by

$$p=-\frac{{\lambda }^2\left(4\alpha +3\beta x_0^2+2\gamma x_0^4\right)+\lambda \left(3\beta x_0^2+4\gamma x_0^4\right)+2\gamma x_0^4}{2{\lambda }^2},$$E96

$$q=-\frac{2\left({\lambda }^2\left(\alpha +\beta x_0^2+\gamma x_0^4\right)+\lambda \left(\beta x_0^2+2\gamma x_0^4\right)+\gamma x_0^4\right)}{\lambda },$$E97

and the value of the quantity $\lambda$ is a solution of the cubic

$$6{\dot{x}}_0^2{\lambda }^3-\left(6\alpha x_0^2+6\beta x_0^4+6\gamma x_0^6\right){\lambda }^2-\left(3\beta x_0^4+6\gamma x_0^6\right)\lambda -2\gamma x_0^6=0.$$E98

The solution to the the i.v.p. (95) is obtained from the formulas in the first section.

b. Suppose that $x_0=0$, in this case, the solution of the i.v.p. (92)–(93) is given by

$$x\left(t\right)=\frac{\sqrt{1+\lambda }v\left(t\right)}{\sqrt{1+\lambda v^2\left(t\right)}},$$E99

where the function $v\equiv v\left(t\right)$ is the solution of the following Duffing equation

$$\left\{\begin{array}{l}\ddot{v}+\mathit{pv}+{\mathit{qv}}^3=0,\\ v\left(0\right)≔v_0=0,\\ v^{\prime }\left(0\right)≔{\dot{v}}_0=\frac{{\dot{x}}_0}{\sqrt{1+\lambda }}.\end{array}\right.$$E100

The values of the coefficients $p$ and $q$ are expressed as

$$p=-\frac{2\gamma +\left(3\beta +4\gamma \right)\lambda +\left(4\alpha +3\beta +2\gamma \right){\lambda }^2}{2{\lambda }^2},$$E101

$$q=-\frac{2\left(\gamma +\left(\beta +2\gamma \right)\lambda +\left(\alpha +\beta +\gamma \right){\lambda }^2\right)}{\lambda },$$E102

and the value of $\lambda$ is a solution of the cubic

$$\left(6\alpha +3\beta +2\gamma -6{\dot{x}}_0^2\right){\lambda }^3+\left(6\alpha +6\beta +6\gamma \right){\lambda }^2+\left(3\beta +6\gamma \right)\lambda +2\gamma =0.$$E103

Note that the solution of the i.v.p. (100) could be obtained from the formulas in the first section.

• Proof: case (a)

Inserting ansatz (94) into Eq. (92) taking the following equation into consideration

$${\dot{v}}^2={\dot{v}}_0^2+{\mathit{pv}}_0^2+\frac{q}{2}{v}_0^4-{\mathit{pv}}^2-\frac{q}{2}{v}^4,$$E104

and using Eq. (100), we have

$$\sum _{j=1}^5{H}_{j}v{\left(t\right)}^j=0,$$E105

with

$$\begin{array}{l}{H}_1=\frac{\left(-6{\mathit{px}}_0^2\lambda -3{\mathit{qx}}_0^2\lambda +2x_0^2\alpha -6{\dot{x}}_0^2{\lambda }^3-12{\dot{x}}_0^2{\lambda }^2-6{\dot{x}}_0^2\lambda -2x_0^2\right)}{2x_0^2},\\ H_3=-\left(-3\mathit{p\lambda }+q-x_0^2\mathit{\beta \lambda }-x_0^2\beta -2\mathit{\alpha \lambda }+\lambda \right),\\ H_5=\frac{1}{2}\left(\mathit{q\lambda }+2x_0^2{\mathit{\beta \lambda }}^2+2x_0^2\mathit{\beta \lambda }+2x_0^4{\mathit{\gamma \lambda }}^2+4x_0^4\mathit{\gamma \lambda }+2x_0^4\gamma +2{\mathit{\alpha \lambda }}^2\right),\end{array}$$

where $j=\mathrm{1,3,5}$.

Equating the coefficients $H_j$ to zero gives an algebraic system: $H_1=0$, $H_3=0$, and $H_5=0$. Solving $H_1=0$ and $H_3=0$ will give the values of $p$ and $q$ that are given in Eqs. (101)–(102). Finally, by inserting the values of $p$ and $q$ into $H_1=0$, we obtain the cubic Eq. (103). Likewise, the case (b) can be proved.

8. Damped Cubic-Quintic Oscillator

Let us define the following i.v.p.

$$\left\{\begin{array}{c}\ddot{x}++2\epsilon \dot{x}+\mathit{px}+{\mathit{qx}}^3+{\mathit{rx}}^5=0,\\ x\left(0\right)=x_0\&x^{\prime }\left(0\right)={\dot{x}}_0.\end{array}\right.$$E106

We seek approximate analytic solution in the ansatz form

$$x\left(t\right)=exp\left(-\mathit{\rho t}\right)y\left(f\left(t\right)\right),$$E107

with

$$f\left(t\right)=\frac{1-exp\left(-2\left(\epsilon -\rho \right)t\right)}{2\left(\epsilon -\rho \right)},$$E108

where $y\equiv y\left(t\right)$ is the exact solution to the i.v.p.

$$\left\{\begin{array}{c}\ddot{y}+\left(p-\mathit{\epsilon \rho }+{\rho }^2\right)y+{\mathit{qy}}^3+{\mathit{ry}}^5=0,\\ y\left(0\right)=x_0\&y^{\prime }\left(0\right)={\dot{x}}_0+x_0\rho .\end{array}\right.$$E109

Define the residual

$$R\left(t\right)=\ddot{x}+2\epsilon \dot{x}+\mathit{px}+{\mathit{qx}}^3+{\mathit{rx}}^5,$$E110

then, the condition $R^{\prime }\left(0\right)=0$ gives

$$2{\rho }^3-6{\mathit{\epsilon \rho }}^2+\left(2p+3{\mathit{qx}}_0^2+4{\mathit{rx}}_0+4{\epsilon }^2\right)\rho -2\mathit{p\epsilon }-2{\mathit{qx}}_0^2\epsilon -2{\mathit{rx}}_0^4\epsilon =0.$$E111

Some real roots of Eq. (111) give the value of $\rho$. For $x_0=0$, the default value of $\rho$ could be chosen as $\rho =\left(2/3\right)\epsilon$.

Example 9.

Let

$$\left\{\begin{array}{c}\ddot{x}+0.05\dot{x}+96.6289x-3.5x^3-0.8x^5=0,\\ x\left(0\right)=1\&x^{\prime }\left(0\right)=0.\end{array}\right.$$E112

The approximate analytical solution of the i.v.p. is given by

$$\begin{array}{l}{x}_{\mathrm{app}}\left(t\right)\\ =\frac{1.00575e^{-0.0257101t}\left(\mathrm{cn}\left(\left.f\left(t\right)\right|m\right)-0.00261773\mathrm{dn}\left(\left.f\left(t\right)\right|m\right)\mathrm{sn}\left(\left.f\left(t\right)\right|m\right)\right)}{\left(1+5.14\times {10}^{-9}\mathrm{sn}{\left(\left.f\left(t\right)\right|m\right)}^2\right)\sqrt{\begin{array}{c}1+\frac{0.0107805{\left(\mathrm{cn}\left(\left.f\left(t\right)\right|m\right)-0.00261773\mathrm{dn}\left(\left.f\left(t\right)\right|m\right)\mathrm{sn}\left(\left.f\left(t\right)\right|m\right)\right)}^2}{{\left(1+5.14\times {10}^{-9}\mathrm{sn}{\left(\left.f\left(t\right)\right|m\right)}^2\right)}^2}+\\ \frac{0.000750033{\left(\mathrm{cn}\left(\left.f\left(t\right)\right|m\right)-0.00261773\mathrm{dn}\left(\left.f\left(t\right)\right|m\right)\mathrm{sn}\left(\left.f\left(t\right)\right|m\right)\right)}^4}{{\left(1+5.14\times {10}^{-9}\mathrm{sn}{\left(\left.f\left(t\right)\right|m\right)}^2\right)}^4}\end{array}}},\end{array}$$E113

where

$$f\left(t\right)=6807.39-6807.39e^{0.00142017t}\&m=-0.000750607.$$E114

The distance error as compared to the RK4 numerical solution reads

$$\underset{0\le t\le 5}{max}\left|x_{\mathrm{app}}\left(t\right)-x_{\mathrm{RK}4}\left(t\right)\right|=0.0561216.$$E115

In Figure 9, we make a comparison between our solution, RK4 solution, and Zuñiga solution given in Ref. [17]. It is clear that the accuracy of our solution is better than the solution of Zuñiga [17].

Figure 9.

A comparison between RK4 solution and (a) Solution (113) and (b) Zuñiga solution [17].

9. Realistic physical applications

The above solutions could be applied to various fields of physics and engineering such as they could be used for describing the behavior of oscillations in RLC electronic circuits, plasma physics etc. In the below section, the above solution will be devoted for studying oscillations in various plasma models.

9.1 Nonlinear oscillations in RLC series circuits with external source

In the RLC series circuits consisting of a linear resistor with resistance $R$ in Ohm unit, a linear inductor with inductance $L$ in Henry unit, and nonlinear capacitor with capacitance $C$ in Farady unit as well as external applied voltage $E$ in voltage unit, the Kirchhoff’s voltage law (KVL) could be written as

$$L{\partial }_t{i}^{\text{'}}\left(t\right)+i\left(t\right)R+\mathit{sq}+{\mathit{aq}}^2=E,$$E116

where the relation between the current the charge is given by $i={\partial }_{t}q\equiv \dot{q}$, $i^{\prime }\equiv {\partial }_{t}i$, the coefficients $\left(a,s\right)$ are related to the nonlinear capacitor, and $E$ represents the voltage of the battery which is constant. By reorganizing Eq. (116), the following constant forced and damped Helmholtz equation could be obtained as

$$\ddot{q}+2\gamma \dot{q}+\mathit{\alpha q}+\beta q^2=F,$$E117

with $\gamma =R/\left(2L\right)$$,\alpha =1/\left(\mathit{LC}\right)$, $\beta =1/\left({\mathit{Cq}}_{0}L\right)$, and $F=E/L$ where $q_0=q\left(t=0\right)$ is the initial charge value at $t=0$, $\ddot{q}\equiv {\partial }_t^{2}q$, and $\dot{q}\equiv {\partial }_{t}q$.

The solution of Eq. (117) can be devoted for interpreting and analyzing the oscillations that can generated in the RLC circuit.

9.2 Duffing-Helmholtz equation for modeling the oscillations in a plasma

For studying the plasma oscillations using fluid theory, the basic equations of plasma particles using the reductive perturbation method (RPM) will be reduced to some evolution equations such as KdV equation and its family [37, 38, 39, 40, 41]. Let us consider a collisionaless and unmagnetized electronegative complex plasma, consisting of inertialess cold positive and negative ion species, inertia non-Maxwellian electrons in addition to stationary negative dust impurities [42]. Thus, the quasi-neutrality condition reads: $n_2^{\left(0\right)}+n_e^{\left(0\right)}=n_1^{\left(0\right)}$ where $n_{s,e}^{\left(0\right)}$ donates the unperturbed number density of the plasma particles (here, the index “$s$” $=$ “1” and “2” point out the positive ion and negative ion, and “$e$” refers to the electron, respectively). It is assumed that the plasma oscillations take place only in $x-$directional which means that the fluid equations of the plasma particles become perturbed only in $x-$directional. If the effect of the ionic kinematic viscosities ${\eta }_s$ for both positive $\left({\eta }_1\right)$ and negative $\left({\eta }_2\right)$ ions are included in the present investigation, as a source of damping/dissipation, in this case we will get a new evolution equation governs the dynamics of damping pulses. The dynamics of plasma oscillations are governed by the following fluid equations: ${\partial }_x{n}_s+{\partial }_t\left(n_s{u}_s\right)=0,$${\partial }_t{u}_s+u_s{\partial }_x{u}_s+\left(\delta /Q_s\right){\partial }_x\varphi -{\eta }_s{\partial }_x^2{u}_s=0,$ and ${\partial }_x^2\varphi -n_d^{\left(0\right)}-n_2-n_e+n_1=0$, where $n_e=\mu \left(1-\mathit{\beta \varphi }+{\mathit{\beta \varphi }}^2\right)exp\varphi$. Here, $n_s$ donates the normalized number density of positive and negative ions, and $u_s$ represents the normalized fluid velocity of positive and negative ions, and $\varphi$ is the normalized electrostatic wave potential. The mass ratio is defined as: $Q_s=m_1/m_s$ (note that $Q_1=m_1/m_1=1$ and $Q_2\equiv Q=m_1/m_2$), where $m_s$ is the ionic mass, $\delta =1\left(-1\right)$ for positive (negative) ion, and $\beta$ illustrates nonthermality parameter. The quasi-neutrality condition in the normalized form reads: $\mu =1-\alpha ,$ where $\alpha =n_2^{\left(0\right)}/n_1^{\left(0\right)}$ gives the the negative ion concentration and $\mu =n_e^{\left(0\right)}/n_1^{\left(0\right)}$ is the electron concentration.

Now, the RPM is introduced to reduce the fluid plasma equations to the evolution equation. According to the RPM, the independent quantities $\left(x,t\right)$ are stretched as: $\xi =\epsilon \left(x-V_{\mathit{ph}}t\right)$, $\tau ={\epsilon }^{3}t$, and ${\eta }_s={\mathit{\epsilon \eta }}_s^{\left(0\right)}$ where $V_{\mathit{ph}}$ is the wave phase velocity of the ion-acoustic waves and $\epsilon$ is a real and small parameter ($0<\epsilon <<1$). The dependent perturbed quantities $\mathrm{\Pi }\left(x,t\right)\equiv {\left(n_1,n_2,u_1,u_2,\varphi \right)}^T$ are expanded as: $\mathrm{\Pi }={\mathrm{\Pi }}^{\left(0\right)}+{\sum }_{j=1}^{\infty }{\epsilon }^j\mathrm{\Pi }{\left(\xi ,\tau \right)}^{\left(j\right)},$ where ${\mathrm{\Pi }}^{\left(0\right)}={\left(1,\alpha \mathrm{,0,0,0}\right)}^T$ and $T$ represents the transpose of the matrix. Inserting both the stretching and expansions of the independent and dependent quantities into the basic fluid equations and after boring but straightforward calculations, the Gardner-Burgers/EKdVB equation is obtained

$${\partial }_{\tau }\phi +\left(P_1\phi +P_2{\phi }^2\right){\partial }_{\xi }\phi +P_3{\partial }_{\xi }^3\phi +P_4{\partial }_{\xi }^2\phi =0,$$E118

with the coefficients of the quadratic nonlinear, cubic nonlinear, dispersion, and dissipation terms $P_1,$$P_2$, $P_3$, and $P_4$, respectively,

$$\begin{array}{l}{P}_1=\frac{3}{2}{P}_3\left[\frac{1}{V_{\mathit{ph}}^4}-\frac{3\alpha }{V_{\mathit{ph}}^4{Q}^2}-\frac{2h_2}{3}\right],P_2=\frac{3}{4}{P}_3\left[\frac{5}{V_{\mathit{ph}}^6}-\frac{5\alpha }{V_{\mathit{ph}}^6{Q}^3}-2h_3\right],\\ P_3=\frac{V_{\mathit{ph}}^{3}Q}{2\left(Q+\alpha \right)},P_4=-P_3\left[\frac{{\eta }_1^{\left(0\right)}}{V_{\mathit{ph}}^4}+\frac{{\mathit{\alpha \eta }}_2^{\left(0\right)}}{{\mathit{QV}}_{\mathit{ph}}^3}\right]\&V_{\mathit{ph}}=\sqrt{\frac{\left(Q+\alpha \right)}{{\mathit{Qh}}_1}},\end{array}$$

where $\phi \equiv {\varphi }^{\left(1\right)}$, $h_1=\mu \left(1-\beta \right)$, $h_2=\mu /2$, and $h_3=\mu \left(1+3\beta \right)/6.$

It is shown that the coefficients $P_1,$$P_2$, $P_3$, and $P_4$, are functions in the physical plasma parameters namely, negative ion concentration $\alpha$, the mass ratio $Q$, and the electron nonthermal parameter $\beta$. It is known that at the critical plasma compositions say ${\beta }_c$ or ${\alpha }_c$ (critical value of negative ion concentration), the coefficient $P_1$ vanishes and in this case Eq. (118) will be reduced to the following mKdVB equation which is used to describe the damped wave dynamics at critical plasma compositions

$${\partial }_{\tau }\phi +P_2{\phi }^2{\partial }_{\xi }\phi +P_3{\partial }_{\xi }^3\phi +P_4{\partial }_{\xi }^2\phi =0,$$E119

To convert EKdVB Eq. (118) to the damped H-D Eq. (4), the traveling wave transformation $\phi \left(\xi ,\tau \right)\to \phi \left(X\right)$ with $X=\left(\xi +\mathit{\lambda \tau }\right)$ should be inserted into Eq. (118) and integrate once over $\eta$, and by applying the boundary conditions: $\left(\phi ,{\phi }^{\prime },{\phi }^{″}\right)\to 0$ as $\mid X\mid \to \infty$, the constant forced damped following constant forced damped Duffing-Helmholtz equation is obtained

$${\phi }^{″}+2\epsilon {\phi }^{\prime }+\mathit{p\phi }+q{\phi }^2+r{\phi }^3+D=0,$$E120

where $\lambda$ represents the reference frame speed, ${\phi }^{\text{'}}$ and ${\phi }^{\text{'}\text{'}}$ denote the first and second ordinary derivative of regarding $X$, $\epsilon =P_4/\left(2P_3\right)$, $p=\lambda /P_3$, $q=P_1/\left(2P_3\right)$, $r=P_2/\left(3P_3\right)$, and $D=C/P_3$.

Note that the coefficient $q$ may be positive or negative according to the values of plasma parameters and for studying oscillations using (120), solution (72) can be devoted for this purpose. In the absence of the ionic kinematic viscosity ($P_4=0$ or $\epsilon =0$), then Eq. (120) reduces to the constant forced undamping Duffing-Helmholtz equation and in this case the solution (63) can be applied for investigating the undamped oscillations in the present plasma model. Also, for $q=0$, Eq. (120) reduces to the constant forced damped Helmholtz equation. Moreover, the constant forced damped Duffing equation can be obtained for $p=0$.

10. Conclusion

The analytical and semi-analytical solutions for nonlinear oscillator integrable and non-integrable equations have been investiagted. First, the standard integrable Duffing equation has been analyized and its solutions have been obtained depending on the sign of its discriminant $\mathrm{\Delta }$. Accordingly, three cases $\left(\mathrm{\Delta }>0,\mathrm{\Delta }<0,\text{and}\mathrm{\Delta }=0\right)$ have been discussed in details and the solutions of each case has be obtained. Second, the analytical and semi-analytical solutions of the integrable Duffing-Helmholtz equation and its non-integrable family including the damped Duffing-Helmholtz equation, forced undamped Duffing-Helmholtz equation, forced damped Duffing-Helmholtz equation, and the damped and trigonometric forced Duffing-Helmholtz equation have been obtained and discussed in details. Third, the solutions to the intgrable cubic-quintic Duffing equation and the non-intgrable damped cubic-quintic Duffing equation have been investigated. Moreover, some realistic applications reaslted to the RLC circuits and physics of plasmas have been introduced and discussed depending on the solutions of the mentioned evolution equations.