Series Solution of the System of Fuzzy Differential Equations

Abstract

The homotopy analysis method (HAM) is proposed to obtain a semianalytical solution of the system of fuzzy differential equations (SFDE). The HAM contains the auxiliary parameter ħ, which provides us with a simple way to adjust and control the convergence region of solution series. Concept of ħ-meshes and contour plots firstly are introduced in this paper which are the generations of traditional h-curves. Convergency of this method for the SFDE has been considered and some examples are given to illustrate the efficiency and power of HAM.

1. Introduction

In many cases of the modeling real world phenomena, information about the behavior of a dynamical system is uncertain. In order to obtain a more realistic model, we have to take into account these uncertainties. Since 1965, when Zadeh published his pioneering paper [1], hundreds of examples have been supplied where the nature of uncertainty in the behavior of a given system processes is fuzzy rather than stochastic nature. The concept of fuzzy derivative was first introduced by Chang and Zadeh in [2]. It was followed up by Dubois and Prade in [3], who defined and used the extension principle. Other methods have been discussed by Puri and Ralescu in [4] and Goetschel and Voxman in [5]. The initial value problem for fuzzy differential equation (FIVP) has been studied by Kaleva in [6, 7] and by Seikkala in [8].

The purpose of this paper is to find the approximate solution of fuzzy differential equations system with the homotopy analysis method (HAM) introduced first by Liao in 1992 [9, 10], that is, analytic approach to get series solutions of various types of linear and nonlinear equations.

Some of numerical methods have been applied to obtain the solution of fuzzy differential equations [11–15]. Sami Bataineh et al. have applied the HAM for systems of ODEs and PDEs in [16, 17]. Also, recently many types of nonlinear problems solved with HAM by others [18–22].

In Section 2, some basic definitions which will be used later in the paper are provided. In Sections 3 and 4, system of fuzzy differential equations and then basic ideas of HAM applied to these types of equations have been reviewed, respectively. Convergency of HAM for SFDE that shows its reliability has been considered in Section 5. The proposed method is illustrated by solving several examples in Section 6, and finally the conclusion is drawn in Section 7.

2. Preliminaries

In this section, the most basic notations used in this paper are introduced.

Definition 1 (See [23]). An arbitrary fuzzy number 𝑢 in parametric form is an ordered pair (𝑢(𝑟),𝑢(𝑟)) of functions 𝑢(𝑟), 𝑢(𝑟); 0≤𝑟≤1, which satisfies the following requirements:(1)𝑢(𝑟) is a bounded left-continuous nondecreasing function over [0,1],(2)𝑢(𝑟) is a bounded left-continuous nonincreasing function over [0,1],(3)𝑢(𝑟)≤𝑢(𝑟), 0≤𝑟≤1.
The set of all such fuzzy numbers is represented by 𝔼¹.

Remark 2. For arbitrary 𝑢=(𝑢(𝑟),𝑢(𝑟)),  𝑣=(𝑣(𝑟),𝑣(𝑟)), and 𝑘∈ℝ, we define addition and multiplication by 𝑘 as (𝑢+𝑘𝑣𝑢)(𝑟)=𝑢(𝑟)+𝑣(𝑟),(𝑢+𝑣)(𝑟)=𝑢(𝑟)+𝑣(𝑟),𝑘𝑢(𝑟)=𝑘𝑢(𝑟),𝑘𝑢(𝑟)=𝑘𝑢(𝑟),if𝑘≥0,(𝑟)=𝑘𝑢(𝑟),𝑘𝑢(𝑟)=𝑘𝑢(𝑟),if𝑘<0.(1)

Definition 3. For arbitrary fuzzy numbers 𝑢, 𝑣∈𝔼¹, we use the distance 𝐷(𝑢,𝑣)=sup_0≤𝑟≤1max||𝑢(𝑟)−𝑣(𝑟)||,||𝑢(𝑟)−𝑣(𝑟)||,(2) and it is shown that (𝔼¹,𝐷) is a complete metric space.

Definition 4. Let 𝑓∶[𝑎,𝑏]→𝔼¹, for each partition 𝑃={𝑡₀,…,𝑡_𝑛} of [𝑎,𝑏] and for arbitrary 𝜉_𝑖∈[𝑡_𝑖−1,𝑡_𝑖], 1≤𝑖≤𝑛 suppose 𝑅_𝑝=^𝑛_𝑖=1Δ∶=max||𝑡_𝑖𝑓𝜉_𝑖𝑡_𝑖−𝑡_𝑖−1,−𝑡_𝑖−1||,𝑖=1,…,𝑛.(3) The definite integral of 𝑓(𝑡) over [𝑎,𝑏] is ^𝑏𝑎𝑓(𝑡)𝑑𝑡=lim_Δ→0𝑅_𝑝,(4) provided that this limit exists in the metric 𝐷. If the fuzzy function 𝑓(𝑡) is continuous in the metric 𝐷, its definite integral exists and also, ^𝑏𝑎𝑓(𝑡;𝑟)𝑑𝑡=^𝑏𝑎𝑓(𝑡;𝑟)𝑑𝑡,^𝑏𝑎𝑓(𝑡;𝑟)𝑑𝑡=^𝑏𝑎𝑓(𝑡;𝑟)𝑑𝑡.(5)

Definition 5. Let 𝑓∶ℝ¹→𝔼¹ be a fuzzy function and let 𝑡₀∈ℝ¹. The derivative 𝑓′(𝑡₀) of 𝑓 at the point 𝑡₀ is defined by 𝑓^𝑡₀=lim_ℎ→0𝑓𝑡₀+ℎ−𝑓𝑡₀ℎ,(6) provided that this limit, taken with respect to the metric 𝐷, exists.

The elements 𝑓(𝑡₀+ℎ), 𝑓(𝑡₀) at the right-hand side of (6) are observed as elements in the Banach space 𝔹=ℂ[0,1]×ℂ[0,1]. Thus, if 𝑓(𝑡₀+ℎ)=(𝑎,𝑎) and 𝑓(𝑡_𝑜)=(𝑏,𝑏), the difference is simply 𝑓(𝑡₀+ℎ)−𝑓(𝑡_𝑜)=(𝑎−𝑏,𝑎−𝑏).

Clearly [𝑓(𝑡₀+ℎ)−𝑓(𝑡₀)]/ℎ may not be a fuzzy number for all ℎ. However, if it approaches 𝑓′(𝑡₀) (in 𝔹) and 𝑓′(𝑡₀) is also a fuzzy number (i.e., in 𝔼¹), this number is the fuzzy derivative of 𝑓(𝑡) at 𝑡₀. In this case, if 𝑓=(𝑓,𝑓), it can be easily shown that 𝑓^𝑡₀=𝑓^𝑡₀,𝑓^𝑡₀,(7) where 𝑓^ and 𝑓^ are the classic derivatives of 𝑓 and 𝑓, respectively.

3. System of Fuzzy Differential Equations

In this section, we will review system of fuzzy differential equations of the forms ^𝑛_𝑖=0𝑌𝐴_𝑖(𝑡)𝑌^(𝑖)(𝑡)=𝐹(𝑡),_(𝑖)(0)=𝐺_𝑖,𝑖=0,…,𝑛−1,(8) where 𝑡 is a scaler and 𝐴₀(𝑡),𝐴₁(𝑡),…,𝐴_𝑛−1(𝑡) are 𝑠×𝑠 matrixes and every component of them is a real function of 𝑡 and 𝐴_𝑛(𝑡)=𝐼_𝑠 denotes the 𝑠×𝑠 identity matrix. 𝑌, 𝐺_𝑖, and 𝐹 are fuzzy 𝑠-dimensional vectors. The 𝑟th component of 𝑌∈𝐸^𝑠 will be denoted by 𝑦_𝑟, so that we may write 𝑦_𝑖𝑌=𝑦₁(𝑡,𝑟),𝑦₂(𝑡,𝑟),…,𝑦_𝑠(𝑡,𝑟)^𝑇,(𝑡,𝑟)𝐹==𝑓𝑦_𝑙1_𝑖(𝑡,𝑟),𝑦_𝑢_𝑖(𝑡,𝑟),𝑖=1,…,𝑠,𝑓_𝑖(𝑡,𝑟),𝑓₂(𝑡,𝑟),…,𝑓_𝑠(𝑡,𝑟)^𝑇,(𝑡,𝑟)=𝐺𝑓_𝑙_𝑖(𝑡,𝑟),𝑓_𝑢_𝑖(𝑡,𝑟),𝑖=1,…,𝑠,_𝑖=𝑔_𝑖1(𝑟),𝑔_𝑖2(𝑟),…,𝑔_𝑖𝑠(𝑟)^𝑇,(9) where 𝑔_𝑖𝑗𝑖(=𝑟)1=,…,𝑔𝑛_𝑙−_𝑖𝑗1(𝑟,)𝑗,=𝑔_𝑢1,_𝑖𝑗…(𝑟,)𝑠.,(10)

The superscript 𝑇 denotes transpose. The component that is in 𝑟th row, 𝑘th column, 1≤𝑟,  𝑘≤𝑠 of matrix 𝐴_𝑝, 0≤𝑝≤𝑛 will be denoted by 𝑎^{𝑝𝑟,𝑘}. Then, (8) can be replaced by the following equivalent system in parametric form: ^𝑛_𝑠𝑖=0_𝑘=1𝑎^{𝑖𝑑,𝑘}(𝑡)𝑦^𝑘(𝑖)(𝑡,𝑟)=𝑓_𝑑(𝑡,𝑟),𝑑=1,…,𝑠,_𝑛_𝑠𝑖=0_𝑘=1𝑦𝑎^{𝑖𝑑,𝑘}(𝑡)𝑦^𝑘(𝑖)(𝑡,𝑟)=𝑓_𝑑(𝑡,𝑟),𝑑=1,…,𝑠,_𝑗(𝑖)(0,𝑟)=𝑔_𝑖𝑗(𝑟),𝑖=0,…,𝑛−1,𝑗=1,…,𝑠,𝑦^𝑗(𝑖)(0,𝑟)=𝑔_𝑖𝑗(𝑟),𝑖=0,…,𝑛−1,𝑗=1,…,𝑠,(11) where 𝑟∈[0,1].

We define functions 𝑏^{𝑖𝑑,𝑘}(𝑡) and 𝑐^{𝑖𝑑,𝑘}(𝑡); 1≤𝑑, 𝑘≤𝑠,  0≤𝑖≤𝑛, in the following form: 𝑏^{𝑖𝑑,𝑘}(𝑡)=0,𝑎^{𝑖𝑑,𝑘}𝑎(𝑡)<0,_{𝑖𝑑,𝑘}𝑐(𝑡),𝑎^{𝑖𝑑,𝑘}(𝑡)≥0,_{𝑖𝑑,𝑘}(𝑡)=𝑎^{𝑖𝑑,𝑘}0,(𝑡),𝑎𝑎^{𝑖𝑑,𝑘}_{𝑖𝑑,𝑘}(𝑡)<0,(𝑡)≥0.(12) Now (11) becomes as follows: ^𝑛_𝑠𝑖=0_𝑘=1𝑏^{𝑖𝑑,𝑘}(𝑡)𝑦_𝑙^𝑘(𝑖)(𝑡,𝑟)+𝑐^{𝑖𝑑,𝑘}(𝑡)𝑑𝑦=_𝑢1^𝑘(𝑖),…(𝑡,𝑟,)𝑠_𝑛=𝑓_𝑙_𝑑(𝑡,𝑟),_𝑠𝑖=0_𝑘=1𝑐^{𝑖𝑑,𝑘}𝑦(𝑡)𝑦_𝑙^𝑘(𝑖)(𝑡,𝑟)+𝑏^{𝑖𝑑,𝑘}(𝑡)𝑑𝑦=_𝑢1,^𝑘(𝑖)…(𝑡,,𝑟𝑠_𝑙)=𝑓_𝑢_𝑑(𝑡,𝑟),𝑦^𝑗(𝑖)(0,𝑟)=𝑔_𝑙_𝑖𝑗(𝑟),𝑖=0,…,𝑛−1,𝑗=1,…,𝑠,_𝑢^𝑗(𝑖)(0,𝑟)=𝑔_𝑢_𝑖𝑗(𝑟),𝑖=0,…,𝑛−1,𝑗=1,…,𝑠.(13)

4. Basic Ideas of HAM

We consider the following differential equations: ^𝑛_𝑖=0𝐴_𝑖(𝑡)𝑍^(𝑖)(𝑡,𝑟)−𝐹(𝑡,𝑟)=0,(14) where 𝐹𝑍((𝑡𝑡,,𝑟𝑟))==𝑓𝑧₁₁(𝑡,𝑟),𝑧₂(𝑡,𝑟),…,𝑧_𝑠(𝑡,𝑟)^𝑇,(𝑡,𝑟),𝑓₂(𝑡,𝑟),…,𝑓_𝑠(𝑡,𝑟)^𝑇,(15) and 𝐴₀(𝑡),𝐴₁(𝑡),…,𝐴_𝑛−1(𝑡) are 𝑠×𝑠 matrixes and every component of them is a real function of 𝑡 and 𝐴_𝑛(𝑡)=𝐼_𝑠, where 𝐼_𝑠 denotes the 𝑠×𝑠 identity matrix.

The above system of equations can be written in following form: 𝒩_𝑖𝑧₁(𝑡,𝑟),𝑧₂(𝑡,𝑟),…,𝑧_𝑠(𝑡,𝑟)=0,𝑖=1,2,…,𝑠,(16) where 𝒩_𝑖 are nonlinear operators that represent the whole equations, 𝑡 and 𝑟 denote the independent variables, and 𝑧_𝑖(𝑡,𝑟),  𝑖=1,2,…,𝑠 are unknown functions, respectively. By means of generalizing the traditional homotopy method constructed the so-called zero-order deformation equations (1−=𝑞)𝑞ℒℏ𝒩_𝑖_𝑖𝜙_𝑖(𝑡,𝑟;𝑞)−𝑧_𝑖,0(𝑡,𝑟)𝜙₁(𝑡,𝑟;𝑞𝑖)=,…1,…,𝜙,_𝑛𝑛(𝑡,,𝑟;𝑞),(17) where 𝑞∈[0,1] is an embedding parameter, ℏ is nonzero auxiliary parameter, ℒ_𝑖,  𝑖=1,2,…,𝑠 are auxiliary linear operators, 𝑧_𝑖,0(𝑡,𝑟) are initial guesses of 𝑧_𝑖(𝑡,𝑟), and 𝜙_𝑖(𝑡,𝑟;𝑞) are unknown functions. It is important to note that one has great freedom to choose auxiliary objects such as ℏ and ℒ_𝑖 in HAM. Obviously, when 𝑞=0 and 𝑞=1, both 𝜙_𝑖(𝑡,𝑟;0)=𝑧_𝑖,0(𝑡,𝑟),𝜙_𝑖(𝑡,𝑟;1)=𝑧_𝑖(𝑡,𝑟)(18) hold. Thus, as 𝑞 increases from 0 to 1, the solutions 𝜙_𝑖(𝑡,𝑟;𝑞) varies from the initial guesses 𝑧_𝑖,0(𝑡,𝑟) to the solutions 𝑧_𝑖(𝑡,𝑟). Expanding 𝜙_𝑖(𝑡,𝑟;𝑞) in Taylor’s series with respect to 𝑞, one has 𝜙_𝑖(𝑡,𝑟;𝑞)=𝑧_𝑖,0(𝑡,𝑟)+^+∞_𝑚₌₁𝑧_𝑖,𝑚(𝑡,𝑟)𝑞^𝑚,(19) where 𝑧_𝑖,𝑚=1𝑚!𝜕^𝑚𝜙_𝑖(𝑡,𝑟;𝑞)𝜕𝑞^𝑚∣_𝑞=0.(20) If auxiliary linear operator, initial guesses, and auxiliary parameter ℏ are properly chosen, then the series equations (19) converges at 𝑞=1 and 𝜙_𝑖(𝑡,𝑟;1)=𝑧_𝑖,0(𝑡,𝑟)+^+∞_𝑚₌₁𝑖=1,…𝑧,_𝑖,𝑚𝑠,(𝑡,𝑟),(21) which must be solution of the original nonlinear equations.

According to (20), the governing equations can be deduced from the zero-order deformation equations (17). Define the vectors ^→𝑧_𝑖,𝑚=𝑧_𝑖,0(𝑡,𝑟),𝑧_𝑖,1(𝑡,𝑟),…,𝑧_𝑖,𝑚(𝑡,𝑟),𝑖=1,…,𝑠.(22) Differentiating (17) 𝑚 time with respect to the embedding parameter 𝑞 and then setting 𝑞=0 and finally dividing them by 𝑚!, we have the so-called 𝑚th order deformation equations ℒ_𝑖𝑧=_𝑖,𝑚ℏ(ℛ𝑡,𝑟_𝑖,𝑚)−𝜒_𝑚𝑧_{𝑖,𝑚−1}(𝑡,𝑟)^→𝑧_1,𝑚−1,…,^→𝑧_{𝑠,𝑚−1}𝑖=1,…,𝑠,(23) where ℛ_𝑖,𝑚^→𝑧_1,𝑚−1,…,^→𝑧_{𝑠,𝑚−1}=1(𝑚−1)!𝜕^𝑚−1𝒩_𝑖𝜙₁(𝑡,𝑟;𝑞),…,𝜙_𝑠(𝑡,𝑟;𝑞)𝜕𝑞^𝑚−1𝜒𝑖=1,…,𝑠,_𝑚∣_𝑞=0,=⎧⎪⎨⎪⎩01,,𝑚𝑚≤>11,.(24)

Now, to simplify and solve (11), we need to know the sign of components of 𝐴_𝑖(𝑡);  0≤𝑖≤𝑛 and use of Remark 2. Also we should be noted that in case that some components of 𝐴₀(𝑡),𝐴₁(𝑡),…,𝐴_𝑛(𝑡) for some 𝑡 are negative, we must divide the defined interval for 𝑡 into small intervals, so that sign of each component of 𝐴₀(𝑡),𝐴₁(𝑡),…,𝐴_𝑛(𝑡) in each small interval be unchanged. Now for each small interval we have a separate system 2𝑠×𝑠(𝑛+1).

5. Convergency of HAM for System of Fuzzy Differential Equations

In this section, we prove a theorem which shows the convergency of approximate HAM solution applied for (8).

First we define ℛ_𝑘,𝑚[(𝑦^→_𝑙)_1,𝑚−1,…,(𝑦^→_𝑙)_{𝑠,𝑚−1}, (𝑦^→_𝑢)_1,𝑚−1,…,(𝑦^→_𝑢)_{𝑠,𝑚−1}]; 1≤𝑘≤𝑠, in the following form:ℛ_𝑘,𝑚=1(𝑚−1)!𝜕^𝑚−1𝒩_𝑘𝑦^𝑙_1,𝑚−1,…,𝑦_𝑙_{𝑠,𝑚−1},𝑦_𝑢_1,𝑚−1,…,𝑦_𝑢_{𝑠,𝑚−1}𝜕𝑞^𝑚−1∣𝑞=0,(25)

if 𝑘th equation of (13) is inclusive 𝑏^{𝑖𝑑,𝑘}(𝑡)(𝑦_𝑙)^𝑘(𝑖)(𝑡,𝑟);  1≤𝑖≤𝑛, 1≤𝑑≤𝑠, and otherwiseℛ_𝑘,𝑚=1(𝑚−1)!𝜕^𝑚−1𝒩_𝑘+𝑠𝑦^𝑙_1,𝑚−1,…,𝑦_𝑙_{𝑠,𝑚−1},𝑦_𝑢_1,𝑚−1,…,𝑦_𝑢_{𝑠,𝑚−1}𝜕𝑞^𝑚−1∣𝑞=0.(26)

Theorem 6. Let the series ∑^∞𝑚=1(𝑦)_𝑘,𝑚(𝑡,𝑟);1≤𝑘≤𝑠 be uniformly converge to (𝑦)_𝑘(𝑡,𝑟); 1≤𝑘≤𝑠. Then ((𝑦_𝑙)_𝑘(𝑡,𝑟),(𝑦_𝑢)_𝑘(𝑡,𝑟)),1≤𝑘≤𝑠 is exact solution of (13).

Proof. Since ∑^∞𝑚=1(𝑦)_𝑘,𝑚(𝑡,𝑟),1≤𝑘≤𝑠 is convergent, we must have lim_𝑚→∞(𝑦)_𝑘,𝑚(𝑡,𝑟)=0,1≤𝑘≤𝑠,0≤𝑡≤𝑇,0≤𝑟≤1.(27) On the other hand, since ℒ_𝑖 are linear operators, thus ^𝑛_𝑚=1=ℒ_𝑘(𝑦)_𝑘,𝑚(𝑡,𝑟)−𝜒_𝑚(𝑦)_{𝑘,𝑚−1}(𝑡,𝑟)_𝑛_𝑚₌₁ℒ_𝑘=ℒ_𝑘(𝑦)_𝑘,𝑚(𝑡,𝑟)−𝜒_𝑚ℒ_𝑘(𝑦)_{𝑘,𝑚−1}(𝑡,𝑟)−ℒ(_𝑘𝑦)_𝑘,1(𝑡,𝑟)+ℒ_𝑘(𝑦)_𝑘,2(𝑡,𝑟)+⋯+(ℒ𝑦)_𝑘,1𝑘(𝑡,𝑟)=ℒ_𝑘(𝑦)_𝑘,𝑛(𝑡,𝑟)−ℒ_𝑘(𝑦)_{𝑘,𝑛−1}(𝑡,𝑟)(𝑦)_𝑘,𝑛(𝑡,𝑟),𝑘=1,2,…,𝑠,(28) then from (27) we can write ^∞_𝑚=1ℒ=_𝑘li(m𝑦_𝑛→∞)_𝑘,𝑚(𝑡,𝑟)−𝜒_𝑚(𝑦)_{𝑘,𝑚−1}(𝑡,𝑟)ℒ_𝑘(𝑦)_𝑘,𝑛(𝑡,𝑟)=ℒ_𝑘lim_𝑛→∞𝑘=1,(2𝑦,)…_𝑘,𝑛,(𝑠𝑡,,𝑟)=0,(29) hence, ℏ^∞_𝑚₌₁ℛ_𝑘,𝑚^→_→𝑦_𝑙_1,𝑚−1,…,𝑦^→_𝑙_{𝑠,𝑚−1},𝑦_𝑢_1,𝑚−1,…,𝑦^→_𝑢_{𝑠,𝑚−1}=_∞_𝑚₌₁ℒ_𝑘(𝑦)_𝑘,𝑚(𝑡,𝑟)−𝜒𝑘_𝑚(=𝑦)1,2,_{𝑘,𝑚−1}…(𝑡,,𝑟𝑠).=0,(30) Since ℏ≠0, thus from the above equations for 𝑘=1,2,…,𝑠, we can write ^∞_𝑚=1ℛ_𝑘,𝑚^→𝑦_𝑙_1,𝑚−1,…,𝑦^→_𝑙_{𝑠,𝑚−1},𝑦^→_𝑢_1,𝑚−1,…,𝑦^→_𝑢_{𝑠,𝑚−1}=0.(31) From uniform convergency we have ^∞_𝑚=1𝜕^𝑖𝜕𝑡^𝑖(𝑦)_{𝑘,𝑚−1}(𝑡,𝑟)=𝜕^𝑖𝜕𝑡^𝑖(𝑦)_𝑘(𝑡𝑘,𝑟=),1,𝑖2,=…0,,1𝑠,,…,𝑛,(32) thus 0=^∞_𝑚₌₁^𝑛_𝑖=0=𝐴_𝑖(𝑡)𝑌^(𝑖)(𝑡)−1−𝜒_𝑚(𝐹)_𝑑(𝑡,𝑟)_𝑛_𝑖=0𝐴_𝑖(𝑡)^∞_𝑚₌₁𝑌^(𝑖)(𝑡)−^∞_𝑚₌₁1−𝜒_𝑚(𝐹)_𝑑=(𝑡,𝑟)_𝑛_𝑖=0𝐴_𝑖(𝑡)^∞_𝑚₌₁𝜕^𝑖𝜕𝑡^𝑖=(𝑌)_{𝑘,𝑚−1}(𝑡,𝑟)−(𝐹)_𝑑(𝑡,𝑟)_𝑛_𝑖=0𝐴_𝑖(𝑡)𝜕^𝑖𝜕𝑡^𝑖(𝑌)_𝑘(𝑡,𝑟)−(𝐹)_𝑑(𝑡,𝑟),𝑑=1,…,𝑠.(33) Therefore, {𝑦_𝑘(𝑡,𝑟)=((𝑦_𝑙)_𝑘(𝑡,𝑟),(𝑦_𝑢)_𝑘(𝑡,𝑟));𝑘=1,2,…𝑠} is the exact solution of (11) and 𝑌=[𝑦₁(𝑡,𝑟),𝑦₂(𝑡,𝑟),…,𝑦_𝑠(𝑡,𝑟)]^𝑇  is the exact solution of (8) and proof is completed.

6. Illustrative Examples

Example 7. Consider the following second-order fuzzy linear differential equation: 𝑌𝑌(0𝑌)^=+(0𝑌.1=𝑟−−𝑡0,.1,𝑡0∈.1[−0,01.]1𝑟),_(0)=(0.088+0.1𝑟,0.288−0.1𝑟).(34) The exact solution is as follows: 𝑌(𝑡,𝑟)=(0.1𝑟−0.1)cos(𝑡)+(1.088+0.1𝑟)sin(𝑡)−𝑡,𝑌(𝑡,𝑟)=(0.1−0.1𝑟)cos(𝑡)+(1.288−0.1𝑟)sin(𝑡)−𝑡.(35) According to (11), we may replace (34) by the following equivalent system: 𝑌^+𝑌=−𝑡,𝑡∈[0,1],𝑌𝑌^+𝑌=𝑡,𝑡∈[0,1],(0)=0.1𝑟−0.1,𝑌^(0)=0.088+0.1𝑟,𝑌(0)=0.1−0.1𝑟,𝑌^(0)=0.288−0.1𝑟.(36) We first construct the zero-order deformation equations (1−𝑞)ℒ_𝑖𝜙_𝑖(𝑡,𝑟;𝑞)−𝑧_𝑖,0(𝑡,𝑟)=𝑖𝑞=ℏ𝒩1,_𝑖2,𝜙_𝑖(𝑡,𝑟;𝑞),(37) subject to the initial conditions 𝑌₀=(0.1𝑟−0.1)+(0.088+0.1𝑟)𝑡,𝑌₀=(0.1−0.1𝑟)+(0.288−0.1𝑟)𝑡,(38) and the linear operator ℒ_𝑖𝜙_𝑖(𝑡,𝑟;𝑞)=𝜎²𝜙_𝑖(𝑡,𝑟;𝑞)𝜎𝑡²,𝑖=1,2,(39) with the property ℒ_𝑖𝑐_𝑖0+𝑡𝑐_𝑖1=0,𝑖=1,2,(40) where 𝑐_𝑖0, 𝑐_𝑖1(𝑖=1,2) are integral constants. Also from (36), we can define 𝒩₁𝜙₁(𝑡,𝑟;𝑞)=𝜎²𝜙₁(𝑡,𝑟;𝑞)𝜎𝑡²𝒩+𝜙₁(𝑡,𝑟;𝑞)+𝑡,₂𝜙₂(𝑡,𝑟;𝑞)=𝜎²𝜙₂(𝑡,𝑟;𝑞)𝜎𝑡²+𝜙₂(𝑡,𝑟;𝑞)+𝑡.(41) Obviously, when 𝑞=0 and 𝑞=1, 𝜙𝜙₁₂(𝑡,𝑟;0)=𝑧_1,0(𝑡,𝑟)=𝑌₀,𝜙₁(𝑡,𝑟;1)=𝑌,(𝑡,𝑟;0)=𝑧_2,0(𝑡,𝑟)=𝑌₀,𝜙₂(𝑡,𝑟;1)=𝑌.(42)
Therefore, when the embedding parameter 𝑞 increases from 0 to 1, the homotopy solutions 𝜙_𝑖(𝑡,𝑟,𝑞) vary from 𝑧_𝑖,0(𝑡,𝑟) to the solutions 𝑧_𝑖(𝑡,𝑟) for 𝑖=1,2. Now, by expanding 𝜙_𝑖(𝑡,𝑟;𝑞) in Taylor’s series with respect to 𝑞, we have 𝜙_𝑖(𝑡,𝑟;𝑞)=𝑧_𝑖,0(𝑡,𝑟)+^+∞_𝑚₌₁𝑧_𝑖,𝑚(𝑡,𝑟)𝑞^𝑚,𝑖=1,2,(43) where 𝑧_𝑖,𝑚(𝑡,𝑟)=1𝑚!𝜕^𝑚𝜙_𝑖(𝑡,𝑟;𝑞)𝜕𝑞^𝑚|_𝑞=0,𝑖=1,2.(44) Assuming that auxiliary parameter ℏ, the initial guesses and the auxiliary linear operator are properly chosen, then the above series is convergent at 𝑞=1, and 𝑌(𝑡,𝑟)=𝑧_1,0(𝑡,𝑟)+^∞_𝑚₌₁𝑧_1,𝑚(𝑡,𝑟),𝑌(𝑡,𝑟)=𝑧_2,0(𝑡,𝑟)+^∞_𝑚₌₁𝑧_2,𝑚(𝑡,𝑟).(45) The 𝑚th order deformation equations are ℒ_𝑖𝑧_𝑖,𝑚(𝑡,𝑟)−𝜒_𝑚𝑧_{𝑖,𝑚−1}(𝑡,𝑟)=ℏℛ_𝑖,𝑚^→𝑧_{𝑖,𝑚−1}𝑖=1,2,,(46) with the initial conditions 𝑧_𝑖,𝑚(0,𝑟)=0,𝑖=1,2,(47) where ℛ_𝑖,𝑚^→𝑧_{𝑖,𝑚−1}=𝜕²𝑧_{𝑖,𝑚−1}𝜕𝑡²+𝑧_{𝑖,𝑚−1}+1−𝜒_𝑚𝑡,𝑖=1,2.(48) Therefore, we recursively obtain 𝑧_1,1(𝑡,𝑟)=−ℏ𝑡²𝑧20+ℏ2𝑟0𝑡²+6387ℏ5𝑡³+ℏ6𝑟0_1,2𝑡³,(𝑡,𝑟)=−ℏ𝑡²20−−ℏ²𝑡²+ℏ20+ℏ2𝑟0₂𝑡²𝑟𝑡²20+6387ℏ5𝑡³+68ℏ²𝑡³+ℏ𝑟𝑡₃37560+ℏ²𝑟𝑡³60−ℏ²𝑡⁴+ℏ240₂𝑟𝑡⁴240+17ℏ²𝑡⁵1875+ℏ²𝑟𝑡⁵𝑧1200_2,1,(𝑡,𝑟)=ℏ𝑡²𝑧20−ℏ2𝑟0𝑡²+16715ℏ0𝑡³−ℏ6𝑟0_2,2𝑡³,(𝑡,𝑟)=ℏ𝑡²20+−ℏ²𝑡²−ℏ20−ℏ2𝑟0₂𝑡²𝑟𝑡²20+16715ℏ0𝑡³+161ℏ²𝑡³−ℏ𝑟𝑡₃75060−ℏ²𝑟𝑡³60+ℏ²𝑡⁴−ℏ240₂𝑟𝑡⁴240+161ℏ²𝑡⁵15000−ℏ²𝑟𝑡⁵⋮.1200,(49) Then the solutions obtained by HAM are as follows: 𝑌(𝑡,𝑟)=+𝑧𝑧_1,0(_1,2𝑡,𝑟)+𝑧_1,1(𝑡,𝑟)(𝑡,𝑟)+𝑧_1,3(𝑡,𝑟)+⋯,𝑌(𝑡,𝑟)=+𝑧𝑧_2,0(_2,2𝑡,𝑟)+𝑧_2,1(𝑡,𝑟)(𝑡,𝑟)+𝑧_2,3(𝑡,𝑟)+⋯.(50)
Figures 1 and 3 show the ℏ-mesh of 𝑌 and 𝑌 to get a proper interval for convergency. ℏ-mesh is a generalization of traditional Contour plots, and for connection between ℏ-meshes and Contour plots we plot Figures 2 and 4. Also to find the best quantity of ℏ that lies in the convergency interval, we use the residual error of norm 2 as follows: Res𝑌_𝑛=¹⁰𝑌_𝑛(𝑡,𝑟)−𝑌(𝑡,𝑟)²𝑑𝑡𝑑𝑟^1/2,(51) which is a function with respect to ℏ. Now, by minimizing the Res[𝑌₁₀] we obtain the best choice for auxiliary parameter to approximate of 𝑌₁₀(𝑡,𝑟) as follows: ℏ=−1.02698,(52) and in this case absolute error for the 10th order approximation by HAM for 𝑌_𝑛(𝑡,𝑟) is plotted in Figure 5.
By minimizing the residual error defined by Res𝑌_𝑛=¹⁰𝑌_𝑛(𝑡,𝑟)−𝑌(𝑡,𝑟)²𝑑𝑡𝑑𝑟^1/2,(53) it is clear that the best choice for auxiliary parameter to approximate 𝑌₁₀(𝑡,𝑟) is ℏ=−1.02683(54) which in this case absolute error for the 10th order approximation by HAM for 𝑌_𝑛(𝑡,𝑟) is plotted in Figure 6.

Figure 1 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 2 

Contour plot of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 3 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 4 

Contour plot of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 5 

Absolute error of 𝑌10(𝑡,𝑟) for ℏ=−1.02698.

Figure 6 

Absolute error of 𝑌10(𝑡,𝑟) for ℏ=−1.02683.

Example 8. Consider the system of fuzzy differential equation 𝑌^+𝑒^2𝑡𝐺^=4−𝑟4²𝑟+₃5𝑟−1𝑒^2𝑡,𝑌_−2𝑟𝑡+2𝑡+12𝑒^2𝑡,+𝐺^=2−𝑟2²𝑟+₃2𝑟𝑒^2𝑡,𝑌(0)=𝑟₂+6𝑒^2𝑡−2𝑟+2,𝑌_+𝑟,3−𝑟³,𝐺((00))==(02,𝑟0²)+,2𝑟,−𝐺2_𝑟³+6,(0)=(𝑟−1,0),𝑡∈[0,1].(55)
The exact solution is as follows: 𝐺𝑌((𝑡𝑡,,𝑟𝑟))==(𝑟𝑟²−+1)𝑟𝑡,𝑒(1^2𝑡,−3𝑟)−𝑡₂𝑟³𝑒^2𝑡,.(56)
We may replace (55) by the equivalent system 𝑌^+𝑒^2𝑡𝐺^=4𝑟²+5𝑟−1𝑒^2𝑡,𝑌^+𝑒𝑌^2𝑡𝐺^_=−4𝑟³−2𝑟𝑡+2𝑡+12𝑒^2𝑡,+𝐺^=2𝑟²+2𝑟𝑒^2𝑡,𝑌𝑌^+𝐺^=−2𝑟³+6𝑒^2𝑡−2𝑟+2,𝑌_(0,𝑟)=𝑟²+𝑟,𝑌(0,𝑟)=3−𝑟³,(0,𝑟)𝐺=2𝑟²+2𝑟,𝑌^(0,𝑟)=−2𝑟³+6,𝐺(0,𝑟)=0,𝐺(0,𝑟)=0,_(0,𝑟)=𝑟−1,𝐺^(0,𝑟)=0.(57) We select the following initial guesses: 𝑌₀=𝑟²+𝑟+𝐺2𝑟²+2𝑟𝑡,𝑌₀=3−𝑟³+−2𝑟³+6𝑡,₀=(𝑟−1)𝑡,𝐺₀=2−𝑟,(58) and the auxiliary linear operators ℒ_𝑖𝜙_𝑖(𝑡,𝑟;𝑞)=𝜕²𝜙_𝑖(𝑡,𝑟;𝑞)𝜕𝑡²,𝑖=1,2,3,4,(59) with the property ℒ_𝑖𝑐_𝑖0+𝑡𝑐_𝑖1=0,𝑖=1,2,3,4,(60) where 𝑐_𝑖0 and 𝑐_𝑖1(𝑖=1,2,3,4) are integral constants. Also the nonlinear operators are 𝒩=₁𝜕𝜙₁(𝑡,𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)₂𝜙₁(𝑡,𝑟;𝑞)𝜕𝑡²+𝑒^2𝑡𝜕𝜙₃(𝑡,𝑟;𝑞)𝜕𝑡−4𝑟²𝒩₂+5𝑟−1𝑒^2𝑡,=𝜕𝜙₁(𝑡,𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)₂𝜙₂(𝑡,𝑟;𝑞)𝜕𝑡²−𝜕𝜙⁴(𝑡,𝑟;𝑞)𝜕𝑡−−4𝑟³𝒩₃−2𝑟𝑡+2𝑡+12𝑒^2𝑡,=𝜕𝜙₁(𝑡,𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)₂𝜙₃(𝑡,𝑟;𝑞)𝜕𝑡²+𝜕𝜙¹(𝑡,𝑟;𝑞)𝜕𝑡−2𝑟²𝒩₄+2𝑟𝑒^2𝑡,=𝜕𝜙₁(𝑡,𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)₂𝜙₄(𝑡,𝑟;𝑞)𝜕𝑡²+𝜕𝜙²(𝑡,𝑟;𝑞)𝜕𝑡−−2𝑟³+6𝑒^2𝑡−2𝑟+2.(61) Then 𝑚th order deformation equations are ℒ_𝑖𝑧_𝑖,𝑚(𝑡,𝑟)−𝜒_𝑚𝑧_{𝑖,𝑚−1}(𝑡,𝑟)=ℏℛ_𝑖,𝑚^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1𝑖=1,2,3,4,,(62) with the initial conditions ℛ𝑧_𝑖,𝑚(0,𝑟)=0,𝑖=1,2,3,4,_1,𝑚^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕₂𝑧_1,𝑚−1𝜕𝑡²+𝑒^2𝑡𝜕𝑧_3,𝑚−1−1−𝜒_𝑚𝜕𝑡ℛ_2,𝑚4𝑟²+5𝑟−1𝑒^2𝑡,^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕₂𝑧_2,𝑚−1𝜕𝑡²−𝜕𝑧^4,𝑚−1−1−𝜒_𝑚𝜕𝑡ℛ_3,𝑚−4𝑟³−2𝑟𝑡+2𝑡+12𝑒^2𝑡,^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕₂𝑧_3,𝑚−1𝜕𝑡²+𝜕𝑧^1,𝑚−1−1−𝜒_𝑚𝜕𝑡ℛ_4,𝑚2𝑟²+2𝑟𝑒^2𝑡,^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕₂𝑧_4,𝑚−1𝜕𝑡²+𝜕𝑧^2,𝑚−1−1−𝜒_𝑚𝜕𝑡−2𝑟³+6𝑒^2𝑡−2𝑟+2.(63) Therefore, we recursively obtain 𝑧𝑧_1,1(𝑡,𝑟)=ℏ𝑟−𝑒^2𝑡ℏ𝑟+ℏ𝑟²−𝑒^2𝑡ℏ𝑟²+2ℏ𝑟𝑡+2ℏ𝑟²𝑡,_1,2(𝑡,𝑟)=ℏ𝑟−𝑒^2𝑡ℏ𝑟+21ℏ²𝑟16−54𝑒^2𝑡ℏ²𝑟−116𝑒^4𝑡+ℏ𝑟₂ℏ²𝑟−𝑒^2𝑡ℏ𝑟²+21ℏ²𝑟²−5164𝑒^2𝑡ℏ²𝑟²−116𝑒^4𝑡+9ℏ²𝑟²+2ℏ𝑟𝑡4ℏ²𝑟𝑡+12𝑒^2𝑡+9ℏ²𝑟𝑡+2ℏ𝑟²𝑡4ℏ²𝑟²𝑡+12𝑒^2𝑡𝑧ℏ²𝑟²𝑡,_2,1(𝑡,𝑟)=5ℏ2−52𝑒^2𝑡ℏ+ℏ𝑟2−12𝑒^2𝑡ℏ𝑟−ℏ𝑟³+𝑒^2𝑡ℏ𝑟³+11ℏ𝑡2−12𝑒^2𝑡+ℏ𝑟𝑡ℏ𝑡2+21𝑒^2𝑡𝑧ℏ𝑟𝑡−2ℏ𝑟³𝑡,_2,2(𝑡,𝑟)=5ℏ2−52𝑒^2𝑡+47ℏ₂ℏ16−141𝑒^2𝑡ℏ²−316𝑒^4𝑡ℏ²+ℏ𝑟−122𝑒^2𝑡ℏ𝑟+ℏ²𝑟−𝑒^2𝑡ℏ²𝑟−ℏ𝑟³+𝑒^2𝑡ℏ𝑟³−21ℏ²𝑟³16+45𝑒^2𝑡ℏ²𝑟³+116𝑒^4𝑡ℏ²𝑟³+11ℏ𝑡2−12𝑒^2𝑡ℏ𝑡+23ℏ²𝑡4+21𝑒^2𝑡+ℏ𝑟𝑡ℏ²𝑡2+12𝑒^2𝑡−2ℏ𝑟₃ℏ𝑟𝑡+ℏ²𝑟𝑡+𝑒^2𝑡ℏ²𝑟𝑡𝑡−94ℏ²𝑟³𝑡−12𝑒^2𝑡𝑧ℏ²𝑟³𝑡,_3,1(𝑡,𝑟)=ℏ𝑟2−12𝑒^2𝑡ℏ𝑟+ℏ𝑟²2−12𝑒^2𝑡+ℏ𝑟₂ℏ𝑟²+ℏ𝑟𝑡𝑧_3,2𝑡+ℏ𝑟𝑡²+ℏ𝑟²𝑡²,(𝑡,𝑟)=ℏ𝑟2−12𝑒^2𝑡+ℏ₂ℏ𝑟𝑟−𝑒^2𝑡ℏ²𝑟+ℏ𝑟²2−12𝑒^2𝑡ℏ𝑟²++ℏℏ𝑟²𝑟²₂−𝑒^2𝑡ℏ²𝑟²+ℏ𝑟𝑡+2ℏ²𝑟𝑡𝑡+2ℏ²𝑟²𝑡+ℏ𝑟𝑡²𝑧_4,1+2ℏ²𝑟𝑡²+ℏ𝑟²𝑡²+2ℏ²𝑟²𝑡²,(𝑡,𝑟)=3ℏ2−32𝑒^2𝑡ℏ−ℏ𝑟³2+12𝑒^2𝑡+2ℏ𝑡₂ℏ𝑟³+3ℏ𝑡−ℏ𝑟³𝑡𝑧_4,2+ℏ𝑟𝑡²−ℏ𝑟³𝑡²,(𝑡,𝑟)=3ℏ2−32𝑒^2𝑡ℏ+21ℏ²−2188𝑒^2𝑡ℏ²+3ℏ²𝑟8−83𝑒^2𝑡ℏ²𝑟−ℏ𝑟³+122𝑒^2𝑡+11ℏ₂ℏ𝑟³−ℏ²𝑟³+𝑒^2𝑡ℏ²𝑟³+3ℏ𝑡𝑡2−41𝑒^2𝑡ℏ²𝑡+12ℏ²𝑟𝑡+14𝑒^2𝑡−ℏ𝑟₃ℏ²𝑟𝑡𝑡−2ℏ²𝑟³𝑡+2ℏ𝑡²+19ℏ²𝑡²+ℏ𝑟𝑡₂4+54ℏ²⋮.𝑟𝑡²−ℏ𝑟³𝑡²−2ℏ²𝑟³𝑡²,(64) Then the solutions obtained by HAM are as follows: 𝑌(𝑡,𝑟)=+𝑧𝑧_1,01,2(𝑡,𝑟)+𝑧_1,1(𝑡,𝑟)(𝑡,𝑟)+𝑧_1,3(𝑡,𝑟)+⋯,𝑌(𝑡,𝑟)=+𝑧𝑧_2,02,2(𝑡,𝑟)+𝑧_2,1(𝑡,𝑟)𝐺(𝑡,𝑟)+𝑧_2,3(𝑡,𝑟)+⋯,(𝑡,𝑟)=+𝑧𝑧_3,03,2(𝑡,𝑟)+𝑧_3,1(𝑡,𝑟)(𝑡,𝑟)+𝑧_3,3(𝑡,𝑟)+⋯,𝐺(𝑡,𝑟)=+𝑧𝑧_4,04,2(𝑡,𝑟)+𝑧_4,1(𝑡,𝑟)(𝑡,𝑟)+𝑧_4,3(𝑡,𝑟)+⋯,(65) and Figures 7, 8, 9, and 10 show the ℏ-meshes of (𝜕³/𝜕𝑡³)𝑌₁₀(0,𝑟),  (𝜕³/𝜕𝑡³)𝑌₁₀(0,𝑟),  (𝜕/𝜕𝑡)𝐺₁₀(0.1,𝑟), and (𝜕/𝜕𝑡)𝐺₁₀(0.1,𝑟) respectively. Also, Figures 11, 12, 13, and 14 exhibit their related Contour plots.

Figure 7 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 8 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 9 

ℏ-mesh of (𝜕/𝜕𝑡)𝐺10(0.1,𝑟).

Figure 10 

ℏ-mesh of (𝜕/𝜕𝑡)𝐺10(0.1,𝑟).

Figure 11 

Contour plot of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 12 

Contour plot of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 13 

Contour plot of (𝜕/𝜕𝑡)𝐺10(0.1,𝑟).

Figure 14 

Contour plot of (𝜕/𝜕𝑡)𝐺10(0.1,𝑟).

By minimizing the residual error, we find that the best choices for auxiliary parameter to approximate of 𝑌₁₅(𝑡,𝑟), 𝑌₁₅(𝑡,𝑟), 𝐺₁₅(𝑡,𝑟), and 𝐺₁₅(𝑡,𝑟) are ℏℏ==−−00..997722966596,,ℏℏ==−−00..99772296654,,(66) respectively. Absolute errors for the 15th order approximation by HAM for 𝑌_𝑛(𝑡,𝑟) and 𝑌_𝑛(𝑡,𝑟), 𝐺_𝑛(𝑡,𝑟), 𝐺_𝑛(𝑡,𝑟) are plotted in Figures 15, 16, 17, and 18, respectively.

Figure 15 

Absolute error of 𝑌15(𝑡,𝑟) for ℏ=−0.972656.

Figure 16 

Absolute error of 𝑌15(𝑡,𝑟) for ℏ=−0.97265.

Figure 17 

Absolute error of 𝐺15(𝑡,𝑟) for ℏ=−0.972969.

Figure 18 

Absolute error of 𝐺15(𝑡,𝑟) for ℏ=−0.972964.

Example 9. Let us consider the first-order system of fuzzy differential equation 𝑌′𝑌−+𝐺𝐺=′=(−(2−𝑡𝑟𝑟𝑡++2𝑡𝑡++32𝑟𝑟++32))𝑒𝑒^2𝑡_2𝑡,(2𝑡𝑟−2𝑡−2𝑟+7)𝑒^2𝑡,𝑌(0)=(2+𝑟,4−𝑡∈𝑟)[,0,1]𝐺.,((𝑟0𝑡)−=𝑡(−𝑟,32𝑟−+𝑟8)),𝑒^2𝑡,(67) The exact solution is as follows: 𝑌(𝑡,𝑟)=(2+𝑟𝐺)𝑒(𝑡^2𝑡,+𝑟)(=1−𝑟𝑒𝑟)_2𝑡𝑡𝑒^2𝑡,(4−𝑟)𝑒^2𝑡+(𝑟−1)𝑡𝑒^2𝑡,,(2−𝑟)𝑒^2𝑡.(68) We may replace (67) by the equivalent system 𝑌^−𝐺=(−2𝑡𝑟+2𝑡+2𝑟+3)𝑒^2𝑡,𝑌𝑌^−𝐺=(2𝑡𝑟−2𝑡−2𝑟+7)𝑒^2𝑡,+𝐺^=(−𝑟𝑡+𝑡+3𝑟+2)𝑒^2𝑡,𝑌𝑌+𝐺^=(𝑟𝑡−𝑡−3𝑟+8)𝑒^2𝑡,𝐺(0,𝑟)=2+𝑟,𝑌(0,𝑟)=4−𝑟,(0,𝑟)=𝑟,𝐺(0,𝑟)=2−𝑟.(69) We choose the initial guesses of system (69) as follows: 𝐺𝑌₀=2+𝑟,𝑌₀=4−𝑟,₀=𝑟,𝐺₀=2−𝑟,(70) and the auxiliary linear operators ℒ_𝑖𝜙_𝑖(𝑡,𝑟;𝑞)=𝜕𝜙^𝑖(𝑡,𝑟;𝑞)𝜕𝑡,𝑖=1,2,3,4,(71) with the properties ℒ_𝑖𝑐_𝑖=0,𝑖=1,2,3,4,(72) where 𝑐_𝑖(𝑖=1,2,3,4) are integral constants, and the nonlinear operators, 𝒩₁=𝜙₁𝜕(𝜙𝑡,₁𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)(𝑡,𝑟;𝑞)𝜕𝑡−𝜙⁴−(−2𝑡𝑟+2𝑡+2𝑟(+𝑡,3𝑟;)𝑞𝑒)_2𝑡𝒩₂,=𝜙₁𝜕(𝜙𝑡,₂𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)(𝑡,𝑟;𝑞)𝜕𝑡−𝜙³−(2𝑡𝑟−2𝑡−2𝑟+(𝑡7,)𝑟𝑒;𝑞_2𝑡)𝒩₃,=𝜙₁𝜕(𝜙𝑡,₃𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)(𝑡,𝑟;𝑞)𝜕𝑡+𝜙¹−(−𝑟𝑡+𝑡+3𝑟+(2𝑡),𝑒𝑟;_2𝑡𝑞)𝒩₄,=𝜙₁𝜕(𝜙𝑡,₄𝑟;𝑞),…,𝜙₄(𝑡,𝑟;𝑞)(𝑡,𝑟;𝑞)𝜕𝑡+𝜙²−(𝑟𝑡−𝑡−3𝑟+8)(𝑡𝑒,𝑟_2𝑡;𝑞).(73) The 𝑚th order deformation equations are ℒ_𝑖𝑧_𝑖,𝑚(𝑡,𝑟)−𝜒_𝑚𝑧_{𝑖,𝑚−1}(𝑡,𝑟)=ℏℛ_𝑖,𝑚^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1𝑖=1,2,3,4,,(74) with the initial conditions 𝑧_𝑖,𝑚(0,𝑟)=0,𝑖=1,2,3,4,(75) where ℛ_1,𝑚^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕𝑧_1,𝑚−1𝜕𝑡−𝑧^4,𝑚−1ℛ_2,𝑚−1−𝜒_𝑚(−2𝑡𝑟+2𝑡+2𝑟+3)𝑒^2𝑡,^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕𝑧_2,𝑚−1𝜕𝑡−𝑧^3,𝑚−1ℛ_3,𝑚−1−𝜒_𝑚(2𝑡𝑟−2𝑡−2𝑟+7)𝑒^2𝑡,^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕𝑧_3,𝑚−1𝜕𝑡+𝑧^1,𝑚−1ℛ_4,𝑚−1−𝜒_𝑚(−𝑟𝑡+𝑡+3𝑟+2)𝑒^2𝑡,^→𝑧_1,𝑚−1,…,^→𝑧_4,𝑚−1=𝜕𝑧_4,𝑚−1𝜕𝑡+𝑧^2,𝑚−1−1−𝜒_𝑚(𝑟𝑡−𝑡−3𝑟+8)𝑒^2𝑡.(76) Therefore, we recursively obtain 𝑧_1,1(𝑡,𝑟)=ℏ−𝑒^2𝑡ℏ+3ℏ𝑟2−3𝑒^2𝑡ℏ𝑟−𝑒_2𝑡2−2ℏ𝑡𝑧ℏ𝑡+ℏ𝑟𝑡+𝑒^2𝑡ℏ𝑟𝑡,_1,2(𝑡,𝑟)=ℏ−𝑒^2𝑡ℏ−5ℏ²4+5𝑒^2𝑡ℏ²+3ℏ𝑟42−3𝑒^2𝑡ℏ𝑟2+5ℏ²𝑟2−5𝑒^2𝑡ℏ²𝑟2−2ℏ𝑡−𝑒^2𝑡−25ℏ₂ℏ𝑡𝑡4−5𝑒^2𝑡ℏ²𝑡4+ℏ𝑟𝑡+𝑒^2𝑡+11ℏ₂ℏ𝑟𝑡𝑟𝑡4+5𝑒^2𝑡ℏ²𝑟𝑡4−2ℏ²𝑡²+ℏ²𝑟𝑡²𝑧_2,12,(𝑡,𝑟)=4ℏ−4𝑒^2𝑡ℏ−3ℏ𝑟2+3𝑒^2𝑡ℏ𝑟2+𝑒^2𝑡−ℏ𝑟𝑡−𝑒_2𝑡ℏ𝑡𝑧ℏ𝑟𝑡,_2,2(𝑡,𝑟)=4ℏ−4𝑒^2𝑡ℏ+15ℏ²4−15𝑒^2𝑡ℏ²−3ℏ𝑟42+3𝑒^2𝑡ℏ𝑟−5ℏ₂2𝑟2+5𝑒^2𝑡ℏ²𝑟2+𝑒^2𝑡ℏ𝑡−3ℏ²𝑡+5𝑒_2𝑡4ℏ²𝑡4−ℏ𝑟𝑡−𝑒^2𝑡−11ℏ₂ℏ𝑟𝑡𝑟𝑡4−5𝑒^2𝑡ℏ²𝑟𝑡4−ℏ²𝑡²−ℏ²𝑟𝑡²𝑧_3,12,(𝑡,𝑟)=3ℏ4−3𝑒^2𝑡ℏ4+74ℏ𝑟−7𝑒^2𝑡ℏ𝑟−𝑒4+2ℏ𝑡_2𝑡ℏ𝑡2+ℏ𝑟𝑡+𝑒^2𝑡ℏ𝑟𝑡𝑧_3,22,(𝑡,𝑟)=3ℏ4−3𝑒^2𝑡ℏ4+ℏ²−𝑒^2𝑡ℏ²+7ℏ𝑟4−7𝑒^2𝑡ℏ𝑟+11ℏ₂4𝑟4−11𝑒^2𝑡ℏ²𝑟4+2ℏ𝑡−𝑒^2𝑡ℏ𝑡+3ℏ₂2𝑡−𝑒^2𝑡ℏ²𝑡+ℏ𝑟𝑡+𝑒^2𝑡ℏ𝑟𝑡+5ℏ₂2𝑟𝑡2+𝑒^2𝑡ℏ²𝑟𝑡−ℏ²𝑡²+ℏ²𝑟𝑡²𝑧_4,12,(𝑡,𝑟)=17ℏ4−17𝑒^2𝑡ℏ4−74ℏ𝑟+7𝑒^2𝑡ℏ𝑟+𝑒4+4ℎ𝑡_2𝑡ℏ𝑡2−ℏ𝑟𝑡−𝑒^2𝑡ℏ𝑟𝑡𝑧_4,22,(𝑡,𝑟)=17ℏ4−17𝑒^2𝑡ℏ+13ℏ₂42−13𝑒^2𝑡ℏ²2−74ℏ𝑟+7𝑒^2𝑡ℏ𝑟−11ℏ₂4𝑟4+11𝑒^2𝑡ℏ²𝑟4+4ℏ𝑡+𝑒^2𝑡ℏ𝑡+8ℏ₂2𝑡+𝑒^2𝑡ℏ²𝑡−ℏ𝑟𝑡−𝑒^2𝑡ℏ𝑟𝑡−5ℏ₂2𝑟𝑡2−𝑒^2𝑡ℏ²𝑟𝑡−ℏ²𝑟𝑡²⋮.2,(77) Then obtained solutions are as follows: 𝑌(𝑡,𝑟)=𝑧_1,0(𝑡,𝑟)+𝑧_1,1(𝑡,𝑟)+𝑧_1,2(𝑡,𝑟)+𝑧_1,3(𝑡,𝑟)+⋯,𝐺𝑌(𝑡,𝑟)=𝑧_2,0(𝑡,𝑟)+𝑧_2,1(𝑡,𝑟)+𝑧_2,2(𝑡,𝑟)+𝑧_2,3(𝑡,𝑟)+⋯,(𝑡,𝑟)=𝑧_3,0(𝑡,𝑟)+𝑧_3,1(𝑡,𝑟)+𝑧_3,2(𝑡,𝑟)+𝑧_3,3(𝑡,𝑟)+⋯,𝐺(𝑡,𝑟)=𝑧_4,0(𝑡,𝑟)+𝑧_4,1(𝑡,𝑟)+𝑧_4,2(𝑡,𝑟)+𝑧_4,3(𝑡,𝑟)+⋯.(78) Figures 19, 20, 21, and 22 show the ℏ-meshes of (𝜕³/𝜕𝑡³)𝑌₁₀(0,𝑟), (𝜕³/𝜕𝑡³)𝑌₁₀(0,𝑟), (𝜕³/𝜕𝑡³)𝐺₁₀(0,𝑟) and (𝜕³/𝜕𝑡³)𝐺₁₀(0,𝑟), respectively. Also, Figures 23, 24, 25, and 26 show the related Contour plots.

Figure 19 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 20 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 21 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝐺10(0,𝑟).

Figure 22 

ℏ-mesh of (𝜕3/𝜕𝑡3)𝐺10(0,𝑟).

Figure 23 

Contour plot of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 24 

Contour plot of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟).

Figure 25 

Contour plot of (𝜕3/𝜕𝑡3)𝐺10(0,𝑟).

Figure 26 

Contour plot of (𝜕3/𝜕𝑡3)𝐺10(0,𝑟).

Best choices for auxiliary parameter ℏ to approximate 𝑌₁₅(𝑡,𝑟),𝑌₁₅(𝑡,𝑟),𝐺₁₅(𝑡,𝑟), and 𝐺₁₅(𝑡,𝑟) from minimizing the residual errors are ℏℏ==−−01..90717755861,,ℏℏ==−−00.9.97688553271,,(79) respectively. Finally, absolute errors for the 15th order approximation by HAM for 𝑌_𝑛(𝑡,𝑟),𝑌_𝑛(𝑡,𝑟),𝐺_𝑛(𝑡,𝑟) and 𝐺_𝑛(𝑡,𝑟), are plotted in Figures 27, 28, 29, and 30.

Figure 27 

Absolute error of 𝑌15(𝑡,𝑟) for ℏ=−0.977581.

Figure 28 

Absolute error of 𝑌15(𝑡,𝑟) for ℏ=−0.968521.

Figure 29 

Absolute error of 𝐺15(𝑡,𝑟) for ℏ=−1.01756.

Figure 30 

Absolute error of 𝐺15(𝑡,𝑟) for ℏ=−0.978537.

7. Conclusion

In this paper, homotopy analysis method has been implemented to derive approximate analytical solutions for the system of fuzzy differential equations. Obtained results show that we can control the convergence region of HAM series solution by the auxiliary parameter ℏ. The concept of traditional ℎ-curves has been generalized to ℏ-meshes and then contour plot firstly has been introduced in HAM. Convergence theorem and given illustrative examples show the efficiency and accuracy of the HAM.

Acknowledgments

The respected anonymous referees have carefully reviewed this paper. As a result of their careful analysis, our paper has been improved. The authors would like to express their thankfulness to them for their helpful constructive comments.