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.
1. Introduction
A Duffing oscillator equation [1–10], which is an ordinary differential nonhomogeneous equation of the second order, can be expressed as where the dampening coefficient , coefficients and in the restoring force, amplitude , and frequency in the external oscillator force are constant. In [11], the constant coefficients , , and 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 [12–15]. The nonlinear equation shown in (1) can be rewritten as a nonlinear dynamic system as follows:
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 , which is the initial condition and period 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. Section 2 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 period T. Changing the variable , where , transforms (2) into or
Denote , where . Thus, the solution for period is equivalent to of the unit period, and (4) is transformed into
An initial value , where 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 depends on , and let for some functions , . According to the assumption that solutions are periodic, ; therefore, (6) becomes Define where denotes the transpose of a vector. Expressed in vector form, (7) can be rewritten as and (9) can be solved using Newton’s iterative method by applying an initial estimate of . The iterative process is executable if the Jacobian matrix is nonsingular, vector can be evaluated, and initial estimate is close to the solution. According to (8) and (6), the Jacobian matrix is where is an identity matrix of dimension . Let
Subsequently, and the preceding equations can be written in matrix form as The matrix is calculated by solving (14) with the initial conditions . The vector can be evaluated using Newton’s iterative process, the RungeKuttaFehlberg (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 . Accordingly, (8) is modified as which have unknowns, , and . But, there are only equations , . Kubíček proposed fixing a direction (e.g., ). The derivation is similar for the prescribed periods and is required to calculate the term Let and This can be abbreviated as With , for . The RKF method is used again to obtain . Therefore,
Denote . When , and is initial estimated appropriately, can be obtained using the RKF. The expression is repeated until .
There exist two numbers and in [] such that is monotonic at each interval [] and [], where the intervals [] and [] are the largest. Define and as the minimum and maximum of the union . The integrals of with the initial value are determined by applying the RKF method from to . Furthermore, where , and can be set as 1/2, which is acceptable if ; otherwise, should be modified as follows.
The interval [] is divided into subintervals, where is monotonic and the minimum and maximum occur at the boundaries. The turning positions and are the minimum and maximum on the monotonic part of with the maximal difference (e.g., []); that is, Therefore, Denote and calculate where , and , yielding
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 () is to fix the specific direction during the whole computation. The corresponding algorithm is given in Algorithm 1.

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 the th direction remains fixed throughout the process. To investigate this possibility, we propose an iterative approach that involves alternating the fixed th direction from to . However, for the fixed th direction, Newton’s method is applied to obtain a convergent solution . 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 Algorithm 2.

2.3. The Postverification Method
In simulations, the computational cost of evaluating and is high. We eliminated the integration and replaced the f value of the th direction by introducing a small random perturbation to . 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, , and should be identical; (ii) the maximal (minimal) values at [] and [] should be equal (similar); and (iii) the time interval between the two maximal (minimal) values should be equal (similar) to . The algorithm for the PVM is detailed as in Algorithm 3.

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 , , , and [20], are identical. For the initial estimate , where , we obtained the solution at the convergence tolerance . 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.
Example 2. The values of the coefficients , , , 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 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 and MK converge to the initial values 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 (leftbottom 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 where the convergence tolerance .
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 . Figure 3 shows the result for the initial value obtained using the PVM with the initial estimate . 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 , 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.
Example 4. This example illustrates an autonomous Duffingtype oscillator [8, 17] expressed as 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 , , and . The system described in (28) then becomes with for each . The coefficients are set as and . The obtained initial conditions for a periodic solution are
Figure 4 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.
4. Applications to PeriodicDoubling Bifurcation Problems
Kubíček and Marek [19] studied the system shown in (2) and replaced the source term depending on the value of parameter . They set the eigenvalue of the monodromy matrix at −1 to identify the periodicdoubling bifurcation points and to extend the solution path with amplitudes from the periodicdoubling bifurcation points. There exist solutions of period where the function has a period equal to [21]. We concentrated on identifying the periodicdoubling bifurcation points and calculating the solution path for a prescribed initial value by using the quasiarclength continuation method proposed in [8]. The other branch solution path of period has been obtained based on the periodicdoubling bifurcation points.
The eigenvalues of the monodromy matrix are called Floquet multipliers, and the numbers are called Floquet exponents. According to Floquet theory, a periodic solution is stable when for each (except ). 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 periodicdoubling bifurcation point exists where the two branches intersect at on the unit circle. Another method for determining the periodicdoubling bifurcation points is detailed as follows. The monodromy matrix has an eigenvalue of −1 at the periodicdoubling bifurcation point; therefore, there exists a nonzero vector such that is fulfilled. Because , where , these parameters satisfy (19) and
The aforementioned system of variables has only equations. Furthermore, each nonzero multiplier of vector is a solution of (19). A normalization approach, enables fixing the component of throughout the computation. Without loss of generality, we set . Let be the unique solution of the system shown in (19) and (20), which satisfy where . The solution of (where ) is the periodicdoubling bifurcation point. Finally, the solution can be determined by applying Newton’s iterative method.
Example 5. The Duffing equation [20, 22] can be expressed as where . Equation (36) then becomes 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 quasiarclength continuation method [8] can be applied to calculate the two curves that intersect at point P. Figure 5 shows the periodicdoubling bifurcation point P .
The remaining of this section is to demonstrate the difficulty on calculation of the unstable periodic solutions. In the left panel of Figure 5, 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 denote the initial values of the solution of 2π, where the initial value A indicates the stable solution and the initial values B represents the unstable solution. Figure 6 depicts the phase plane and the unstable periodic solution requires a further study.
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 with variables and 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to anonymous referees for their valuable comments. ChienChang Yen was supported in part by the Ministry of Science and Technology, Taiwan, under the Grant NSC 1022115M030003 and the other FJU Project of A0502004.