Table of Contents Author Guidelines Submit a Manuscript
Advances in Fuzzy Systems
Volume 2012 (2012), Article ID 407647, 16 pages
Research Article

Series Solution of the System of Fuzzy Differential Equations

1Dipartimento di Mathematica e Informatica, Università degli Stueli di Perugia, 06123, Italy
2Department of Mathematics and Computer Science, Engineering and Technical Faculty, University of Bonab, Bonab 55517, Iran

Received 9 April 2012; Accepted 28 June 2012

Academic Editor: Rustom M. Mamlook

Copyright © 2012 M. S. Hashemi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 [1115]. 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 [1822].

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 𝔼1.

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

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

Definition 4. Let 𝑓[𝑎,𝑏]𝔼1, for each partition 𝑃={𝑡0,,𝑡𝑛} of [𝑎,𝑏] and for arbitrary 𝜉𝑖[𝑡𝑖1,𝑡𝑖], 1𝑖𝑛 suppose 𝑅𝑝=𝑛𝑖=1𝑓𝜉𝑖𝑡𝑖𝑡𝑖1,||𝑡Δ=max𝑖𝑡𝑖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 𝑓1𝔼1 be a fuzzy function and let 𝑡01. The derivative 𝑓(𝑡0) of 𝑓 at the point 𝑡0 is defined by 𝑓𝑡0=lim0𝑓𝑡0𝑡+𝑓0,(6) provided that this limit, taken with respect to the metric 𝐷, exists.

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

Clearly [𝑓(𝑡0+)𝑓(𝑡0)]/ may not be a fuzzy number for all . However, if it approaches 𝑓(𝑡0) (in 𝔹) and 𝑓(𝑡0) is also a fuzzy number (i.e., in 𝔼1), this number is the fuzzy derivative of 𝑓(𝑡) at 𝑡0. In this case, if 𝑓=(𝑓,𝑓), it can be easily shown that 𝑓𝑡0=𝑓𝑡0,𝑓𝑡0,(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 𝐴0(𝑡),𝐴1(𝑡),,𝐴𝑛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(𝑡,𝑟),𝑦2(𝑡,𝑟),,𝑦𝑠(𝑡,𝑟)𝑇,𝑦𝑖𝑦(𝑡,𝑟)=𝑙𝑖𝑦(𝑡,𝑟),𝑢𝑖𝑓(𝑡,𝑟),𝑖=1,,𝑠,𝐹=1(𝑡,𝑟),𝑓2(𝑡,𝑟),,𝑓𝑠(𝑡,𝑟)𝑇,𝑓𝑖(𝑓𝑡,𝑟)=𝑙𝑖(𝑓𝑡,𝑟),𝑢𝑖(𝐺𝑡,𝑟),𝑖=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 𝑧𝑍(𝑡,𝑟)=1(𝑡,𝑟),𝑧2(𝑡,𝑟),,𝑧𝑠(𝑡,𝑟)𝑇,𝑓𝐹(𝑡,𝑟)=1(𝑡,𝑟),𝑓2(𝑡,𝑟),,𝑓𝑠(𝑡,𝑟)𝑇,(15) and 𝐴0(𝑡),𝐴1(𝑡),,𝐴𝑛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: 𝒩𝑖𝑧1(𝑡,𝑟),𝑧2(𝑡,𝑟),,𝑧𝑠(𝑡,𝑟)=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(𝑡,𝑟;𝑞),,𝜙𝑛(,𝑡,𝑟;𝑞)𝑖=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(𝑡,𝑟)++𝑚=1𝑧𝑖,𝑚(𝑡,𝑟)𝑞𝑚,(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𝑧𝑖,𝑚(𝑡,𝑟),𝑖=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,𝜒𝑖=1,,𝑠,𝑚=0,𝑚1,1,𝑚>1.(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 𝐴0(𝑡),𝐴1(𝑡),,𝐴𝑛(𝑡) for some 𝑡 are negative, we must divide the defined interval for 𝑡 into small intervals, so that sign of each component of 𝐴0(𝑡),𝐴1(𝑡),,𝐴𝑛(𝑡) 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(𝑡,𝑟)=𝑘(𝑦)𝑘,1+(𝑡,𝑟)𝑘(𝑦)𝑘,2(𝑡,𝑟)𝑘(𝑦)𝑘,1(𝑡,𝑟)++𝑘(𝑦)𝑘,𝑛(𝑡,𝑟)𝑘(𝑦)𝑘,𝑛1(𝑡,𝑟)=𝑘(𝑦)𝑘,𝑛(𝑡,𝑟),𝑘=1,2,,𝑠,(28) then from (27) we can write 𝑚=1𝑘(𝑦)𝑘,𝑚(𝑡,𝑟)𝜒𝑚(𝑦)𝑘,𝑚1(𝑡,𝑟)=lim𝑛𝑘(𝑦)𝑘,𝑛(𝑡,𝑟)=𝑘lim𝑛(𝑦)𝑘,𝑛(𝑡,𝑟)=0,𝑘=1,2,,𝑠,(29) hence, 𝑚=1𝑘,𝑚𝑦𝑙1,𝑚1,,𝑦𝑙𝑠,𝑚1,𝑦𝑢1,𝑚1,,𝑦𝑢𝑠,𝑚1=𝑚=1𝑘(𝑦)𝑘,𝑚(𝑡,𝑟)𝜒𝑚(𝑦)𝑘,𝑚1(𝑡,𝑟)=0,𝑘=1,2,,𝑠.(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(𝜕𝑡,𝑟)=𝑖𝜕𝑡𝑖(𝑦)𝑘(𝑡,𝑟),𝑖=0,1,,𝑛,𝑘=1,2,,𝑠,(32) thus 0=𝑚=1𝑛𝑖=0𝐴𝑖(𝑡)𝑌(𝑖)(𝑡)1𝜒𝑚(𝐹)𝑑=(𝑡,𝑟)𝑛𝑖=0𝐴𝑖(𝑡)𝑚=1𝑌(𝑖)(𝑡)𝑚=11𝜒𝑚(𝐹)𝑑(=𝑡,𝑟)𝑛𝑖=0𝐴𝑖(𝑡)𝑚=1𝜕𝑖𝜕𝑡𝑖(𝑌)𝑘,𝑚1(𝑡,𝑟)(𝐹)𝑑=(𝑡,𝑟)𝑛𝑖=0𝐴𝑖𝜕(𝑡)𝑖𝜕𝑡𝑖(𝑌)𝑘(𝑡,𝑟)(𝐹)𝑑(𝑡,𝑟),𝑑=1,,𝑠.(33) Therefore, {𝑦𝑘(𝑡,𝑟)=((𝑦𝑙)𝑘(𝑡,𝑟),(𝑦𝑢)𝑘(𝑡,𝑟));𝑘=1,2,𝑠} is the exact solution of (11) and 𝑌=[𝑦1(𝑡,𝑟),𝑦2(𝑡,𝑟),,𝑦𝑠(𝑡,𝑟)]𝑇  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,1𝑌(0)=(0.1𝑟0.1,0.10.1𝑟),(0)=(0.088+0.1𝑟,0.2880.1𝑟).(34) The exact solution is as follows: 𝑌(𝑡,𝑟)=(0.1𝑟0.1)cos(𝑡)+(1.088+0.1𝑟)sin(𝑡)𝑡,𝑌(𝑡,𝑟)=(0.10.1𝑟)cos(𝑡)+(1.2880.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.10.1𝑟,𝑌(0)=0.2880.1𝑟.(36) We first construct the zero-order deformation equations (1𝑞)𝑖𝜙𝑖(𝑡,𝑟;𝑞)𝑧𝑖,0(𝑡,𝑟)=𝑞𝒩𝑖𝜙𝑖,(𝑡,𝑟;𝑞)𝑖=1,2,(37) subject to the initial conditions 𝑌0=(0.1𝑟0.1)+(0.088+0.1𝑟)𝑡,𝑌0=(0.10.1𝑟)+(0.2880.1𝑟)𝑡,(38) and the linear operator 𝑖𝜙𝑖=𝜎(𝑡,𝑟;𝑞)2𝜙𝑖(𝑡,𝑟;𝑞)𝜎𝑡2,𝑖=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 𝒩1𝜙1=𝜎(𝑡,𝑟;𝑞)2𝜙1(𝑡,𝑟;𝑞)𝜎𝑡2+𝜙1𝒩(𝑡,𝑟;𝑞)+𝑡,2𝜙2=𝜎(𝑡,𝑟;𝑞)2𝜙2(𝑡,𝑟;𝑞)𝜎𝑡2+𝜙2(𝑡,𝑟;𝑞)+𝑡.(41) Obviously, when 𝑞=0 and 𝑞=1, 𝜙1(𝑡,𝑟;0)=𝑧1,0(𝑡,𝑟)=𝑌0,𝜙1(𝑡,𝑟;1)=𝑌,𝜙2(𝑡,𝑟;0)=𝑧2,0(𝑡,𝑟)=𝑌0,𝜙2(𝑡,𝑟;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𝑧𝑖,𝑚(𝑡,𝑟)𝑞𝑚,𝑖=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𝑧1,𝑚(𝑡,𝑟),𝑌(𝑡,𝑟)=𝑧2,0(𝑡,𝑟)+𝑚=1𝑧2,𝑚(𝑡,𝑟).(45) The 𝑚th order deformation equations are 𝑖𝑧𝑖,𝑚(𝑡,𝑟)𝜒𝑚𝑧𝑖,𝑚1(𝑡,𝑟)=𝑖,𝑚𝑧𝑖,𝑚1,𝑖=1,2,(46) with the initial conditions 𝑧𝑖,𝑚(0,𝑟)=0,𝑖=1,2,(47) where 𝑖,𝑚𝑧𝑖,𝑚1=𝜕2𝑧𝑖,𝑚1𝜕𝑡2+𝑧𝑖,𝑚1+1𝜒𝑚𝑡,𝑖=1,2.(48) Therefore, we recursively obtain 𝑧1,1(𝑡,𝑟)=𝑡2+20𝑟𝑡2+2068𝑡3+375𝑟𝑡3,𝑧601,2(𝑡,𝑟)=𝑡2202𝑡2+20𝑟𝑡2+202𝑟𝑡2+2068𝑡3+375682𝑡3+375𝑟𝑡3+602𝑟𝑡3602𝑡4+2402𝑟𝑡4+240172𝑡5+18752𝑟𝑡5,𝑧12002,1(𝑡,𝑟)=𝑡220𝑟𝑡2+20161𝑡3750𝑟𝑡3,𝑧602,2(𝑡,𝑟)=𝑡2+202𝑡220𝑟𝑡2202𝑟𝑡2+20161𝑡3+7501612𝑡3750𝑟𝑡3602𝑟𝑡3+602𝑡42402𝑟𝑡4+2401612𝑡5150002𝑟𝑡5,1200.(49) Then the solutions obtained by HAM are as follows: 𝑌(𝑡,𝑟)=𝑧1,0(𝑡,𝑟)+𝑧1,1(𝑡,𝑟)+𝑧1,2(𝑡,𝑟)+𝑧1,3(𝑡,𝑟)+,𝑌(𝑡,𝑟)=𝑧2,0(𝑡,𝑟)+𝑧2,1(𝑡,𝑟)+𝑧2,2(𝑡,𝑟)+𝑧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𝑛=10𝑌𝑛(𝑡,𝑟)𝑌(𝑡,𝑟)2𝑑𝑡𝑑𝑟1/2,(51) which is a function with respect to . Now, by minimizing the Res[𝑌10] we obtain the best choice for auxiliary parameter to approximate of 𝑌10(𝑡,𝑟) 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𝑌𝑛=10𝑌𝑛(𝑡,𝑟)𝑌(𝑡,𝑟)2𝑑𝑡𝑑𝑟1/2,(53) it is clear that the best choice for auxiliary parameter to approximate 𝑌10(𝑡,𝑟) 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𝑟2𝑒+5𝑟12𝑡,4𝑟3𝑒2𝑟𝑡+2𝑡+122𝑡,𝑌+𝐺=2𝑟2𝑒+2𝑟2𝑡,2𝑟3𝑒+62𝑡,𝑟2𝑟+2𝑌(0)=2+𝑟,3𝑟3,𝑌(0)=2𝑟2+2𝑟,2𝑟3,+6𝐺(0)=(0,0),𝐺([].0)=(𝑟1,0),𝑡0,1(55)
The exact solution is as follows: 𝑟𝑌(𝑡,𝑟)=2𝑒+𝑟2𝑡,3𝑟3𝑒2𝑡,𝐺(𝑡,𝑟)=(𝑟1)𝑡,(1𝑟)𝑡2.(56)
We may replace (55) by the equivalent system 𝑌+𝑒2𝑡𝐺=4𝑟2𝑒+5𝑟12𝑡,𝑌+𝑒2𝑡𝐺=4𝑟3𝑒2𝑟𝑡+2𝑡+122𝑡,𝑌+𝐺=2𝑟2𝑒+2𝑟2𝑡,𝑌+𝐺=2𝑟3𝑒+62𝑡𝑌2𝑟+2,(0,𝑟)=𝑟2+𝑟,𝑌(0,𝑟)=3𝑟3,𝑌(0,𝑟)=2𝑟2+2𝑟,𝑌(0,𝑟)=2𝑟3𝐺+6,(0,𝑟)=0,𝐺𝐺(0,𝑟)=0,(0,𝑟)=𝑟1,𝐺(0,𝑟)=0.(57) We select the following initial guesses: 𝑌0=𝑟2+𝑟+2𝑟2+2𝑟𝑡,𝑌0=3𝑟3+2𝑟3𝐺+6𝑡,0=(𝑟1)𝑡,𝐺0=2𝑟,(58) and the auxiliary linear operators 𝑖𝜙𝑖=𝜕(𝑡,𝑟;𝑞)2𝜙𝑖(𝑡,𝑟;𝑞)𝜕𝑡2,𝑖=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 𝒩1𝜙1(𝑡,𝑟;𝑞),,𝜙4=𝜕(𝑡,𝑟;𝑞)2𝜙1(𝑡,𝑟;𝑞)𝜕𝑡2+𝑒2𝑡𝜕𝜙3(𝑡,𝑟;𝑞)𝜕𝑡4𝑟2𝑒+5𝑟12𝑡,𝒩2𝜙1(𝑡,𝑟;𝑞),,𝜙4=𝜕(𝑡,𝑟;𝑞)2𝜙2(𝑡,𝑟;𝑞)𝜕𝑡2𝜕𝜙4(𝑡,𝑟;𝑞)𝜕𝑡4𝑟3𝑒2𝑟𝑡+2𝑡+122𝑡,𝒩3𝜙1(𝑡,𝑟;𝑞),,𝜙4=𝜕(𝑡,𝑟;𝑞)2𝜙3(𝑡,𝑟;𝑞)𝜕𝑡2+𝜕𝜙1(𝑡,𝑟;𝑞)𝜕𝑡2𝑟2𝑒+2𝑟2𝑡,𝒩4𝜙1(𝑡,𝑟;𝑞),,𝜙4=𝜕(𝑡,𝑟;𝑞)2𝜙4(𝑡,𝑟;𝑞)𝜕𝑡2+𝜕𝜙2(𝑡,𝑟;𝑞)𝜕𝑡2𝑟3𝑒+62𝑡.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=𝜕2𝑧1,𝑚1𝜕𝑡2+𝑒2𝑡𝜕𝑧3,𝑚1𝜕𝑡1𝜒𝑚4𝑟2𝑒+5𝑟12𝑡,2,𝑚𝑧1,𝑚1,,𝑧4,𝑚1=𝜕2𝑧2,𝑚1𝜕𝑡2𝜕𝑧4,𝑚1𝜕𝑡1𝜒𝑚4𝑟3𝑒2𝑟𝑡+2𝑡+122𝑡,3,𝑚𝑧1,𝑚1,,𝑧4,𝑚1=𝜕2𝑧3,𝑚1𝜕𝑡2+𝜕𝑧1,𝑚1𝜕𝑡1𝜒𝑚2𝑟2𝑒+2𝑟2𝑡,4,𝑚𝑧1,𝑚1,,𝑧4,𝑚1=𝜕2𝑧4,𝑚1𝜕𝑡2+𝜕𝑧2,𝑚1𝜕𝑡1𝜒𝑚2𝑟3𝑒+62𝑡.2𝑟+2(63) Therefore, we recursively obtain 𝑧1,1(𝑡,𝑟)=𝑟𝑒2𝑡𝑟+𝑟2𝑒2𝑡𝑟2+2𝑟𝑡+2𝑟2𝑧𝑡,1,2(𝑡,𝑟)=𝑟𝑒2𝑡𝑟+212𝑟5164𝑒2𝑡21𝑟𝑒164𝑡2𝑟+𝑟2𝑒2𝑡𝑟2+212𝑟25164𝑒2𝑡2𝑟21𝑒164𝑡2𝑟2+9+2𝑟𝑡421𝑟𝑡+2𝑒2𝑡2𝑟𝑡+2𝑟2𝑡+942𝑟21𝑡+2𝑒2𝑡2𝑟2𝑧𝑡,2,1(𝑡,𝑟)=5252𝑒2𝑡+𝑟212𝑒2𝑡𝑟𝑟3+𝑒2𝑡𝑟3+11𝑡212𝑒2𝑡+𝑡𝑟𝑡2+12𝑒2𝑡𝑟𝑡2𝑟3𝑧𝑡,2,2(𝑡,𝑟)=5252𝑒2𝑡+47216114𝑒2𝑡23𝑒164𝑡2+𝑟212𝑒2𝑡𝑟+2𝑟𝑒2𝑡2𝑟𝑟3+𝑒2𝑡𝑟3212𝑟3+5164𝑒2𝑡2𝑟3+1𝑒164𝑡2𝑟3+11𝑡212𝑒2𝑡𝑡+232𝑡4+12𝑒2𝑡2𝑡+𝑟𝑡2+12𝑒2𝑡𝑟𝑡+2𝑟𝑡+𝑒2𝑡2𝑟𝑡2𝑟39𝑡42𝑟31𝑡2𝑒2𝑡2𝑟3𝑧𝑡,3,1(𝑡,𝑟)=𝑟212𝑒2𝑡𝑟+𝑟2212𝑒2𝑡𝑟2+𝑟𝑡+𝑟2𝑡+𝑟𝑡2+𝑟2𝑡2,𝑧3,2(𝑡,𝑟)=𝑟212𝑒2𝑡𝑟+2𝑟𝑒2𝑡2𝑟+𝑟2212𝑒2𝑡𝑟2+2𝑟2𝑒2𝑡2𝑟2+𝑟𝑡+22𝑟𝑡+𝑟2𝑡+22𝑟2𝑡+𝑟𝑡2+22𝑟𝑡2+𝑟2𝑡2+22𝑟2𝑡2,𝑧4,1(𝑡,𝑟)=3232𝑒2𝑡𝑟32+12𝑒2𝑡𝑟3+3𝑡𝑟3𝑡+2𝑡2+𝑟𝑡2𝑟3𝑡2,𝑧4,2(𝑡,𝑟)=3232𝑒2𝑡+2128218𝑒2𝑡2+32𝑟838𝑒2𝑡2𝑟𝑟32+12𝑒2𝑡𝑟32𝑟3+𝑒2𝑡2𝑟3++3𝑡112𝑡214𝑒2𝑡21𝑡+221𝑟𝑡+4𝑒2𝑡2𝑟𝑡𝑟3𝑡22𝑟3𝑡+2𝑡2+192𝑡24+𝑟𝑡2+542𝑟𝑡2𝑟3𝑡222𝑟3𝑡2,.(64) Then the solutions obtained by HAM 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(𝑡,𝑟)+,(65) and Figures 7, 8, 9, and 10 show the -meshes of (𝜕3/𝜕𝑡3)𝑌10(0,𝑟),  (𝜕3/𝜕𝑡3)𝑌10(0,𝑟),  (𝜕/𝜕𝑡)𝐺10(0.1,𝑟), and (𝜕/𝜕𝑡)𝐺10(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 𝑌15(𝑡,𝑟), 𝑌15(𝑡,𝑟), 𝐺15(𝑡,𝑟), and 𝐺15(𝑡,𝑟) are =0.972656,=0.97265,=0.972969,=0.972964,(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𝑡+2𝑟+3)𝑒2𝑡,(2𝑡𝑟2𝑡2𝑟+7)𝑒2𝑡,𝑌+𝐺=(𝑟𝑡+𝑡+3𝑟+2)𝑒2𝑡,(𝑟𝑡𝑡3𝑟+8)𝑒2𝑡,[].𝑌(0)=(2+𝑟,4𝑟),𝐺(0)=(𝑟,2𝑟),𝑡0,1(67) The exact solution is as follows: 𝑌(𝑡,𝑟)=(2+𝑟)𝑒2𝑡+(1𝑟)𝑡𝑒2𝑡,(4𝑟)𝑒2𝑡+(𝑟1)𝑡𝑒2𝑡,𝐺(𝑡,𝑟)=𝑟𝑒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: 𝑌0=2+𝑟,𝑌0𝐺=4𝑟,0=𝑟,𝐺0=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, 𝒩1𝜙1(𝑡,𝑟;𝑞),,𝜙4=(𝑡,𝑟;𝑞)𝜕𝜙1(𝑡,𝑟;𝑞)𝜕𝑡𝜙4(𝑡,𝑟;𝑞)(2𝑡𝑟+2𝑡+2𝑟+3)𝑒2𝑡,𝒩2𝜙1(𝑡,𝑟;𝑞),,𝜙4=(𝑡,𝑟;𝑞)𝜕𝜙2(𝑡,𝑟;𝑞)𝜕𝑡𝜙3(𝑡,𝑟;𝑞)(2𝑡𝑟2𝑡2𝑟+7)𝑒2𝑡,𝒩3𝜙1(𝑡,𝑟;𝑞),,𝜙4=(𝑡,𝑟;𝑞)𝜕𝜙3(𝑡,𝑟;𝑞)𝜕𝑡+𝜙1(𝑡,𝑟;𝑞)(𝑟𝑡+𝑡+3𝑟+2)𝑒2𝑡,𝒩4𝜙1(𝑡,𝑟;𝑞),,𝜙4=(𝑡,𝑟;𝑞)𝜕𝜙4(𝑡,𝑟;𝑞)𝜕𝑡+𝜙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,𝑚11𝜒𝑚(2𝑡𝑟+2𝑡+2𝑟+3)𝑒2𝑡,2,𝑚𝑧1,𝑚1,,𝑧4,𝑚1=𝜕𝑧2,𝑚1𝜕𝑡𝑧3,𝑚11𝜒𝑚(2𝑡𝑟2𝑡2𝑟+7)𝑒2𝑡,3,𝑚𝑧1,𝑚1,,𝑧4,𝑚1=𝜕𝑧3,𝑚1𝜕𝑡+𝑧1,𝑚11𝜒𝑚(𝑟𝑡+𝑡+3𝑟+2)𝑒2𝑡,4,𝑚𝑧1,𝑚1,,𝑧4,𝑚1=𝜕𝑧4,𝑚1𝜕𝑡+𝑧2,𝑚11𝜒𝑚(𝑟𝑡𝑡3𝑟+8)𝑒2𝑡.(76) Therefore, we recursively obtain 𝑧1,1(𝑡,𝑟)=𝑒2𝑡+3𝑟23𝑒2𝑡𝑟22𝑡𝑒2𝑡𝑡+𝑟𝑡+𝑒2𝑡𝑧𝑟𝑡,1,2(𝑡,𝑟)=𝑒2𝑡524+5𝑒2𝑡24+3𝑟23𝑒2𝑡𝑟2+52𝑟25𝑒2𝑡2𝑟22𝑡𝑒2𝑡𝑡252𝑡45𝑒2𝑡2𝑡4+𝑟𝑡+𝑒2𝑡+𝑟𝑡112𝑟𝑡4+5𝑒2𝑡2𝑟𝑡422𝑡2+2𝑟𝑡22,𝑧2,1(𝑡,𝑟)=44𝑒2𝑡3𝑟2+3𝑒2𝑡𝑟2+𝑒2𝑡𝑡𝑟𝑡𝑒2𝑡𝑧𝑟𝑡,2,2(𝑡,𝑟)=44𝑒2𝑡+152415𝑒2𝑡243𝑟2+3𝑒2𝑡𝑟252𝑟2+5𝑒2𝑡2𝑟2+𝑒2𝑡𝑡32𝑡4+5𝑒2𝑡2𝑡4𝑟𝑡𝑒2𝑡𝑟𝑡112𝑟𝑡45𝑒2𝑡2𝑟𝑡42𝑡22𝑟𝑡22,𝑧3,1(𝑡,𝑟)=343𝑒2𝑡4+7𝑟47𝑒2𝑡𝑟4𝑒+2𝑡2𝑡𝑡2𝑒+𝑟𝑡+2𝑡𝑟𝑡2,𝑧3,2(𝑡,𝑟)=343𝑒2𝑡4+2𝑒2𝑡2+7𝑟47𝑒2𝑡𝑟4+112𝑟411𝑒2𝑡2𝑟4𝑒+2𝑡2𝑡𝑡2+32𝑡𝑒2𝑡2𝑒𝑡+𝑟𝑡+2𝑡𝑟𝑡2+52𝑟𝑡2+𝑒2𝑡2𝑟𝑡2𝑡2+2𝑟𝑡22,𝑧4,1(𝑡,𝑟)=17417𝑒2𝑡47𝑟4+7𝑒2𝑡𝑟4+𝑒+4𝑡2𝑡𝑡2𝑒𝑟𝑡2𝑡𝑟𝑡2,𝑧4,2(𝑡,𝑟)=17417𝑒2𝑡4+132213𝑒2𝑡227𝑟4+7𝑒2𝑡𝑟4112𝑟4+11𝑒2𝑡2𝑟4𝑒+4𝑡+2𝑡𝑡2+82𝑡+𝑒2𝑡2𝑒𝑡𝑟𝑡2𝑡𝑟𝑡252𝑟𝑡2𝑒2𝑡2𝑟𝑡2𝑟𝑡22,.(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 (𝜕3/𝜕𝑡3)𝑌10(0,𝑟), (𝜕3/𝜕𝑡3)𝑌10(0,𝑟), (𝜕3/𝜕𝑡3)𝐺10(0,𝑟) and (𝜕3/𝜕𝑡3)𝐺10(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 𝑌15(𝑡,𝑟),𝑌15(𝑡,𝑟),𝐺15(𝑡,𝑟), and 𝐺15(𝑡,𝑟) from minimizing the residual errors are =0.977581,=0.968521,=1.01756,=0.978537,(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.


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.


  1. L. A. Zadeh, “Fuzzy sets,” Information and Control, vol. 8, no. 3, pp. 338–353, 1965. View at Google Scholar · View at Scopus
  2. S. L. Chang and L. A. Zadeh, “On fuzzy mapping and control,” IEEE Transactions on Systems, Man and Cybernetics, vol. 2, no. 1, pp. 30–34, 1972. View at Google Scholar · View at Scopus
  3. D. Dubois and H. Prade, “Towards fuzzy differential calculus—part 3: differentiation,” Fuzzy Sets and Systems, vol. 8, no. 3, pp. 225–233, 1982. View at Google Scholar · View at Scopus
  4. M. L. Puri and D. Ralescu, “Differentials of fuzzy functions,” Journal of Mathematical Analysis and Applications, vol. 91, pp. 321–325, 1983. View at Google Scholar
  5. R. Goetschel and W. Voxman, “Elementary fuzzy calculus,” Fuzzy Sets and Systems, vol. 18, no. 1, pp. 31–43, 1986. View at Google Scholar · View at Scopus
  6. O. Kaleva, “Fuzzy differential equations,” Fuzzy Sets and Systems, vol. 24, no. 3, pp. 301–317, 1987. View at Google Scholar · View at Scopus
  7. O. Kaleva, “The cauchy problem for fuzzy differential equations,” Fuzzy Sets and Systems, vol. 35, no. 3, pp. 389–396, 1990. View at Google Scholar · View at Scopus
  8. S. Seikkala, “On the fuzzy initial value problem,” Fuzzy Sets and Systems, vol. 24, no. 3, pp. 319–330, 1987. View at Google Scholar ·