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

Authors: Muhammad Ahsar Karim; Agus Yodi Gunawan
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–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–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–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.


## 2. Materials

Some concepts related to our discussion are here mentioned.


### 2.1. 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:(1)y¨+y+2εμẏ+εy3=ωcosβt,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:
(2)y=acos1/3βt−γ+2Λcosβt+Oε,with Λ=ω/21−β2 and a,γ are numerically obtained by the Runge-Kutta method from the system:
(3)a′=−μa−,γ′=σ−9Λ2−98a2−94aΛcosγ.


### 2.2. Fuzzy Concepts

Some concepts of the fuzzy theory are given as follows.Definition 1.
A fuzzy subsetA 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 subsetA can be expressed by a set of ordered pairs consisting of the elements x∈X and a certain degree of the membership function μAx of the form:
(4)A=x,μAx:x∈X.
(2)
The alpha-cut ofA, 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:
(5)Aα=x∈X:μAx≥α,α∈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
(1)
A is normal, that is ∃x∈ℝμAx=1
(2)
All alpha-cuts ofA are closed intervals of ℝ
(3)
The Support ofA, that is SuppA=x∈R:μAx>0, is boundedThe collection of all fuzzy subsets ofℝ is denoted by Fℝ, and the alpha-cut of A∈Fℝ is shortened by Aα=aα−,aα+, with aα−,aα+=inf,supx∈ℝ:μAx≥α.Fuzzy arithmetic for fuzzy numbers based on extension principle is given as follows:Definition 3.
LetA and B be fuzzy numbers with alpha-cuts Aα=aα−,aα+ and Bα=bα−,bα+, respectively, and δ be a real number.
(1)
The sum and the difference ofAα and Bα:
(6)A+Bα=Aα+Bα=aα−+bα−,aα++bα+,A−Bα=Aα−Bα=aα−−bα+,aα+−bα−.
(2)
The multiplication ofAα by δ:
(7)δAα=δAα=δaα−,aα+=δaα−,δaα+,δ≥0,δaα+,δaα−,δ<0.
(3)
The multiplication ofAα and Bα:
(8)A·Bα=Aα·Bα=minP,maxP,with P=aαs·bαr;s,r∈−,+.
(4)
The division ofAα by Bα, if 0∉SuppB:
(9)ABα=AαBα=Aα·1Bα,with 1/Bα=1/bα+,1/bα−.


### 2.3. Fuzzy Differential Equations (FDEs)

Some basics concepts of FDEs, i.e., fuzzy function, fuzzy differences, Seikkala differential, and Hukuhara and generalized Hukuhara differentials will be presented below.Definition 4.
LetF:a,b⟶Fℝ; a,b⊆ℝ be a fuzzy function with Fxα=fxα−,fxα+. Then, F′xα=f′xα−,f′xα+; ∀α∈0,1 and F′x∈Fℝ, which is called by Seikkala derivative of F. The fuzzy function F is called by Seikkala differentiable.Lemma 5 (see [11–13]).
LetF,G:a,b⟶Fℝ; a,b⊆ℝ be fuzzy functions. If F and G are the Seikkala differentiable, then F+G′=F′+G′ and kF′=kF′; ∀k∈ℝ.Definition 6.
LetIℝ⊂Fℝ be the family of all the fuzzy numbers and F:a,b⟶Iℝ. If the limits of some pair:
(1)
limh→0+Fx0+h−HFx0/h and limh→0+Fx0−HFx0−h/h, or
(2)
limh→0+Fx0−HFx0+h/−h and limh→0+Fx0−h−HFx0/−hexist and are equal to some element F′x0∈Iℝ, then F is strongly generalized differentiable at x0 and F′x0 is the strongly generalized derivative of F at x0.The Hukuhara difference−H has rules: if A,B∈Iℝ, then A−HB=C⟺A=B+C, provided C exists. Here, the +sign is the standard addition operation on fuzzy numbers. In terms of alpha-cut,
(10)A−HBα=aα−−bα−,aα+−bα+,with Aα=aα−,aα+ and Bα=bα−,bα+.If functionF satisfies Definition 6, (1) then F is called a Hukuhara differentiable (Hd), and if its satisfies Definition 6, (2) then F is called a generalized Hukuhara differentiable (gHd).Lemma 7 (see [12]).
LetF:a,b⟶Iℝ. For Fx∈Iℝ with Fxα=fα−x,fα+x; ∀x∈a,b then
(1)
ifF is Hd, then F′xα=fα−′x,fα+′x
(2)
ifF is gHd, then F′xα=fα+′x,fα−′x


### 2.4. Fuzzy Differential Inclusion (FDI)

The differential inclusion can be expressed in the general form:(11)y′t∈Ft,yt,y0∈Γ,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=ft,yt, y0=y0∈Γ, where f is a selection of F depending on y0=y0∈Γ.An FDI is a generalization of a differential inclusion that is defined by [19, 20]:
(12)y′t∈Ft,yt,y0∈Y0,and is interpreted as the family of differential inclusions
(13)y′t∈Ft,ytα,y0∈Y0α,for all α∈0,1. Here, Fα:0,T×ℝn⟶Ωn and Y0α∈Ω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 y0=y0∈Y0α. The set of all solutions of Equation (6) is denoted by ∑αy0,T and the attainable set at t∈0,T by Λαy0,t. Diamond [19] has proved that the sets ∑αy0,T are the alpha-cut of the fuzzy solution of Problem (5). Gomez et al. [20] has guaranteed that if F is continuous and bounded, then ∑αy0,T are defined.


### 2.5. Classical and Extended Runge-Kutta (RK) Methods

Let the system of ordinary differential equations be(14)y′t=ft,yt.The general form of the classical RK method is given by [28]:
(15)yi+1=yi+a1k1+a2k2+⋯+anknhwith the evaluation functions ki are
(16)k1=fti,yi,k2=fti+p1h,yi+q11k1h,k3=fti+p2h,yi+q21k1h+q22k2h,⋮kn=fti+pn−1h,yi+qn−1,1k1h+qn−1,2k2h+⋯+qn−1,n−1kn−1h,where ai, pi and qi,i are constants.The general form of the extended RK method is given by (Wu and Xia [25]):
(17)yn+1=yn+∑i=1mhbiki1+h2ciki2with the evaluation functions ki1 and ki2 being
(18)ki1=ftn+c¯ih,yn+h∑s=1i−1aisks1ki2=f′tn+c¯ih,yn+h∑s=1i−1aisks1where bi,ci,c¯i, and ais 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′yn≈fyn−fyn−1/h.


## 2.1. 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:(1)y¨+y+2εμẏ+εy3=ωcosβt,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:
(2)y=acos1/3βt−γ+2Λcosβt+Oε,with Λ=ω/21−β2 and a,γ are numerically obtained by the Runge-Kutta method from the system:
(3)a′=−μa−,γ′=σ−9Λ2−98a2−94aΛcosγ.


## 2.2. Fuzzy Concepts

Some concepts of the fuzzy theory are given as follows.Definition 1.
A fuzzy subsetA 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 subsetA can be expressed by a set of ordered pairs consisting of the elements x∈X and a certain degree of the membership function μAx of the form:
(4)A=x,μAx:x∈X.
(2)
The alpha-cut ofA, 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:
(5)Aα=x∈X:μAx≥α,α∈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
(1)
A is normal, that is ∃x∈ℝμAx=1
(2)
All alpha-cuts ofA are closed intervals of ℝ
(3)
The Support ofA, that is SuppA=x∈R:μAx>0, is boundedThe collection of all fuzzy subsets ofℝ is denoted by Fℝ, and the alpha-cut of A∈Fℝ is shortened by Aα=aα−,aα+, with aα−,aα+=inf,supx∈ℝ:μAx≥α.Fuzzy arithmetic for fuzzy numbers based on extension principle is given as follows:Definition 3.
LetA and B be fuzzy numbers with alpha-cuts Aα=aα−,aα+ and Bα=bα−,bα+, respectively, and δ be a real number.
(1)
The sum and the difference ofAα and Bα:
(6)A+Bα=Aα+Bα=aα−+bα−,aα++bα+,A−Bα=Aα−Bα=aα−−bα+,aα+−bα−.
(2)
The multiplication ofAα by δ:
(7)δAα=δAα=δaα−,aα+=δaα−,δaα+,δ≥0,δaα+,δaα−,δ<0.
(3)
The multiplication ofAα and Bα:
(8)A·Bα=Aα·Bα=minP,maxP,with P=aαs·bαr;s,r∈−,+.
(4)
The division ofAα by Bα, if 0∉SuppB:
(9)ABα=AαBα=Aα·1Bα,with 1/Bα=1/bα+,1/bα−.


## 2.3. Fuzzy Differential Equations (FDEs)

Some basics concepts of FDEs, i.e., fuzzy function, fuzzy differences, Seikkala differential, and Hukuhara and generalized Hukuhara differentials will be presented below.Definition 4.
LetF:a,b⟶Fℝ; a,b⊆ℝ be a fuzzy function with Fxα=fxα−,fxα+. Then, F′xα=f′xα−,f′xα+; ∀α∈0,1 and F′x∈Fℝ, which is called by Seikkala derivative of F. The fuzzy function F is called by Seikkala differentiable.Lemma 5 (see [11–13]).
LetF,G:a,b⟶Fℝ; a,b⊆ℝ be fuzzy functions. If F and G are the Seikkala differentiable, then F+G′=F′+G′ and kF′=kF′; ∀k∈ℝ.Definition 6.
LetIℝ⊂Fℝ be the family of all the fuzzy numbers and F:a,b⟶Iℝ. If the limits of some pair:
(1)
limh→0+Fx0+h−HFx0/h and limh→0+Fx0−HFx0−h/h, or
(2)
limh→0+Fx0−HFx0+h/−h and limh→0+Fx0−h−HFx0/−hexist and are equal to some element F′x0∈Iℝ, then F is strongly generalized differentiable at x0 and F′x0 is the strongly generalized derivative of F at x0.The Hukuhara difference−H has rules: if A,B∈Iℝ, then A−HB=C⟺A=B+C, provided C exists. Here, the +sign is the standard addition operation on fuzzy numbers. In terms of alpha-cut,
(10)A−HBα=aα−−bα−,aα+−bα+,with Aα=aα−,aα+ and Bα=bα−,bα+.If functionF satisfies Definition 6, (1) then F is called a Hukuhara differentiable (Hd), and if its satisfies Definition 6, (2) then F is called a generalized Hukuhara differentiable (gHd).Lemma 7 (see [12]).
LetF:a,b⟶Iℝ. For Fx∈Iℝ with Fxα=fα−x,fα+x; ∀x∈a,b then
(1)
ifF is Hd, then F′xα=fα−′x,fα+′x
(2)
ifF is gHd, then F′xα=fα+′x,fα−′x


## 2.4. Fuzzy Differential Inclusion (FDI)

The differential inclusion can be expressed in the general form:(11)y′t∈Ft,yt,y0∈Γ,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=ft,yt, y0=y0∈Γ, where f is a selection of F depending on y0=y0∈Γ.An FDI is a generalization of a differential inclusion that is defined by [19, 20]:
(12)y′t∈Ft,yt,y0∈Y0,and is interpreted as the family of differential inclusions
(13)y′t∈Ft,ytα,y0∈Y0α,for all α∈0,1. Here, Fα:0,T×ℝn⟶Ωn and Y0α∈Ω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 y0=y0∈Y0α. The set of all solutions of Equation (6) is denoted by ∑αy0,T and the attainable set at t∈0,T by Λαy0,t. Diamond [19] has proved that the sets ∑αy0,T are the alpha-cut of the fuzzy solution of Problem (5). Gomez et al. [20] has guaranteed that if F is continuous and bounded, then ∑αy0,T are defined.


## 2.5. Classical and Extended Runge-Kutta (RK) Methods

Let the system of ordinary differential equations be(14)y′t=ft,yt.The general form of the classical RK method is given by [28]:
(15)yi+1=yi+a1k1+a2k2+⋯+anknhwith the evaluation functions ki are
(16)k1=fti,yi,k2=fti+p1h,yi+q11k1h,k3=fti+p2h,yi+q21k1h+q22k2h,⋮kn=fti+pn−1h,yi+qn−1,1k1h+qn−1,2k2h+⋯+qn−1,n−1kn−1h,where ai, pi and qi,i are constants.The general form of the extended RK method is given by (Wu and Xia [25]):
(17)yn+1=yn+∑i=1mhbiki1+h2ciki2with the evaluation functions ki1 and ki2 being
(18)ki1=ftn+c¯ih,yn+h∑s=1i−1aisks1ki2=f′tn+c¯ih,yn+h∑s=1i−1aisks1where bi,ci,c¯i, and ais 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′yn≈fyn−fyn−1/h.


## 3. Main Results

### 3.1.α−cut Deterministic Systems

In the form of initial value problem, fory1=y and y1′=y2, Equation (1) can be expressed as follows:
(19)y1′=y2,y10=y1_0,y2′=−y1−2εμy2−εy13+ωcosβt,y20=y2_0,with y1_0,y2_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:
(20)y~1α′=y~2α,y~10α=y~1_0α,y~2α′=−y~1α−2εμy~2α−εy~1α3+ωcosβt,y~20α=y~2_0α,with y~1_0,y~2_0∈Fℝ. The problem of Equation (11) is called the fuzzy forced Duffing equation.Lety~α=yα−,yα+. Then, we get
(21)y1α−,y1α+′=y2α−,y2α+,y2α−,y2α+′=−y1α−,y1α+−2εμy2α−,y2α+−εy1α−,y1α+3+ωcosβt,with the initial conditions:
(22)y10α−,y10α+=y10α−,y10α+,y20α−,y20α+=y20α−,y20α+.Using the Hd concept in Lemma7 (1), then, we obtain the α−cut deterministic system of Equation (2):
(23)y1α−′=y2α−,y1α+′=y2α+,y2α−′=−y1α+−2εμy2α+−εmaxPα+ωcosβt,y2α+′=−y1α−−2εμy2α−−εminPα+ωcosβtand by using the gHd concept in Lemma 7 (2), then, we obtain
(24)y1α−′=y2α+,y1α+′=y2α−,y2α−′=−y1α−−2εμy2α−−εminPα+ωcosβt,y2α+′=−y1α+−2εμy2α+−εmaxPα+ωcosβtwith
(25)Pα=y1α−3,y1α−2y1α+,y1α−y1α+2,y1α+3.and the initial conditions:
(26)y10α−,y10α+=y10α−,y10α+,y20α−,y20α+=y20α−,y20α+.For our simulations, we take [4]:
(27)μ=0.1,ε=0.1,σ=1,β=3+εσ,ω=−4.87and the initial values:
(28)y10α−,y10α+=1.14,1.98,y2_0α−,y2_0α+=0.10,1.00.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.Figure 1
Solutions of Equation (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 a0=γ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α+.
(a)
Solution
(b)
Phase PlaneFigure 2
Solutions of Equation (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) and the initial values which are in Equation (28). For (a) and (b), similar information as in Figures 1(a) and 1(b) are applied, except in (b), the black-marked square states coordinate point yα−4.5,yα+4.5.From Figure1(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 Figure2, 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
(29)y1′=y2,y10=y1_0,y2′=−y1−2εμy2−εy13+ωcosβt,y20=y2_0,for all y1_0∈y1_0α−,y1_0α+,y2_0∈y2_0α−,y2_0α+, and α∈0,1. The fuzzy solution of System (17) is obtained as
(30)y~1α=y1α−,y1α+=miny1,maxy1,y~2α=y2α−,y2α+=miny2,maxy2,as long as the system is continuous and is bounded for given y~1_0α,y~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α+=y1α−,y1α+ 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.Figure 3
Graphs of Equation (30), i.e., the α−cut deterministic system which is obtained from applying the FDI concept. It is performed by solving Equation (29) using the extended RK method with the parameters which are in Equation (27) and the initial values which are in Equation (28). For (a) and (b), similar information as in Figures 1(a) and 1(b) are applied, except in (b), the dotted line states equation yα−=yα+.


### 3.2. Parameter Estimation

From Section3.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,
(31)ydatatα=ydataα−t,ydataα+t.Equation (31) is simulated from the approximation solution in Equation (2), with the parameters which are in Equation (27) and the initial values are a0=γ0=1). The α−cut of data fuzzy simulation is given in Figure 4.Figure 4
Theα−cut of data simulation ydataα−t,ydataα+t, with the left-end ydataα−t denoted by the unmarked square and the right one ydataα+t by the black-marked square. The solid line states the approximation solution by a multiple-scale method.To perform parameter estimation by using the least square nonlinear method, objective functionV to be optimized is given as follows:
(32)minε∗,ω∗Vε∗,ω∗22=minε∗,ω∗1M∑t=0N−1y1α−−ydataα−2+∑t=0N−1y1α+−ydataα+2,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, y1α− and y1α+ are the fuzzy solutions in Equation (30), and ydataα− and ydataα+ are the data fuzzy simulations as in Figure 4.Optimization by using objective functionV 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:
(33)ε∗=9.467×10−2,ω∗=−7.731By using the parameters in Equation (33) and the initial values y1_0α−,y1_0α+=ydataα−0,ydataα+0 and y2_0α−,y2_0α+ as in Equation (28), graphs of the solutions of Equation (29) which are performed by using the extended Runge-Kutta method are given in Figure 5. Addition graphs with the initial values y1_0α−=y1_0α+=1.51 (y2_0α−,y2_0α+ still remains) and the phase plane of the solutions of Equation (29) are given in Figure 6.Figure 5
Theα−cut of the data simulation ydataα−t,ydataα+t and the α−cut of solution of Problem (29), i.e., yα−,yα+, with the parameters as in Equation (33).Figure 6
(a) Theα−cut of solutions of Problem (29) as in Figure 5, i.e., yα−,yα+ and the ycα−,ycα+, with the parameters as in Equation (33). (b) Similar information as in Figure 3(b) are applied with the addition of a thick dashed line that represents the graph of yα−,yα+ for t∈1,4.


## 3.1.α−cut Deterministic Systems

In the form of initial value problem, fory1=y and y1′=y2, Equation (1) can be expressed as follows:
(19)y1′=y2,y10=y1_0,y2′=−y1−2εμy2−εy13+ωcosβt,y20=y2_0,with y1_0,y2_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:
(20)y~1α′=y~2α,y~10α=y~1_0α,y~2α′=−y~1α−2εμy~2α−εy~1α3+ωcosβt,y~20α=y~2_0α,with y~1_0,y~2_0∈Fℝ. The problem of Equation (11) is called the fuzzy forced Duffing equation.Lety~α=yα−,yα+. Then, we get
(21)y1α−,y1α+′=y2α−,y2α+,y2α−,y2α+′=−y1α−,y1α+−2εμy2α−,y2α+−εy1α−,y1α+3+ωcosβt,with the initial conditions:
(22)y10α−,y10α+=y10α−,y10α+,y20α−,y20α+=y20α−,y20α+.Using the Hd concept in Lemma7 (1), then, we obtain the α−cut deterministic system of Equation (2):
(23)y1α−′=y2α−,y1α+′=y2α+,y2α−′=−y1α+−2εμy2α+−εmaxPα+ωcosβt,y2α+′=−y1α−−2εμy2α−−εminPα+ωcosβtand by using the gHd concept in Lemma 7 (2), then, we obtain
(24)y1α−′=y2α+,y1α+′=y2α−,y2α−′=−y1α−−2εμy2α−−εminPα+ωcosβt,y2α+′=−y1α+−2εμy2α+−εmaxPα+ωcosβtwith
(25)Pα=y1α−3,y1α−2y1α+,y1α−y1α+2,y1α+3.and the initial conditions:
(26)y10α−,y10α+=y10α−,y10α+,y20α−,y20α+=y20α−,y20α+.For our simulations, we take [4]:
(27)μ=0.1,ε=0.1,σ=1,β=3+εσ,ω=−4.87and the initial values:
(28)y10α−,y10α+=1.14,1.98,y2_0α−,y2_0α+=0.10,1.00.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.Figure 1
Solutions of Equation (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 a0=γ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α+.
(a)
Solution
(b)
Phase PlaneFigure 2
Solutions of Equation (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) and the initial values which are in Equation (28). For (a) and (b), similar information as in Figures 1(a) and 1(b) are applied, except in (b), the black-marked square states coordinate point yα−4.5,yα+4.5.From Figure1(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 Figure2, 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
(29)y1′=y2,y10=y1_0,y2′=−y1−2εμy2−εy13+ωcosβt,y20=y2_0,for all y1_0∈y1_0α−,y1_0α+,y2_0∈y2_0α−,y2_0α+, and α∈0,1. The fuzzy solution of System (17) is obtained as
(30)y~1α=y1α−,y1α+=miny1,maxy1,y~2α=y2α−,y2α+=miny2,maxy2,as long as the system is continuous and is bounded for given y~1_0α,y~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α+=y1α−,y1α+ 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.Figure 3
Graphs of Equation (30), i.e., the α−cut deterministic system which is obtained from applying the FDI concept. It is performed by solving Equation (29) using the extended RK method with the parameters which are in Equation (27) and the initial values which are in Equation (28). For (a) and (b), similar information as in Figures 1(a) and 1(b) are applied, except in (b), the dotted line states equation yα−=yα+.


## 3.2. Parameter Estimation

From Section3.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,
(31)ydatatα=ydataα−t,ydataα+t.Equation (31) is simulated from the approximation solution in Equation (2), with the parameters which are in Equation (27) and the initial values are a0=γ0=1). The α−cut of data fuzzy simulation is given in Figure 4.Figure 4
Theα−cut of data simulation ydataα−t,ydataα+t, with the left-end ydataα−t denoted by the unmarked square and the right one ydataα+t by the black-marked square. The solid line states the approximation solution by a multiple-scale method.To perform parameter estimation by using the least square nonlinear method, objective functionV to be optimized is given as follows:
(32)minε∗,ω∗Vε∗,ω∗22=minε∗,ω∗1M∑t=0N−1y1α−−ydataα−2+∑t=0N−1y1α+−ydataα+2,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, y1α− and y1α+ are the fuzzy solutions in Equation (30), and ydataα− and ydataα+ are the data fuzzy simulations as in Figure 4.Optimization by using objective functionV 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:
(33)ε∗=9.467×10−2,ω∗=−7.731By using the parameters in Equation (33) and the initial values y1_0α−,y1_0α+=ydataα−0,ydataα+0 and y2_0α−,y2_0α+ as in Equation (28), graphs of the solutions of Equation (29) which are performed by using the extended Runge-Kutta method are given in Figure 5. Addition graphs with the initial values y1_0α−=y1_0α+=1.51 (y2_0α−,y2_0α+ still remains) and the phase plane of the solutions of Equation (29) are given in Figure 6.Figure 5
Theα−cut of the data simulation ydataα−t,ydataα+t and the α−cut of solution of Problem (29), i.e., yα−,yα+, with the parameters as in Equation (33).Figure 6
(a) Theα−cut of solutions of Problem (29) as in Figure 5, i.e., yα−,yα+ and the ycα−,ycα+, with the parameters as in Equation (33). (b) Similar information as in Figure 3(b) are applied with the addition of a thick dashed line that represents the graph of yα−,yα+ for t∈1,4.


## 4. 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 multiple-scale method.



---

*Source: 6179591-2020-04-16.xml*