# An Elementary Solution to a Duffing Equation

Authors: Alvaro H. Salas
Journal: The Scientific World Journal
 (2022)

---

## Abstract

In this work, we study the Duffing equation. Analytical solution for undamped and unforced case is provided for any given arbitrary initial conditions. An approximate analytical solution is given for the damped or trigonometrically forced Duffing equation for arbitrary initial conditions. The analytical solutions are expressed in terms of elementary trigonometric functions as well as in terms of the Jacobian elliptic functions. Examples are added to illustrate the obtained results. We also introduce new functions for approximating the Jacobian and Weierstrass elliptic functions in terms of the trigonometric functions sine and cosine. Results are high accurate.


---

## Body

## 1. Introduction

Many physical phenomena are modeled by nonlinear systems of ordinary differential equations. The Duffing equation is an externally forced and damped oscillator equation that exhibits a range of interesting dynamic behavior in its solutions. The Duffing oscillator is an important model of nonlinear and chaotic dynamics. It was introduced by Germanic engineer Duffing in 1918 [1]. The Duffing oscillator is described by the differential equation:(1)x¨+rx˙+ω02x+βx3=Fcosωt.It differs from the classical forced and damped harmonic oscillator only by the nonlinear termβx3, which changes the dynamics of the system drastically. Motivated by potential applications in physics, engineering, biology, and communication theory, the damped Duffing equation(2)x¨+rx˙+ω02x+βx3=0,is considered. Equation (2) is a ubiquitous model arising in many branches of physics and engineering, such as the study of oscillations of a rigid pendulum undergoing with moderately large amplitude motion [2, 3], vibrations of a buckled beam, and so on [3–5].It has provided a useful paradigm for studying nonlinear oscillations and chaotic dynamical systems, dating back to the development of approximate analytical methods based on perturbative ideas [2], and continuing with the advent of fast numerical integration by the computer, to be used as an archetypal illustration of chaos [2, 5–7]. Various methods for studying the damped Duffing equation and the forced Duffing equation (1) in feedback control, strange attractor, stability, periodic solutions, and numerical simulations have been proposed, and a vast number of profound results have been established [2].The Duffing equation has been studied extensively in the literature. However, only few works are devoted to the study of its analytical solutions not using perturbation methods [8, 9]. Our aim is to avoid using such perturbation methods. This study is organized as follows. In the first section, we give exact analytical solution for the undamped and unforced Duffing equation for any given arbitrary initial conditions. In the second section, we provide formulas for obtaining a good approximate analytical solution using a new ansatz. The problems are solved for any arbitrary initial conditions. Finally, in the last section, we give approximate analytical solution to (1) and we compare it with Runge–Kutta numerical solution. Other useful methods are the homotopy perturbation method (HPM) [10–17], the Lindstedt–Poincaré method, and the Krylov–Bogoliubov–Mitropolsky method. The importance of numerical solution of differential equations in different fields of science and engineering is given in [18, 19].


## 2. Undamped and Unforced Duffing Equation

This is the equation:(3)x¨+px+qx3=0,x=xt,and given the initial conditions,(4)x0=x0andx′0=x˙0.The general solution to equation (3) may be written in terms of any of the twelve Jacobian elliptic functions [20]. Let, for example,(5)xt=c1cnωt+c2,m.Then,(6)x¨+px+qx3=c13q−2c1mωcn3+2c1mω+c1p−c1ωcn,wherecn=cnωt+c2,m.Equating to zero, the coefficients of cnj to zero gives an algebraic system whose solution is(7)ω=p+qc12andm=qc122p+qc12.Thus, the general solution to the Duffing equation is(8)xt=cnp+qc12t+c2,qc122p+qc12.The values for the constantsc1 and c2 are determined from the initial conditions.Definition 1.
The numberΔ=p+qx022+2qx˙02 is called the discriminant for the Duffing equations (3) and (4).
We will distinguish three cases depending on the sign of the discriminant [20].


### 2.1. First Case:Δ>0

The solution to the i.v.p. (3) and (4) is given by(9)xt=Δ−pqcnΔt4−signx˙0cn−1qΔ−px0,12−p2Δ,12−p2Δ.Making use of the addition formula,(10)cnx+y,m=cnx,mcny,m+snx,mdnx,msny,mdny,m1−msnx,msny,m.The solution (9) may be expressed as(11)xt=x0cnΔt4|m+x˙0/Δ4snΔt4|mdnΔt4|m1+p+qx02/2Δ−1/2snΔt4|m2,where(12)m=121−pΔ.Solution (11) is a periodic solution with period(13)T=4K1/21−p/ΔΔ4.Example 1.
Let us consider the i.v.p.(14)x″t+xt+x3t=0,x0=1&x′0=−1.
Using formula (9), the exact solution to (14) is given by(15)xt=6−1cn6t4−cn−11−1+6,12−126|12−126.
According to the relations (11) and (12), the exact solution to the i.v.p. (14) may also be written as(16)xt=cn6t4|1/21−1/6+6dn46t4|1/21−1/6sn6t4|1/21−1/61+2−6/26sn6t4|1/21−1/62.
The period is given by(17)T=28434K121−16≈4.37417.
In Figure1, the comparison between the exact analytical solution (??) and the approximate numerical RK4 solution is presented. Full compatibility between the two analytical and numerical solutions is observed.Figure 1
Comparison between the exact solution and the numerical solution for Example1.


### 2.2. Second Case:Δ<0

In this case,q<0. Define(18)δ=−2px02+qx04+2x˙02q.Observe that(19)δ=p2−Δq2>0.Let(20)xt=c−2c1+yt,where y=yt is a solution to Duffing equation(21)y′′t+ayt+by3twith initial conditions(22)y0=y0=2cx˙0c−x02andy′0=y˙0=c+x0c−x0.Inserting ansatz (21) into the ode x′′t+pxt+qx3t=0 and taking into account the relation,(23)y′t2=y˙0+ay02+b2y04−ay2t−b2y4t,we get(24)cyt−c−x0yt−x04c−x04yt+14x¨+px+qx3=8ac5−8ac4x0−8ac3x02+8ac2x03+4bc5+12bc4x0+12bc3x02+4bc2x03+c7q−5c6qx0+2c5p+11c5qx02−10c4px0−15c4qx03+20c3px02+15c3qx04+30c3x˙02−20c2px03−11c2qx05−22c2x0x˙02+10cpx04+5cqx06+10cx02x˙02−2px05−qx07−2x03x˙02+c−x08ac4−16ac3x0+8ac2x02+4bc4+8bc3x0+4bc2x02−3c6q+10c5qx0+2c4p−11c4qx02−12c3px0+24c2px02+11c2qx04+22c2x˙02−20cpx03−10cqx05−20cx0x˙02+6px04+3qx06+6x02x˙02yt+c−x024bc3+4bc2x0+3c5q−5c4qx0−2c3p−2c2px0+10cpx02+5cqx04+10cx˙02−6px03−3qx05−6x0x˙02y2t+c−x034bc2−qc4−2c2p+2px02+qx04+2x˙02y3t.Equating to zero, the coefficients ofyjt give an algebraic system. A solution to this system is(25)a=12−p−3δq,b=12p−δq,c=δ4.Observe that the Duffing equations (21) and (22) have a positive discriminant given by(26)a+by022+2by˙02=δc−x042c4x˙02+δc2+x0224c8x02.Then, the problem reduces to the first case.Example 2.
Let us assume the following i.v.p.:(27)x′′t+2xt−x3t=0,x0=−1andx′0=1.
The solution of i.v.p. (27) according to the relation (??) reads(28)xt=1.49535−2.99071+0.198509cn2.08627t|1.00018−0.230219dn2.08627t|1.00018sn2.08627t|1.00018/1−1.00072sn2.08627t|1.000182.
The period of solution (28) is given by(29)T=4K1/mwm=10.9034form=1.00018andw=2.08627.
Comparison between the exact solution and numerical solution is shown in Figure2.Figure 2
Comparison between the exact solution and the numerical solution for Example2.


### 2.3. Third Case:Δ=0 and p≠0

If the discriminant vanishesΔ=0, then q<0, and the only solution to problem (??) with(30)x′02=x˙02=p+qy022−2qreads(31)xt=−pqtanhp2t±tanh−1x0−qp,which may be verified by direct computation.


### 2.4. Fourth Case:Δ=0 and p=0

The solution is given by(32)xt=2x02+−2qx0t.Remark 1.
The solution to the i.v.p.(33)x¨+px+qx3=0,x0=x0andx′0=0,is(34)xt=x0cnp+qx02,qx022p+qx02.Remark 2.
Letp+p2+2qx˙02>0. Then, the solution to the i.v.p.(35)x¨+px+qx3=0,x0=0andx′0=x˙0,is(36)xt=2x˙0p2+2qx˙02+psnp+p2+2qx˙022t,−p2+qx˙02−p2+2qx˙02pqx˙02.Remark 3.
Using the identity(37)cnωt,m=1−64m+11+12/4m+1ω℘t;1/1216m2−16m+1ω2,1/2162m−132m2−32m−1ω3,the solution to the Duffing equations (3) and (4) may be written in terms of the Weierstrass elliptic function ℘. More precisely, if Δ>0, then(38)xt=A−A2+4p/3A2q+p1+12/3A2q+p℘t+t0;1/12−3q2A4−6pqA2+p2,1/216p9q2A4+18pqA2+p2,where(39)t0=℘−13A3q+3A2qx0+5Ap+px012A−x0;112−3A4q2−6A2pq+p2,1216p9A4q2+18A2pq+p2,(40)A=−p±p+qx022+2qx˙02q=±−p±Δq.
The solution (38) is periodic with period(41)T=2∫ρ+∞dx4x3−g2x−g3,where ρ is the greatest real root to the cubic 4x3−g2x−g3=0 and(42)g2=112−3q2A4−6pqA2+p2andg3=1216p9A4q2+18A2pq+p2.
On the other hand,(43)℘t;g2,g3=3g216m2−16m+11−cn23g2/16m2−16m+14t,m−g24m+1248m2−48m+3,where m is a root to the sextic(44)4096g23−27g32z6−12288g23−27g32z5+13056g23−27g32z4−5632g23−27g32z3+1259g23−1836g32z2+125g23+108g32z+g23−27g32=0.
Thus,(45)∫ρ+∞dx4x3−g2x−g3=216m2−16m+13g24Kmform<1.


## 2.1. First Case:Δ>0

The solution to the i.v.p. (3) and (4) is given by(9)xt=Δ−pqcnΔt4−signx˙0cn−1qΔ−px0,12−p2Δ,12−p2Δ.Making use of the addition formula,(10)cnx+y,m=cnx,mcny,m+snx,mdnx,msny,mdny,m1−msnx,msny,m.The solution (9) may be expressed as(11)xt=x0cnΔt4|m+x˙0/Δ4snΔt4|mdnΔt4|m1+p+qx02/2Δ−1/2snΔt4|m2,where(12)m=121−pΔ.Solution (11) is a periodic solution with period(13)T=4K1/21−p/ΔΔ4.Example 1.
Let us consider the i.v.p.(14)x″t+xt+x3t=0,x0=1&x′0=−1.
Using formula (9), the exact solution to (14) is given by(15)xt=6−1cn6t4−cn−11−1+6,12−126|12−126.
According to the relations (11) and (12), the exact solution to the i.v.p. (14) may also be written as(16)xt=cn6t4|1/21−1/6+6dn46t4|1/21−1/6sn6t4|1/21−1/61+2−6/26sn6t4|1/21−1/62.
The period is given by(17)T=28434K121−16≈4.37417.
In Figure1, the comparison between the exact analytical solution (??) and the approximate numerical RK4 solution is presented. Full compatibility between the two analytical and numerical solutions is observed.Figure 1
Comparison between the exact solution and the numerical solution for Example1.


## 2.2. Second Case:Δ<0

In this case,q<0. Define(18)δ=−2px02+qx04+2x˙02q.Observe that(19)δ=p2−Δq2>0.Let(20)xt=c−2c1+yt,where y=yt is a solution to Duffing equation(21)y′′t+ayt+by3twith initial conditions(22)y0=y0=2cx˙0c−x02andy′0=y˙0=c+x0c−x0.Inserting ansatz (21) into the ode x′′t+pxt+qx3t=0 and taking into account the relation,(23)y′t2=y˙0+ay02+b2y04−ay2t−b2y4t,we get(24)cyt−c−x0yt−x04c−x04yt+14x¨+px+qx3=8ac5−8ac4x0−8ac3x02+8ac2x03+4bc5+12bc4x0+12bc3x02+4bc2x03+c7q−5c6qx0+2c5p+11c5qx02−10c4px0−15c4qx03+20c3px02+15c3qx04+30c3x˙02−20c2px03−11c2qx05−22c2x0x˙02+10cpx04+5cqx06+10cx02x˙02−2px05−qx07−2x03x˙02+c−x08ac4−16ac3x0+8ac2x02+4bc4+8bc3x0+4bc2x02−3c6q+10c5qx0+2c4p−11c4qx02−12c3px0+24c2px02+11c2qx04+22c2x˙02−20cpx03−10cqx05−20cx0x˙02+6px04+3qx06+6x02x˙02yt+c−x024bc3+4bc2x0+3c5q−5c4qx0−2c3p−2c2px0+10cpx02+5cqx04+10cx˙02−6px03−3qx05−6x0x˙02y2t+c−x034bc2−qc4−2c2p+2px02+qx04+2x˙02y3t.Equating to zero, the coefficients ofyjt give an algebraic system. A solution to this system is(25)a=12−p−3δq,b=12p−δq,c=δ4.Observe that the Duffing equations (21) and (22) have a positive discriminant given by(26)a+by022+2by˙02=δc−x042c4x˙02+δc2+x0224c8x02.Then, the problem reduces to the first case.Example 2.
Let us assume the following i.v.p.:(27)x′′t+2xt−x3t=0,x0=−1andx′0=1.
The solution of i.v.p. (27) according to the relation (??) reads(28)xt=1.49535−2.99071+0.198509cn2.08627t|1.00018−0.230219dn2.08627t|1.00018sn2.08627t|1.00018/1−1.00072sn2.08627t|1.000182.
The period of solution (28) is given by(29)T=4K1/mwm=10.9034form=1.00018andw=2.08627.
Comparison between the exact solution and numerical solution is shown in Figure2.Figure 2
Comparison between the exact solution and the numerical solution for Example2.


## 2.3. Third Case:Δ=0 and p≠0

If the discriminant vanishesΔ=0, then q<0, and the only solution to problem (??) with(30)x′02=x˙02=p+qy022−2qreads(31)xt=−pqtanhp2t±tanh−1x0−qp,which may be verified by direct computation.


## 2.4. Fourth Case:Δ=0 and p=0

The solution is given by(32)xt=2x02+−2qx0t.Remark 1.
The solution to the i.v.p.(33)x¨+px+qx3=0,x0=x0andx′0=0,is(34)xt=x0cnp+qx02,qx022p+qx02.Remark 2.
Letp+p2+2qx˙02>0. Then, the solution to the i.v.p.(35)x¨+px+qx3=0,x0=0andx′0=x˙0,is(36)xt=2x˙0p2+2qx˙02+psnp+p2+2qx˙022t,−p2+qx˙02−p2+2qx˙02pqx˙02.Remark 3.
Using the identity(37)cnωt,m=1−64m+11+12/4m+1ω℘t;1/1216m2−16m+1ω2,1/2162m−132m2−32m−1ω3,the solution to the Duffing equations (3) and (4) may be written in terms of the Weierstrass elliptic function ℘. More precisely, if Δ>0, then(38)xt=A−A2+4p/3A2q+p1+12/3A2q+p℘t+t0;1/12−3q2A4−6pqA2+p2,1/216p9q2A4+18pqA2+p2,where(39)t0=℘−13A3q+3A2qx0+5Ap+px012A−x0;112−3A4q2−6A2pq+p2,1216p9A4q2+18A2pq+p2,(40)A=−p±p+qx022+2qx˙02q=±−p±Δq.
The solution (38) is periodic with period(41)T=2∫ρ+∞dx4x3−g2x−g3,where ρ is the greatest real root to the cubic 4x3−g2x−g3=0 and(42)g2=112−3q2A4−6pqA2+p2andg3=1216p9A4q2+18A2pq+p2.
On the other hand,(43)℘t;g2,g3=3g216m2−16m+11−cn23g2/16m2−16m+14t,m−g24m+1248m2−48m+3,where m is a root to the sextic(44)4096g23−27g32z6−12288g23−27g32z5+13056g23−27g32z4−5632g23−27g32z3+1259g23−1836g32z2+125g23+108g32z+g23−27g32=0.
Thus,(45)∫ρ+∞dx4x3−g2x−g3=216m2−16m+13g24Kmform<1.


## 3. Approximate Analytical Solution Using Elementary Functions

We define the generalized sine and cosine functions as follows:(46)sinmt=sin1+κt1+κcos21+κt,(47)cosmt=1+κcos1+κt1+κcos21+κt,(48)κ=114m−12+144−144m+m2.These functions are good approximations to the Jacobian elliptic functions sn and cn for−1≤m≤1/2. For example, let(49)T=2Km≈π1+κ.Then,(50)Max−T≤t≤T,−1≤m≤1/2snt,m−sinmt=0.00290582,(51)Max−T≤t≤T,−1≤m≤1/2cnt,m−cosmt=0.00541969.Table1 provides the errors for different values of m.Table 1
Errors for approximations (46) and (47).
mMax−T≤t≤Tsnt,m−sinmtMax−T≤t≤Tcnt,m−cosmt0000.050.00001861420.00003324840.10.00007756410.0001390510.150.0001821020.0003277220.20.0003384160.0006115340.250.0005538610.001005220.30.0008372560.001526650.350.00119930.002197750.40.001653120.003045750.450.002215060.004105020.50.002905820.005419690.550.003752090.007047580.60.004789160.0.00906640.660.006064990.01158380.70.007647180.01475460.750.009635460.01881290.80.01218650.02413640.850.0559620.03139490.90.02031370.04196350.950.02776810.0595356More accurate approximations are obtained by letting(52)sinmt=sin1+κt1+κcos21+κt+μcos41+κt,(53)cosmt=1+κ+μcos1+κt1+κcos21+κt+μcos41+κt,(54)sinmt=sint−κ−μ−1/μ−11−μcos2t−κ−μ−1/μ−1μcos4t−κ−μ−1/μ−1+κcos2t−κ−μ−1/μ−1+1.(55)κ=5184m2−14256m+15633μ4+−24192m2+92016m−89856μ3+42048m2−211392m+193536μ2+−32256m2+207360m−184320μ+9216m2−73728m+65536+−72m−9μ2+168m−960μ−96m−256483μ+10(56)μ=4m21061m4−7768m3+24128m2−32768m+1638410521m6−136752m5+786336m4−2345984m3+3792896m2−3145728m+1048576.For these new approximations, we will have(57)Max−T≤t≤T,−1≤m≤0.9snt,m−sinmt=0.0607706,(58)Max−T≤t≤T,−1≤m≤0.9cnt,m−cosmt=0.030797.Table2 provides the errors for different values of m.Table 2
Errors for approximations (52) and (53).
mMax−T≤t≤Tsnt,m−sinmtMax−T≤t≤Tcnt,m−cosmt0000.058.066e-96.548e-90.11.4551e-71.158e-70.158035e-76.51e-70.23.0125e-62.23e-60.258.45e-66.29e-60.30.000020.0000150.350.000040.0000310.40.0000890.0000610.450.000170.000110.50.000320.000210.550.000580.000370.60.00110.000650.650.00190.00120.70.00360.00210.750.00710.00390.80.01460.00790.850.03270.01730.90.08320.0430.950.2650.136From Tables1 and 2, it is seen that for the values 0.8<m<1, the approximations (46) and (47) are better than (53) and (54). Thus, we have(59)snt,m≈sinmt,cnt,m≈cosmtanddnt,m≈1−msn2t,m=dnmtfor−1≤m<1.We may write approximate elementary solution to Duffing equations (3) and (4) as follows:(60)xt=x0cosmΔt4+x˙0/Δ4sinmΔt4dnmΔt41+p+qx02/2Δ−1/2sinmΔt42,where(61)m=121−pΔ,p>0,Δ=p+qx022+2qx˙02>0.(62)xt=A−2A1+Bb0cosmωt+b1sinmωtdnmωt/1+b2sinmωt2,Δ<0.The values for the constants in (62) are the same as in (??).Remark 4.
In the case whenm>1, we use the approximations:(63)cnt,m≈dn1/mmt,(64)snt,m≈1msin1/mmt,(65)dnt,m≈cos1/mmt.
Then, formula (60) takes the following form when p<0 and m>1:(66)xt=x0dn1/mΔ4mt+x˙0/Δ41/msin1/mmΔt4cos1/mΔ4mt1+1/mp+qx02/2Δ−1/2sin1/m2Δ4mt,Δ>0.
Similar formula for (62) when Δ<0.
In the case when0.9<m≤1.1, we may use the following approximations:(67)cnt,m≈18m−1sinh2t−2ttanhtsecht+secht,sn≈14m−1tsech2t−m−5tanht,dn≈secht−14m−1tanhtsinht+tsecht.Example 3.
Let us return to Example1.(68)x′′t+xt+x3t=0,x0=1&x′0=−1.
The approximate elementary analytical solution is(69)xappt=3.4319cos1.43552t/14.−2.22208cos21.43552t−0.638943sin1.43552t1+4.14226sin21.43552t/2.22208cos21.43552t−14./1.−0.15872cos21.43552t1−0.0917517sin21.43552t/1.−0.15872cos21.43552t.
The exact period is given byT=4.37417. The approximate period is that of (59), and it is given by(70)Tapp=2π1.4355163606=4.37695.
This value differs from the exact value by 0.00278457. The error of the approximate solution comparted with exact solution is(71)max−T/2≤t≤T/2xappt−xexactt=0.001847.
Comparison between the exact solution and the approximate analytical solution is shown in Figure3.Figure 3
Comparison between the exact solution and the approximate analytical solution for Example3.Example 4.
Let us return to Example2. Let (Figure 4)(72)x′′t+2xt−x3t=0,x0=−1andx0=1.Figure 4
Comparison between the exact solution and the approximate analytical solution for Example4.Remark 5.
From (47)–(44) or (53)–(44), we obtain the following approximate expression for the Weierstrass elliptic function by means of the cosine function:(73)℘t;g2,g3=3g216m2−16m+11−cosm23g2/16m2−16m+14t−g24m+1248m2−48m+3,where m is a root to the sextic (44).Example 5.
Letg2=2 and g3=1. For this choice, m=0.0119056. We have(74)℘t;2,1≈2.718671−3.73048cos2.32484t/14.−0.0835474cos22.32484t−0.474691.
The period isT=2.70262, and the error on the interval −T/2≤t≤T/2 in the sup norm is E=8.52×10−7. Comparison between the approximate analytical solution and the numerical solution is shown in Figure 5. The reciprocals of the two functions are plotted on the interval −T/2≤t≤T/2.Figure 5
Comparison between the approximate analytical solution and the numerical solution for Example5.


## 4. Analytical Solution to a Generalized Duffing Equation

Let us consider the i.v.p. [21]:(75)u¨+αu+βu3=F,given that(76)u0=u0andu′0=u˙0.We will say that (75) is a constantly forced Duffing equation. When F=0, that becomes an undamped and unforced Duffing equation, and we already know how to solve it for arbitrary initial conditions. Let(77)xt=λ+μ1+℘ρt+t0;g2,g3,where λ,μ, w,g2, g3, and t0 are some constants to be determined. Plugging ansatz (77) into (75) gives(78)u¨+αu+βu3−F=121+℘3−2F−αλ−βλ3−2μρ2−23F−3αλ−αμ−3βλ3−3βλ2μ+6μρ2℘+−6F−3g2μρ2+6αλ+4αμ+6βλ3+12βλ2μ+6βλμ2℘2+−2F+g2μρ2−4g3μρ2+2αλ+2αμ+2βλ3+6βλ2μ+6βλμ2+2βμ3℘3,where ℘=℘ρt+t0;g2,g3. Equating to zero, the coefficients of ℘jj=0,1,2,3 in the right-hand side of (78) gives an algebraic system. A nontrivial solution to this system is(79)μ=6F−λα+βλ2α+3βλ2,ρ=12α3+βλ2,g2=12−144βλ−F+αλ+βλ3α+3βλ22,g3=827F2β+9αβλβλ3−4F+α3+18α2βλ2α+3βλ23.Now, to find the values oft0 and λ, we make use of the addition formula:(80)℘w+z;g2,g3=14℘′w;g2,g3−℘′z;g2,g3℘w;g2,g3−℘z;g2,g32−℘w;g2,g3−℘z;g2,g3.We then find that(81)t0=℘−1u0−λ−μλ−u0;g2,g3.The numberλ must be a solution to the quartic(82)4Fu0−βu04−2u02α−2u˙02−4Fλ+2αλ2+βλ4=0.Using (75), we also may obtain an approximate analytical solution in terms of the cosine function.Remark 6.
The i.v.p.(83)x¨+n+px+qx2+rx3=0,x=xt,given the initial conditions(84)x0=x0andx′0=x˙0,is a particular case of (75) and (76). Indeed, let xt=ut−q/3r. Then, problems (83) and (84) reduce to the problem(85)n+utp−q23r−pq3r+2q327r2+rut3+u′′t=0,(86)u0=x0+q3randu′0=x˙0.Example 6.
Let(87)u¨+ut+u3t=1∧u0=1∧u′0=1.
The exact solution is(88)uexactt=0.29637−2.457921+℘0.370192−0.709561t;−0.657491,1.34218,with period(89)Texact=4.24726758.
An approximate analytical solution is(90)uappt=1.29637−2.457920.474364+2.28411/1−3.64963cos0.771764−1.47927t/14.−0.68022cos20.771764−1.47927t,with period(91)Tapp=4.24748881.
The error is(92)max−T/2≤t≤T/2uexactt−uappt=0.000105178.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure6.Figure 6
Comparison between the approximate analytical solution and the numerical solution for Example6.


## 5. Damped and Unforced Duffing Equation

Let us consider the i.v.p.(93)u¨+2εu˙+αu+βu3=0,ε>0,given that(94)u0=u0andu′0=u˙0.We will suppose thatlimt⟶∞ut=0. Define the residual as(95)Rt=u′′t+2εu′t+put+qu3t.


### 5.1. First Case:u0≠0

Assume the ansatz(96)ut=exp−ρtu0cnft,mt+b1snft,mtdnft,mt1+b2sn2ft,mt.Then, from results in [22],(97)f′t2=α−2ερ+ρ2+βu02e−2tρandmt=βu02e−2tρ2f′t2,so that(98)ft=1ρβu02+μ−βu02e−2tρ+μ+μtanh−11+βu02e−2tρμ−tanh−11+βu02μ,(99)mt=1/21+μe2tρ/βu02,being(100)μ=α−2ερ+ρ2.The numbersb1, b2, and ρ are obtained from the following conditions:(101)u′0=u˙0,R0=0,R′0=0.Solving the two equations in (88) gives(102)b1=u0ρ+u˙0u02β+μ.b2=u03β3α−4ερ−3μ+u˙0u02β2ε−3ρ+2ε−ρα+ρρ−2ε+2u0ρε−ρα+ρρ−2ε2u0u02β+μ2.The numberρ is a root to the septic(103)−2εu02β+α2u02α+u04β+2u˙0u0ε+3u˙02+u02β+α2u02α2+5u04αβ+8u02αε2−4u˙0u0αε+6u˙02α+3u06β2+8u04βε2+2u˙0u03βε+9u˙02u02β+16u˙0u0ε3+24u˙02ε2z−28u02α2ε−4u˙0u0α2+15u04αβε−9u˙0u03αβ+4u02αε3−4u˙0u0αε2+18u˙02αε+7u06β2ε−3u˙0u05β2+4u04βε3+2u˙0u03βε2+21u˙02u02βε+8u˙0u0ε4+12u˙02ε3z2+8u02α2+17u04αβ+32u02αε2−40u˙0u0αε+12u˙02α+5u06β2+24u04βε2−38u˙0u03βε+15u˙02u02β+48u˙02ε2z3−217u02αε−8u˙0u0α+15u04βε−9u˙0u03β+8u02ε3−22u˙0u0ε2+15u˙02εz4+25u02α+5u04β+16u02ε2−18u˙0u0ε+3u˙02z5−4u05u0ε−2u˙0z6+4u02z7=0.R′0=0. To avoid solving this, the seventh-degree equation, we may set the default value ρ=ε. Taking this value for ρ, we get the following simplified expressions:(104)b1=u0ε+u˙0u02β+α−ε2,b2=−u0βεu0ε+u˙02u02β+α−ε22.Remark 7.
In the integrable case, we haveρ=2ε/3 and then α=8/9ε2 From (99) and (100), μ=0 and m=1/2. Thus, our approach covers the only integrable case for the damped Duffing equation.Example 7.
Let(105)u¨+0.4u˙+u+u3=0,u0=14andu′0=−15.
The error of the approximate analytical solution compared with numerical solution is(106)max0≤t≤20uappt−uRuinge−Kuttat=0.00141579.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure7.Figure 7
Comparison between the approximate analytical solution and the numerical solution for Example7.Example 8.
Let(107)u¨+0.02u˙+u−u3=0,u0=0.1andu′0=0.
The error of the approximate analytical solution compared with numerical solution is(108)max0≤t≤150uappt−uRuinge−Kuttat=0.00135922.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure8.Figure 8
Comparison between the approximate analytical solution and the numerical solution for Example8.Remark 8.
Letρ=ε. An approximate analytical solution to the i.v.p.(109)u¨+2εu˙+αu+βu3=0,u0=u0andu′0=0,is given by(110)ut=u01−u02βε2/2u02β+α−ε22snft|mt2cnft|mt+εu02β+α−ε2dnft|mtsnft|mte−tε,where ft and mt are given by (98) and (99).


### 5.2. Second Case:u0=0

Let(111)u¨+2εu˙+αu+βu3=0,u0=0andu′0=u˙0.Assuming the ansatz [22],(112)ut=−λ1−m0texp−ρtsdf0t,m0t,we will have(113)f0′t2=12βλ2e−2tρ+κ,(114)m0t=−βλ22f0′t2e−2tρ.Then,(115)ft=12ρ4κ+2βλ2−4κ+2βλ2e−2tρ+2κtanh−11+qλ22κe−2tρ−2κtanh−11+qλ22κ,(116)mt=−11+2κ/βλ2e2tρ,κ=p−2ρε+ρ2.The numberλ is found form the initial condition u′0=u˙0, and its value reads(117)λ=−4x˙02β+α−2ερ+ρ22−α+2ερ−ρ2/β2.The numberρ is a solution to some decic equation. Default value is ρ=ε.Example 9.
Let(118)u¨+0.1u˙+u+u3=0,u0=0andu′0=−0.1.
The error of the approximate analytical solution compared with numerical solution is(119)max0≤t≤150uappt−uRuinge−Kuttat=0.000303296.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure9.Figure 9
Comparison between the approximate analytical solution and the numerical solution for Example9.


## 5.1. First Case:u0≠0

Assume the ansatz(96)ut=exp−ρtu0cnft,mt+b1snft,mtdnft,mt1+b2sn2ft,mt.Then, from results in [22],(97)f′t2=α−2ερ+ρ2+βu02e−2tρandmt=βu02e−2tρ2f′t2,so that(98)ft=1ρβu02+μ−βu02e−2tρ+μ+μtanh−11+βu02e−2tρμ−tanh−11+βu02μ,(99)mt=1/21+μe2tρ/βu02,being(100)μ=α−2ερ+ρ2.The numbersb1, b2, and ρ are obtained from the following conditions:(101)u′0=u˙0,R0=0,R′0=0.Solving the two equations in (88) gives(102)b1=u0ρ+u˙0u02β+μ.b2=u03β3α−4ερ−3μ+u˙0u02β2ε−3ρ+2ε−ρα+ρρ−2ε+2u0ρε−ρα+ρρ−2ε2u0u02β+μ2.The numberρ is a root to the septic(103)−2εu02β+α2u02α+u04β+2u˙0u0ε+3u˙02+u02β+α2u02α2+5u04αβ+8u02αε2−4u˙0u0αε+6u˙02α+3u06β2+8u04βε2+2u˙0u03βε+9u˙02u02β+16u˙0u0ε3+24u˙02ε2z−28u02α2ε−4u˙0u0α2+15u04αβε−9u˙0u03αβ+4u02αε3−4u˙0u0αε2+18u˙02αε+7u06β2ε−3u˙0u05β2+4u04βε3+2u˙0u03βε2+21u˙02u02βε+8u˙0u0ε4+12u˙02ε3z2+8u02α2+17u04αβ+32u02αε2−40u˙0u0αε+12u˙02α+5u06β2+24u04βε2−38u˙0u03βε+15u˙02u02β+48u˙02ε2z3−217u02αε−8u˙0u0α+15u04βε−9u˙0u03β+8u02ε3−22u˙0u0ε2+15u˙02εz4+25u02α+5u04β+16u02ε2−18u˙0u0ε+3u˙02z5−4u05u0ε−2u˙0z6+4u02z7=0.R′0=0. To avoid solving this, the seventh-degree equation, we may set the default value ρ=ε. Taking this value for ρ, we get the following simplified expressions:(104)b1=u0ε+u˙0u02β+α−ε2,b2=−u0βεu0ε+u˙02u02β+α−ε22.Remark 7.
In the integrable case, we haveρ=2ε/3 and then α=8/9ε2 From (99) and (100), μ=0 and m=1/2. Thus, our approach covers the only integrable case for the damped Duffing equation.Example 7.
Let(105)u¨+0.4u˙+u+u3=0,u0=14andu′0=−15.
The error of the approximate analytical solution compared with numerical solution is(106)max0≤t≤20uappt−uRuinge−Kuttat=0.00141579.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure7.Figure 7
Comparison between the approximate analytical solution and the numerical solution for Example7.Example 8.
Let(107)u¨+0.02u˙+u−u3=0,u0=0.1andu′0=0.
The error of the approximate analytical solution compared with numerical solution is(108)max0≤t≤150uappt−uRuinge−Kuttat=0.00135922.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure8.Figure 8
Comparison between the approximate analytical solution and the numerical solution for Example8.Remark 8.
Letρ=ε. An approximate analytical solution to the i.v.p.(109)u¨+2εu˙+αu+βu3=0,u0=u0andu′0=0,is given by(110)ut=u01−u02βε2/2u02β+α−ε22snft|mt2cnft|mt+εu02β+α−ε2dnft|mtsnft|mte−tε,where ft and mt are given by (98) and (99).


## 5.2. Second Case:u0=0

Let(111)u¨+2εu˙+αu+βu3=0,u0=0andu′0=u˙0.Assuming the ansatz [22],(112)ut=−λ1−m0texp−ρtsdf0t,m0t,we will have(113)f0′t2=12βλ2e−2tρ+κ,(114)m0t=−βλ22f0′t2e−2tρ.Then,(115)ft=12ρ4κ+2βλ2−4κ+2βλ2e−2tρ+2κtanh−11+qλ22κe−2tρ−2κtanh−11+qλ22κ,(116)mt=−11+2κ/βλ2e2tρ,κ=p−2ρε+ρ2.The numberλ is found form the initial condition u′0=u˙0, and its value reads(117)λ=−4x˙02β+α−2ερ+ρ22−α+2ερ−ρ2/β2.The numberρ is a solution to some decic equation. Default value is ρ=ε.Example 9.
Let(118)u¨+0.1u˙+u+u3=0,u0=0andu′0=−0.1.
The error of the approximate analytical solution compared with numerical solution is(119)max0≤t≤150uappt−uRuinge−Kuttat=0.000303296.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure9.Figure 9
Comparison between the approximate analytical solution and the numerical solution for Example9.


## 6. Damped and Forced Duffing Equation

Let us consider the Duffing equation as originally was introduced by Georg Duffing:(120)x¨+rx˙+ω02x+βx3=Fcosωt,given the initial conditions(121)x0=x0andx′0=x˙0.Let(122)xt=ut+c1cosωt+c2sinωt.We will suppose that the functionu=ut is a solution to the Duffing equation(123)u¨+ru˙+αu+βu3=0,u0=x0−c1andu′0=x˙0−c2ω,where(124)α=123c12β+3c22β+2ω02.The numbersc1 and c2 are chosen, so that(125)432F2β3c13+1152Fr2β2ω2c12−192β3F2βω2−3F2βω02−4r4ω4−4r2ω6+8r2ω02ω4−4r2ω04ω2c1−192Fβ3F2β−4r2ω4+4r2ω02ω2=0(126)432F2β3c23−1152c22Frβ2ωω2−ω02c22+768r2βω2r2ω2+ω4−2ω02ω2+ω04c2−768Fr3βω3=0.Example 10.
Let(127)x¨+0.1x˙+x+x3=0.1cos0.4t,x0=0andx′0=0.
The error of the approximate analytical solution compared with numerical solution is(128)max0≤t≤150xappt−xRuinge−Kuttat=0.00287382.
The approximate solution is(129)xt=e−0.05t−0.117339cn−0.00791022dnsn1+1.411×10−7sn2+0.00551957sin0.4t+0.117339cos0.4t,where(130)cn=cnft,mt,sn=snft,mt,dn=dnft,mt,ft=−20.13431.0182+0.0137685e−0.0993331t−5.07916log1−0.9910231.0182+0.0137685e−0.0993331t2+10.1583log0.9910231.0182+0.0137685e−0.0993331t+1−37.4127,mt=1147.903e0.0993331t+2.
Comparison between the approximate analytical solution and the numerical solution is shown in Figure10.
Finally, let us compare the accuracy of the obtained results in comparation with the homotopy perturbation method (HPM). This method gives the approximate solution:(131)xHPMt=e−1/2tr+κ2κ2r2ω2+ω2−ω022,(132)Fr2ω2−ω02etκ+1−rκω2+ω02etκ−1+4ω02ω02−ω2etκ+1+2Fκ2e12tr+κrωsintω+ω02−ω2costω+r2ω2+ω2−ω022x0rr+κetκ+r−κ−4ω02etκ+1+2x˙0κetκ−1,(133)κ=r2−4ω02.
The error of this approximation compared with numerical solution is(134)Max0≤t≤150xHPMt−xRuinge−Kuttat=0.0136613.
Comparison between the approximate analytical solution and the homotopy solution is shown in Figure11.Figure 10
Comparison between the approximate analytical solution and the numerical solution for Example10.Figure 11
Comparison between the approximate analytical solution and the homotopy solution for Example10.


## 7. Analysis and Discussion

We have solved the undamped and constantly forced Duffing equation exactly. Trigonometric approximant was also provided. For the damped or forced case, we derived approximate analytical solution. As far as we know, the Duffing equation (1) has not been solved using the tools we employed in this work. For the damped unforced case, author in [8] obtained approximate analytical solution using generalized Jacobian elliptic functions. More exactly, author considered the following equation:(135)x¨+2βx˙+αx−εx3=0,x0=x0andx′0=x˙0.The obtained solution in [8] has the form(136)xt=c0exp−βtsnωt+c1,mt,where(137)ωt=12β22α−2β2−c02ε−2α−2β2−c02εe−2βt+2c0εβ2−α1−2α−β2/c02εcsch−1c0ε/2β2−2α2α−2β2−c02ε+eβt2α−2β2−c02εe−2βtcsch−1c0εeβ−t/2β2−2αc02ε2β2−αe2βt/c02ε+1,(138)mt=−c02εc02ε+2β2−αe2βt.The constantsc0 and c1 are determined from the initial conditions as follows:(139)c1=sn−1x0/c0,m0,(140)c0=±αε−β2ε±α2−2αβ2+β4−2αεx02−4βεx0x˙0+ε2x04−2εx˙02ε.This approach is different from the method we used in this work. Let us compare the solution (136) with the solution we obtained in Example 8:(141)u¨+0.02u˙+u−u3=0,u0=0.1andu′0=0for0≤t≤150.Using formula gives the approximate analytical solution(142)xt=0.100005e−0.01tsn50.21.41767−1.9998+0.010001e−0.02t−0.233.4329−99.9901e0.01t1.9998+0.010001e−0.02tcsch−10.0707177e−0.01t1+199.96e0.02t+1.55892,0.010001−0.010001−1.9998e0.02t.The error of this solution compared with the Runge–Kutta numerical method equals 0.00142148. The error obtained in our method equals 0.00135922. We may try a simpler ansatz in the form(143)xt=c0exp−ρtsinwt,with(144)wt=−4α+4ρρ−2β−3c02ε+4α+4ρρ−2β−3c02εe−2ρt+−4α+4ρρ−2β−3c02ε+4α+4ρρ−2β−3c02εe−2ρt+2α+ρρ−2βcoth−12α+ρρ−2β/4α+4ρρ−2β−3c02ε−coth−12α+ρρ−2β/4α+4ρρ−2β−3c02εe−2ρt2ρ.The numbersc0 and ρ are determined from the system(145)−8αβ+8αρ+16β2ρ−24βρ2+6βc02ε−9c02ερ+8ρ3=0,−4αc02+8βc02ρ+3βc04ε−4c02ρ2−3c02εx02+4αx02+8βρx02+8ρ2x02+8ρx0x˙o+4x˙02=0.Using this ansatz, we obtain the approximate analytical solution:(146)xtrigot=0.100005e−0.00996291tsin1.56087+50.18613.9996+0.030003e−0.0199258t+1.9999−coth−11.99993.9996+0.030003e−0.0199258t+3.14134−1.5708i−2.00739The error of the trigonometric solution (146) compared with the numerical solution using the Runge–Kutta method equals 0.00195254, so that the trigonometric solution is good as well.


## 8. Conclusions and Future Work

The methods employed here may be useful to study other nonlinear oscillators of the form(147)x¨+2δx˙+fx=Fcosωt,x0=x0andx′0=x˙0,where the function f is odd: f−x=−fx. To this end, we approximate this function on some interval −A,A by means of Chebyshev polynomials in the form(148)fx≈px+qx3.Then, the i.v.p. is replaced with the i.v.p.(149)x¨+2δx˙+px+qx3=Fcosωt,x0=x0andx′0=x˙0.On the other hand, we may study the following cubic quintic Duffing oscillator:(150)x¨+2δx˙+px+qx3+rx5=Fcosωt,x0=x0andx′0=x˙0.For the unforced and damped cases, we may try the ansatz(151)xt=c0exp−ρtsinft+arccosx0c0,where(152)ft=122∫0t8p−16δρ+8ρ2+6c02e−2τρq+5c04e−4τρrdτ,and(153)8px02−2δρx02+2ρ2x02+2ρx˙0x0+x˙02−2−8δρ+4p+4ρ2−3qx02c02+5rx02−6qc04−5rc06=0.The numberρ is a free parameter that is chosen in order to minimize the residual error



---

*Source: 2357258-2022-05-17.xml*