Next Article in Journal
Supersymmetric Partners of the One-Dimensional Infinite Square Well Hamiltonian: Special Cases
Previous Article in Journal
Impact of Buoyancy and Stagnation-Point Flow of Water Conveying Ag-MgO Hybrid Nanoparticles in a Vertical Contracting/Expanding Riga Wedge
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Periodic Orbits of Nonlinear Ordinary Differential Equations Computed by a Boundary Shape Function Method

1
Center of Excellence for Ocean Engineering, National Taiwan Ocean University, Keelung 202301, Taiwan
2
Center of Excellence for the Oceans, National Taiwan Ocean University, Keelung 202301, Taiwan
3
Department of Mechanical Engineering, National United University, Miaoli 360302, Taiwan
4
Department of Marine Engineering, National Taiwan Ocean University, Keelung 202301, Taiwan
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(7), 1313; https://doi.org/10.3390/sym14071313
Submission received: 6 May 2022 / Revised: 28 May 2022 / Accepted: 22 June 2022 / Published: 24 June 2022
(This article belongs to the Section Mathematics)

Abstract

:
In the paper, we determine the period of an n-dimensional nonlinear dynamical system by using a derived formula in an (n + 1)-dimensional augmented space. To form a periodic motion, the periodic conditions in the state space and nonlinear first-order differential equations constitute a special periodic problem within a time interval with an unknown length. Two periodic problems are considered: (a) boundary values are given and (b) boundary values are unknown. By using the shape functions, a boundary shape function method (BSFM) is devised to obtain an initial value problem with the initial values of the new variables given. The unknown terminal values of the new variables and period are determined by two iterative algorithms for the case (a) and one iterative algorithm for the case (b). The periodic solutions obtained from the BSFM satisfy the periodic conditions automatically. For the numerical example, the computed order of convergence displays the merit of the BSFM. For the sake of comparison, the iterative algorithms based on the shooting method for cases (a) and (b) were developed by directly implementing the Poincaré map into the fictitious time-integration method to determine the period. The BSFM is better than the shooting method in terms of convergence speed, accuracy, and stability.

1. Introduction

The study of nonlinear ordinary differential equations (ODEs) is highly important for many practical applications, since the most dynamical phenomena in our real world are described by nonlinear ODEs. In the last few decades, considerable attention was focused on the periodic solutions of nonlinear oscillators [1], which are ubiquitous in every area of science concerned with oscillatory phenomena, not only in the areas of mechanics and physics, but also in other disciplines involving engineering applications.
Computational methods have been developed for determining the periodic solutions of nonlinear ODEs, such as the linearized Lindstedt–Poincaré method [2], the harmonic balance method [3,4,5], the variational iteration method [6,7], the homotopy perturbation method [8,9,10], the parameter-expanding method [11], the exp-function method [12], the differential transform method [13], and the optimal scale polynomial interpolation technique [14]. Dai et al. [15,16] employed a simple time-domain collocation method to solve nonlinear oscillatory problems. Khan et al. [17] solved the periodic solutions of a nonlinear Duffing oscillator quickly by using the modified harmonic balance method with optimal scales and the optimally scaled polynomial-Fourier-series method. Viswanath [18] developed a method to compute high-precision periodic orbits of nonlinear dynamical systems based on the linearized Lindstedt–Poincaré method.
The periodic problems of nonlinear dynamical systems involve complicated situations. Not all nonlinear dynamical systems allow periodic motions; even if periodic motion is present in a nonlinear dynamical system, not all the initial points necessarily possess a periodic orbit to emanate from it. In practice, the periodic conditions are the basic requirements of the numerical method for determining the unknown period of a system, but for many nonlinear dynamical systems, the periodic conditions are highly unsuitable, since most of the initial points do not lead to periodic orbits in the state space. For most nonlinear dynamical systems, the periods are unknown, and the determination of these periods and periodic solutions is a difficult issue. Although AUTO, COCO, or MATCONT can be used to solve the periodic problems of nonlinear ODEs, new and efficient methods are still important to solve this complex periodic problem.
The idea of the boundary shape function (BSF) was first introduced by Liu and Chang [19] to find the periodic solutions of nonlinear jerk equations, and by Liu et al. [20] to solve the optimal control problems of nonlinear Duffing oscillators. Since the methods of the BSF are versatile, they have been successfully applied to solving boundary value problems. However, the current problems with unknown period and periodic boundary values in the n-dimensional space are more difficult to solve.
This paper is structured as follows. In Section 2, a novel formula is derived in an (n + 1)-dimensional augmented space for n-dimensional first-order ODEs to compute the period. For the two cases considered in Section 3, we derive two sets of shape functions. In Section 4, we introduce new variables whose initial values are given, such that the periodic problem is framed as an initial value problem, and three iterative algorithms, based on the boundary shape function method (BSFM), to account for different periodic problems. Seven numerical examples are given in Section 5 to test the proposed new iterative algorithms’ ability to solve the periodic problems of nonlinear ODEs. The periodic problem for the non-autonomous system of nonlinear ODEs is discussed in Section 6, where the periodic solutions of the forced Duffing oscillator are given. Finally, we draw some conclusions in Section 7.

2. Nonlinear First-Order ODEs and Period

We can transform any n-th order ODE into a system of n first-order ODEs in the state space. Here, we consider the periodic problem of n first-order ODEs:
x ˙ j ( t ) = f j ( x 1 , , x n ) ,       j = 1 , , n .
The extension to a non-autonomous system of n first-order ODEs is discussed in Section 6. In terms of τ = t / T , it follows from Equation (1) that
x j ( τ ) = T f j ( x 1 , , x n ) ,       j = 1 , , n .
When T is an unknown constant, there are two periodic boundary conditions to be considered:
( a )   x j ( 0 ) = x j 0 = x j ( 1 ) ,   j = 1 , , n ,
( b )   x j ( 0 ) = x j ( 1 ) ,   j = 1 , , n .
For the case (a), x j 0 ,   j = 1 , , n are given, such that there is only one unknown period T to be determined. For the case (b), x j 0 ,   j = 1 , , n are not given, such that n + 1 unknown values x j 0 ,   j = 1 , , n and T are to be determined. In (a) and (b), the common parts x j ( 0 ) = x j ( 1 ) ,   j = 1 , , n are termed periodic conditions in the state space. Equations (2) and (3) or (4) constitute periodic problems in the state space. It is apparent that the periodic problem (b) is more difficult than the periodic problem (a).
It must be emphasized that periodic solutions may not exist for some ODEs, and that there may be infinitely more periodic solutions for others, such as the Lorenz dynamical system. Therefore, the periodic problem we consider here must be present in at least one periodic motion in the state space, but it may be that not all initial points permit the periodic motion. Hence, the periodic problem with case (b) is proposed to seek such a periodic motion in the state space, with both period and initial point being unknown. Some ODEs permit periodic motions for most initial points, and the period depends on the given initial point. Hence, the periodic problem with case (a) is proposed for seeking the unknown period of such a periodic motion in the state space. Notice that the periodicity condition in Equations (3) and (4) is, in essence, the Poincaré map of the sought periodic orbit, also known as the first return map.
In order to evaluate the period T in Equation (2), we introduce a supplementary variable x n + 1 ( τ ) by considering its evolution:
x n + 1 ( τ ) = x I ( τ ) ,   x n + 1 ( 0 ) = x n + 1 0 ,   I { 1 , , n } ,
where x n + 1 0 is a given initial value and the index I is one integer among { 1 , , n } . Taking the product of x n + 1 to the I-th equation, and using Equation (5), Equation (2) with j = I and the following formula:
( x n + 1 x I ) x n + 1 x I = x n + 1 x I ,
we can obtain
( x n + 1 x I ) x I 2 = T x n + 1 f I ( x 1 , , x n ) ,
which is integrated by using Equations (3) and (4) with j = I ,
[ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) 0 1 x I 2 ( τ ) d τ = T 0 1 x n + 1 ( τ ) f I ( x 1 ( τ ) , , x n ( τ ) ) d τ .
We can choose a suitable index I { 1 , , n } and some suitable values of x n + 1 0 to cause
{ [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) 0 1 x I 2 ( τ ) d τ } { 0 1 x n + 1 ( τ ) f I ( x 1 ( τ ) , , x n ( τ ) ) d τ } > 0 ,
such that the period T can be computed from Equation (8) by
T = [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) 0 1 x I 2 ( τ ) d τ 0 1 x n + 1 ( τ ) f I ( x 1 ( τ ) , , x n ( τ ) ) d τ ,
which expresses the period T in an (n + 1)-dimensional augmented space ( x 1 , , x n ) .

3. Shape Functions

According to [19], two sets of shape functions are derived to solve the periodic problem of Equation (1) for the two cases (a) and (b) in Equations (3) and (4).
Theorem 1.
For case (a), we take two shape functions by
p 1 ( τ ) = 1 τ ,   p 2 ( τ ) = τ .
For any function y j ( τ ) C 1 [ 0 ,   1 ] ,
x j ( τ ) = y j ( τ ) p 1 ( τ ) [ y j ( 0 ) x j 0 ] p 2 ( τ ) [ y j ( 1 ) x j 0 ] ,   j = 1 , , n
Satisfy the periodic conditions in Equation (3).
Proof. 
Inserting τ = 0 into Equation (12), leads to
x j ( 0 ) = y j ( 0 ) p 1 ( 0 ) [ y j ( 0 ) x j 0 ] p 2 ( 0 ) [ y j ( 1 ) x j 0 ] ,
Upon using p 1 ( 0 ) = 1 and p 2 ( 0 ) = 0 from Equation (11), which becomes
x j ( 0 ) = y j ( 0 ) [ y j ( 0 ) x j 0 ] = x j 0 ,
Thus, we prove the first part of Equation (3). Similarly, inserting τ = 1 into Equation (12) and using p 1 ( 1 ) = 0 and p 2 ( 1 ) = 1 from Equation (11) yields
x j ( 1 ) = y j ( 1 ) p 1 ( 1 ) [ y j ( 0 ) x j 0 ] p 2 ( 1 ) [ y j ( 1 ) x j 0 ] = y j ( 1 ) [ y j ( 1 ) x j 0 ] = x j 0 .
Thus, the second part of Equation (3) is proven. □
Theorem 2.
For the case (b), we let the shape functions be
s j ( τ ) = s j 0 τ ,   j = 1 , , n ,
where s j 0 are constant parameters. For any function y j ( τ ) C 1 [ 0 ,   1 ] ,
x j ( τ ) = y j ( τ ) s ( τ ) [ y j ( 0 ) y j ( 1 ) ] ,   j = 1 , , n
Satisfy the periodic conditions in Equation (4).
Proof. 
It follows from Equation (14) after inserting τ = 0 and τ = 1 that
x j ( 0 ) = y j ( 0 ) s j ( 0 ) [ y j ( 0 ) y j ( 1 ) ] ,
x j ( 1 ) = y j ( 1 ) s j ( 1 ) [ y j ( 0 ) y j ( 1 ) ] .
Subtracting the first by the second yields
x j ( 0 ) x j ( 1 ) = y j ( 0 ) y j ( 1 ) [ s j ( 0 ) s j ( 1 ) ] [ y j ( 0 ) y j ( 1 ) ] ,
which, with the aid of Equation (13), becomes
x j ( 0 ) x j ( 1 ) = y j ( 0 ) y j ( 1 ) y j ( 0 ) + y j ( 1 ) = 0 .
The proof is ended. □

4. Iterative Algorithms

4.1. First BSFM Iterative Algorithm

From Theorem 1, for case (a), we transform the periodic problem into the initial value problem with the help of the shape functions. For this iterative algorithm, we can consider the variable transformation from x ( τ ) to a new variable y ( τ ) by
y j ( τ ) = x j ( τ ) + G j ( τ ) ,   j = 1 , , n ,
where
G j ( τ ) = ( 1 τ ) [ y j ( 0 ) x j 0 ] + τ [ y j ( 1 ) x j 0 ] .
Equation (2) in terms of y j ( τ ) is presented by a new system of first-order ODEs:
y j ( τ ) = G j + T f j ( y 1 G 1 , , y n G n ) ,   j = 1 , , n ,
where
G j = y j ( 1 ) y j ( 0 ) ,   j = 1 , , n ,
are constants. Upon knowing T and giving y j ( 0 ) ,   j = 1 , , n as the initial values, Equation (17) can be integrated to obtain y j ( τ ) ,   j = 1 , , n .
In terms of y j . in Equation (15), T in Equation (10) reads as
T = [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) 0 1 [ y I ( ξ ) G I ( ξ ) ] 2 d ξ 0 1 x n + 1 ( ξ ) f I ( y 1 ( ξ ) G 1 ( ξ ) , , y n ( ξ ) G n ( ξ ) ) d ξ ,
which by letting
y n + 1 ( τ ) 0 τ [ y I ( ξ ) G I ( ξ ) ] 2 d ξ ,
y n + 2 ( τ ) 0 τ x n + 1 ( ξ ) f I ( y 1 ( ξ ) G 1 ( ξ ) , , y n ( ξ ) G n ( ξ ) ) d ξ ,
yields
T = [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) y n + 1 ( 1 ) y n + 2 ( 1 ) .
Next, our task is to use Equations (20)–(22) to determine the unknown period T .
Here, y j ( 1 ) in the function G j ( τ ) are unknown constants, which are denoted by c j . The first BSFM iterative algorithm for Equation (1) of the case (a) is given as follows: (i) Give x j 0 , y j ( 0 ) , j = 1 , , n , I { 1 , , n } , x n + 1 0 , T 0 , c j 0 , ϵ and Δ τ = 1 / N , and (ii) for k = 1 ,   2 ,   , with N steps, we integrate n + 3 first-order ODEs involving Equation (17) and
x n + 1 ( τ ) = y I ( τ ) G I ( τ ) ,   x n + 1 ( 0 ) = x n + 1 0 ,
y n + 1 ( τ ) = [ y I ( τ ) G I ( τ ) ] 2 ,   y n + 1 ( 0 ) = 0 ,
y n + 2 ( τ ) = x n + 1 ( τ ) f I ( y 1 ( τ ) G 1 ( τ ) , , y n ( τ ) G n ( τ ) ) ,   y n + 2 ( 0 ) = 0 ,
where the first equation is derived from Equations (5) and (15), and the last two equations are derived from Equations (20) and (21), with initial values obtained by inserting τ = 0 . Take
c j k + 1 = y j ( 1 ) ,     j = 1 , , n ,   T k + 1 = [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) y n + 1 ( 1 ) y n + 2 ( 1 ) ,
where x I ( 0 ) = x I 0 is given. The iteration is terminated until
r k j = 1 n ( c j k + 1 c j k ) 2 + ( T k + 1 T k ) 2 < ϵ .

4.2. Second BSFM Iterative Algorithm

Next, we come to a more difficult periodic problem of Equations (2) and (4). The second BSFM iterative algorithm for solving Equation (1) of case (b) is similar to that in Section 4.1, merely replacing G j by
G j ( τ ) = s j ( τ ) [ y j ( 0 ) y j ( 1 ) ] ,     s j ( τ ) = s j 0 τ .
We can obtain
x j ( 0 ) = y j ( 0 ) s j 0 [ y j ( 0 ) y j ( 1 ) ] ,     s j 0 0 ,
where y j ( 1 ) are the convergent values of the sequence c j k , k = 1 ,   2 ,   , and Equation (25) is obtained from Equation (14) by inserting τ = 0 . and using s j ( 0 ) = s j 0 .
For case (b), we can insert Equation (25) with j = I for x I ( 0 ) into Equation (23) to accelerate the speed of convergence of the computation of T k + 1 .

4.3. Third BSFM Iterative Algorithm

Theorem 1 is modified to the following for case (a) in Equation (3):
Theorem 3.
For any functions y j ( τ ) C 1 [ 0 ,   1 ] , x j ( τ ) given by Equation (14) satisfy the periodic conditions (3) for the case (a), if y j ( 0 ) are given by
y j ( 0 ) = x j 0 s j 0 y j ( 1 ) 1 s j 0 ,
where s j 0 0 and s j 0 1 . Moreover, we have
y j ( 1 ) = 1 s j 0 [ x j 0 ( 1 s j 0 ) y j ( 0 ) ] .
Proof. 
In Theorem 2, x j ( 0 ) = x j ( 1 ) has been proven for x j ( τ ) given by Equation (14); hence, we merely need to prove that x j ( τ ) given by Equation (14) satisfy x j ( 0 ) = x j 0 , if the condition (26) holds. Inserting τ = 0 into Equation (14), leads to
x j ( 0 ) = y j ( 0 ) s j ( 0 ) [ y j ( 0 ) y j ( 1 ) ] ,
Using Equation (13), which becomes
x j ( 0 ) = y j ( 0 ) s j 0 [ y j ( 0 ) y j ( 1 ) ] .
Inserting Equation (26) for y j ( 0 ) , it changes to
x j ( 0 ) = ( 1 s j 0 ) y j ( 0 ) s j 0 y j ( 1 ) = ( 1 s j 0 ) x j 0 s j 0 y j ( 1 ) 1 s j 0 + s j 0 y j ( 1 ) = x j 0 .
Hence, we prove x j ( 1 ) = x j ( 0 ) = x j 0 in Equation (3). Equation (27) follows from Equation (26) directly.
From Theorem 3, we can either specify y j ( 0 ) by Equation (26) or specify y j ( 1 ) by Equation (27). For the latter, y j ( 1 ) are obtained upon giving x j 0 and y j ( 0 ) , and we can reduce the n + 1 unknown constants from ( y j ( 1 ) , , y n ( 1 ) ,   T ) to T since it is the only unknown constant to be determined. Accordingly, we can develop the third BSFM iterative algorithm for Equation (1) of case (a): (i) Give x j 0 , y j ( 0 ) , s j 0 , j = 1 , , n , I { 1 , , n } , x n + 1 0 , initial guess of T 0 , ϵ , and Δ τ = 1 / N , and compute y j ( 1 ) , j = 1 , , n by Equation (27). (ii) For k = 1 ,   2 ,   , integrate Equation (17) with G j given by Equation (24) and
x n + 1 ( τ ) = y I ( τ ) G I ( τ ) ,   x n + 1 ( 0 ) = x n + 1 0 ,
y n + 1 ( τ ) = [ y I ( τ ) G I ( τ ) ] 2 ,   y n + 1 ( 0 ) = 0 ,
y n + 2 ( τ ) = x n + 1 ( τ ) f I ( y 1 ( τ ) G 1 ( τ ) , , y n ( τ ) G n ( τ ) ) ,   y n + 2 ( 0 ) = 0 .
Take
T k + 1 = [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) y n + 1 ( 1 ) y n + 2 ( 1 ) .
If T k + 1 converges according to
r k | T k + 1 T k | < ϵ ,
then the iteration is terminated. □
In the iterative algorithms of the above three BSFMs, the initial guesses of c j 0 , j = 1 , , n are required in the variable transformations from x j to y j ; however, the iterative algorithms are not sensitive to c j 0 , j = 1 , , n , of which c j 0 = 0 , j = 1 , , n are taken for all testing examples.

4.4. Shooting Method

Let x = ( x 1 , , x n ) be a vector. If x * is located on a periodic orbit with the minimal period T , then the Poincaré map is given by [21]:
x ( T ,   x * ) = x * ,
where x ( T ,   x * ) is a flow of the dynamical system emanating from the initial point x * . The shooting method has been described in [21], wherein the Poincaré map was solved by the Newton–Raphson iterative scheme. However, the initial guesses of x * and T must be close to the actual values and the Jacobian must be known or approximated.
We can observe that the vector periodicity condition in Equation (28) is equivalent to a scalar equation:
x ( T ,   x * ) x * = 0 ,
where x ( T ,   x * ) x * is the Euclidean norm. Instead of solving Equation (28), we can solve Equation (29) in the newly developed shooting method.
The shooting method to solve Equations (2) and (3) requires adjustment of an unknown period T , such that after integrating Equation (2) with the given initial conditions x j ( 0 ) = x j 0 , the integrated results of x j ( 1 ) can match x j ( 1 ) = x j 0 in Equation (3). Let
F ( T ) j = 1 n [ x j ( 1 ) x j 0 ] 2
be an implicit function of T and we solve
F ( T ) = 0 ,
To obtain the period T .
By introducing a fictitious time variable, ξ , Equation (31) is recast to a first-order ODE [22]:
T ( ξ ) = ν 1 + ξ F ( T ( ξ ) ) ,
Integrating it for solving F ( T ) = 0 , yields
T k + 1 = T k + ν Δ ξ 1 + ξ k F ( T k ) ,   k = 1 ,   2 ,   ,
where Δ ξ is a fictitious time increment and ξ k = k Δ ξ . The iterations in Equation (33) can be terminated if | T k + 1 T k | < ϵ 0 . The above method is termed the fictitious time-integration method [22].
For the periodic problem in Equations (2) and (4), we can develop the following shooting method: (i) Give initial guesses of T 0 and c j 0 : = x j 0 , j = 1 , , n , ν , Δ ξ , ϵ 0 and Δ τ = 1 / N . (ii) For k = 1 ,   2 ,   , integrate Equation (2) and take
c j k + 1 = x j ( 1 ) ,   j = 1 , , n ,   F ( T k ) j = 1 n [ c j k + 1 c j k ] 2 ,   T k + 1 = T k + ν Δ ξ 1 + ξ k F ( T k ) .
If c j k + 1 converge with
F ( T k ) < ϵ 0 ,
then stop. For the periodic problem (a), c j k in Equation (34) is replaced by the given x j 0 .
Remark 1.
In the shooting method, x j ( τ ) in Equation (2) or x j ( t ) in Equation (1) are obtained with a guessed value of the period T . In the first BSFM, y j ( τ ) in Equation (17) is obtained with definite initial values y j ( 0 ) being given. For the shooting method, it is not guaranteed that x j ( τ ) in Equation (2) satisfy the right boundary conditions; however, when we insert y j ( τ ) into Equation (12), Theorem 1 guarantees that x j ( τ ) satisfy the right boundary conditions x j ( τ ) = x j 0 .

5. Testing through Examples

In summary, we have developed three iterative algorithms based on the BSFM to determine the period and periodic solution of n-dimensional nonlinear dynamical system. For case (a), we developed two algorithms, namely the first (Theorem 1) and third BSFM (Theorem 3) iterative algorithms, based on either two shape functions in Equation (11) or n shape functions in Equation (13). For case (b), we developed one iterative algorithm, namely the second BSFM (Theorem 2), based on the shape functions in Equation (13). Below are examples through which the performance of these iterative algorithms is evaluated.
For solving a scalar equation f ( x ) = 0 iteratively, the numerically computed order of convergence (COC) is approximated by [23,24]:
COC : = ln   | ( x n + 1 r ) / ( x n r ) | ln   | ( x n r ) / ( x n 1 r ) | ,
where r is a solution of f ( x ) = 0 and the sequence x n are generated from an iterative scheme. In the computation of COC, we store the values of x n and take n k 0 1 , where k 0 is the number of iterations for the convergence. Indeed, COC consists of a sequence of numbers, where x n 1 r 0 , x n r 0 , and x n + 1 r 0 .
For the first BSFM to solve the periodic problem, we can generate the sequence ( T k ,   c j k ,   j = 1 , , n ) . We define
E ( k ) ( T k T e ) 2 + j = 1 n ( c j k c j e ) 2 ,
COC ( k ) ln E ( k + 1 ) / E ( k ) ln E ( k ) / E ( k 1 ) ,   k = 1 , , k 0 1 ,
ACOC 1 k 0 1 k = 1 k 0 1 COC ( k ) ,
where k 0 is the number of iterations for satisfying the convergence criterion, and we take T e = T k 0 , c j e = c j k 0 instead of the exact values, which are in general unknown. For the other two BSFMs, the processes to compute COC(k) and ACOC are the same.

5.1. Example 1

Next, the first BSFM is used to solve a two-dimensional predator–prey equation of Lotka–Volterra type [25]:
x ˙ = x + x y ,   y ˙ = y x y ,
where x is the population of predator and y is the population of prey. The system is conserved by
H ( x , y ) H ( x 0 , y 0 ) = ln   x x + ln   y y ( ln   x 0 x 0 + ln   y 0 y 0 ) = 0 ,
whose level set is a periodic orbit. For any given initial point ( x 0 ,   y 0 ) with x 0 > 0 and y 0 > 0 , a periodic orbit is passed through ( x 0 ,   y 0 ) .
Under the parameters I = 2 , x 0 = y 0 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = 0 , T 0 = 10 , x 3 0 = 1.5 , ϵ = 10 7 and N = 200 , the first BSFM is convergent with 7 iterations. T = 6.6939297 is obtained and the errors of the periodic conditions are | x ( T ) 0.5 | = 1.34 × 10 9 and | y ( T ) 0.5 | = 5.49 × 10 10 . In Table 1, the sequences of COC(n) are listed, whose average is ACOC = 1.6638. This reveals that the current iterative algorithm is rapidly, with high performance.
For the sake of comparison, the shooting method in Section 4.4 is applied with x 0 = y 0 = 0.5 , T 0 = 10 , ϵ 0 = 10 7 , N = 200 , Δ ξ = 0.005 , and ν = 300 , which is convergent with 23 iterations. T = 6.6939299 is obtained and the errors of the periodic conditions are | x ( T ) 0.5 | = 6.58 × 10 8 and | y ( T ) 0.5 | = 6.66 × 10 8 . The shooting method with ACOC = 0.9679 is more slowly convergent than the first BSFM; moreover, the first BSFM is more accurate than the shooting method to preserve the periodic conditions.
Under the parameters I = 1 , x 0 = y 0 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = 0 , T 0 = 5 , x 3 0 = 0.6 , ϵ = 10 7 , and N = 200 , the first BSFM is convergent with 16 iterations, as shown in Figure 1a by solid line. T = 6.6939298 is obtained, and | x ( T ) 0.5 | = 3.82 × 10 8 , and | y ( T ) 0.5 | = 3.9 × 10 8 are shown in Figure 1b by solid line. In this case, ACOC = 1.468 is reduced, which shows that the parameter I = 2 is better than I = 1 . If I = 1 , T 0 = 10 , x 3 0 = 0.6 , and other parameters are unchanged, the first BSFM is convergent with 76 iterations and the error in the order 10 8 . ACOC = 14.85 is large than 2. That is, if the T 0 increases, it is not easy to search for the solution under a large step size. From the above results, we can state that the BSFM obtained through the different initial guess values, I, x 3 0 , and T 0 can approximate the same periodic T. This means that the BSFM can avoid the multi-solutions problem.
If we change some parameters to x 0 = y 0 = 0.6 and x 3 0 = 0.8 and other parameters are unchanged, which is convergent with 21 iterations, as shown in Figure 1a by dashed line, T = 6.517386507 is obtained with the error in the order 10 7 , as shown in Figure 1b by dashed line. Rothe [26] developed a formula to compute the period of Lotka–Volterra system; however, the result is complicated and is not as accurate.

5.2. Example 2

It is known that
x ˙ = y ,         y ˙ = z ,         z ˙ = k y ε x 3
possess a periodic solution for arbitrary nonzero initial values with a period 2 π / k when ε = 0 .
Under the parameters I = 1 , k = 1 , ε = 0 , x 0 = 1 , y 0 = 0.5 , z 0 = 1 , y 1 ( 0 ) = y 2 ( 0 ) = y 3 ( 0 ) = 0 , T 0 = 6 , x 4 0 = 0 , ϵ = 10 5 , and N = 100 , the first BSFM is convergent with 28 iterations. T = 6.2831861 is obtained with the error in the order of 10 7 compared to the exact period 2 π . As shown in Figure 1c by dashed line, the orbit constructed is good, with the error 10 8 of the periodic conditions.
Next, we consider a perturbed example with ε = 0.001 , which is solved by the second BSFM with I = 3 , y 1 ( 0 ) = 1 , y 2 ( 0 ) = 1 , y 3 ( 0 ) = 1 , s 10 = s 20 = s 30 = 0.01 , T 0 = 6 , x 4 0 = 0.1 , and N = 200 . The number of iterations is seven under ϵ = 10 2 . T = 6.27824126 , x 0 = 1.01345 , y 0 = 1.00677 , and z 0 = 1.01347 are obtained with the errors of | x ( T ) x 0 | = 4.3 × 10 5 , | y ( T ) y 0 | = 9.89 × 10 3 , and | z ( T ) z 0 | = 1.06 × 10 4 . The periodic orbit is shown in Figure 1c by solid line. For this problem, COC(n) = 1.649 is obtained, which reveals that the second BSFM is rapidly convergent.
The shooting method in Section 4.4 is used to solve this problem with case (b), where we take the initial guesses to be x 0 = 0.5 , y 0 = 0.5 , z 0 = 0 , and T 0 = 6 , and the other parameters are ϵ 0 = 10 2 , N = 500 , Δ ξ = 0.001 , and ν = 100 . The shooting method is convergent with 45 iterations. T = 6.268137713 is obtained with the errors | x ( T ) x ( 0 ) | = 3.03 × 10 3 , | y ( T ) y ( 0 ) | = 7.17 × 10 3 , and | z ( T ) z ( 0 ) | = 4.72 × 10 3 . The shooting method is more slowly convergent than the second BSFM; moreover, the second BSFM is more accurate than the shooting method at preserving the periodic conditions. We find that the shooting method is unstable in its sensitivity to the initial guesses. For instance, with the initial guesses of x 0 = 0.5 , y 0 = 0.5 , and z 0 = 1 , the shooting method is divergent.

5.3. Example 3

We consider a jerk equation [27,28]:
x ˙ = y ,         y ˙ = z ,         z ˙ = x ( 1 ε ) ( y + z ) ε ( y + z ) ( x 2 + y 2 + z 2 ) ,
which has a periodic solution for certain nonzero initial values.
Under the parameters I = 3 , ε = 0.1 , y 1 ( 0 ) = 0 , y 2 ( 0 ) = 0.5 , y 3 ( 0 ) = 0 , s 10 = 0 , s 20 = 0.001 , s 30 = 0 , T 0 = 6 , x 4 0 = 0 , and N = 500 , the second BSFM obtains T = 6.3359174155 with the absolute error in the order 10 5 compared to the period 2 π / [ 1 ε / 12 ] obtained by the harmonic balance method [28]. x 0 = 0 , y 0 = 0.8150 and z 0 = 0 are obtained and the periodic conditions are satisfied well with the error in the order 10 3 as shown in Figure 2a by solid line.
With ε = 1 , I = 3 , s 10 = s 20 = s 30 = 0.1 , T 0 = 7 , x 4 0 = 0.01 , and N = 500 , the third BSFM based on Equation (26) is convergent with 25 iterations under ϵ = 10 3 . T = 6.9480378865 is obtained with the error in the order 10 4 , as shown in Figure 2a by dashed line.

5.4. Example 4

We consider Chen’s system [29]:
x ˙ = a ( y x ) ,     y ˙ = ( c a ) x + c y x z ,     z ˙ = x y b z .
Under the parameters I = 2 , a = 45 , b = 3 , c = 28 , x 0 = 1 , y 0 = 1.133 , z 0 = 6.648 , y 1 ( 0 ) = y 2 ( 0 ) = y 3 ( 0 ) = 0 , T 0 = 2 , x 4 0 = 0.2 , and N = 500 , the first BSFM is convergent with 13 iterations under ϵ = 10 2 . T = 2.66182489 is obtained and, as shown in Figure 2b, the orbit to preserve the periodic conditions is good, with error of 10 3 .

5.5. Example 5

We consider Lorenz’s system [30,31]:
x ˙ = σ ( y x ) , y ˙ = r x y x z , z ˙ = x y b z .
As mentioned in [31], for non-standard parameters with ( σ ,   r ,   b ) = ( 10 ,   250 ,   8 / 3 ) the Lorenz system no longer displays chaotic dynamics. Especially for the initial values of ( x 0 ,   y 0 ,   z 0 ) = ( 16.21325444114593 ,   55.78140243373939 ,   249 ) , Tucker [31] found a periodic orbit with T = 0.4600941506 .
Under the parameters I = 1 , y 1 ( 0 ) = y 2 ( 0 ) = y 3 ( 0 ) = 0 , T 0 = 0.45 , x 4 0 = 1 , and N = 500 , the first BSFM is convergent with 5 iterations under ϵ = 10 5 . T = 0.4600941568 is obtained with the error in the order of 10 6 , as shown in Figure 2c. ACOC = 1.5764 is obtained to show that the first BSFM is very rapidly convergent.

5.6. Example 6

We consider the Euler equation of a rigid body motion [25,32]:
x ˙ = ( e 3 e 2 ) y z ,     y ˙ = ( e 1 e 3 ) x z ,     z ˙ = ( e 2 e 1 ) x y .
Under the parameters e 1 = 3 , e 2 = 2 , and e 3 = 1 , and the initial values x 0 = 1 , y 0 = 2 , and z 0 = 3 , the period is computed from the elliptic integral, which is T e x = 1.4482931752118 .
With I = 1 , s 10 = s 20 = s 30 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = y 3 ( 0 ) = 0 , T 0 = 1.5 , x 4 0 = 0.1 , and N = 500 , the third BSFM based on Equation (27) is convergent with 34 iterations under ϵ = 10 15 , as shown in Figure 3a. T = 1.4482931752143 is obtained, of which the error is 10 10 , as shown in Figure 3b.

5.7. Example 7

We consider a more difficult four-dimensional ODE [33]:
x ˙ 1 = [ μ 1 x 1 + ν 1 x 2 + Λ 1 ( x 2 x 3 x 1 x 4 ) ] ,
x ˙ 2 = [ μ 1 x 2 ν 1 x 1 + Λ 1 ( x 1 x 3 + x 2 x 4 ) ] ,
x ˙ 3 = [ μ 2 x 3 + ν 2 x 4 2 Λ 2 x 1 x 2 ] ,
x ˙ 4 = [ μ 2 x 2 ν 2 x 3 + Λ 2 ( x 1 2 x 2 2 ) F 2 ] .
Under the parameters σ 1 = 1.131 , σ 2 = 0.457 . or σ 2 = 0.38 , μ 1 = 0.09 , μ 2 = 0.22 , Λ 1 = 294.732 , Λ 2 = 2213.729 , ν 1 = ( σ 1 + σ 2 ) / 2 , ν 2 = σ 2 , and F 2 = 0.003936 , the orbit tends to be periodic in the steady state after t > 5000 , as shown in Figure 4a by dashed line, where x Λ 1 Λ 2 x 1 and y Λ 1 Λ 2 x 3 . The transient trajectory emanating from the initial point x 1 ( 0 ) = 0.00059 , x 2 ( 0 ) = 0.00119 , x 3 ( 0 ) = 0.003 , and x 4 ( 0 ) = 0.0035 is not plotted in the figure.
With I = 1 , s 10 = s 20 = s 30 = s 40 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = y 3 ( 0 ) = y 4 ( 0 ) = 0 , T 0 = 6 , x 5 0 = 0.001 and N = 500 , the third BSFM to solve this problem with σ 2 = 0.457 is convergent with 58 iterations under ϵ = 10 8 . T = 8.7830547366 is obtained with the errors of the periodic conditions being 10 10 . The periodic orbit, shown by solid line in Figure 4a, coincides with that obtained by the RK4.
For σ 2 = 0.38 , we change some parameters to I = 2 , T 0 = 9 , and x 5 0 = 0.0001 , and the number of iterations reduces to 14. T = 9.043227973 is obtained and the periodic conditions feature an error in the order 10 9 . As shown in Figure 4b by solid line, the periodic orbit obtained from the third BSFM almost coincides with that obtained by the RK4, as shown by dashed line.

6. Non-Autonomous Dynamical System

For the n-dimensional non-autonomous dynamical system, we can consider
x j ( τ ) = T f j ( x 1 , ,   x n ,   T τ ) ,   j = 1 , , n ,
x n + 1 ( τ ) = x I ( τ ) , x n + 1 ( 0 ) = x n + 1 0 ,   I { 1 , , n } ,  
where the period T is determined by
T = [ x n + 1 ( 1 ) x n + 1 0 ] x I ( 0 ) 0 1 x I 2 ( τ ) d τ 0 1 x n + 1 ( τ ) f I ( x 1 ( τ ) , , x n ( τ ) ,   T τ ) d τ .
For the two cases in Equations (3) and (4), the numerical processes of these three BSFM iterative algorithms are the same.
Finally, we consider a more difficult forced Duffing equation:
x ¨ ( t ) + γ x ˙ ( t ) + α x ( t ) + β x 3 ( t ) = f 0 cos ( ω t ) ,
which can be recast to
x ( τ ) = T y ( τ ) ,   y ( τ ) = T [ f 0 cos ( ω T τ ) γ y ( τ ) α x ( τ ) β x 3 ( τ ) ] ,
where τ = t / T , and T , x ( 0 ) , and y ( 0 ) are to be determined. Under certain parameters, there are periodic motions in the phase plane [34].
Under the parameters f 0 = 0.2 , γ = 0.3 , α = 1 , β = 1, and ω = 1.2 , the orbit tends to be periodic in the steady state after t > 990 , as shown in Figure 5a by solid line. The transient trajectory emanating from the initial point x ( 0 ) = y ( 0 ) = 0 is not plotted in the figure. With f 0 = 0.2 , I = 2 , s 10 = s 20 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = 0 , T 0 = 4 , x 3 0 = 0 , and N = 200 , the third BSFM is used to solve this problem, with its initial values obtained from the RK4 to a final time t f = 1000 . T = 5.09856306156 is obtained with the periodic conditions with errors in the order 10 5 .
With f 0 = 0.28 , the orbit tends to be 1 / 2 -subharmonic in the steady state after t > 5000 , as shown in Figure 5b by solid line. The transient trajectory emanating from the initial point x ( 0 ) = y ( 0 ) = 0 is not plotted in the figure. With I = 2 , s 10 = s 20 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = 0 , T 0 = 10 , x 3 0 = 0 , and N = 200 , we apply the third BSFM to solve this problem, with its initial values obtained from the RK4 to a final time t f = 5100 . T = 10.70055549431 is obtained with the error 10 5 of the periodic conditions. On the other hand, with I = 2 , s 10 = s 20 = 0.5 , y 1 ( 0 ) = y 2 ( 0 ) = 0 , x 3 0 = 0 , ϵ = 10 7 , and N = 500 , the second BSFM solves the problem for case (b). In Table 2, we list T , x ( 0 ) , y ( 0 ) , Err 1 | x ( T ) x ( 0 ) | , Err 2 | y ( T ) y ( 0 ) | , and the number of iterations (NIs) for different initial guess of T 0 .
In Figure 6, we show the periodic orbits in the phase plane from shorter to longer periods listed in Table 2. It can be seen that 1/2-subharmonic motions with different periods in (a) T = 3.045354512 , (b) T = 7.895660986 , (c) T = 9.684802491 , and (d) T = 10.06086361 , 1/3-subharmonic motion in (e) T = 13.82756607 , and 1/4-subharmonic motion in (f) T = 18.63694088 .

7. Conclusions

We normalized a nonlinear dynamical system with an unknown period that appeared in ODEs and derived an extra formula to compute the period in the augmented space, with its evolution being one variable in the nonlinear ODEs of the dynamical system. The novelties of the present paper include the methods based on the derived shape functions for the two considered cases, which transformed the periodic problem of the nonlinear dynamical system to the initial value problem for new variables. The period, the derived formula, and the terminal values of the new variables were determined iteratively. Based on the BSFM, three novel iterative algorithms were developed, and eight numerical examples were tested to confirm the high performance of the BSFMs, in which, for most of the examples, the COC was much larger than one for a few iterations. The numerical result shows that the accuracy and convergence speed is better than with the shooting method. For solving numerical procedures, the presented BSFM is easy to implement and easy to program when determining the period and periodic solution of the nonlinear dynamical system, which features no initial values. Again, we transformed the n-dimensional periodic conditions into an equivalent nonlinear scalar equation and derived an iteration formula to determine the period based on the fictitious time-integration method (FTIM). The BSFM is better than the FTIM in terms of accuracy and convergence speed. Hence, the proposed algorithm is accurate and convergent, and better than the traditional harmonic balance method, shooting method, and FTIM for nonlinear dynamic systems.

Author Contributions

C.-S.L. contributed to the conception and supervision of the work (conceptualization, resources, methodology, writing—original draft), collected and analyzed the data, and interpreted the results. C.-W.C. contributed to the design and supervision of the work (data collection, investigation, project administration, software), and validated and visualized the data. Y.-W.C. contributed to the validation of the work (writing—review and editing, validating and visualizing the data) and the funding acquisition. Y.-S.C. contributed to the writing of the work (review and editing, software, project administration). All authors have read and agreed to the published version of the manuscript.

Funding

The corresponding authors would like to thank the Ministry of Science and Technology, Taiwan for their financial support (grant number MOST 110-2221-E-019-044).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare that there are no conflict of interest regarding the publication of this paper.

References

  1. Farkas, M. Periodic Motions; Springer: New York, NY, USA, 1994. [Google Scholar] [CrossRef]
  2. Liu, C.S.; Chen, Y.W. A simplified Lindstedt-Poincaré method for saving computational cost to determine higher order nonlinear free vibrations. Mathematics 2021, 9, 3070. [Google Scholar] [CrossRef]
  3. Donescu, P.; Virgin, L.N.; Wu, J.J. Periodic solutions of an unsymmetric oscillator including a comprehensive study of their stability characteristics. J. Sound Vib. 1996, 192, 959–976. [Google Scholar] [CrossRef]
  4. Wu, B.S.; Sun, W.P.; Lim, C.W. An analytical approximate technique for a class of strongly non-linear oscillators. Int. J. Non-Linear Mech. 2006, 41, 766–774. [Google Scholar] [CrossRef]
  5. Liu, L.; Thomas, J.P.; Dowell, E.H.; Attar, P.; Hall, K.C. A comparison of classical and high dimension harmonic balance approaches for a Duffing oscillator. J. Comput. Phys. 2006, 215, 298–320. [Google Scholar] [CrossRef]
  6. He, J.H. Variational iteration method—A kind of non-linear analytic technique: Some examples. Int. J. Non-Linear Mech. 1999, 34, 699–708. [Google Scholar] [CrossRef]
  7. Öziş, T.; Yildirim, A. A study of nonlinear oscillators with u1/3 force by He’s variational iteration method. J. Sound Vib. 2007, 306, 372–376. [Google Scholar] [CrossRef]
  8. He, J.H. A coupling method of a homotopy technique and a perturbation technique for non-linear problems. Int. J. Non-Linear Mech. 2000, 35, 37–43. [Google Scholar] [CrossRef]
  9. Shou, D.H. The homotopy perturbation method for nonlinear oscillators. Comput. Math. Appl. 2009, 58, 2456–2459. [Google Scholar] [CrossRef] [Green Version]
  10. Liu, C.S. Linearized homotopy perturbation method for two nonlinear problems of Duffing equations. J. Math. Research 2021, 13, 10–19. [Google Scholar] [CrossRef]
  11. Köroğlu, C.; Öziş, T. Applications of parameter-expanding method to nonlinear oscillators in which the restoring force is inversely proportional to the dependent variable or in form of rational function of dependent variable. Comput. Model. Eng. Sci. 2011, 75, 223–234. [Google Scholar] [CrossRef]
  12. He, J.H.; Abdou, A. New periodic solutions for nonlinear evolution equations using exp-function method. Chaos Soliton Fract. 2007, 34, 1421–1429. [Google Scholar] [CrossRef]
  13. Chu, H.P.; Lo, C.Y. Application of the differential transform method for solving periodic solutions of strongly non-linear oscillators. Comput. Model. Eng. Sci. 2011, 77, 161–172. [Google Scholar] [CrossRef]
  14. Yue, X.K.; Dai, H.H.; Liu, C.S. Optimal scale polynomial interpolation technique for obtaining periodic solutions to the Duffing oscillator. Nonlinear Dyn. 2014; 77, 1455–1468. [Google Scholar] [CrossRef]
  15. Dai, H.H.; Schnoor, M.; Atluri, S.N. A simple collocation scheme for obtaining the periodic solutions of the duffing equation, and its equivalence to the high dimensional harmonic balance method: Subharmonic oscillations. Comput. Model. Eng. Sci. 2012, 84, 459–497. [Google Scholar]
  16. Dai, H.H.; Yue, X.K.; Yuan, J.P. A time domain collocation method for obtaining the third superharmonic solutions to the duffing oscillator. Nonlinear Dyn. 2013, 73, 593–609. [Google Scholar] [CrossRef]
  17. Khan, N.A.; Liu, C.S.; Riaz, F. An optimally scaled polynomial-Fourier-series method for the numerical solution of the Duffing oscillator. Int. J. Appl. Nonlinear Sci. 2016, 2, 290–310. [Google Scholar] [CrossRef]
  18. Viswanath, D. The Lindstedt-Poincaré technique as an algorithm for computing periodic orbits. SIAM Rev. 2001, 43, 478–495. [Google Scholar] [CrossRef] [Green Version]
  19. Liu, C.S.; Chang, J.R. The periods and periodic solutions of nonlinear jerk equations solved by an iterative algorithm based on a shape function method. Appl. Math. Lett. 2020, 102, 106151. [Google Scholar] [CrossRef]
  20. Liu, C.S.; Kuo, C.L.; Chang, J.R. Solving the optimal control problems of nonlinear Duffing oscillators by using an iterative shape functions method. Comput. Model. Eng. Sci. 2020, 122, 33–48. [Google Scholar] [CrossRef]
  21. Parker, T.S.; Chua, L.O. Practical Numerical Algorithms for Chaotic Systems; Springer: New York, NY, USA, 1989. [Google Scholar] [CrossRef]
  22. Liu, C.S.; Atluri, S.N. A novel time integration method for solving a large system of non-linear algebraic equations. Comput. Model. Eng. Sci. 2008, 31, 71–83. [Google Scholar] [CrossRef]
  23. Weerakoon, S.; Fernando, T.G.I. A variant of Newton’s method with accelerated third-order convergence. Appl. Math. Lett. 2000, 13, 87–93. [Google Scholar] [CrossRef]
  24. Liu, C.S.; Hong, H.K.; Lee, T.L. A splitting method to solve a single nonlinear equation with derivative-free iterative schemes. Math. Comput. Simul. 2021, 190, 837–847. [Google Scholar] [CrossRef]
  25. Liu, C.S. Preserving constraints of differential equations by numerical methods based on integrating factors. Comput. Model. Eng. Sci. 2006, 12, 83–107. [Google Scholar] [CrossRef]
  26. Rothe, F. The periods of the Volterra-Lotka system. J. Reine Angew Math. 1985, 355, 129–138. [Google Scholar] [CrossRef]
  27. Mulholland, R.J. Non-linear oscillations of a third order differential equation. Int. J. Non-linear Mech. 1971, 6, 279–294. [Google Scholar] [CrossRef]
  28. Gottlieb, H.P.W. Harmonic balance approach to limit cycles for nonlinear jerk equations. J. Sound Vib. 2006, 297, 243–250. [Google Scholar] [CrossRef] [Green Version]
  29. Chen, G.; Ueta, T. Yet another chaotic attractor. Int. J. Bifur. Chaos 1999, 9, 1465–1466. [Google Scholar] [CrossRef]
  30. Thompson, J.M.T.; Stewart, H.B. Nonlinear Dynamics and Chaos; John Wiley and Sons: Hoboken, NJ, USA, 1986. [Google Scholar]
  31. Tucker, W. Computing accurate Poincaré maps. Physica D 2002, 171, 127–137. [Google Scholar] [CrossRef]
  32. Calvo, M.; Laburta, M.P.; Montijano, J.I.; Rández, L. Error growth in the numerical integration of periodic orbits. Math. Comput. Simul. 2011, 81, 2646–2661. [Google Scholar] [CrossRef]
  33. Nayfeh, A.H.; Balachandran, B. Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods; John Wiley & Sons: Hoboken, NJ, USA, 1995. [Google Scholar] [CrossRef]
  34. Liu, C.S. Cone of non-linear dynamical system and group preserving schemes. Int. J. Non-Linear Mech. 2001, 36, 1047–1068. [Google Scholar] [CrossRef]
Figure 1. Showing (a) convergence rates and (b) periodic orbits for example 1 and (c) two periodic orbits for example 2.
Figure 1. Showing (a) convergence rates and (b) periodic orbits for example 1 and (c) two periodic orbits for example 2.
Symmetry 14 01313 g001
Figure 2. For example 3: (a) two periodic orbits. For example 4: (b) a periodic orbit. For example 5: (c) a periodic orbit.
Figure 2. For example 3: (a) two periodic orbits. For example 4: (b) a periodic orbit. For example 5: (c) a periodic orbit.
Symmetry 14 01313 g002
Figure 3. For example 6: (a) convergence rate and (b) a periodic orbit.
Figure 3. For example 6: (a) convergence rate and (b) a periodic orbit.
Symmetry 14 01313 g003
Figure 4. For example 7 of the four-dimensional ODEs comparing the periodic solution in the steady state obtained by the RK4 and the one obtained from the third BSFM: (a) σ 2 = 0.457 and (b) σ 2 = 0.38 .
Figure 4. For example 7 of the four-dimensional ODEs comparing the periodic solution in the steady state obtained by the RK4 and the one obtained from the third BSFM: (a) σ 2 = 0.457 and (b) σ 2 = 0.38 .
Symmetry 14 01313 g004
Figure 5. For the forced Duffing oscillators solved by the third BSFM, comparing results in the phase plane by dashed line with that computed by RK4 to the steady state as shown by solid line: (a) periodic orbit and (b) 1 / 2 -subharmonic motion.
Figure 5. For the forced Duffing oscillators solved by the third BSFM, comparing results in the phase plane by dashed line with that computed by RK4 to the steady state as shown by solid line: (a) periodic orbit and (b) 1 / 2 -subharmonic motion.
Symmetry 14 01313 g005
Figure 6. For the forced Duffing oscillators solved by the second BSFM, showing periodic orbits in the phase plane from shorter to longer periods. The first four in (ad) are 1/2-subharmonic motions, (e) is a 1/3-subharmonic motion, and (f) is a 1/4-subharmonic motion.
Figure 6. For the forced Duffing oscillators solved by the second BSFM, showing periodic orbits in the phase plane from shorter to longer periods. The first four in (ad) are 1/2-subharmonic motions, (e) is a 1/3-subharmonic motion, and (f) is a 1/4-subharmonic motion.
Symmetry 14 01313 g006aSymmetry 14 01313 g006b
Table 1. The sequence of COC(n) when applying the first BSFM to example 1.
Table 1. The sequence of COC(n) when applying the first BSFM to example 1.
n123456
COC(n)4.85730.76811.25511.09291.00881.0276
Table 2. For the forced Duffing equation, the periodic data obtained from the second BSFM for different T 0 .
Table 2. For the forced Duffing equation, the periodic data obtained from the second BSFM for different T 0 .
T 0 T x 0 y 0 Err1Err2NI
19.6848024910.4694063590.137037918 3.96 × 10 8 4.36 × 10 9 22
510.060863610.4839624690.264694491 4.98 × 10 8 8.51 × 10 9 21
63.0453545120.9515416217−0.180806635 1.49 × 10 8 7.51 × 10 8 22
813.827566070.5068769586−0.194995296 4.94 × 10 8 3.27 × 10 8 28
1018.636940880.5893259280−0.354623798 4.25 × 10 9 9.48 × 10 9 26
127.8956609860.8312980124−0.3218575774 1.43 × 10 8 7.54 × 10 8 24
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, C.-S.; Chang, C.-W.; Chen, Y.-W.; Chang, Y.-S. Periodic Orbits of Nonlinear Ordinary Differential Equations Computed by a Boundary Shape Function Method. Symmetry 2022, 14, 1313. https://doi.org/10.3390/sym14071313

AMA Style

Liu C-S, Chang C-W, Chen Y-W, Chang Y-S. Periodic Orbits of Nonlinear Ordinary Differential Equations Computed by a Boundary Shape Function Method. Symmetry. 2022; 14(7):1313. https://doi.org/10.3390/sym14071313

Chicago/Turabian Style

Liu, Chein-Shan, Chih-Wen Chang, Yung-Wei Chen, and Yen-Shen Chang. 2022. "Periodic Orbits of Nonlinear Ordinary Differential Equations Computed by a Boundary Shape Function Method" Symmetry 14, no. 7: 1313. https://doi.org/10.3390/sym14071313

APA Style

Liu, C. -S., Chang, C. -W., Chen, Y. -W., & Chang, Y. -S. (2022). Periodic Orbits of Nonlinear Ordinary Differential Equations Computed by a Boundary Shape Function Method. Symmetry, 14(7), 1313. https://doi.org/10.3390/sym14071313

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop