# A Postverification Method for Solving Forced Duffing Oscillator Problems without Prescribed Periods

Authors: Hong-Yen Lin; Chien-Chang Yen; Kuo-Ching Jen; Kang C. Jea
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper proposes a postverification method (PVM) for solving forced Duffing oscillator problems without prescribed periods. Comprising a postverification procedure and small random perturbation, the proposed PVM improves the sensitivity of the convergence of Newton’s iteration. Numerical simulations revealed that the PVM is more accurate and robust than Kubíček’s approach. We applied the PVM to previous research on bifurcation problems.


---

## Body

## 1. Introduction

A Duffing oscillator equation [1–10], which is an ordinary differential nonhomogeneous equation of the second order, can be expressed as
(1)u′′+au′+bu+cu3=αcos⁡ωt,
where the dampening coefficient a, coefficients b and c in the restoring force, amplitude α, and frequency ω in the external oscillator force are constant. In [11], the constant coefficients a, b, and c were generalized to variable coefficients, and the existence of positive solutions was shown. However, researchers are interested in periodic and bifurcation solutions with the prescribed period T=2π/ω [12–15]. The nonlinear equation shown in (1) can be rewritten as a nonlinear dynamic system as follows:
(2)dxidt=fi(x1,x2,…,xn,t),i=1,2,…,n,t≥0.There are several applications for nonlinear dynamic systems with ambiguous or unknown force oscillators [16–18]. However, further research is necessary to solve (1) where the frequency of the external oscillator force is unknown. Specifically, it is necessary to identify a state where I0=(x1(0),x2(0),…,xn(0),T), which is the initial condition and period T for (2).A continuation method involving Newton’s iterative method has frequently been applied to solve the branch of bifurcation points. The approach proposed by Kubíček and Marek [19] for solving the aforementioned problem does not guarantee that the obtained numerical solution is always periodic, potentially because fixing the specific direction or approach is sensitive to the convergence of Newton’s iteration. First, we modified Kubíček’s approach by applying an alternative method to eliminate the drawback of fixing a specific direction; hence, we term the proposed approach the modified Kubíček (MK) approach. A postverification process, called the postverification method (PVM), is then introduced into the MK approach to ensure periodicity. The PVM, which is based on the stability of solutions, has the following three characteristics: (i) an alternative direction for the initial state is fixed; (ii) a randomly small perturbation approach is used to confirm the stability of the solution; and (iii) a postverification test of periodicity is applied.Forced Duffing oscillators, which lack nonlinear restoring force, were simulated using Kubíček’s approach, the MK approach, and proposed PVM to evaluate the performances of these methods. The numerical results indicated that the proposed approaches can be used to obtain periodic numerical solutions with greater stability. Furthermore, the PVM was applied to autonomous equations to obtain the periodic solutions.The remainder of this paper is organized as follows. Section2 introduces the proposed method, and the numerical results are reported in Section 3. Section 4 describes the application of the obtained results to a bifurcation calculation. Finally, Section 5 provides our conclusion.


## 2. The Postverification Method

Based on the assumption that the solution is stable, we developed the MK approach [19] and PVM. The proposed method is feasible because of the stable criterion and postverification process. The proposed method involves starting from the system shown in (2) by using a prescribed periodT. Changing the variable t=Ts, where 0≤s≤1, transforms (2) into
(3)dxidt=dxiTds=fi(x1,x2,…,xn,Ts),i=1,2,…,n
or
(4)dxids=Tfi(x1,x2,…,xn,Ts),i=1,2,…,n.Denoteyi(s)=xi(Ts), where 0≤s≤1. Thus, the solution xi for period T is equivalent to yi of the unit period, and (4) is transformed into
(5)dyids=Tfi(y1,y2,…,yn,Ts),i=1,2,…,n.An initial valueyi(0)=ηi, where i=1,2,…,n such that the solutions that are both periodic and have bifurcation points, can be identified because a periodic solution of (5) exists. Moreover, the value yi(1) depends on η1,η2,…,ηn, and let
(6)yi(1)≡ϕi(η1,η2,…,ηn)
for some functions ϕi, i=1,…,n. According to the assumption that solutions are periodic, yi(0)=yi(1); therefore, (6) becomes
(7)ηi-ϕi(η1,η2,…,ηn)=0,i=1,2,…,n.
Define
(8)Fi=ηi-ϕi(η1,η2,…,ηn),i=1,2,…,n,η=(η1,η2,…,ηn)t,y=(y1,y2,…,yn)t,ϕ=(ϕ1,ϕ2,…,ϕn)t,F=(F1,F2,…,Fn)t,
where (·)t denotes the transpose of a vector. Expressed in vector form, (7) can be rewritten as
(9)F(η)=0,
and (9) can be solved using Newton’s iterative method
(10)η(l+1)=η(l)-[∂F(η(l))∂η]-1F(η(l)),l=0,1,2,…
by applying an initial estimate of η(0). The iterative process is executable if the Jacobian matrix Fη(η(l))=∂F(η(l))/∂η is nonsingular, vector F(η(l)) can be evaluated, and initial estimate η(0) is close to the solution. According to (8) and (6), the Jacobian matrix is
(11)Fη=In×n-∂ϕ∂ηwith∂ϕi∂ηj=∂yi(1)∂ηj,
where In×n is an identity matrix of dimension n×n. Let
(12)Pij(s)=∂yi(s)∂ηj,0≤s≤1,i,j=1,2,…,n.Subsequently,(13)ddsPij(s)=dds∂yi(s)∂ηj=∂∂ηj(dyi(s)ds)=∂∂ηj(Tfi(y1,y2,…,yn,Ts))=T∑m=1n∂fi∂ym∂ym∂ηj,
and the preceding n×n equations can be written in matrix form as
(14)dPds=T∂f∂yP.
The matrix ∂ϕ/∂η=P(1) is calculated by solving (14) with the initial conditions P(0)=In×n. The vector ϕ=y(1) can be evaluated using Newton’s iterative process, the Runge-Kutta-Fehlberg (RKF) method, and (5) and (6). We introduce both Kubíček’s method and the proposed MK and PVM to solve (2) without the prescribed period T. Accordingly, (8) is modified as
(15)Fi(η1,η2,…,ηn,T)=ηi-ϕi(η1,η2,…,ηn,T),iiiiiiiiiiiiiiiiiii=1,2,…,n,
which have n+1 unknowns, η1,η2,…,ηn, and T. But, there are only n equations Fi=0, i=1,2,…,n. Kubíček proposed fixing a direction (e.g., k). The derivation is similar for the prescribed periods and is required to calculate the term
(16)∂Fi∂T=-∂ϕi∂T=-∂yi(1)∂T.
Let ri(s)=∂yi(s)/∂T and
(17)ddsri(s)=dds∂yi∂T=∂∂Tdyi(s)ds=∂∂T(Tfi(y1,y2,…,yn,Ts))=fi(y1,y2,…,yn,Ts)+T(∂fi(y1,y2,…,yn,Ts)∂T+∑m=1n∂fi∂ym∂ym∂T).
This can be abbreviated as
(18)r′=f+T∂f∂T+T∂f∂yr.
With ri(0)=∂yi(0)/∂T=∂ηi/∂T=0, for i=1,2,…,n. The RKF method is used again to obtain ri(1)=∂Fi/∂T=∂ϕi/∂T. Therefore,
(19)∂ϕi∂ηj={Pij(1),i≤j≤k-1,Pij+1(1),k≤j≤n-1,ri(1),j=n.DenoteU=(y1,…,yk-1,yk+1,…,yn,T). When Fi(U)=Fi(y1,y2,…,yn,T)=ηi-ϕi(η1,η2,…,ηn)=ηi-yi(1),
(20)dyids=Tfi(y1,y2,…,yn,Ts),
and yi(0)=ηi is initial estimated appropriately, y(1) can be obtained using the RKF. The expression
(21)U(l+1)=U(l)-(FU(U(l)))-1(y(1)-η(l))
is repeated until ∥U(l+1)-U(l)∥<ε.There exist two numberst1 and t2 in [0,T] such that yk(t) is monotonic at each interval [0,t1] and [t2,T], where the intervals [0,t1] and [t2,T] are the largest. Define yk- and yk+ as the minimum and maximum of the union [0,t1]∪[t2,T]. The integrals of
(22)y′=fi(y1,y2,…,yn,T),fori=1,2,…,n.
with the initial value y(0)=(η1,η2,…,ηn) are determined by applying the RKF method from t=0 to T. Furthermore,
(23)ηk-=12((1+ω)yk-+(1-ω)yk+),ηk+=12((1-ω)yk-+(1+ω)yk+),
where 0<ω<1, and ω can be set as 1/2, which is acceptable if ηk-<ηk<ηk+; otherwise, ηk should be modified as follows.The interval [0,1] is divided into subintervals, where yk is monotonic and the minimum and maximum occur at the boundaries. The turning positions yk-- and yk++ are the minimum and maximum on the monotonic part of yk(t) with the maximal difference (e.g., [ti,ti+1]); that is,
(24)yk--=min⁡t∈[ti,ti+1]yk(t),yk++=max⁡t∈[ti,ti+1]yk(t).
Therefore,
(25)ηk=12(yk++-yk--).
Denote ηn+1=T and calculate
(26)FU=[∂Fi∂ηj]=[∂ϕi∂ηj]-[∂ηi∂ηj],
where i=1,2,…,n, and j=1,2,…,k-1,k+1,…,n+1, yielding
(27)∂ηi∂ηj={δij,1≤i,j≤k-1,0,i=k,1≤j≤n-1,1,j=i-1,k+1≤i≤n-1,0,j=n,1≤i≤n.The aforementioned argument can be summarized as the following algorithm, which is referred to as Kubíček’s approach [12, 19].


### 2.1. Kubíček’s Approach

The approach proposed by Kubíček and Marek [19] for solving (q) is to fix the specific direction during the whole computation. The corresponding algorithm is given in Algorithm 1.Algorithm 1:Kubíček method [19].
Letε be given and fix  k
Set an initial estimateU(0)=(η1,…,ηk-1,ηk+1,…,ηn,T)t.
Repeat
U(l+1)=U(l)-(FU(U(l)))-1F(U(l)),whereF=(F1,…,Fn)t.
Modifyηk by applying (25)
until∥U(l+1)-U(l)∥<ε.


### 2.2. The Modified Kubíček Approach

When using Kubíček’s approach, a periodic solution is not always obtainable, even when the criterion of convergence is satisfied, potentially because thekth direction remains fixed throughout the process. To investigate this possibility, we propose an iterative approach that involves alternating the fixed kth direction from k=1 to n. However, for the fixed kth direction, Newton’s method is applied to obtain a convergent solution (η1,…,ηk-1,ηk+1,…ηn,T). The same procedure is applied for each direction until the criterion for convergence is satisfied.The algorithm for the modified Kubíček(MK) approach is given in Algorithm2.Algorithm 2:The modified Kubíček method.
Letε be a given value.
Set the initial estimate(η1,…,ηn,T).
Repeat
Fork=1 to n // let ηk be fixed
Let U=(η1,…,ηk-1,ηk+1,…,ηn,T)t.
Repeat
δ←[FU(U)]-1F(U).
U←U-δ.
until  ∥δ∥<ε.
Modify ηk by applying (25).
Nextk
until(η1,…,ηn,T) becomes constant.


### 2.3. The Postverification Method

In simulations, the computational cost of evaluatingyk-- and yk++ is high. We eliminated the integration and replaced the f value of the kth direction by introducing a small random perturbation to ηk. In addition, we investigated the sensitivity of the convergence of Newton’s iteration. The proposed PVM was employed to overcome the sensitivity and, thus, confirm the following three conditions: (i) the values at 0, T, and 2T should be identical; (ii) the maximal (minimal) values at [0,T] and [T,2T] should be equal (similar); and (iii) the time interval between the two maximal (minimal) values should be equal (similar) to T. The algorithm for the PVM is detailed as in Algorithm 3.Algorithm 3:The postverification method.
Letε be given.
Set the initial estimate(η1,…,ηn,T).
Repeat
Fork=1 to n // let ηk be fixed.
Let U=(η1,…,ηk-1,ηk+1,…,ηn,T)t.
Repeat
δ←[FU(U)]-1F(U).
U←U-δ.
until ∥δ∥<ε.
ηk←ηk+randomnumber
Nextk
until(η1,…,ηn,T) remains constant and the postverification is successful.


## 2.1. Kubíček’s Approach

The approach proposed by Kubíček and Marek [19] for solving (q) is to fix the specific direction during the whole computation. The corresponding algorithm is given in Algorithm 1.Algorithm 1:Kubíček method [19].
Letε be given and fix  k
Set an initial estimateU(0)=(η1,…,ηk-1,ηk+1,…,ηn,T)t.
Repeat
U(l+1)=U(l)-(FU(U(l)))-1F(U(l)),whereF=(F1,…,Fn)t.
Modifyηk by applying (25)
until∥U(l+1)-U(l)∥<ε.


## 2.2. The Modified Kubíček Approach

When using Kubíček’s approach, a periodic solution is not always obtainable, even when the criterion of convergence is satisfied, potentially because thekth direction remains fixed throughout the process. To investigate this possibility, we propose an iterative approach that involves alternating the fixed kth direction from k=1 to n. However, for the fixed kth direction, Newton’s method is applied to obtain a convergent solution (η1,…,ηk-1,ηk+1,…ηn,T). The same procedure is applied for each direction until the criterion for convergence is satisfied.The algorithm for the modified Kubíček(MK) approach is given in Algorithm2.Algorithm 2:The modified Kubíček method.
Letε be a given value.
Set the initial estimate(η1,…,ηn,T).
Repeat
Fork=1 to n // let ηk be fixed
Let U=(η1,…,ηk-1,ηk+1,…,ηn,T)t.
Repeat
δ←[FU(U)]-1F(U).
U←U-δ.
until  ∥δ∥<ε.
Modify ηk by applying (25).
Nextk
until(η1,…,ηn,T) becomes constant.


## 2.3. The Postverification Method

In simulations, the computational cost of evaluatingyk-- and yk++ is high. We eliminated the integration and replaced the f value of the kth direction by introducing a small random perturbation to ηk. In addition, we investigated the sensitivity of the convergence of Newton’s iteration. The proposed PVM was employed to overcome the sensitivity and, thus, confirm the following three conditions: (i) the values at 0, T, and 2T should be identical; (ii) the maximal (minimal) values at [0,T] and [T,2T] should be equal (similar); and (iii) the time interval between the two maximal (minimal) values should be equal (similar) to T. The algorithm for the PVM is detailed as in Algorithm 3.Algorithm 3:The postverification method.
Letε be given.
Set the initial estimate(η1,…,ηn,T).
Repeat
Fork=1 to n // let ηk be fixed.
Let U=(η1,…,ηk-1,ηk+1,…,ηn,T)t.
Repeat
δ←[FU(U)]-1F(U).
U←U-δ.
until ∥δ∥<ε.
ηk←ηk+randomnumber
Nextk
until(η1,…,ηn,T) remains constant and the postverification is successful.


## 3. Numerical Results

The performance of the proposed methods was evaluated by simulating the Duffing oscillator equation with a nonlinear restoring force [20] and autonomous equations representing the case of the nonlinear restoring force of an unknown period.Example 1.
This example illustrates that the results obtained by using Kubíček’s approach, the MK approach, and proposed PVM to solve (1), where a=b=0.1, c=1, α=5.0259, and ω=1 [20], are identical. For the initial estimate (η1(0),η2(0),T(0))=(3.0,1.2,5.6), where k=2, we obtained the solution (η1,η2,T)=(2.6885114,1.2,6.28134528) at the convergence tolerance ε=1E-6. Figure 1 shows the numerical solutions and trajectories in the phase plane, which were simulated using Kubíček’s approach, the MK approach, and proposed PVM. The solutions of the two periods were obtained quickly at the maximum, and a slow variation occurred near the minimum. The trajectories of the solutions in the phase plane show that the periodicity was preserved. The results of these three methods were identical.The numerical solution (top panels) and the trajectories in the phase planes (bottom panels) simulated by Kubíček’s approach, MK, and PVM are shown from left to right panels, respectively. These results are identical for these three methods.Example 2.
The values of the coefficientsa, b, c, and ω were reset to those used in Example 1, butα was changed from 5.0259 to 5.0258827. This example illustrates that Kubíček’s approach is sensitive to the parameterα. Figure 2 shows the profiles of the numerical solutions Figure 2 (top panels) and phase plane Figure 2 (bottom panels) for the initial estimate (η1(0),η2(0),T(0))=(3.0,1.67,5.6) when simulated using Kubíček’s approach (left panels), the MK approach (middle panels), and the proposed PVM (right panels). Kubíček’s approach with fixed k=2 and MK converge to the initial values (η1,η2,T)=(2.6630419,1.4719480,2.418418) and (−2.3248997, 1.67, and 34.8045258), respectively. Unexpectedly, the solutions were not periodic because the trajectory of the solution in the phase plane (left-bottom panel) was not closed, potentially because the solutions are sensitive to the convergence of Newton’s iteration. However, the proposed PVM yielded a periodic solution with the initial value (η1,η2,T)=(2.69035153,1.1774813,6.2842755) where the convergence tolerance ε=1E-6.These results show the numerical solutions on [0,4π] (top panels) and the trajectories in the phase planes (bottom panels) simulated by Kubíček’s approach, MK, and PVM are shown from left to right panels, respectively.Example 3.
Because the period is unprescribed, the problem (1) allows alternative solutions for various periods. The parameters were set the same as they were set in Example 1, except α=1. Figure 3 shows the result for the initial value (η1,η2,T)=(2.70149522,1.00802403,12.56646993) obtained using the PVM with the initial estimate (η1(0),η2(0),T(0))=(3.0,1.0,5.6). The numerical solutions and trajectories in the phase plane are shown in the left and right panels of the figure, respectively. Compared with Example 2, the period of the solution was 12.56646993≈4π, which can be confirmed according to the trajectory in the phase plane. According to this simulation, the unprescribed period enables the period to differ from the forced oscillator. However, the periodicity of the PVM was still robust.These results show the numerical solution on [0,8π] (left) and the trajectories in the phase plane (right) for Example 3. The periodicity can be observed by the profiles of the trajectories in the phase planes.Example 4.
This example illustrates an autonomous Duffing-type oscillator [8, 17] expressed as
(28)x′′-bx′-x+x3+κz=0,z′=ω(x′-z),
where ω is the characteristic rate and κ is the feedback coefficient. This example shows that the external force is determined implicitly. Our goal was to calculate a periodic solution that differs from the motivation of [9]. Let x1=x, x2=x′, and x3=z. The system described in (28) then becomes
(29)dx1dt=x2,dx2dt=x1-x13+bx2-κx3,dx3dt=ω(x2-x3)
with xi(t+T)=xi(t) for each i. The coefficients are set as b=κ=0.5 and ω=1. The obtained initial conditions for a periodic solution are
(30)x1(0)=-1.3241577E+7,x2(0)=-6.9441420E+13,x3(0)=-1.3241578E+7,T=5.23147E-7.Figure4 shows the solution of two periods and the corresponding trajectory in three dimensions. Although the period was small and the amplitudes of the solutions were large, the solutions were stable. This stability enables identifying the periodic solutions. In this simulation, it was necessary to narrow the initial estimate to a value that was close to the solution.These results show the numerical solutionx1 (left) and the trajectories in 3D (right) for Example 4.


## 4. Applications to Periodic-Doubling Bifurcation Problems

Kubíček and Marek [19] studied the system shown in (2) and replaced the source term fi(x1,x2,…,xn,α) depending on the value of parameter α. They set the eigenvalue of the monodromy matrix
(31)M≡∂ϕi∂ηj
at −1 to identify the periodic-doubling bifurcation points and to extend the solution path with amplitudes from the periodic-doubling bifurcation points. There exist solutions yi(t) of period T where the function fi has a period equal to T [21]. We concentrated on identifying the periodic-doubling bifurcation points and calculating the solution path for a prescribed initial value by using the quasi-arclength continuation method proposed in [8]. The other branch solution path of period 2T has been obtained based on the periodic-doubling bifurcation points.The eigenvaluesρj of the monodromy matrix M are called Floquet multipliers, and the numbers γi=log⁡ρi are called Floquet exponents. According to Floquet theory, a periodic solution is stable when |ρj|<1 for each ρj (except ρj=1). The stability of the periodic solutions on one branch may change according to variation in α; one or two multipliers may be located outside the unit disk. Once a multiplier passes through the unit circle at −1 on the real axis, the periodic stable solution becomes unstable, and another branch of periodic solutions can be identified at this point. In other words, the periodic-doubling bifurcation point exists where the two branches intersect at (-1,0) on the unit circle. Another method for determining the periodic-doubling bifurcation points is detailed as follows. The monodromy matrix M has an eigenvalue of −1 at the periodic-doubling bifurcation point; therefore, there exists a nonzero vector ν=(ν1,…,νn)t≠0 such that
(32)(M+I)ν=0
is fulfilled. Because yi(0)=yi(1), where i=1,2,…,n, these parameters satisfy (19) and
(33)Fi(η1,η2,…,ηn,α)=ϕi(η1,η2,…,ηn,α)-ηi=0,iiiiiiiiiiiiiiiiiiiiiiii=1,2,…,n.The aforementioned system of2n+1 variables (η1,…,ηn,α,ν1,…,νn)t has only 2n equations. Furthermore, each nonzero multiplier of vector ν is a solution of (19). A normalization approach,
(34)ν=νk(ν1νk,ν2νk,…,νnνk),forsomeνk≠0,
enables fixing the component νk of ν=1 throughout the computation. Without loss of generality, we set k=1. Let P=(P1,…,Pn,Pn+1,…,P2n)t be the unique solution of the system shown in (19) and (20), which satisfy
(35)Pi=Fi=0,Pn+i=∑m=1n∂ϕi∂ηmνm+νi=0,
where i=1,2,…,n. The solution w of P(w)=0 (where w1=1) is the periodic-doubling bifurcation point. Finally, the solution w can be determined by applying Newton’s iterative method.Example 5.
The Duffing equation [20, 22] can be expressed as
(36)u′′+0.1u′+0.1u+u3=αcos⁡t,
where y1=u,y2=u’. Equation (36) then becomes
(37)dy1dt=y2,dy2dt=-0.1y2-0.1y1-y13+αcos⁡t,y1(0)=η1,y2(0)=η2.
We apply the PVM to solving (37) without prescribed periods, but for fixed parameter α. The obtained initial values represent the solutions of periods 2π and 4π, similarly for Examples 2 and 3. Furthermore, the quasi-arclength continuation method [8] can be applied to calculate the two curves that intersect at point P. Figure 5 shows the periodic-doubling bifurcation point P(η1,η2,α)=(2.68964196,1.18618911,5.02588278).The calculation of bifurcation points depicts the periods of2π (blue line) and 4π (red line) in the left panel and illustrates the stable (blue line) and unstable (red line) curves in the right panel.The remaining of this section is to demonstrate the difficulty on calculation of the unstable periodic solutions. In the left panel of Figure5, each point on the curve represents an initial value determined by (37). The blue curve denotes the initial values, which can be used to ascertain the solution of period 2π, and the red curve represents the initial values, which can be used to determine the solution of period 4π Figure 5 (left panel). The Floquet multipliers are used to evaluate the stability. According to the definition of stability, the initial values for the solution of periods 2π and 4π can be divided into stable (blue line) and unstable (red line) regions. In the right panel of Figure 5, the markers
(38)A(y1(0),y2(0),α)=(1.42420721847224,-0.39906340734434,iiiiiiii4.10971356401897),B(y1(0),y2(0),α)=(3.41661920410986,-0.61423445258396,iiiiiiii12.69779142606720)
denote the initial values of the solution of 2π, where the initial valueA indicates the stable solution and the initial valuesB represents the unstable solution. Figure 6 depicts the phase plane and the unstable periodic solution requires a further study.It illustrates the phase plane with the initial value pointA (stable) in the left panel and shows the phase plane with the initial value pointB (unstable) expressed in (37) in the right panel.


## 5. Conclusion

This study proposes a PVM for obtaining the periodic solutions of forced Duffing equations without a prescribed period. Specifically, to solve a system withn+1 variables and n equations, we alternatively fixed a variable to improve Kubíček's method and introduced a random small perturbation to verify the stability of the solutions. The postverification ensures that specific conditions are met whenever a solution is obtained. The proposed PVM is more robust on periodicity based on alternative convergence criterion, as demonstrated by our numerical results. We conclude that the most crucial application of the proposed PVM is determining the periodicity of solutions.The application of the bifurcation problem demonstrated the branching of other periodic solutions that can be traced to the bifurcation points. However, the simulation and autonomous case revealed that the convergent initial values must be near to the initial estimate. Thus, calculating unstable periodic solutions remains challenging.



---

*Source: 317460-2014-08-26.xml*