Parameter Estimations of Fuzzy Forced Duffing Equation: Numerical Performances by the Extended Runge-Kutta Method

In this work, the forced Duffing equation with secondary resonance will be considered our subject by assuming that the initial values has uncertainty in terms of a fuzzy number. The resulted fuzzy models will be studied by three fuzzy differential approaches, namely, Hukuhara differential and its generalization and fuzzy differential inclusion. Applications of fuzzy arithmetics to the models lead to a set of alpha-cut deterministic systems with some additional equations. These systems are then solved by the extended Runge-Kutta method. The extended Runge-Kutta method is chosen as our numerical approach in order to enhance the order of accuracy of the solutions by including both function and its first derivative values in calculations. Among our fuzzy approaches, our simulations show that the fuzzy differential inclusion is the most appropriate approach to capture oscillation behaviors of the model. Using the aforementioned fuzzy approach, we then demonstrate how to estimate parameters to our generated fuzzy simulation data.


Introduction
The dynamic behavior of the complex systems in the real world has long been widely studied by researchers through mathematical modeling, by assuming the variables and parameters are the set of real numbers. This is of course too strict (crisp) to be used as variables or parameters sourced from the data obtained through measurements that contain uncertainty. To accommodate these uncertainties in the modeling, intensive studies are needed to describe the structure of mathematical models, develop methodologies to determine solutions of the model, and make procedures for estimating parameters of the model.
Many interesting behaviors can be observed in a system, such as nonlinear oscillation behavior. This behavior may show complex dynamics, depending on initial values and parameters. One of the mathematical models that illustrates this behavior is the Duffing equation that was first introduced by Georg Duffing in 1918. The equation is widely applied in physics and also in biology [1], disease predictions [2], and population dynamics problems [3]. The Duffing equation produces a useful model for examining nonlinear oscillations and chaotic dynamical systems. Another aspect that is interesting to observe is the presence of external force which leads to resonance phenomena, either primary or secondary resonances [4,5]. This attracted much attention to many researchers to study further in determining the solution of the model, both analytical and numerical approaches [6][7][8].
Besides the appearance of oscillation phenomena in a system, the involvement of uncertainties in the system has to be taken into account in the model. It can be caused by several factors, including limitations of available data, complexity of a system, or changes in the environment or demographics beyond the control of researchers when conducting experiments. The model which can describe uncertainties has been known in the last few decades, the so-called a fuzzy differential equation. This concept was first introduced by Chang and Zadeh [9] and currently has been developed by many other researchers with several extensions. The first proposal was given by Hukuhara [10], which is based on an interval-valued function, referred to as the Hukuhara differential. Furthermore, Seikkala proposed a fuzzy differential based on the alpha-cut concept, known as the Seikkala differential [11]. Then, Kaleva [12,13] proved that the Hukuhara differential solution is equivalent to the Seikkala differential solution and its derivatives are the same. Furthermore, the concept of Hukuhara differential was later extended to what is called generalized Hukuhara differential [14]. Later, Baidosov [15,16] used the generalization of the concept of differential inclusion to produce a new concept, known as fuzzy differential inclusion [17][18][19][20]. In principle, all the aforementioned concepts will transform the fuzzy models into what is called the alpha-cut deterministic equations, by using the fuzzy arithmetic method [21,22].
To get insight into oscillation phenomena due to an uncertainty factor, in the present study, we take the forced Duffing equation with secondary resonance as a model that represents an oscillation system having uncertainties in the initial values. We choose the secondary resonance type to provide a different oscillation behavior from our previous studies, i.e., oscillation behaviors with dumping [23] and primary resonance [24] types. Such uncertainties can be classified as fuzzy numbers, so that the equation is then called by fuzzy forced Duffing equations. We will examine the comparison of solutions from the fuzzy forced Duffing equation using these three types of fuzzy differentials, i.e., Hukuhara differential, generalized Hukuhara differential, and fuzzy differential inclusion. The alpha-cut deterministic equations generated from these three types will then be solved using the fourth-order extended Runge-Kutta method [25][26][27][28]. In contrast to the standard fourth-order Runge-Kutta method, this extended method uses new parameters to improve the accuracy of the solution by adding the first derivative of the main function to be evaluated in the calculation. This is chosen because it has been proven to be close to the exact solution than the standard method, in several types of system behavior (crisp), such as growth, logistics, and periodic models [29]. Finally, we demonstrate how to estimate parameters using the least square nonlinear method, by choosing the right fuzzy differential type and using simulated fuzzy data. The fuzzy data will be determined through an approximate solution of the multiple scale method of the forced Duffing equation with random noises.

Materials
Some concepts related to our discussion are here mentioned.

Forced Duffing Equations.
The Duffing equation is a mathematical model that illustrates nonlinear oscillation behavior and chaotic dynamical system. One of the interesting aspects to observe is the presence of external force which leads to resonance phenomena, either primary or secondary resonances. The forced Duffing equation with secondary resonance is given as follows: € y + y + 2εμ_ y + εy 3 where the y is the variable of coordinate position, which is a function of time t; the ω is the amplitude (or displacement) of the wave function; the β is the angular velocity; the ε is the damping forced (strength); and the μ is the damping controller.
The approximation solution of Equation (1) by the multiple-scale method is given as follows: with Λ = ðω/ð2ð1 − β 2 ÞÞÞ and a, γ are numerically obtained by the Runge-Kutta method from the system: 2.2. Fuzzy Concepts. Some concepts of the fuzzy theory are given as follows.

Definition 1.
A fuzzy subset A of universe X is characterized by a function μ A : X ⟶ ½0, 1, called a membership function of A, that represents the degree of membership of element in fuzzy subset A.
(1) The fuzzy subset A can be expressed by a set of ordered pairs consisting of the elements x ∈ X and a certain degree of the membership function μ A ðxÞ of the form: (2) The alpha-cut of A, denoted by ½A α , is the crisp set of all elements in X that belong to a fuzzy set A at least to the degree α ∈ ½0, 1: Definition 2. Let ℝ be a set of all real numbers and A be a fuzzy subset of ℝ. The fuzzy subset A is called by fuzzy number when The collection of all fuzzy subsets of ℝ is denoted by F ℝ , and the alpha- Fuzzy arithmetic for fuzzy numbers based on extension principle is given as follows: Definition 3. Let A and B be fuzzy numbers with alpha-cuts ½A α = ½a − α , a + α and ½B α = ½b − α , b + α , respectively, and δ be a real number.
2.4. Fuzzy Differential Inclusion (FDI). The differential inclusion can be expressed in the general form: where F : ℝ × ℝ n ⟶ Pðℝ n Þ; Pðℝ n Þ is the family of all subsets of ℝ n and Γ ⊆ ℝ n . A solution of Equation (11) is obtained by solving the differential equation y ′ ðtÞ = f ðt, y ðtÞÞ, yð0Þ = y 0 ∈ Γ, where f is a selection of F depending on yð0Þ = y 0 ∈ Γ. An FDI is a generalization of a differential inclusion that is defined by [19,20]: and is interpreted as the family of differential inclusions for all α ∈ ½0, 1. Here, ½F α : ½0, T × ℝ n ⟶ Ω n and ½Y 0 α ∈ Ω n , where Ω n is the family of all nonempty compact and convex subsets of ℝ n . A solution of Equation (6) is an absolutely continuous function y : ½0, T ⟶ ℝ n that satisfies the inclusion in ½0, T and yð0Þ = y 0 ∈ ½Y 0 α . The set of all solutions of Equation (6) is denoted by ∑ α ðy 0 , TÞ and the attainable set at t ∈ ½0, T by Λ α ðy 0 , tÞ. Diamond [19] has proved that the sets ∑ α ðy 0 , TÞ are the alpha-cut of the fuzzy solution of Problem (5). Gomez 3 Abstract and Applied Analysis et al. [20] has guaranteed that if F is continuous and bounded, then ∑ α ðy 0 , TÞ are defined.

Classical and Extended
Runge-Kutta (RK) Methods. Let the system of ordinary differential equations be The general form of the classical RK method is given by [28]: with the evaluation functions k i are where a i , p i and q i,i are constants. The general form of the extended RK method is given by (Wu and Xia [25]): with the evaluation functions k i1 and k i2 being where b i , c i , c i , and a is are constants. The f ′ is approximated by the forward difference method. Specifically, f ′ is embedded in f , i.e., f ′ is approximated by a difference quotient of past and current evaluations of f [27], f ′ ðy n Þ ≈ ð f ðy n Þ − f ðy n−1 ÞÞ/h.

Main Results
3.1. α − cut Deterministic Systems. In the form of initial value problem, for y 1 = y and y 1 ′ = y 2 , Equation (1) can be expressed as follows: with y 1 0 , y 2 0 ∈ ℝ. By assuming that initial values are fuzzy numbers, from Equation (1), we obtain the fuzzy initial value problem in the form of α − cut as follows: withỹ 1 0 ,ỹ 2 0 ∈ F ℝ . The problem of Equation (11) is called the fuzzy forced Duffing equation. Let ½ỹ α = ½y − α , y + α . Then, we get with the initial conditions: Using the Hd concept in Lemma 7 (1), then, we obtain the α − cut deterministic system of Equation (2): and by using the gHd concept in Lemma 7 (2), then, we obtain with and the initial conditions: 4 Abstract and Applied Analysis For our simulations, we take [4]: The graphs of the fuzzy solutions and its phase plane of Equations (23) and (24) by using the extended RK method are given in Figures 1 and 2, respectively.
From Figure 1(a), the graph of y − α decreases from the beginning, while the graph of y + α decreases for a moment and then rises; then, both move away from the graph of approximation solution y quickly. The α − cut of the solutions resulted by the Hd concept did not show an oscillatory behavior. It means that the uncertainty of solutions increases even since the beginning of the evolution. This is also shown through the phase plane in Figure 1(b). The coordinates of the phase plane curve starts at the origin ðy − α ð0Þ, y + α ð0ÞÞ, then the curve of y − α decreases, whereas conversely, the curve of y + α decreases for a moment then increases. In Figure 2, almost the same thing with that of Figure 1 happened to the α − cut of the solutions resulted by the gHd concept. The difference is that there is once oscillation before the uncertainty increases. This result does not indicate a switch point, which is different from the our previous results [23,24], where both of them also examine the oscillation behavior of the nonlinear fuzzy model. In Karim et al. [23], the solution of the fuzzy harmonic oscillator model that uses the gHd concept produces several oscillations with the appearance of the switch point, whereas in Karim et al. [24], the solution of the fuzzy Goodwin model that uses the same differential type produces oscillations throughout evolution with the appearance of switch points in each period.
On the other hand, an FDI concept of Equation (2) is the family of all differential inclusions of    (23), the α − cut deterministic system which is obtained from applying the Hd concept. It is performed by using the extended RK method with the parameters which are in Equation (27) and the initial values which are in Equation (28). (a) The α − cut ½y − α , y + α , with the left-end y − α is denoted by the no marked square, and the right one y + α by the black-marked square. In addition, the solid line states the approximation solution y of Equation (2), with að0Þ = γð0Þ = 1 [4]. (b) Phase plane ðy − α , y + α Þ, with y − α placed horizontally and y + α vertically. In addition, the black-marked cycle states coordinate point ðy − α ð0Þ, y + α ð0ÞÞ, the black-marked square states coordinate point ðy − α ð2:5Þ, y + α ð2:5ÞÞ, and the dotted line states equation y − α = −y + α .

Abstract and Applied Analysis
The fuzzy solution of System (17) is obtained as as long as the system is continuous and is bounded for given ½ỹ 1 0 α , ½ỹ 2 0 α . By solving Equation (29) using the extended RK method, with the parameters in Equation (23) and the initial values in Equation (24), the graphs of ½y − α , y + α = ½y 1 − α , y 1 + α in Equation (18) are obtained as in Figure 3. By using this concept, we may capture the oscillatory behavior of the approximation solution of Equation (2), with the uncertainty which also oscillates.

Parameter Estimation.
From Section 3.1, we find that the concept of FDI is able to capture oscillatory behavior and maintain the uncertainty of the solution of fuzzy forced Duffing equations. This leads us to choose Equation (29) as the basis for estimating the parameters of the model. Parameter estimation is performed by using a nonlinear least square (lsqnonlin) method.
To illustrate the process, we set the α − cut of data fuzzy simulation, namely, Equation (31) is simulated from the approximation solution in Equation (2), with the parameters which are in Equa-tion (27) and the initial values are að0Þ = γð0Þ = 1). The α − cut of data fuzzy simulation is given in Figure 4.
To perform parameter estimation by using the least square nonlinear method, objective function V to be optimized is given as follows: where M = 2N ; N is the number of the data series of t, ε * , and ω * are the amplitude of the wave function and the damping forced (respectively) to be estimated, y 1 − α and y 1 + α are the fuzzy solutions in Equation (30), and y data − α and y data + α are the data fuzzy simulations as in Figure 4.
Optimization by using objective function V in Equation (32) requires an initial value, and also, optional lower and upper bounds lb and ub on the parameters are to be estimated. Therefore, the parameters ε and ω in Equation (27), which are used in the approximation solution in Equation (2), are chosen as initial parameters, namely, ε 0 = 0:1 and ω 0 = −4:87. Then, by choosing lb = 0, ub = 1 for the ε * and also lb = −10, ub = 0 for the ω * , optimization by using Equation (32) produces parameters:  (24), the α − cut deterministic system which is obtained from applying the gHd concept. It is performed by using the extended RK method with the parameters which are in Equation (27)

Concluding Remarks
Three fuzzy differential concepts were examined by using the extended Runge-Kutta method to capture the oscillatory behavior of the forced Duffing equation with secondary resonance. Neither the Hd nor gHd concepts can capture the oscillation. Conversely, the concept of FDI was able to capture the oscillation of the equation and   maintained the uncertainty of the fuzzy forced Duffing equation. This prompted us to apply the concept of FDI to estimate the parameters of the fuzzy equation to a set data fuzzy simulation, which was performed by using the nonlinear least square method. The data fuzzy was simulated from the approximation solution by the multiplescale method.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no conflicts of interest.   Figure 3(b) are applied with the addition of a thick dashed line that represents the graph of ðy − α , y + α Þ for t ∈ ð1, 4Þ.