# Approximate Periodic Solutions for Oscillatory Phenomena Modelled by Nonlinear Differential Equations

Authors: Constantin Bota; Bogdan Căruntu; Olivia Bundău
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We apply the Fourier-least squares method (FLSM) which allows us to find approximate periodic solutions for a very general class of nonlinear differential equations modelling oscillatory phenomena. We illustrate the accuracy of the method by using several significant examples of nonlinear problems including the cubic Duffing oscillator, the Van der Pol oscillator, and the Jerk equations. The results are compared to those obtained by other methods.


---

## Body

## 1. Introduction

Oscillatory phenomena are frequently encountered in various fields of science such as, for example, physics, molecular biology, and many branches of engineering. These oscillatory phenomena can be modelled using nonlinear differential equations. Nonlinear differential equations are one of the most important mathematical tools required for understanding these oscillatory phenomena present in everyday life. As is known, finding exact solutions of nonlinear differential equations is possible only in some particular cases. This justifies the need to resort to approximate methods for the computation of approximate periodic solutions, which in turn could provide important information about the phenomena studied.In the present paper we apply the Fourier-least squares method (FLSM) for the computation of approximate periodic solutions for oscillatory phenomena modelled by nonlinear differential equations of the type(1)x(n)(t)=F(x(n-1)(t),x(n-2)(t),…,x(1)(t),x(t),t)
with the initial conditions
(2)x(k)(0)=Ak,k∈0,1,2,…,n-1,
where F is a nonlinear continuous function, t∈ℝ, Ak∈ℝ.Equation (1) is a very general one, being able to model a large class of oscillatory phenomena. As a consequence, since our method can be applied for (1), it follows that it can be considered a powerful and useful method.We remark that recently there has been much interest in finding approximate periodic solutions of nonlinear differential equations of type (1) and conditions (2) and many approximate methods were proposed for the computation of approximate periodic solutions of equations of this type.Among the methods used to compute such approximate periodic solutions we mention the following.The homotopy perturbation method (see, [1–5]), the variational formulation method (see, [5–7]), the harmonic balance methods (see, [8–13]), the quasilinearization technique (see, [14]), the reproducing kernel space method (see, [15]), the Adomian decomposition method (see, [16]), the parameter-expansion method (see, [17, 18]), the variational iteration methods (see, [19–21]), the energy balance method (see, [20, 22]), the amplitude-frequency formulation (see, [4, 18]), the homotopy analysis method (see, [23–25]), the max-min approach (see, [18]), the optimal homotopy asymptotic method (see, [26]), the residue harmonic balance method (see, [27]), the enhanced cubication method (see, [28]), the linearisation method (see, [29]), perturbation methods (see, [30]), and numerical methods (see, [31–35]).If the problem consisting of (1) and conditions (2) admits a periodic solution, FLSM allows us to determine an accurate approximate solution of this problem. In order to test the accuracy of the method, we apply it to several well-known examples of nonlinear equations and compare the approximate solutions obtained with this method with approximate solutions obtained by other methods.


## 2. The Fourier-Least Squares Method

We suppose that the problem consisting of (1) and conditions (2) admit a periodic solution with the period T and corresponding frequency ω=(2·π)/T. We consider the operator
(3)D(x)=x(n)(t)-F(x(n-1)(t),x(n-2)(t),…,x(1)(t),x(t),t).Ifx~ is an approximate periodic solution of (1), we evaluate the error obtained by replacing the exact solution x with the approximate one x~ as the remainder:
(4)R(t,x~)=D(x~(t)),t∈ℝ.Definition 1.
One calls anFourier-function a function of the form
(5)f(t)=∑k=0m‍[ak·cos⁡(k·ω·t)+bk·sin(k·ω·t)],
where ak,bk,ω∈ℝ, m∈ℕ.We will find approximate Fourier-solutionsx~ of the problem consisting of (1) and conditions (2) on ℝ which satisfy the following conditions:
(6)|R(t,x~)|<ϵ,(7)x~(0)=A0,x~(1)(0)=A1,
where
(8)x~=∑k=0p‍[a~k·cos⁡(k·ω~·t)+b~k·sin(k·ω~·t)].Definition 2.
One calls aϵ-approximate Fourier-solution of the problem consisting of (1) and conditions (2) a Fourier-function x~ satisfying the relations (6) and (7).Definition 3.
One calls aweakδ-approximate Fourier-solution of the problem consisting of (1) and conditions (2) a Fourier-function x~ satisfying the relation
(9)∫02π/ω~R2(t,x~)dt≤δ‍
together with the initial conditions (7).Definition 4.
One calls aFourier-sequence a sequence of Fourier-functions {sp(t)}p∈ℕ:
(10)sp(t)=∑k=0p‍[apk·cos⁡(k·ωp·t)+bpk·sin(k·ωp·t)],
where apk,bpk∈ℝ, bp0=0.Definition 5.
One considers a Fourier-sequence satisfying the conditions(11)sp(0)=A0,sp(1)(0)=A1.
We call the Fourier-sequence {sp(t)}p∈ℕ convergent to the solution of the problem consisting of (1) and conditions (2) if
(12)lim⁡p→∞R(t,sp(t))=0.
We will find a weakϵ-approximate Fourier-solution of the type (8), where p≥1 and the constants a~k, b~k, ω~ are calculated using the following steps.(i)
By substituting the approximate solution (8) in (1) we obtain the following expression:
(13)ℜ(t,a~0,a~1,…,a~p,b~1,…,b~p,ω~)=R(t,x~)=x~(n)(t)=R(t,x~)-F(x~(n-1)(t),x~(n-2)(t),…,x~(1)(t),x~(t),t).(ii)
We attach to the problem consisting of (1) and conditions (2) the following real functional:
(14)J(a~1,…,a~p,b~2,…,b~p,ω~)=∫02π/ω~ℜ2(t,a~0,a~1,…,a~p,b~1,…,b~p,ω~)dt,‍wherea~0, b~1 are computed as functions of a~1,…,a~p, b~2,…,b~p by using the initial conditions (7).(iii)
We compute the values ofa10,…,ap0, b20,…,bp0,ω0 as the values which give the minimum of the functional (14) and the values of a00, b10 again as functions of a10,…,ap0, b20,…,bp0 by using the initial conditions (7).(iv)
Using the constantsa00,a10,…,ap0, b10,…,bp0,ω0 thus determined, we consider the Fourier-sequence
(15)sp(t)=∑k=0p‍[apk·cos⁡(k·ωp·t)+bpk·sin(k·ωp·t)]withapk=ak0 for k=0,…,p, bpk=bk0 for k=1,…,p, and ωp=ω0 (where b00=0).The following convergence theorem holds.Theorem 6.
If the problem consisting of (1) and conditions (2) admits a periodic solution, then the Fourier-sequence sp(t) from (15) satisfies the property
(16)lim⁡p→∞∫02π/ωp‍R2(t,sp(t))dt=0.
Moreover, for all ϵ>0, ∃p0∈ℕ such that for all p∈ℕ, p>p0 it follows that sp(t) is a weak ϵ-approximate Fourier-solution of the problem consisting of (1) and conditions (2).Proof.
From the fact that the problem consisting of (1) and conditions (2) admits a periodic solution it follows that the series
(17)∑k=0∞‍[ak·cos⁡(k·ω·t)+bk·sin(k·ω·t)]
exists and its sequence of partial sums fp(18)fp(t)=∑k=0p‍[ak·cos⁡(k·ω·t)+bk·sin(k·ω·t)]
converges to the solution of the problem consisting of (1) and conditions (2); that is,
(19)lim⁡p→∞R(t,fp(t))=0.
Based on the way the Fourier-function sp(t) is computed and taking into account the relations (8)–(15), the following inequality holds:
(20)0≤∫02π/ωp‍R2(t,sp(t))dt≤∫02π/ωp‍R2(t,fp(t))dt,∀p∈ℕ.
It follows that
(21)0≤lim⁡p→∞∫02π/ωp‍R2(t,sp(t))dt≤lim⁡p→∞∫02π/ωp‍R2(t,fp(t))dt=0,∀p∈ℕ.
We obtain
(22)lim⁡p→∞∫02π/ωp‍R2(t,sp(t))dt=0.
From this limit we obtain that for all ϵ>0, ∃p0∈ℕ such that for all p∈ℕ, p>p0 it follows that sp(t) is a weak ϵ-approximate Fourier-solution of the problem consisting of (1) and conditions (2).Remark 7.
Anyϵ-approximate Fourier-solution of the problem consisting of (1) and conditions (2) is also a weak approximate Fourier-solution, but the opposite is not always true. It follows that the set of weak approximate Fourier-solutions of the problem consisting of (1) and conditions (2) also contains the approximate Fourier-solutions of the problem.Taking into account the above remark, in order to findϵ-approximate Fourier-solutions of the problem consisting of (1) and conditions (2) by the Fourier-Least Squares Method, we will first determine weak approximate Fourier-solutions, x~. If |R(t,x~)|<ϵ, then x~ is also an ϵ-approximate Fourier-solution of the problem.


## 3. Applications

The test problems included this section are the Duffing oscillator (two cases, an autonomous one and one involving integral forcing terms) and the Jerk equations.These problems were extensively studied over the years, and various solutions, both approximate analytical ones and numerical ones, were proposed.The qualitative properties of these oscillators were also extensively studied. Stability and bifurcation studies for the Duffing oscillators include [36–38] among many others. For Jerk-type equations a comprehensive bifurcation study can be found in [39] and a study of the limit cycles can be found in [40]. Following the computations presented in these papers, corresponding conclusions can be drawn for the problems studied in the Sections 3.1–3.3. For example, in the case of the autonomous Duffing oscillator (23), for the values of the parameters considered in the computations (ϵ=1.25 and A=2), a quick computation similar to the one in [36] indicates that the only equilibrium point is the origin, which is a center; similar computations can be performed for the other problems.In the following we compute approximate solutions for the Duffing oscillator and the Jerk equations and compare our results with similar analytical approximations previously computed by using other methods.


### 3.1. Application 1: The Autonomous Duffing Oscillator

Our first test problem is the autonomous Duffing oscillator:(23)x(2)(t)+x(t)+ϵ·x3(t)=0,x(0)=A,x(1)(0)=0.The Duffing oscillator is extensively studied in literature and some relatively recent results are presented in [12, 19]. In [12] approximate solutions are computed using the rational harmonic balance method (RHB) and in [19] approximate solutions of (23) are computed using a variational iteration procedure (VI). The approximate frequency of the oscillations obtained by using these methods is compared with the exact period known in the literature.In the following, in order to obtain our approximation of the frequency and of the solution, we will perform the steps described in the previous section. The computations were performed using the SAGE open source software (version 5.5, available athttp://www.sagemath.org/).We will perform in detail our computations for the valuesϵ=1.25 and A=2.Since the computations of the minimum of the functional (14) are relatively difficult for large values of p in (8), we will actually use an iterative procedure, starting with p=1 and increasing the value until we achieve the desired precision.


#### 3.1.1. Approximate Solution forp=1

Taking into account these considerations, first we choose the approximate solution (8) of the form
(24)x~1(t)=a~0+a~1·cos⁡(ω~1·t)+b~1·sin(ω~1·t).In Step 1, the expression (13) becomes
(25)ℜ(t,a~0,a~1,b1,ω~1)=a~1ω~12cos⁡(ω~1·t)-b~1ω~12sin(ω~1·t)+54(a~1cos⁡(ω~1·t)+b~1sin(ω~1·t)+a~0)3+a~1cos⁡(ω~1·t)+b~1sin(ω~1·t)+a~0.Taking into account the initial conditionsx~(0)=A, x~(1)(t)=0 we obtain the relations
(26)a~0=-a~1+2,b~1=0.Replacing these values, the corresponding functional (14) from Step 2 is
(27)J(a~1,ω~1)=(36864π+5775πa~16+128πa~12ω~14-37800πa~15+107800πa~14)·1128ω~1+(-172800πa~13-16×(75πa~14-240πa~13+256πa~12)ω~12+167424πa~12-98304πa~1)·1128ω~1.In Step 3 we must compute the minimum ofJ with respect to a~1 and ω~1. For relatively simple problems such as this it is possible to compute directly the critical points of J and subsequently select the value corresponding to the minimum.In general, the critical points corresponding to the functional (14) J(a~1,…,a~p,b~2,…,b~p,ω~) are the solution of the system:
(28)∂J∂a~1=0,⋮∂J∂a~p=0,∂J∂b~2=0,⋮∂J∂b~p=0,∂J∂ω~1=0.Forp=1 the system becomes
(29)-(49152π-17325πa~15-128πa~1ω~4+94500πa~14-215600πa~13)·164ω~1+(32(75πa~13-180πa~12+128πa~1)ω~2+259200πa~12-167424πa~1)·164ω~1=0,(30)(16πa~12w3-(75πa~14-240πa~13+256πa~12)ω~)·14ω~1-(36864π+5775πa~16+128πa~12ω~4-37800πa~15+107800πa~14-172800πa~13)·1128ω~12+(-16(75πa~14-240πa~13+256πa~12)ω~2+167424πa~12-98304πa~1)·1128ω~12=0.This system can be solved directly. We used the “solve” command in SAGE and, after we excluded the complex solutions, we found the critical points(31)a~1=1.95861253449,ω~1=2.16432748538,a~1=1.95861253449,ω~1=-2.16432748538,a~1=0.181528334876,ω~1=-6.92142088267,a~1=0.181528334876,ω~1=6.92142088267.In order to find the minimum, we use the second partial derivative test, which is easy enough to implement in SAGE, and find thata~1=a10=1.95861253449, ω~1=ω10=2.16432748538 is the local minima. Using again the relations a00=-a10+2,b10=0, we replace all the values in the expression of the approximate solution (8) and we obtain our first approximation:
(32)x~1(t)=1.95861253449cos⁡(2.16432748538t)+0.0413874655104.In Figure1 we present the comparison between our x~1 approximate solution (solid line) and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 1
Comparison betweenx~1 (solid line) and numerical solution (dotted line).As we already saw, the approximate value of the frequency obtained here isω~1=ω10=2.16432748538 and it is already close to the exact value which in this case is ωex=2.150416169536.We observe that, generally speaking, if the value ofp in the expression of x~(t) is larger than 1, then the solution of the corresponding system of equations which gives the critical points (28) cannot be found directly. In the particular case of SAGE, the command “solve” fails to find the solutions, exiting with some kind of error message.In this situation it is still possible to find good approximations of the solutions of the problem solving the system (28) by means of a numerical method. More precisely, we can find approximate solutions for the given problem (23) solving (28) by means of a SAGE implementation of the well-known Newton method.


#### 3.1.2.  Approximate Solution forp>1

As the following results will show, the Newton method is able to find approximate solutions of (28) which can lead to highly accurate approximate solutions of the problem (23).Forp=2 the approximate solution (8) has the form
(33)x~2(t)=a~0+a~1·cos⁡(ω~2·t)+b~1·sin(ω~2·t)+a~2·cos⁡(2·ω~2·t)+b~2·sin(2·ω~2·t).After we compute the corresponding expressions ofℜ(t,a~0,a~1,a~2,b~1,b~2,ω~2) and J(a~1,a~2,b~2,ω~2) and of the system (28) (all too large to insert here), we apply Newton’s method taking as the starting point of the iteration (a~1,s,a~2,s,b~2,s,ω~2,s), where a~1,s and ω~2,s are the values of a10 and ω10 computed for the previous approximation p=1, namely, a~1,s=1.95861253449 and ω~2,s=2.16432748538. In order for the sequence of approximations given by Newton’s method to converge to the solution(s) of the system (28), a~2,s and b~2,s will take successively values on a given grid of the type G=Ia~2×Ib~2, where Ii is a division of a symmetric interval centered in zero.For all the test problems included in this paper and for all the values ofp tested, a grid of the form G={-1,0.9,-0.8,…,1}×{-1,0.9,-0.8,…,1} (i.e., from −1 to 1 with step size 0.1) is large enough in the sense that if the starting point (s0,s1) scans G we can obtain using Newton's method the desired solutions of (28).In the particular case of the problem (23) it was actually sufficient to choose a~2,s=0 and b~2,s=0 and the sequence converged to the minimum of J.Forp=3 we repeat the same procedure: we compute the corresponding expressions of ℜ(t,a~0,a~1,a~2,a~3,b~1,b~2,b~3,ω~3) and J(a~1,a~2,a~3,b~2,b~3,ω~3) and of the system (28) and we apply Newton's method taking as the starting point of the iteration (a~1,s,a~2,s,a~3,s,b~2,s,b~3,s,ω~3,s), where (a~1,s,a~2,s,b~3,s,ω~3,s) are the minimum values (a10,a20,b20,ω20) computed for the previous approximation p=2 and a3,s=b3,s=0.The process can be carried on for increased values ofp until the desired accuracy is reached.While at a first look the computations may seem long and tedious, by using a software program such as SAGE, we were actually able to perform them easy and quick, obtaining a very good accuracy. Thus, forp=7, we obtained for the approximate frequency the value ω70=2.15041853722 and the corresponding approximate solution (8) has the form
(34)x~7(t)=1.0618674103207715·10-08+1.930547907025286·cos⁡(2.150418537217718·t)+1.5538232860698776·10-08·cos⁡(4.300837074435436·t)+0.0671185574473875·cos⁡(6.4512556116531545·t)+1.4239017145731592·10-09·cos⁡(8.601674148870872·t)+0.0022576839638690305·cos⁡(10.75209268608859·t)+1.3333451751464932·10-10·cos⁡(12.902511223306309·t)+7.582384931426247·10-05·cos⁡(15.052929760524027·t).In Figure2 we present the comparison between the x~7 approximate solution (solid line) and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 2
Comparison betweenx~7 (solid line) and numerical solution (dotted line).In [19] approximate solutions for the problem (23) were computed for the cases ϵ=1.25,A=2 (studied above) and ϵ=250, A=2.For the caseϵ=1.25, A=2, the approximate solution from [19] is
(35)xVI⁡(t)=1.989351555·cos⁡(2.15031·t)+0.008334403·cos⁡(6.45093·t)+0.002206664·cos⁡(10.75155·t)+0.000053689·cos⁡(15.05217·t)+(2.6·10-7)·cos⁡(19.35279·t).Table1 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xVI⁡ (see, [19]) and x~7 for the case ϵ=1.25, A=2.Table 1
Comparison of the absolute errors of the approximate solutions for problem (23) in the case ϵ=1.25,A=2.
x
VI
⁡
x
~
7
t
=
0
5.34289999998
·
1
0
-
5
0.00000000000
t
=
0.4
8.81935946956
·
1
0
-
2
3.92745813649
·
1
0
-
7
t
=
0.8
3.40733004552
·
1
0
-
2
1.66866686269
·
1
0
-
6
t
=
1.2
5.63129918651
·
1
0
-
2
5.23059939095
·
1
0
-
6
t
=
1.6
1.96411171977
·
1
0
-
2
1.89375331883
·
1
0
-
6
t
=
2
7.93642203903
·
1
0
-
2
6.81658032309
·
1
0
-
6
t
=
2.4
8.23170554343
·
1
0
-
2
1.36265840209
·
1
0
-
5
t
=
2.8
1.50931697081
·
1
0
-
2
9.46274851943
·
1
0
-
6For the caseϵ=250, A=2, the approximate solution from [19] is
(36)xVI⁡(t)=1.977507645·cos⁡(26.8026·t)+0.018844215·cos⁡(80.4078·t)+0.003619986·cos⁡(134.013·t)+0.000027438·cos⁡(187.6182·t)+0.000164633·cos⁡(241.2234·t).For this case, the approximate solution computed using our method is(37)x~7(t)=3.22758664474·10-08+1.91015330733·cos⁡(26.810846231·t)+(6.9616769357·10-08)·cos⁡(53.6216924619·t)+0.0859766881274·cos⁡(80.4325386929·t)+(7.78739552572·10-09)·cos⁡(107.243384924·t)+0.00371019879647·cos⁡(134.054231155·t)+(7.82719602479·10-10)·cos⁡(160.865077386·t)+0.0001596952849·cos⁡(187.675923617·t).Table2 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xVI⁡ (see, [19]) and x~7 for the case ϵ=250, A=2.Table 2
Comparison of the absolute errors of the approximate solutions for problem (23) in the case ϵ=250,A=2.
x
VI
⁡
x
~
7
t
=
0
1.63917000000
·
1
0
-
4
0.00000000000
t
=
0.04
9.94430509655
·
1
0
-
2
1.59432943392
·
1
0
-
6
t
=
0.08
1.01602511979
·
1
0
-
1
2.50754265524
·
1
0
-
5
t
=
0.12
1.84880609089
·
1
0
-
3
1.44364392085
·
1
0
-
6
t
=
0.16
9.38684852713
·
1
0
-
2
2.52297827337
·
1
0
-
5
t
=
0.2
1.00379880594
·
1
0
-
1
4.54516346458
·
1
0
-
5Finally, in Table3, we compare the approximate values for the frequency ω70 computed using our method with approximate values computed in [19] (ωVI⁡) and [12] (ωRHB). The comparison is made by means of the percentage error, which for a given approximate error ωapprox is defined as eωapprox=100·(|ωexact-ωapprox|/ωexact), where ωexact is the corresponding exact error. It is easy to see that our approximations are far better than the ones previously computed and they remain accurate even for the case of a very strong nonlinearity.Table 3
Comparison of approximate errors for the frequencyω for problem (23).
ϵ
·
A
2
5
100
1000
10000
ω
exact
2.1504161695
8.5335861889
26.8107384581
84.7274799361
ω
VI
⁡
2.15031
x
26.8026
84.7015
e
ω
VI
⁡ (%)
0.004937
x
0.0303552
0.0306629
ω
RHB
2.15045
8.53402
x
x
e
ω
RHB (%)
0.0015
0.0051
x
x
ω
7
2.15041853722
8.53361825153
26.810846231
84.7278228489
e
ω
7 (%)
0.0001101034
0.0003757217
0.0004019763
0.0004047244


### 3.2. Application 2: The Duffing Equation Involving Integral Forcing Terms

Our second test problem is(38)x(2)(t)+x(t)·x(1)(t)+∫0t‍tsx2(s)ds=f(t),hhhhhhhhhhhhhhhhhhhhhhhhhhh0<t<1,x(0)-x(1)(0)=0,x(1)+x(1)(1)=0,
where f(t)=-3t-3t2+5t3/2+2t4/3-t5/4-2t6/5+t7/6.The problem (38) (see, [14, 15]) is a version of the well-known Duffing equation involving both integral and nonintegral forcing terms with separated boundary conditions. This equation has been studied in a series of recent papers including [14, 15].In [14], the authors applied a generalized quasilinearization technique to prove the existence and uniqueness of the solution of Duffing equation involving both integral and nonintegral forcing terms. They showed that there are sequences of approximate solutions converging monotonically and quadratically to the unique solution of the problem.In [15], the authors gave a representation of exact solution and approximate solution of Duffing equation involving both integral and nonintegral forcing terms in the reproducing kernel space (RKS). They represented the exact solution in the form of a series and they showed that the n-term approximation of the exact solution converges to the exact solution.Next we present our results for (38) using FLSM. Also, we will compare these results with those obtained in [15].Thus, forp=3, we obtained the fact that the approximate periodic solution (8) has the form
(39)x~3(t)=1.45453sin(t)-0.26409sin(2t)+0.0245502sin(3t)+2.6625cos⁡(t)-0.169569cos⁡(2t)+0.00174065cos⁡(3t)-1.49467.Since in [15] only the numerical results are presented while the expression of approximate solution expression is not, we cannot perform a direct graphical comparison of our approximate solution with the corresponding solution from [15].Therefore, in Table4, we present the comparison of several values of the absolute errors (computed as the difference in absolute value between the exact solution and the approximate solution) corresponding to the approximate solutions xRKS from [15] and to our approximate solutions for t=0.1,0.2,…,1, as given in [15].Table 4
Comparison of the absolute errors of the approximate solutions for problem (38).
x
RKS
x
~
3
t
=
0
0
0
t
=
0.1
9.46600
·
1
0
-
6
9.85192
·
1
0
-
8
t
=
0.2
9.77708
·
1
0
-
6
9.2974
·
1
0
-
8
t
=
0.3
9.81938
·
1
0
-
6
1.58651
·
1
0
-
7
t
=
0.4
9.66575
·
1
0
-
6
5.76015
·
1
0
-
8
t
=
0.5
9.37204
·
1
0
-
6
2.07316
·
1
0
-
7
t
=
0.6
8.97794
·
1
0
-
6
5.9859
·
1
0
-
8
t
=
0.7
8.50917
·
1
0
-
6
1.55722
·
1
0
-
7
t
=
0.8
7.97984
·
1
0
-
6
9.17107
·
1
0
-
8
t
=
0.9
7.39469
·
1
0
-
6
9.65862
·
1
0
-
8
t
=
1
6.75105
·
1
0
-
6
7.84328
·
1
0
-
10


### 3.3. Application 3: The Jerk Equation Containing Velocity-Cubed and Velocity Times Displacement-Squared

Our last test is a Jerk nonlinear equation, which describes several physical problems using mechanical oscillations of the third order. The most general form of the Jerk nonlinear equations, which contains the third temporal derivative of displacement, is(40)x(3)+α·x·x(1)·x(2)+β·x(1)·x(2)2+δ·x2·x(1)+ɛ·x(1)3+γ·x(1)=0,
where the parameters α,β,δ,ɛ, and γ are constants.Nonlinear Jerk equations (40) are intensely studied by several authors in the literature and some recent results are presented in [3, 8, 9, 27].We consider the following Jerk equation:(41)x(3)+x2·x(1)+x(1)3=0,x(0)=0,x(1)(0)=B,x(2)(0)=0.Recently, Ma et al. in [3], using the Homotopy Perturbation Method, obtained high-order analytic approximate periods and periodic solutions of the Jerk Equation (41). In [8, 9], Gottlieb used the lowest-order harmonic balance method to determine analytical approximations to the periodic solution of the Jerk equations. Also, Leung and Guo in [27] obtained approximations for the angular frequency and the limit cycle for (41) based on the residue harmonic balance approach.ForB=0.5, the approximate periodic solution xHPM from [3] is
(42)xHPM(t)=0.8708027·sin(0.61924·t)-0.0106486·sin(1.85772·t).Applying FLSM we computed an approximate periodic solution of the problem (41) of the same order (containing terms up to sin(3·ω·t)):
(43)x~3(t)=0.847907380098·sin(0.615691681067·t)-(3.59275929188×10-05)·sin(1.23138336213·t)-0.0119135817912·sin(1.8470750432·t).
The approximate frequency and period are ω30=0.615691681067, T30=10.2050839736674 with an error of eω30=-0.0555984645279972.In Figure3, we can visualize and compare our approximate solution (solid line), the approximate solution xHPM from [3] (dashed line), and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 3
Comparison between our solutionx~3(solid line), solution xHPM in [3] (dashed line), and numerical solution (dotted line).Table5 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solution xHPM from [3] and our approximate solution x~3 for the case B=0.5.Table 5
Comparison of the absolute errors corresponding to the solutionxHPM from [3] and to our solution x~3.
x
HPM
x
~
3
t
=
0
0
0
t
=
1.0
1.72354316028
·
1
0
-
2
2.11130244099
·
1
0
-
4
t
=
2.0
2.3307287671
·
1
0
-
2
1.31078376172
·
1
0
-
4
t
=
3.0
1.79768181549
·
1
0
-
2
3.46325784747
·
1
0
-
4
t
=
4.0
4.72752471124
·
1
0
-
3
8.73925906836
·
1
0
-
4
t
=
5.0
1.46107699703
·
1
0
-
2
1.42164922819
·
1
0
-
3
t
=
6.0
3.05145610435
·
1
0
-
2
1.51326363254
·
1
0
-
3
t
=
7.0
3.04290934289
·
1
0
-
2
8.83072635498
·
1
0
-
4
t
=
8.0
1.45339383549
·
1
0
-
2
7.00497125303
·
1
0
-
4
t
=
9.0
6.86925713324
·
1
0
-
3
2.00167777983
·
1
0
-
3
t
=
10.0
2.91276307352
·
1
0
-
2
2.84571545126
·
1
0
-
3Also, forB=0.5, the approximate periodic solution xRHB from [27] is given by
(44)xRHB(t)=0.8473604·sin(0.6151734·t)-0.01168397·sin(1.8455202·t)+0.000094103·sin(3.075867·t).Applying FLSM we computed an approximate periodic solution of the problem (41) of the same order (containing terms up to sin(5·ω·t)):
(45)x~5(t)=0.847966627859·sin(0.61535099365·t)+(2.70675259884·10-08)·sin⁡(1.2307019873·t)-0.0119562711384·sin⁡(1.84605298095·t)-(9.51418832279·10-11)·sin(2.4614039746·t)+(8.93050847544·10-05)·sin(3.07675496825·t).The approximate frequency and period areω50=0.61535099365, T50=10.2107339908715, with an error of -0.000264516312410058.Table6 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xRHB from [27] and our approximate solutions x~5 for the case B=0.5. In this case we omitted the graphical representation of the approximate solutions xRHB and x~5 since they are both very close to the numerical solution.Table 6
Comparison of the absolute errors corresponding to the solutionxRHB from [27] and to our solution x~5.
x
RHB
x
~
5
t
=
0
0
0
t
=
1.0
2.07336138332
·
1
0
-
4
4.76705901847
·
1
0
-
6
t
=
2.0
8.18802985278
·
1
0
-
4
7.01416299886
·
1
0
-
6
t
=
3.0
6.29650984827
·
1
0
-
4
8.81003685005
·
1
0
-
7
t
=
4.0
3.27159888572
·
1
0
-
4
1.19800236493
·
1
0
-
5
t
=
5.0
7.15072428261
·
1
0
-
4
2.0505386494
·
1
0
-
5
t
=
6.0
7.95019646471
·
1
0
-
4
2.11112229619
·
1
0
-
5
t
=
7.0
1.09803461819
·
1
0
-
3
1.31470477691
·
1
0
-
5
t
=
8.0
5.35184978947
·
1
0
-
4
5.37333608552
·
1
0
-
6
t
=
9.0
8.0197576264
·
1
0
-
4
2.75784654626
·
1
0
-
5
t
=
10.0
1.42381620457
·
1
0
-
3
3.9140208578
·
1
0
-
5It is easy to see from the Tables5 and 6 that the solutions obtained by using FLSM are more accurate than the ones computed by using other methods. This fact is emphasized by Table 7 which presents a comparison of approximate periods computed in several papers, as presented in [27]. In this table, Tk denotes an approximate period obtained by means of an approximate solution containing terms up to sin(k·ω·t).Table 7
Comparison of approximate periods.
B
=
0.5
Period
Error
T
exact
10.210761
0
FLSMT1
10.1976548730449
-
0.128356025129356
FLSMT3
10.2050839736674
-
0.0555984645279972
FLSMT5
10.2107339908715
-
0.000264516312410058
Wu et al. [41] T2
10.1884
-
0.220
Wu et al. [41] T3
10.2107
-
0.001
Leung and Guo [27] T2
10.213682
0.0286
Feng and Chen [23] T2
10.23452
0.232
Feng and Chen [23] T3
10.1935
-
0.169
Feng and Chen [23] T4
10.2029
-
0.077
Hu [42] T2
10.247682
0.362
Hu [42] T3
10.226464
0.154
Hu et al. [43] T2
10.244058
0.326
Ramos [44] T2
10.80674
5.84
Ramos [44] T3
10.55034
3.33
Ramos [44] T2′
10.24328
0.318
Ramos [44] T2′′
10.26079
0.49
Ma et al. [3] T2
10.146603
-
2.27


## 3.1. Application 1: The Autonomous Duffing Oscillator

Our first test problem is the autonomous Duffing oscillator:(23)x(2)(t)+x(t)+ϵ·x3(t)=0,x(0)=A,x(1)(0)=0.The Duffing oscillator is extensively studied in literature and some relatively recent results are presented in [12, 19]. In [12] approximate solutions are computed using the rational harmonic balance method (RHB) and in [19] approximate solutions of (23) are computed using a variational iteration procedure (VI). The approximate frequency of the oscillations obtained by using these methods is compared with the exact period known in the literature.In the following, in order to obtain our approximation of the frequency and of the solution, we will perform the steps described in the previous section. The computations were performed using the SAGE open source software (version 5.5, available athttp://www.sagemath.org/).We will perform in detail our computations for the valuesϵ=1.25 and A=2.Since the computations of the minimum of the functional (14) are relatively difficult for large values of p in (8), we will actually use an iterative procedure, starting with p=1 and increasing the value until we achieve the desired precision.


### 3.1.1. Approximate Solution forp=1

Taking into account these considerations, first we choose the approximate solution (8) of the form
(24)x~1(t)=a~0+a~1·cos⁡(ω~1·t)+b~1·sin(ω~1·t).In Step 1, the expression (13) becomes
(25)ℜ(t,a~0,a~1,b1,ω~1)=a~1ω~12cos⁡(ω~1·t)-b~1ω~12sin(ω~1·t)+54(a~1cos⁡(ω~1·t)+b~1sin(ω~1·t)+a~0)3+a~1cos⁡(ω~1·t)+b~1sin(ω~1·t)+a~0.Taking into account the initial conditionsx~(0)=A, x~(1)(t)=0 we obtain the relations
(26)a~0=-a~1+2,b~1=0.Replacing these values, the corresponding functional (14) from Step 2 is
(27)J(a~1,ω~1)=(36864π+5775πa~16+128πa~12ω~14-37800πa~15+107800πa~14)·1128ω~1+(-172800πa~13-16×(75πa~14-240πa~13+256πa~12)ω~12+167424πa~12-98304πa~1)·1128ω~1.In Step 3 we must compute the minimum ofJ with respect to a~1 and ω~1. For relatively simple problems such as this it is possible to compute directly the critical points of J and subsequently select the value corresponding to the minimum.In general, the critical points corresponding to the functional (14) J(a~1,…,a~p,b~2,…,b~p,ω~) are the solution of the system:
(28)∂J∂a~1=0,⋮∂J∂a~p=0,∂J∂b~2=0,⋮∂J∂b~p=0,∂J∂ω~1=0.Forp=1 the system becomes
(29)-(49152π-17325πa~15-128πa~1ω~4+94500πa~14-215600πa~13)·164ω~1+(32(75πa~13-180πa~12+128πa~1)ω~2+259200πa~12-167424πa~1)·164ω~1=0,(30)(16πa~12w3-(75πa~14-240πa~13+256πa~12)ω~)·14ω~1-(36864π+5775πa~16+128πa~12ω~4-37800πa~15+107800πa~14-172800πa~13)·1128ω~12+(-16(75πa~14-240πa~13+256πa~12)ω~2+167424πa~12-98304πa~1)·1128ω~12=0.This system can be solved directly. We used the “solve” command in SAGE and, after we excluded the complex solutions, we found the critical points(31)a~1=1.95861253449,ω~1=2.16432748538,a~1=1.95861253449,ω~1=-2.16432748538,a~1=0.181528334876,ω~1=-6.92142088267,a~1=0.181528334876,ω~1=6.92142088267.In order to find the minimum, we use the second partial derivative test, which is easy enough to implement in SAGE, and find thata~1=a10=1.95861253449, ω~1=ω10=2.16432748538 is the local minima. Using again the relations a00=-a10+2,b10=0, we replace all the values in the expression of the approximate solution (8) and we obtain our first approximation:
(32)x~1(t)=1.95861253449cos⁡(2.16432748538t)+0.0413874655104.In Figure1 we present the comparison between our x~1 approximate solution (solid line) and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 1
Comparison betweenx~1 (solid line) and numerical solution (dotted line).As we already saw, the approximate value of the frequency obtained here isω~1=ω10=2.16432748538 and it is already close to the exact value which in this case is ωex=2.150416169536.We observe that, generally speaking, if the value ofp in the expression of x~(t) is larger than 1, then the solution of the corresponding system of equations which gives the critical points (28) cannot be found directly. In the particular case of SAGE, the command “solve” fails to find the solutions, exiting with some kind of error message.In this situation it is still possible to find good approximations of the solutions of the problem solving the system (28) by means of a numerical method. More precisely, we can find approximate solutions for the given problem (23) solving (28) by means of a SAGE implementation of the well-known Newton method.


### 3.1.2.  Approximate Solution forp>1

As the following results will show, the Newton method is able to find approximate solutions of (28) which can lead to highly accurate approximate solutions of the problem (23).Forp=2 the approximate solution (8) has the form
(33)x~2(t)=a~0+a~1·cos⁡(ω~2·t)+b~1·sin(ω~2·t)+a~2·cos⁡(2·ω~2·t)+b~2·sin(2·ω~2·t).After we compute the corresponding expressions ofℜ(t,a~0,a~1,a~2,b~1,b~2,ω~2) and J(a~1,a~2,b~2,ω~2) and of the system (28) (all too large to insert here), we apply Newton’s method taking as the starting point of the iteration (a~1,s,a~2,s,b~2,s,ω~2,s), where a~1,s and ω~2,s are the values of a10 and ω10 computed for the previous approximation p=1, namely, a~1,s=1.95861253449 and ω~2,s=2.16432748538. In order for the sequence of approximations given by Newton’s method to converge to the solution(s) of the system (28), a~2,s and b~2,s will take successively values on a given grid of the type G=Ia~2×Ib~2, where Ii is a division of a symmetric interval centered in zero.For all the test problems included in this paper and for all the values ofp tested, a grid of the form G={-1,0.9,-0.8,…,1}×{-1,0.9,-0.8,…,1} (i.e., from −1 to 1 with step size 0.1) is large enough in the sense that if the starting point (s0,s1) scans G we can obtain using Newton's method the desired solutions of (28).In the particular case of the problem (23) it was actually sufficient to choose a~2,s=0 and b~2,s=0 and the sequence converged to the minimum of J.Forp=3 we repeat the same procedure: we compute the corresponding expressions of ℜ(t,a~0,a~1,a~2,a~3,b~1,b~2,b~3,ω~3) and J(a~1,a~2,a~3,b~2,b~3,ω~3) and of the system (28) and we apply Newton's method taking as the starting point of the iteration (a~1,s,a~2,s,a~3,s,b~2,s,b~3,s,ω~3,s), where (a~1,s,a~2,s,b~3,s,ω~3,s) are the minimum values (a10,a20,b20,ω20) computed for the previous approximation p=2 and a3,s=b3,s=0.The process can be carried on for increased values ofp until the desired accuracy is reached.While at a first look the computations may seem long and tedious, by using a software program such as SAGE, we were actually able to perform them easy and quick, obtaining a very good accuracy. Thus, forp=7, we obtained for the approximate frequency the value ω70=2.15041853722 and the corresponding approximate solution (8) has the form
(34)x~7(t)=1.0618674103207715·10-08+1.930547907025286·cos⁡(2.150418537217718·t)+1.5538232860698776·10-08·cos⁡(4.300837074435436·t)+0.0671185574473875·cos⁡(6.4512556116531545·t)+1.4239017145731592·10-09·cos⁡(8.601674148870872·t)+0.0022576839638690305·cos⁡(10.75209268608859·t)+1.3333451751464932·10-10·cos⁡(12.902511223306309·t)+7.582384931426247·10-05·cos⁡(15.052929760524027·t).In Figure2 we present the comparison between the x~7 approximate solution (solid line) and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 2
Comparison betweenx~7 (solid line) and numerical solution (dotted line).In [19] approximate solutions for the problem (23) were computed for the cases ϵ=1.25,A=2 (studied above) and ϵ=250, A=2.For the caseϵ=1.25, A=2, the approximate solution from [19] is
(35)xVI⁡(t)=1.989351555·cos⁡(2.15031·t)+0.008334403·cos⁡(6.45093·t)+0.002206664·cos⁡(10.75155·t)+0.000053689·cos⁡(15.05217·t)+(2.6·10-7)·cos⁡(19.35279·t).Table1 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xVI⁡ (see, [19]) and x~7 for the case ϵ=1.25, A=2.Table 1
Comparison of the absolute errors of the approximate solutions for problem (23) in the case ϵ=1.25,A=2.
x
VI
⁡
x
~
7
t
=
0
5.34289999998
·
1
0
-
5
0.00000000000
t
=
0.4
8.81935946956
·
1
0
-
2
3.92745813649
·
1
0
-
7
t
=
0.8
3.40733004552
·
1
0
-
2
1.66866686269
·
1
0
-
6
t
=
1.2
5.63129918651
·
1
0
-
2
5.23059939095
·
1
0
-
6
t
=
1.6
1.96411171977
·
1
0
-
2
1.89375331883
·
1
0
-
6
t
=
2
7.93642203903
·
1
0
-
2
6.81658032309
·
1
0
-
6
t
=
2.4
8.23170554343
·
1
0
-
2
1.36265840209
·
1
0
-
5
t
=
2.8
1.50931697081
·
1
0
-
2
9.46274851943
·
1
0
-
6For the caseϵ=250, A=2, the approximate solution from [19] is
(36)xVI⁡(t)=1.977507645·cos⁡(26.8026·t)+0.018844215·cos⁡(80.4078·t)+0.003619986·cos⁡(134.013·t)+0.000027438·cos⁡(187.6182·t)+0.000164633·cos⁡(241.2234·t).For this case, the approximate solution computed using our method is(37)x~7(t)=3.22758664474·10-08+1.91015330733·cos⁡(26.810846231·t)+(6.9616769357·10-08)·cos⁡(53.6216924619·t)+0.0859766881274·cos⁡(80.4325386929·t)+(7.78739552572·10-09)·cos⁡(107.243384924·t)+0.00371019879647·cos⁡(134.054231155·t)+(7.82719602479·10-10)·cos⁡(160.865077386·t)+0.0001596952849·cos⁡(187.675923617·t).Table2 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xVI⁡ (see, [19]) and x~7 for the case ϵ=250, A=2.Table 2
Comparison of the absolute errors of the approximate solutions for problem (23) in the case ϵ=250,A=2.
x
VI
⁡
x
~
7
t
=
0
1.63917000000
·
1
0
-
4
0.00000000000
t
=
0.04
9.94430509655
·
1
0
-
2
1.59432943392
·
1
0
-
6
t
=
0.08
1.01602511979
·
1
0
-
1
2.50754265524
·
1
0
-
5
t
=
0.12
1.84880609089
·
1
0
-
3
1.44364392085
·
1
0
-
6
t
=
0.16
9.38684852713
·
1
0
-
2
2.52297827337
·
1
0
-
5
t
=
0.2
1.00379880594
·
1
0
-
1
4.54516346458
·
1
0
-
5Finally, in Table3, we compare the approximate values for the frequency ω70 computed using our method with approximate values computed in [19] (ωVI⁡) and [12] (ωRHB). The comparison is made by means of the percentage error, which for a given approximate error ωapprox is defined as eωapprox=100·(|ωexact-ωapprox|/ωexact), where ωexact is the corresponding exact error. It is easy to see that our approximations are far better than the ones previously computed and they remain accurate even for the case of a very strong nonlinearity.Table 3
Comparison of approximate errors for the frequencyω for problem (23).
ϵ
·
A
2
5
100
1000
10000
ω
exact
2.1504161695
8.5335861889
26.8107384581
84.7274799361
ω
VI
⁡
2.15031
x
26.8026
84.7015
e
ω
VI
⁡ (%)
0.004937
x
0.0303552
0.0306629
ω
RHB
2.15045
8.53402
x
x
e
ω
RHB (%)
0.0015
0.0051
x
x
ω
7
2.15041853722
8.53361825153
26.810846231
84.7278228489
e
ω
7 (%)
0.0001101034
0.0003757217
0.0004019763
0.0004047244


## 3.1.1. Approximate Solution forp=1

Taking into account these considerations, first we choose the approximate solution (8) of the form
(24)x~1(t)=a~0+a~1·cos⁡(ω~1·t)+b~1·sin(ω~1·t).In Step 1, the expression (13) becomes
(25)ℜ(t,a~0,a~1,b1,ω~1)=a~1ω~12cos⁡(ω~1·t)-b~1ω~12sin(ω~1·t)+54(a~1cos⁡(ω~1·t)+b~1sin(ω~1·t)+a~0)3+a~1cos⁡(ω~1·t)+b~1sin(ω~1·t)+a~0.Taking into account the initial conditionsx~(0)=A, x~(1)(t)=0 we obtain the relations
(26)a~0=-a~1+2,b~1=0.Replacing these values, the corresponding functional (14) from Step 2 is
(27)J(a~1,ω~1)=(36864π+5775πa~16+128πa~12ω~14-37800πa~15+107800πa~14)·1128ω~1+(-172800πa~13-16×(75πa~14-240πa~13+256πa~12)ω~12+167424πa~12-98304πa~1)·1128ω~1.In Step 3 we must compute the minimum ofJ with respect to a~1 and ω~1. For relatively simple problems such as this it is possible to compute directly the critical points of J and subsequently select the value corresponding to the minimum.In general, the critical points corresponding to the functional (14) J(a~1,…,a~p,b~2,…,b~p,ω~) are the solution of the system:
(28)∂J∂a~1=0,⋮∂J∂a~p=0,∂J∂b~2=0,⋮∂J∂b~p=0,∂J∂ω~1=0.Forp=1 the system becomes
(29)-(49152π-17325πa~15-128πa~1ω~4+94500πa~14-215600πa~13)·164ω~1+(32(75πa~13-180πa~12+128πa~1)ω~2+259200πa~12-167424πa~1)·164ω~1=0,(30)(16πa~12w3-(75πa~14-240πa~13+256πa~12)ω~)·14ω~1-(36864π+5775πa~16+128πa~12ω~4-37800πa~15+107800πa~14-172800πa~13)·1128ω~12+(-16(75πa~14-240πa~13+256πa~12)ω~2+167424πa~12-98304πa~1)·1128ω~12=0.This system can be solved directly. We used the “solve” command in SAGE and, after we excluded the complex solutions, we found the critical points(31)a~1=1.95861253449,ω~1=2.16432748538,a~1=1.95861253449,ω~1=-2.16432748538,a~1=0.181528334876,ω~1=-6.92142088267,a~1=0.181528334876,ω~1=6.92142088267.In order to find the minimum, we use the second partial derivative test, which is easy enough to implement in SAGE, and find thata~1=a10=1.95861253449, ω~1=ω10=2.16432748538 is the local minima. Using again the relations a00=-a10+2,b10=0, we replace all the values in the expression of the approximate solution (8) and we obtain our first approximation:
(32)x~1(t)=1.95861253449cos⁡(2.16432748538t)+0.0413874655104.In Figure1 we present the comparison between our x~1 approximate solution (solid line) and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 1
Comparison betweenx~1 (solid line) and numerical solution (dotted line).As we already saw, the approximate value of the frequency obtained here isω~1=ω10=2.16432748538 and it is already close to the exact value which in this case is ωex=2.150416169536.We observe that, generally speaking, if the value ofp in the expression of x~(t) is larger than 1, then the solution of the corresponding system of equations which gives the critical points (28) cannot be found directly. In the particular case of SAGE, the command “solve” fails to find the solutions, exiting with some kind of error message.In this situation it is still possible to find good approximations of the solutions of the problem solving the system (28) by means of a numerical method. More precisely, we can find approximate solutions for the given problem (23) solving (28) by means of a SAGE implementation of the well-known Newton method.


## 3.1.2.  Approximate Solution forp>1

As the following results will show, the Newton method is able to find approximate solutions of (28) which can lead to highly accurate approximate solutions of the problem (23).Forp=2 the approximate solution (8) has the form
(33)x~2(t)=a~0+a~1·cos⁡(ω~2·t)+b~1·sin(ω~2·t)+a~2·cos⁡(2·ω~2·t)+b~2·sin(2·ω~2·t).After we compute the corresponding expressions ofℜ(t,a~0,a~1,a~2,b~1,b~2,ω~2) and J(a~1,a~2,b~2,ω~2) and of the system (28) (all too large to insert here), we apply Newton’s method taking as the starting point of the iteration (a~1,s,a~2,s,b~2,s,ω~2,s), where a~1,s and ω~2,s are the values of a10 and ω10 computed for the previous approximation p=1, namely, a~1,s=1.95861253449 and ω~2,s=2.16432748538. In order for the sequence of approximations given by Newton’s method to converge to the solution(s) of the system (28), a~2,s and b~2,s will take successively values on a given grid of the type G=Ia~2×Ib~2, where Ii is a division of a symmetric interval centered in zero.For all the test problems included in this paper and for all the values ofp tested, a grid of the form G={-1,0.9,-0.8,…,1}×{-1,0.9,-0.8,…,1} (i.e., from −1 to 1 with step size 0.1) is large enough in the sense that if the starting point (s0,s1) scans G we can obtain using Newton's method the desired solutions of (28).In the particular case of the problem (23) it was actually sufficient to choose a~2,s=0 and b~2,s=0 and the sequence converged to the minimum of J.Forp=3 we repeat the same procedure: we compute the corresponding expressions of ℜ(t,a~0,a~1,a~2,a~3,b~1,b~2,b~3,ω~3) and J(a~1,a~2,a~3,b~2,b~3,ω~3) and of the system (28) and we apply Newton's method taking as the starting point of the iteration (a~1,s,a~2,s,a~3,s,b~2,s,b~3,s,ω~3,s), where (a~1,s,a~2,s,b~3,s,ω~3,s) are the minimum values (a10,a20,b20,ω20) computed for the previous approximation p=2 and a3,s=b3,s=0.The process can be carried on for increased values ofp until the desired accuracy is reached.While at a first look the computations may seem long and tedious, by using a software program such as SAGE, we were actually able to perform them easy and quick, obtaining a very good accuracy. Thus, forp=7, we obtained for the approximate frequency the value ω70=2.15041853722 and the corresponding approximate solution (8) has the form
(34)x~7(t)=1.0618674103207715·10-08+1.930547907025286·cos⁡(2.150418537217718·t)+1.5538232860698776·10-08·cos⁡(4.300837074435436·t)+0.0671185574473875·cos⁡(6.4512556116531545·t)+1.4239017145731592·10-09·cos⁡(8.601674148870872·t)+0.0022576839638690305·cos⁡(10.75209268608859·t)+1.3333451751464932·10-10·cos⁡(12.902511223306309·t)+7.582384931426247·10-05·cos⁡(15.052929760524027·t).In Figure2 we present the comparison between the x~7 approximate solution (solid line) and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 2
Comparison betweenx~7 (solid line) and numerical solution (dotted line).In [19] approximate solutions for the problem (23) were computed for the cases ϵ=1.25,A=2 (studied above) and ϵ=250, A=2.For the caseϵ=1.25, A=2, the approximate solution from [19] is
(35)xVI⁡(t)=1.989351555·cos⁡(2.15031·t)+0.008334403·cos⁡(6.45093·t)+0.002206664·cos⁡(10.75155·t)+0.000053689·cos⁡(15.05217·t)+(2.6·10-7)·cos⁡(19.35279·t).Table1 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xVI⁡ (see, [19]) and x~7 for the case ϵ=1.25, A=2.Table 1
Comparison of the absolute errors of the approximate solutions for problem (23) in the case ϵ=1.25,A=2.
x
VI
⁡
x
~
7
t
=
0
5.34289999998
·
1
0
-
5
0.00000000000
t
=
0.4
8.81935946956
·
1
0
-
2
3.92745813649
·
1
0
-
7
t
=
0.8
3.40733004552
·
1
0
-
2
1.66866686269
·
1
0
-
6
t
=
1.2
5.63129918651
·
1
0
-
2
5.23059939095
·
1
0
-
6
t
=
1.6
1.96411171977
·
1
0
-
2
1.89375331883
·
1
0
-
6
t
=
2
7.93642203903
·
1
0
-
2
6.81658032309
·
1
0
-
6
t
=
2.4
8.23170554343
·
1
0
-
2
1.36265840209
·
1
0
-
5
t
=
2.8
1.50931697081
·
1
0
-
2
9.46274851943
·
1
0
-
6For the caseϵ=250, A=2, the approximate solution from [19] is
(36)xVI⁡(t)=1.977507645·cos⁡(26.8026·t)+0.018844215·cos⁡(80.4078·t)+0.003619986·cos⁡(134.013·t)+0.000027438·cos⁡(187.6182·t)+0.000164633·cos⁡(241.2234·t).For this case, the approximate solution computed using our method is(37)x~7(t)=3.22758664474·10-08+1.91015330733·cos⁡(26.810846231·t)+(6.9616769357·10-08)·cos⁡(53.6216924619·t)+0.0859766881274·cos⁡(80.4325386929·t)+(7.78739552572·10-09)·cos⁡(107.243384924·t)+0.00371019879647·cos⁡(134.054231155·t)+(7.82719602479·10-10)·cos⁡(160.865077386·t)+0.0001596952849·cos⁡(187.675923617·t).Table2 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xVI⁡ (see, [19]) and x~7 for the case ϵ=250, A=2.Table 2
Comparison of the absolute errors of the approximate solutions for problem (23) in the case ϵ=250,A=2.
x
VI
⁡
x
~
7
t
=
0
1.63917000000
·
1
0
-
4
0.00000000000
t
=
0.04
9.94430509655
·
1
0
-
2
1.59432943392
·
1
0
-
6
t
=
0.08
1.01602511979
·
1
0
-
1
2.50754265524
·
1
0
-
5
t
=
0.12
1.84880609089
·
1
0
-
3
1.44364392085
·
1
0
-
6
t
=
0.16
9.38684852713
·
1
0
-
2
2.52297827337
·
1
0
-
5
t
=
0.2
1.00379880594
·
1
0
-
1
4.54516346458
·
1
0
-
5Finally, in Table3, we compare the approximate values for the frequency ω70 computed using our method with approximate values computed in [19] (ωVI⁡) and [12] (ωRHB). The comparison is made by means of the percentage error, which for a given approximate error ωapprox is defined as eωapprox=100·(|ωexact-ωapprox|/ωexact), where ωexact is the corresponding exact error. It is easy to see that our approximations are far better than the ones previously computed and they remain accurate even for the case of a very strong nonlinearity.Table 3
Comparison of approximate errors for the frequencyω for problem (23).
ϵ
·
A
2
5
100
1000
10000
ω
exact
2.1504161695
8.5335861889
26.8107384581
84.7274799361
ω
VI
⁡
2.15031
x
26.8026
84.7015
e
ω
VI
⁡ (%)
0.004937
x
0.0303552
0.0306629
ω
RHB
2.15045
8.53402
x
x
e
ω
RHB (%)
0.0015
0.0051
x
x
ω
7
2.15041853722
8.53361825153
26.810846231
84.7278228489
e
ω
7 (%)
0.0001101034
0.0003757217
0.0004019763
0.0004047244


## 3.2. Application 2: The Duffing Equation Involving Integral Forcing Terms

Our second test problem is(38)x(2)(t)+x(t)·x(1)(t)+∫0t‍tsx2(s)ds=f(t),hhhhhhhhhhhhhhhhhhhhhhhhhhh0<t<1,x(0)-x(1)(0)=0,x(1)+x(1)(1)=0,
where f(t)=-3t-3t2+5t3/2+2t4/3-t5/4-2t6/5+t7/6.The problem (38) (see, [14, 15]) is a version of the well-known Duffing equation involving both integral and nonintegral forcing terms with separated boundary conditions. This equation has been studied in a series of recent papers including [14, 15].In [14], the authors applied a generalized quasilinearization technique to prove the existence and uniqueness of the solution of Duffing equation involving both integral and nonintegral forcing terms. They showed that there are sequences of approximate solutions converging monotonically and quadratically to the unique solution of the problem.In [15], the authors gave a representation of exact solution and approximate solution of Duffing equation involving both integral and nonintegral forcing terms in the reproducing kernel space (RKS). They represented the exact solution in the form of a series and they showed that the n-term approximation of the exact solution converges to the exact solution.Next we present our results for (38) using FLSM. Also, we will compare these results with those obtained in [15].Thus, forp=3, we obtained the fact that the approximate periodic solution (8) has the form
(39)x~3(t)=1.45453sin(t)-0.26409sin(2t)+0.0245502sin(3t)+2.6625cos⁡(t)-0.169569cos⁡(2t)+0.00174065cos⁡(3t)-1.49467.Since in [15] only the numerical results are presented while the expression of approximate solution expression is not, we cannot perform a direct graphical comparison of our approximate solution with the corresponding solution from [15].Therefore, in Table4, we present the comparison of several values of the absolute errors (computed as the difference in absolute value between the exact solution and the approximate solution) corresponding to the approximate solutions xRKS from [15] and to our approximate solutions for t=0.1,0.2,…,1, as given in [15].Table 4
Comparison of the absolute errors of the approximate solutions for problem (38).
x
RKS
x
~
3
t
=
0
0
0
t
=
0.1
9.46600
·
1
0
-
6
9.85192
·
1
0
-
8
t
=
0.2
9.77708
·
1
0
-
6
9.2974
·
1
0
-
8
t
=
0.3
9.81938
·
1
0
-
6
1.58651
·
1
0
-
7
t
=
0.4
9.66575
·
1
0
-
6
5.76015
·
1
0
-
8
t
=
0.5
9.37204
·
1
0
-
6
2.07316
·
1
0
-
7
t
=
0.6
8.97794
·
1
0
-
6
5.9859
·
1
0
-
8
t
=
0.7
8.50917
·
1
0
-
6
1.55722
·
1
0
-
7
t
=
0.8
7.97984
·
1
0
-
6
9.17107
·
1
0
-
8
t
=
0.9
7.39469
·
1
0
-
6
9.65862
·
1
0
-
8
t
=
1
6.75105
·
1
0
-
6
7.84328
·
1
0
-
10


## 3.3. Application 3: The Jerk Equation Containing Velocity-Cubed and Velocity Times Displacement-Squared

Our last test is a Jerk nonlinear equation, which describes several physical problems using mechanical oscillations of the third order. The most general form of the Jerk nonlinear equations, which contains the third temporal derivative of displacement, is(40)x(3)+α·x·x(1)·x(2)+β·x(1)·x(2)2+δ·x2·x(1)+ɛ·x(1)3+γ·x(1)=0,
where the parameters α,β,δ,ɛ, and γ are constants.Nonlinear Jerk equations (40) are intensely studied by several authors in the literature and some recent results are presented in [3, 8, 9, 27].We consider the following Jerk equation:(41)x(3)+x2·x(1)+x(1)3=0,x(0)=0,x(1)(0)=B,x(2)(0)=0.Recently, Ma et al. in [3], using the Homotopy Perturbation Method, obtained high-order analytic approximate periods and periodic solutions of the Jerk Equation (41). In [8, 9], Gottlieb used the lowest-order harmonic balance method to determine analytical approximations to the periodic solution of the Jerk equations. Also, Leung and Guo in [27] obtained approximations for the angular frequency and the limit cycle for (41) based on the residue harmonic balance approach.ForB=0.5, the approximate periodic solution xHPM from [3] is
(42)xHPM(t)=0.8708027·sin(0.61924·t)-0.0106486·sin(1.85772·t).Applying FLSM we computed an approximate periodic solution of the problem (41) of the same order (containing terms up to sin(3·ω·t)):
(43)x~3(t)=0.847907380098·sin(0.615691681067·t)-(3.59275929188×10-05)·sin(1.23138336213·t)-0.0119135817912·sin(1.8470750432·t).
The approximate frequency and period are ω30=0.615691681067, T30=10.2050839736674 with an error of eω30=-0.0555984645279972.In Figure3, we can visualize and compare our approximate solution (solid line), the approximate solution xHPM from [3] (dashed line), and the numerical solution obtained by using a fourth order Runge-Kutta method (dotted line).Figure 3
Comparison between our solutionx~3(solid line), solution xHPM in [3] (dashed line), and numerical solution (dotted line).Table5 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solution xHPM from [3] and our approximate solution x~3 for the case B=0.5.Table 5
Comparison of the absolute errors corresponding to the solutionxHPM from [3] and to our solution x~3.
x
HPM
x
~
3
t
=
0
0
0
t
=
1.0
1.72354316028
·
1
0
-
2
2.11130244099
·
1
0
-
4
t
=
2.0
2.3307287671
·
1
0
-
2
1.31078376172
·
1
0
-
4
t
=
3.0
1.79768181549
·
1
0
-
2
3.46325784747
·
1
0
-
4
t
=
4.0
4.72752471124
·
1
0
-
3
8.73925906836
·
1
0
-
4
t
=
5.0
1.46107699703
·
1
0
-
2
1.42164922819
·
1
0
-
3
t
=
6.0
3.05145610435
·
1
0
-
2
1.51326363254
·
1
0
-
3
t
=
7.0
3.04290934289
·
1
0
-
2
8.83072635498
·
1
0
-
4
t
=
8.0
1.45339383549
·
1
0
-
2
7.00497125303
·
1
0
-
4
t
=
9.0
6.86925713324
·
1
0
-
3
2.00167777983
·
1
0
-
3
t
=
10.0
2.91276307352
·
1
0
-
2
2.84571545126
·
1
0
-
3Also, forB=0.5, the approximate periodic solution xRHB from [27] is given by
(44)xRHB(t)=0.8473604·sin(0.6151734·t)-0.01168397·sin(1.8455202·t)+0.000094103·sin(3.075867·t).Applying FLSM we computed an approximate periodic solution of the problem (41) of the same order (containing terms up to sin(5·ω·t)):
(45)x~5(t)=0.847966627859·sin(0.61535099365·t)+(2.70675259884·10-08)·sin⁡(1.2307019873·t)-0.0119562711384·sin⁡(1.84605298095·t)-(9.51418832279·10-11)·sin(2.4614039746·t)+(8.93050847544·10-05)·sin(3.07675496825·t).The approximate frequency and period areω50=0.61535099365, T50=10.2107339908715, with an error of -0.000264516312410058.Table6 presents the comparison of the absolute errors (computed as the difference in absolute value between the approximate solution and the corresponding numerical solution given by the Runge-Kutta method) corresponding to the approximate solutions xRHB from [27] and our approximate solutions x~5 for the case B=0.5. In this case we omitted the graphical representation of the approximate solutions xRHB and x~5 since they are both very close to the numerical solution.Table 6
Comparison of the absolute errors corresponding to the solutionxRHB from [27] and to our solution x~5.
x
RHB
x
~
5
t
=
0
0
0
t
=
1.0
2.07336138332
·
1
0
-
4
4.76705901847
·
1
0
-
6
t
=
2.0
8.18802985278
·
1
0
-
4
7.01416299886
·
1
0
-
6
t
=
3.0
6.29650984827
·
1
0
-
4
8.81003685005
·
1
0
-
7
t
=
4.0
3.27159888572
·
1
0
-
4
1.19800236493
·
1
0
-
5
t
=
5.0
7.15072428261
·
1
0
-
4
2.0505386494
·
1
0
-
5
t
=
6.0
7.95019646471
·
1
0
-
4
2.11112229619
·
1
0
-
5
t
=
7.0
1.09803461819
·
1
0
-
3
1.31470477691
·
1
0
-
5
t
=
8.0
5.35184978947
·
1
0
-
4
5.37333608552
·
1
0
-
6
t
=
9.0
8.0197576264
·
1
0
-
4
2.75784654626
·
1
0
-
5
t
=
10.0
1.42381620457
·
1
0
-
3
3.9140208578
·
1
0
-
5It is easy to see from the Tables5 and 6 that the solutions obtained by using FLSM are more accurate than the ones computed by using other methods. This fact is emphasized by Table 7 which presents a comparison of approximate periods computed in several papers, as presented in [27]. In this table, Tk denotes an approximate period obtained by means of an approximate solution containing terms up to sin(k·ω·t).Table 7
Comparison of approximate periods.
B
=
0.5
Period
Error
T
exact
10.210761
0
FLSMT1
10.1976548730449
-
0.128356025129356
FLSMT3
10.2050839736674
-
0.0555984645279972
FLSMT5
10.2107339908715
-
0.000264516312410058
Wu et al. [41] T2
10.1884
-
0.220
Wu et al. [41] T3
10.2107
-
0.001
Leung and Guo [27] T2
10.213682
0.0286
Feng and Chen [23] T2
10.23452
0.232
Feng and Chen [23] T3
10.1935
-
0.169
Feng and Chen [23] T4
10.2029
-
0.077
Hu [42] T2
10.247682
0.362
Hu [42] T3
10.226464
0.154
Hu et al. [43] T2
10.244058
0.326
Ramos [44] T2
10.80674
5.84
Ramos [44] T3
10.55034
3.33
Ramos [44] T2′
10.24328
0.318
Ramos [44] T2′′
10.26079
0.49
Ma et al. [3] T2
10.146603
-
2.27


## 4. Conclusions

In the present paper the Fourier-least squares method (FLSM) is introduced as a straightforward and efficient method to compute approximate periodic solutions for a very general class of nonlinear differential equations modeling oscillatory phenomena. Since (1) is a very general one, being able to model a large class of oscillatory phenomena, FLSM can be considered a powerful and useful method.The test problems include the cubic Duffing oscillator, the Van der Pol oscillator, and the Jerk equation. The computation of approximate solutions by FLSM clearly illustrates the accuracy of the method by comparison with approximate solutions previously computed by using other methods.



---

*Source: 513473-2014-04-23.xml*