About this Journal Submit a Manuscript Table of Contents
Advances in Fuzzy Systems
Volume 2012 (2012), Article ID 407647, 16 pages
doi:10.1155/2012/407647
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.

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 [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 ( 𝑢 + 𝑣 ) ( 𝑟 ) = 𝑢 ( 𝑟 ) + 𝑣 ( 𝑟 ) , ( 𝑢 + 𝑣 ) ( 𝑟 ) = 𝑢 ( 𝑟 ) + 𝑣 ( 𝑟 ) , 𝑘 𝑢 ( 𝑟 ) = 𝑘 𝑢 ( 𝑟 ) , 𝑘 𝑢 ( 𝑟 ) = 𝑘 𝑢 ( 𝑟 ) , i f 𝑘 0 , 𝑘 𝑢 ( 𝑟 ) = 𝑘 𝑢 ( 𝑟 ) , 𝑘 𝑢 ( 𝑟 ) = 𝑘 𝑢 ( 𝑟 ) , i f 𝑘 < 0 . ( 1 )

Definition 3. For arbitrary fuzzy numbers 𝑢 , 𝑣 𝔼 1 , we use the distance 𝐷 ( 𝑢 , 𝑣 ) = s u p 0 𝑟 1 | | m a x 𝑢 ( 𝑟 ) 𝑣 | | , | | 𝑢 ( 𝑟 ) ( 𝑟 ) 𝑣 | | , ( 𝑟 ) ( 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 , | | 𝑡 Δ = m a x 𝑖 𝑡 𝑖 1 | | . , 𝑖 = 1 , , 𝑛 ( 3 ) The definite integral of 𝑓 ( 𝑡 ) over [ 𝑎 , 𝑏 ] is 𝑏 𝑎 𝑓 ( 𝑡 ) 𝑑 𝑡 = l i m Δ 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 𝑡 0 1 . The derivative 𝑓 ( 𝑡 0 ) of 𝑓 at the point 𝑡 0 is defined by 𝑓 𝑡 0 = l i m 0 𝑓 𝑡 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 , , 𝑠 . ( 1 0 )

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 , , 𝑠 , ( 1 1 ) where 𝑟 [ 0 , 1 ] .

We define functions 𝑏 𝑖 𝑑 , 𝑘 ( 𝑡 ) and 𝑐 𝑖 𝑑 , 𝑘 ( 𝑡 ) ; 1 𝑑 , 𝑘 𝑠 ,   0 𝑖 𝑛 , in the following form: 𝑏 𝑖 𝑑 , 𝑘 ( 𝑡 ) = 0 , 𝑎 𝑖 𝑑 , 𝑘 𝑎 ( 𝑡 ) < 0 , 𝑖 𝑑 , 𝑘 ( 𝑡 ) , 𝑎 𝑖 𝑑 , 𝑘 ( 𝑐 𝑡 ) 0 , 𝑖 𝑑 , 𝑘 𝑎 ( 𝑡 ) = 𝑖 𝑑 , 𝑘 ( 𝑡 ) , 𝑎 𝑖 𝑑 , 𝑘 ( 𝑡 ) < 0 , 0 , 𝑎 𝑖 𝑑 , 𝑘 ( 𝑡 ) 0 . ( 1 2 ) Now (11) becomes as follows: 𝑛 𝑠 𝑖 = 0 𝑘 = 1 𝑏 𝑖 𝑑 , 𝑘 𝑦 ( 𝑡 ) 𝑙 𝑘 ( 𝑖 ) ( 𝑡 , 𝑟 ) + 𝑐 𝑖 𝑑 , 𝑘 𝑦 ( 𝑡 ) 𝑢 𝑘 ( 𝑖 ) = 𝑓 ( 𝑡 , 𝑟 ) 𝑙 𝑑 ( 𝑡 , 𝑟 ) , 𝑑 = 1 , , 𝑠 𝑛 𝑠 𝑖 = 0 𝑘 = 1 𝑐 𝑖 𝑑 , 𝑘 𝑦 ( 𝑡 ) 𝑙 𝑘 ( 𝑖 ) ( 𝑡 , 𝑟 ) + 𝑏 𝑖 𝑑 , 𝑘 𝑦 ( 𝑡 ) 𝑢 𝑘 ( 𝑖 ) = 𝑓 ( 𝑡 , 𝑟 ) 𝑢 𝑑 𝑦 ( 𝑡 , 𝑟 ) , 𝑑 = 1 , , 𝑠 𝑙 𝑗 ( 𝑖 ) 𝑔 ( 0 , 𝑟 ) = 𝑙 𝑖 𝑗 𝑦 ( 𝑟 ) , 𝑖 = 0 , , 𝑛 1 , 𝑗 = 1 , , 𝑠 , 𝑢 𝑗 ( 𝑖 ) 𝑔 ( 0 , 𝑟 ) = 𝑢 𝑖 𝑗 ( 𝑟 ) , 𝑖 = 0 , , 𝑛 1 , 𝑗 = 1 , , 𝑠 . ( 1 3 )

4. Basic Ideas of HAM

We consider the following differential equations: 𝑛 𝑖 = 0 𝐴 𝑖 ( 𝑡 ) 𝑍 ( 𝑖 ) ( 𝑡 , 𝑟 ) 𝐹 ( 𝑡 , 𝑟 ) = 0 , ( 1 4 ) where 𝑧 𝑍 ( 𝑡 , 𝑟 ) = 1 ( 𝑡 , 𝑟 ) , 𝑧 2 ( 𝑡 , 𝑟 ) , , 𝑧 𝑠 ( 𝑡 , 𝑟 ) 𝑇 , 𝑓 𝐹 ( 𝑡 , 𝑟 ) = 1 ( 𝑡 , 𝑟 ) , 𝑓 2 ( 𝑡 , 𝑟 ) , , 𝑓 𝑠 ( 𝑡 , 𝑟 ) 𝑇 , ( 1 5 ) 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 , , 𝑠 , ( 1 6 ) 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 , , 𝑛 , ( 1 7 ) 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 ) = 𝑧 𝑖 ( 𝑡 , 𝑟 ) ( 1 8 ) 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 𝑧 𝑖 , 𝑚 ( 𝑡 , 𝑟 ) 𝑞 𝑚 , ( 1 9 ) where 𝑧 𝑖 , 𝑚 = 1 𝜕 𝑚 ! 𝑚 𝜙 𝑖 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑞 𝑚 𝑞 = 0 . ( 2 0 ) 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 , , 𝑠 , ( 2 1 ) 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 , , 𝑠 . ( 2 2 ) 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 , , 𝑠 , ( 2 3 ) where 𝑖 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 𝑠 , 𝑚 1 = 1 ( 𝜕 𝑚 1 ) ! 𝑚 1 𝒩 𝑖 𝜙 1 ( 𝑡 , 𝑟 ; 𝑞 ) , , 𝜙 𝑠 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑞 𝑚 1 𝑞 = 0 , 𝜒 𝑖 = 1 , , 𝑠 , 𝑚 = 0 , 𝑚 1 , 1 , 𝑚 > 1 . ( 2 4 )

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 , ( 2 5 )

if 𝑘 th equation of (13) is inclusive 𝑏 𝑖 𝑑 , 𝑘 ( 𝑡 ) ( 𝑦 𝑙 ) 𝑘 ( 𝑖 ) ( 𝑡 , 𝑟 ) ;   1 𝑖 𝑛 ,  1 𝑑 𝑠 , and otherwise 𝑘 , 𝑚 = 1 𝜕 ( 𝑚 1 ) ! 𝑚 1 𝒩 𝑘 + 𝑠 𝑦 𝑙 1 , 𝑚 1 𝑦 , , 𝑙 𝑠 , 𝑚 1 , 𝑦 𝑢 1 , 𝑚 1 𝑦 , , 𝑢 𝑠 , 𝑚 1 𝜕 𝑞 𝑚 1 𝑞 = 0 . ( 2 6 )

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 l i m 𝑚 ( 𝑦 ) 𝑘 , 𝑚 ( 𝑡 , 𝑟 ) = 0 , 1 𝑘 𝑠 , 0 𝑡 𝑇 , 0 𝑟 1 . ( 2 7 ) On the other hand, since 𝑖 are linear operators, thus 𝑛 𝑚 = 1 𝑘 ( 𝑦 ) 𝑘 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 ( 𝑦 ) 𝑘 , 𝑚 1 = ( 𝑡 , 𝑟 ) 𝑛 𝑚 = 1 𝑘 ( 𝑦 ) 𝑘 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 𝑘 ( 𝑦 ) 𝑘 , 𝑚 1 ( 𝑡 , 𝑟 ) = 𝑘 ( 𝑦 ) 𝑘 , 1 + ( 𝑡 , 𝑟 ) 𝑘 ( 𝑦 ) 𝑘 , 2 ( 𝑡 , 𝑟 ) 𝑘 ( 𝑦 ) 𝑘 , 1 ( 𝑡 , 𝑟 ) + + 𝑘 ( 𝑦 ) 𝑘 , 𝑛 ( 𝑡 , 𝑟 ) 𝑘 ( 𝑦 ) 𝑘 , 𝑛 1 ( 𝑡 , 𝑟 ) = 𝑘 ( 𝑦 ) 𝑘 , 𝑛 ( 𝑡 , 𝑟 ) , 𝑘 = 1 , 2 , , 𝑠 , ( 2 8 ) then from (27) we can write 𝑚 = 1 𝑘 ( 𝑦 ) 𝑘 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 ( 𝑦 ) 𝑘 , 𝑚 1 ( 𝑡 , 𝑟 ) = l i m 𝑛 𝑘 ( 𝑦 ) 𝑘 , 𝑛 ( 𝑡 , 𝑟 ) = 𝑘 l i m 𝑛 ( 𝑦 ) 𝑘 , 𝑛 ( 𝑡 , 𝑟 ) = 0 , 𝑘 = 1 , 2 , , 𝑠 , ( 2 9 ) hence, 𝑚 = 1 𝑘 , 𝑚 𝑦 𝑙 1 , 𝑚 1 , , 𝑦 𝑙 𝑠 , 𝑚 1 , 𝑦 𝑢 1 , 𝑚 1 , , 𝑦 𝑢 𝑠 , 𝑚 1 = 𝑚 = 1 𝑘 ( 𝑦 ) 𝑘 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 ( 𝑦 ) 𝑘 , 𝑚 1 ( 𝑡 , 𝑟 ) = 0 , 𝑘 = 1 , 2 , , 𝑠 . ( 3 0 ) Since 0 , thus from the above equations for 𝑘 = 1 , 2 , , 𝑠 , we can write 𝑚 = 1 𝑘 , 𝑚 𝑦 𝑙 1 , 𝑚 1 , , 𝑦 𝑙 𝑠 , 𝑚 1 , 𝑦 𝑢 1 , 𝑚 1 , , 𝑦 𝑢 𝑠 , 𝑚 1 = 0 . ( 3 1 ) From uniform convergency we have 𝑚 = 1 𝜕 𝑖 𝜕 𝑡 𝑖 ( 𝑦 ) 𝑘 , 𝑚 1 ( 𝜕 𝑡 , 𝑟 ) = 𝑖 𝜕 𝑡 𝑖 ( 𝑦 ) 𝑘 ( 𝑡 , 𝑟 ) , 𝑖 = 0 , 1 , , 𝑛 , 𝑘 = 1 , 2 , , 𝑠 , ( 3 2 ) thus 0 = 𝑚 = 1 𝑛 𝑖 = 0 𝐴 𝑖 ( 𝑡 ) 𝑌 ( 𝑖 ) ( 𝑡 ) 1 𝜒 𝑚 ( 𝐹 ) 𝑑 = ( 𝑡 , 𝑟 ) 𝑛 𝑖 = 0 𝐴 𝑖 ( 𝑡 ) 𝑚 = 1 𝑌 ( 𝑖 ) ( 𝑡 ) 𝑚 = 1 1 𝜒 𝑚 ( 𝐹 ) 𝑑 ( = 𝑡 , 𝑟 ) 𝑛 𝑖 = 0 𝐴 𝑖 ( 𝑡 ) 𝑚 = 1 𝜕 𝑖 𝜕 𝑡 𝑖 ( 𝑌 ) 𝑘 , 𝑚 1 ( 𝑡 , 𝑟 ) ( 𝐹 ) 𝑑 = ( 𝑡 , 𝑟 ) 𝑛 𝑖 = 0 𝐴 𝑖 𝜕 ( 𝑡 ) 𝑖 𝜕 𝑡 𝑖 ( 𝑌 ) 𝑘 ( 𝑡 , 𝑟 ) ( 𝐹 ) 𝑑 ( 𝑡 , 𝑟 ) , 𝑑 = 1 , , 𝑠 . ( 3 3 ) 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 . 1 0 . 1 𝑟 ) , ( 0 ) = ( 0 . 0 8 8 + 0 . 1 𝑟 , 0 . 2 8 8 0 . 1 𝑟 ) . ( 3 4 ) The exact solution is as follows: 𝑌 ( 𝑡 , 𝑟 ) = ( 0 . 1 𝑟 0 . 1 ) c o s ( 𝑡 ) + ( 1 . 0 8 8 + 0 . 1 𝑟 ) s i n ( 𝑡 ) 𝑡 , 𝑌 ( 𝑡 , 𝑟 ) = ( 0 . 1 0 . 1 𝑟 ) c o s ( 𝑡 ) + ( 1 . 2 8 8 0 . 1 𝑟 ) s i n ( 𝑡 ) 𝑡 . ( 3 5 ) According to (11), we may replace (34) by the following equivalent system: 𝑌 + 𝑌 [ ] , = 𝑡 , 𝑡 0 , 1 𝑌 + [ ] , 𝑌 𝑌 = 𝑡 , 𝑡 0 , 1 ( 0 ) = 0 . 1 𝑟 0 . 1 , 𝑌 ( 0 ) = 0 . 0 8 8 + 0 . 1 𝑟 , 𝑌 ( 0 ) = 0 . 1 0 . 1 𝑟 , 𝑌 ( 0 ) = 0 . 2 8 8 0 . 1 𝑟 . ( 3 6 ) We first construct the zero-order deformation equations ( 1 𝑞 ) 𝑖 𝜙 𝑖 ( 𝑡 , 𝑟 ; 𝑞 ) 𝑧 𝑖 , 0 ( 𝑡 , 𝑟 ) = 𝑞 𝒩 𝑖 𝜙 𝑖 , ( 𝑡 , 𝑟 ; 𝑞 ) 𝑖 = 1 , 2 , ( 3 7 ) subject to the initial conditions 𝑌 0 = ( 0 . 1 𝑟 0 . 1 ) + ( 0 . 0 8 8 + 0 . 1 𝑟 ) 𝑡 , 𝑌 0 = ( 0 . 1 0 . 1 𝑟 ) + ( 0 . 2 8 8 0 . 1 𝑟 ) 𝑡 , ( 3 8 ) and the linear operator 𝑖 𝜙 𝑖 = 𝜎 ( 𝑡 , 𝑟 ; 𝑞 ) 2 𝜙 𝑖 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜎 𝑡 2 , 𝑖 = 1 , 2 , ( 3 9 ) with the property 𝑖 𝑐 𝑖 0 + 𝑡 𝑐 𝑖 1 = 0 , 𝑖 = 1 , 2 , ( 4 0 ) 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 ( 𝑡 , 𝑟 ; 𝑞 ) + 𝑡 . ( 4 1 ) Obviously, when 𝑞 = 0 and 𝑞 = 1 , 𝜙 1 ( 𝑡 , 𝑟 ; 0 ) = 𝑧 1 , 0 ( 𝑡 , 𝑟 ) = 𝑌 0 , 𝜙 1 ( 𝑡 , 𝑟 ; 1 ) = 𝑌 , 𝜙 2 ( 𝑡 , 𝑟 ; 0 ) = 𝑧 2 , 0 ( 𝑡 , 𝑟 ) = 𝑌 0 , 𝜙 2 ( 𝑡 , 𝑟 ; 1 ) = 𝑌 . ( 4 2 )
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 , ( 4 3 ) where 𝑧 𝑖 , 𝑚 1 ( 𝑡 , 𝑟 ) = 𝜕 𝑚 ! 𝑚 𝜙 𝑖 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑞 𝑚 | 𝑞 = 0 , 𝑖 = 1 , 2 . ( 4 4 ) 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 , 𝑚 ( 𝑡 , 𝑟 ) . ( 4 5 ) The 𝑚 th order deformation equations are 𝑖 𝑧 𝑖 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 𝑧 𝑖 , 𝑚 1 ( 𝑡 , 𝑟 ) = 𝑖 , 𝑚 𝑧 𝑖 , 𝑚 1 , 𝑖 = 1 , 2 , ( 4 6 ) with the initial conditions 𝑧 𝑖 , 𝑚 ( 0 , 𝑟 ) = 0 , 𝑖 = 1 , 2 , ( 4 7 ) where 𝑖 , 𝑚 𝑧 𝑖 , 𝑚 1 = 𝜕 2 𝑧 𝑖 , 𝑚 1 𝜕 𝑡 2 + 𝑧 𝑖 , 𝑚 1 + 1 𝜒 𝑚 𝑡 , 𝑖 = 1 , 2 . ( 4 8 ) Therefore, we recursively obtain 𝑧 1 , 1 ( 𝑡 , 𝑟 ) = 𝑡 2 + 2 0 𝑟 𝑡 2 + 2 0 6 8 𝑡 3 + 3 7 5 𝑟 𝑡 3 , 𝑧 6 0 1 , 2 ( 𝑡 , 𝑟 ) = 𝑡 2 2 0 2 𝑡 2 + 2 0 𝑟 𝑡 2 + 2 0 2 𝑟 𝑡 2 + 2 0 6 8 𝑡 3 + 3 7 5 6 8 2 𝑡 3 + 3 7 5 𝑟 𝑡 3 + 6 0 2 𝑟 𝑡 3 6 0 2 𝑡 4 + 2 4 0 2 𝑟 𝑡 4 + 2 4 0 1 7 2 𝑡 5 + 1 8 7 5 2 𝑟 𝑡 5 , 𝑧 1 2 0 0 2 , 1 ( 𝑡 , 𝑟 ) = 𝑡 2 2 0 𝑟 𝑡 2 + 2 0 1 6 1 𝑡 3 7 5 0 𝑟 𝑡 3 , 𝑧 6 0 2 , 2 ( 𝑡 , 𝑟 ) = 𝑡 2 + 2 0 2 𝑡 2 2 0 𝑟 𝑡 2 2 0 2 𝑟 𝑡 2 + 2 0 1 6 1 𝑡 3 + 7 5 0 1 6 1 2 𝑡 3 7 5 0 𝑟 𝑡 3 6 0 2 𝑟 𝑡 3 + 6 0 2 𝑡 4 2 4 0 2 𝑟 𝑡 4 + 2 4 0 1 6 1 2 𝑡 5 1 5 0 0 0 2 𝑟 𝑡 5 , 1 2 0 0 . ( 4 9 ) 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 ( 𝑡 , 𝑟 ) + . ( 5 0 )
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: 𝑌 R e s 𝑛 = 1 0 𝑌 𝑛 ( 𝑡 , 𝑟 ) 𝑌 ( 𝑡 , 𝑟 ) 2 𝑑 𝑡 𝑑 𝑟 1 / 2 , ( 5 1 ) which is a function with respect to . Now, by minimizing the R e s [ 𝑌 1 0 ] we obtain the best choice for auxiliary parameter to approximate of 𝑌 1 0 ( 𝑡 , 𝑟 ) as follows: = 1 . 0 2 6 9 8 , ( 5 2 ) 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 R e s 𝑌 𝑛 = 1 0 𝑌 𝑛 ( 𝑡 , 𝑟 ) 𝑌 ( 𝑡 , 𝑟 ) 2 𝑑 𝑡 𝑑 𝑟 1 / 2 , ( 5 3 ) it is clear that the best choice for auxiliary parameter to approximate 𝑌 1 0 ( 𝑡 , 𝑟 ) is = 1 . 0 2 6 8 3 ( 5 4 ) which in this case absolute error for the 10th order approximation by HAM for 𝑌 𝑛 ( 𝑡 , 𝑟 ) is plotted in Figure 6.

407647.fig.001
Figure 1: -mesh of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.002
Figure 2: Contour plot of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.003
Figure 3: -mesh of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.004
Figure 4: Contour plot of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.005
Figure 5: Absolute error of 𝑌 1 0 ( 𝑡 , 𝑟 ) for = 1 . 0 2 6 9 8 .
407647.fig.006
Figure 6: Absolute error of 𝑌 1 0 ( 𝑡 , 𝑟 ) for = 1 . 0 2 6 8 3 .

Example 8. Consider the system of fuzzy differential equation 𝑌 + 𝑒 2 𝑡 𝐺 = 4 𝑟 2 𝑒 + 5 𝑟 1 2 𝑡 , 4 𝑟 3 𝑒 2 𝑟 𝑡 + 2 𝑡 + 1 2 2 𝑡 , 𝑌 + 𝐺 = 2 𝑟 2 𝑒 + 2 𝑟 2 𝑡 , 2 𝑟 3 𝑒 + 6 2 𝑡 , 𝑟 2 𝑟 + 2 𝑌 ( 0 ) = 2 + 𝑟 , 3 𝑟 3 , 𝑌 ( 0 ) = 2 𝑟 2 + 2 𝑟 , 2 𝑟 3 , + 6 𝐺 ( 0 ) = ( 0 , 0 ) , 𝐺 ( [ ] . 0 ) = ( 𝑟 1 , 0 ) , 𝑡 0 , 1 ( 5 5 )
The exact solution is as follows: 𝑟 𝑌 ( 𝑡 , 𝑟 ) = 2 𝑒 + 𝑟 2 𝑡 , 3 𝑟 3 𝑒 2 𝑡 , 𝐺 ( 𝑡 , 𝑟 ) = ( 𝑟 1 ) 𝑡 , ( 1 𝑟 ) 𝑡 2 . ( 5 6 )
We may replace (55) by the equivalent system 𝑌 + 𝑒 2 𝑡 𝐺 = 4 𝑟 2 𝑒 + 5 𝑟 1 2 𝑡 , 𝑌 + 𝑒 2 𝑡 𝐺 = 4 𝑟 3 𝑒 2 𝑟 𝑡 + 2 𝑡 + 1 2 2 𝑡 , 𝑌 + 𝐺 = 2 𝑟 2 𝑒 + 2 𝑟 2 𝑡 , 𝑌 + 𝐺 = 2 𝑟 3 𝑒 + 6 2 𝑡 𝑌 2 𝑟 + 2 , ( 0 , 𝑟 ) = 𝑟 2 + 𝑟 , 𝑌 ( 0 , 𝑟 ) = 3 𝑟 3 , 𝑌 ( 0 , 𝑟 ) = 2 𝑟 2 + 2 𝑟 , 𝑌 ( 0 , 𝑟 ) = 2 𝑟 3 𝐺 + 6 , ( 0 , 𝑟 ) = 0 , 𝐺 𝐺 ( 0 , 𝑟 ) = 0 , ( 0 , 𝑟 ) = 𝑟 1 , 𝐺 ( 0 , 𝑟 ) = 0 . ( 5 7 ) We select the following initial guesses: 𝑌 0 = 𝑟 2 + 𝑟 + 2 𝑟 2 + 2 𝑟 𝑡 , 𝑌 0 = 3 𝑟 3 + 2 𝑟 3 𝐺 + 6 𝑡 , 0 = ( 𝑟 1 ) 𝑡 , 𝐺 0 = 2 𝑟 , ( 5 8 ) and the auxiliary linear operators 𝑖 𝜙 𝑖 = 𝜕 ( 𝑡 , 𝑟 ; 𝑞 ) 2 𝜙 𝑖 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 2 , 𝑖 = 1 , 2 , 3 , 4 , ( 5 9 ) with the property 𝑖 𝑐 𝑖 0 + 𝑡 𝑐 𝑖 1 = 0 , 𝑖 = 1 , 2 , 3 , 4 , ( 6 0 ) 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 𝑟 1 2 𝑡 , 𝒩 2 𝜙 1 ( 𝑡 , 𝑟 ; 𝑞 ) , , 𝜙 4 = 𝜕 ( 𝑡 , 𝑟 ; 𝑞 ) 2 𝜙 2 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 2 𝜕 𝜙 4 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 4 𝑟 3 𝑒 2 𝑟 𝑡 + 2 𝑡 + 1 2 2 𝑡 , 𝒩 3 𝜙 1 ( 𝑡 , 𝑟 ; 𝑞 ) , , 𝜙 4 = 𝜕 ( 𝑡 , 𝑟 ; 𝑞 ) 2 𝜙 3 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 2 + 𝜕 𝜙 1 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 2 𝑟 2 𝑒 + 2 𝑟 2 𝑡 , 𝒩 4 𝜙 1 ( 𝑡 , 𝑟 ; 𝑞 ) , , 𝜙 4 = 𝜕 ( 𝑡 , 𝑟 ; 𝑞 ) 2 𝜙 4 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 2 + 𝜕 𝜙 2 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 2 𝑟 3 𝑒 + 6 2 𝑡 . 2 𝑟 + 2 ( 6 1 ) Then 𝑚 th order deformation equations are 𝑖 𝑧 𝑖 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 𝑧 𝑖 , 𝑚 1 ( 𝑡 , 𝑟 ) = 𝑖 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 , 𝑖 = 1 , 2 , 3 , 4 , ( 6 2 ) 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 𝑟 1 2 𝑡 , 2 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 = 𝜕 2 𝑧 2 , 𝑚 1 𝜕 𝑡 2 𝜕 𝑧 4 , 𝑚 1 𝜕 𝑡 1 𝜒 𝑚 4 𝑟 3 𝑒 2 𝑟 𝑡 + 2 𝑡 + 1 2 2 𝑡 , 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 𝑒 + 6 2 𝑡 . 2 𝑟 + 2 ( 6 3 ) Therefore, we recursively obtain 𝑧 1 , 1 ( 𝑡 , 𝑟 ) = 𝑟 𝑒 2 𝑡 𝑟 + 𝑟 2 𝑒 2 𝑡 𝑟 2 + 2 𝑟 𝑡 + 2 𝑟 2 𝑧 𝑡 , 1 , 2 ( 𝑡 , 𝑟 ) = 𝑟 𝑒 2 𝑡 𝑟 + 2 1 2 𝑟 5 1 6 4 𝑒 2 𝑡 2 1 𝑟 𝑒 1 6 4 𝑡 2 𝑟 + 𝑟 2 𝑒 2 𝑡 𝑟 2 + 2 1 2 𝑟 2 5 1 6 4 𝑒 2 𝑡 2 𝑟 2 1 𝑒 1 6 4 𝑡 2 𝑟 2 + 9 + 2 𝑟 𝑡 4 2 1 𝑟 𝑡 + 2 𝑒 2 𝑡 2 𝑟 𝑡 + 2 𝑟 2 𝑡 + 9 4 2 𝑟 2 1 𝑡 + 2 𝑒 2 𝑡 2 𝑟 2 𝑧 𝑡 , 2 , 1 ( 𝑡 , 𝑟 ) = 5 2 5 2 𝑒 2 𝑡 + 𝑟 2 1 2 𝑒 2 𝑡 𝑟 𝑟 3 + 𝑒 2 𝑡 𝑟 3 + 1 1 𝑡 2 1 2 𝑒 2 𝑡 + 𝑡 𝑟 𝑡 2 + 1 2 𝑒 2 𝑡 𝑟 𝑡 2 𝑟 3 𝑧 𝑡 , 2 , 2 ( 𝑡 , 𝑟 ) = 5 2 5 2 𝑒 2 𝑡 + 4 7 2 1 6 1 1 4 𝑒 2 𝑡 2 3 𝑒 1 6 4 𝑡 2 + 𝑟 2 1 2 𝑒 2 𝑡 𝑟 + 2 𝑟 𝑒 2 𝑡 2 𝑟 𝑟 3 + 𝑒 2 𝑡 𝑟 3 2 1 2 𝑟 3 + 5 1 6 4 𝑒 2 𝑡 2 𝑟 3 + 1 𝑒 1 6 4 𝑡 2 𝑟 3 + 1 1 𝑡 2 1 2 𝑒 2 𝑡 𝑡 + 2 3 2 𝑡 4 + 1 2 𝑒 2 𝑡 2 𝑡 + 𝑟 𝑡 2 + 1 2 𝑒 2 𝑡 𝑟 𝑡 + 2 𝑟 𝑡 + 𝑒 2 𝑡 2 𝑟 𝑡 2 𝑟 3 9 𝑡 4 2 𝑟 3 1 𝑡 2 𝑒 2 𝑡 2 𝑟 3 𝑧 𝑡 , 3 , 1 ( 𝑡 , 𝑟 ) = 𝑟 2 1 2 𝑒 2 𝑡 𝑟 + 𝑟 2 2 1 2 𝑒 2 𝑡 𝑟 2 + 𝑟 𝑡 + 𝑟 2 𝑡 + 𝑟 𝑡 2 + 𝑟 2 𝑡 2 , 𝑧 3 , 2 ( 𝑡 , 𝑟 ) = 𝑟 2 1 2 𝑒 2 𝑡 𝑟 + 2 𝑟 𝑒 2 𝑡 2 𝑟 + 𝑟 2 2 1 2 𝑒 2 𝑡 𝑟 2 + 2 𝑟 2 𝑒 2 𝑡 2 𝑟 2 + 𝑟 𝑡 + 2 2 𝑟 𝑡 + 𝑟 2 𝑡 + 2 2 𝑟 2 𝑡 + 𝑟 𝑡 2 + 2 2 𝑟 𝑡 2 + 𝑟 2 𝑡 2 + 2 2 𝑟 2 𝑡 2 , 𝑧 4 , 1 ( 𝑡 , 𝑟 ) = 3 2 3 2 𝑒 2 𝑡 𝑟 3 2 + 1 2 𝑒 2 𝑡 𝑟 3 + 3 𝑡 𝑟 3 𝑡 + 2 𝑡 2 + 𝑟 𝑡 2 𝑟 3 𝑡 2 , 𝑧 4 , 2 ( 𝑡 , 𝑟 ) = 3 2 3 2 𝑒 2 𝑡 + 2 1 2 8 2 1 8 𝑒 2 𝑡 2 + 3 2 𝑟 8 3 8 𝑒 2 𝑡 2 𝑟 𝑟 3 2 + 1 2 𝑒 2 𝑡 𝑟 3 2 𝑟 3 + 𝑒 2 𝑡 2 𝑟 3 + + 3 𝑡 1 1 2 𝑡 2 1 4 𝑒 2 𝑡 2 1 𝑡 + 2 2 1 𝑟 𝑡 + 4 𝑒 2 𝑡 2 𝑟 𝑡 𝑟 3 𝑡 2 2 𝑟 3 𝑡 + 2 𝑡 2 + 1 9 2 𝑡 2 4 + 𝑟 𝑡 2 + 5 4 2 𝑟 𝑡 2 𝑟 3 𝑡 2 2 2 𝑟 3 𝑡 2 , . ( 6 4 ) 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 ( 𝑡 , 𝑟 ) + , ( 6 5 ) and Figures 7, 8, 9, and 10 show the -meshes of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) ,   ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) ,   ( 𝜕 / 𝜕 𝑡 ) 𝐺 1 0 ( 0 . 1 , 𝑟 ) , and ( 𝜕 / 𝜕 𝑡 ) 𝐺 1 0 ( 0 . 1 , 𝑟 ) respectively. Also, Figures 11, 12, 13, and 14 exhibit their related Contour plots.

407647.fig.007
Figure 7: -mesh of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.008
Figure 8: -mesh of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.009
Figure 9: -mesh of ( 𝜕 / 𝜕 𝑡 ) 𝐺 1 0 ( 0 . 1 , 𝑟 ) .
407647.fig.0010
Figure 10: -mesh of ( 𝜕 / 𝜕 𝑡 ) 𝐺 1 0 ( 0 . 1 , 𝑟 ) .
407647.fig.0011
Figure 11: Contour plot of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.0012
Figure 12: Contour plot of ( 𝜕 3 / 𝜕 𝑡 3 ) 𝑌 1 0 ( 0 , 𝑟 ) .
407647.fig.0013
Figure 13: Contour plot of ( 𝜕 / 𝜕 𝑡 ) 𝐺 1 0 ( 0 . 1 , 𝑟 ) .
407647.fig.0014
Figure 14: Contour plot of ( 𝜕 / 𝜕 𝑡 ) 𝐺 1 0 ( 0 . 1 , 𝑟 ) .

By minimizing the residual error, we find that the best choices for auxiliary parameter to approximate of 𝑌 1 5 ( 𝑡 , 𝑟 ) , 𝑌 1 5 ( 𝑡 , 𝑟 ) , 𝐺 1 5 ( 𝑡 , 𝑟 ) , and 𝐺 1 5 ( 𝑡 , 𝑟 ) are = 0 . 9 7 2 6 5 6 , = 0 . 9 7 2 6 5 , = 0 . 9 7 2 9 6 9 , = 0 . 9 7 2 9 6 4 , ( 6 6 ) respectively. Absolute errors for the 15th order approximation by HAM for 𝑌 𝑛 ( 𝑡 , 𝑟 ) and 𝑌 𝑛 ( 𝑡 , 𝑟 ) , 𝐺 𝑛 ( 𝑡 , 𝑟 ) , 𝐺 𝑛 ( 𝑡 , 𝑟 ) are plotted in Figures 15, 16, 17, and 18, respectively.

407647.fig.0015
Figure 15: Absolute error of 𝑌 1 5 ( 𝑡 , 𝑟 ) for = 0 . 9 7 2 6 5 6 .
407647.fig.0016
Figure 16: Absolute error of 𝑌 1 5 ( 𝑡 , 𝑟 ) for = 0 . 9 7 2 6 5 .
407647.fig.0017
Figure 17: Absolute error of 𝐺 1 5 ( 𝑡 , 𝑟 ) for = 0 . 9 7 2 9 6 9 .
407647.fig.0018
Figure 18: Absolute error of 𝐺 1 5 ( 𝑡 , 𝑟 ) for = 0 . 9 7 2 9 6 4 .

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 ( 6 7 ) The exact solution is as follows: 𝑌 ( 𝑡 , 𝑟 ) = ( 2 + 𝑟 ) 𝑒 2 𝑡 + ( 1 𝑟 ) 𝑡 𝑒 2 𝑡 , ( 4 𝑟 ) 𝑒 2 𝑡 + ( 𝑟 1 ) 𝑡 𝑒 2 𝑡 , 𝐺 ( 𝑡 , 𝑟 ) = 𝑟 𝑒 2 𝑡 , ( 2 𝑟 ) 𝑒 2 𝑡 . ( 6 8 ) 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 𝑟 . ( 6 9 ) We choose the initial guesses of system (69) as follows: 𝑌 0 = 2 + 𝑟 , 𝑌 0 𝐺 = 4 𝑟 , 0 = 𝑟 , 𝐺 0 = 2 𝑟 , ( 7 0 ) and the auxiliary linear operators 𝑖 𝜙 𝑖 = ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝜙 𝑖 ( 𝑡 , 𝑟 ; 𝑞 ) 𝜕 𝑡 , 𝑖 = 1 , 2 , 3 , 4 , ( 7 1 ) with the properties 𝑖 𝑐 𝑖 = 0 , 𝑖 = 1 , 2 , 3 , 4 , ( 7 2 ) 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 𝑡 . ( 7 3 ) The 𝑚 th order deformation equations are 𝑖 𝑧 𝑖 , 𝑚 ( 𝑡 , 𝑟 ) 𝜒 𝑚 𝑧 𝑖 , 𝑚 1 ( 𝑡 , 𝑟 ) = 𝑖 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 , 𝑖 = 1 , 2 , 3 , 4 , ( 7 4 ) with the initial conditions 𝑧 𝑖 , 𝑚 ( 0 , 𝑟 ) = 0 , 𝑖 = 1 , 2 , 3 , 4 , ( 7 5 ) where 1 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 = 𝜕 𝑧 1 , 𝑚 1 𝜕 𝑡 𝑧 4 , 𝑚 1 1 𝜒 𝑚 ( 2 𝑡 𝑟 + 2 𝑡 + 2 𝑟 + 3 ) 𝑒 2 𝑡 , 2 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 = 𝜕 𝑧 2 , 𝑚 1 𝜕 𝑡 𝑧 3 , 𝑚 1 1 𝜒 𝑚 ( 2 𝑡 𝑟 2 𝑡 2 𝑟 + 7 ) 𝑒 2 𝑡 , 3 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 = 𝜕 𝑧 3 , 𝑚 1 𝜕 𝑡 + 𝑧 1 , 𝑚 1 1 𝜒 𝑚 ( 𝑟 𝑡 + 𝑡 + 3 𝑟 + 2 ) 𝑒 2 𝑡 , 4 , 𝑚 𝑧 1 , 𝑚 1 , , 𝑧 4 , 𝑚 1 = 𝜕 𝑧 4 , 𝑚 1 𝜕 𝑡 + 𝑧 2 , 𝑚 1 1 𝜒 𝑚 ( 𝑟 𝑡 𝑡 3 𝑟 + 8 ) 𝑒 2 𝑡 . ( 7 6 ) Therefore, we recursively obtain 𝑧 1 , 1 ( 𝑡 , 𝑟 ) = 𝑒 2 𝑡 + 3 𝑟 2 3 𝑒 2 𝑡 𝑟 2 2 𝑡 𝑒 2 𝑡 𝑡 + 𝑟 𝑡 + 𝑒 2 𝑡 𝑧 𝑟 𝑡 , 1 , 2 ( 𝑡 , 𝑟 ) = 𝑒 2 𝑡 5 2 4 + 5 𝑒 2 𝑡 2 4 + 3 𝑟 2 3 𝑒 2 𝑡 𝑟 2 + 5 2 𝑟 2 5 𝑒 2 𝑡 2 𝑟 2 2 𝑡 𝑒 2 𝑡 𝑡 2 5 2 𝑡 4 5 𝑒 2 𝑡 2 𝑡 4 + 𝑟 𝑡 + 𝑒 2 𝑡 + 𝑟 𝑡 1 1 2 𝑟 𝑡 4 + 5 𝑒 2 𝑡 2 𝑟 𝑡 4 2 2 𝑡 2 + 2 𝑟 𝑡 2 2 , 𝑧 2 , 1 ( 𝑡 , 𝑟 ) = 4 4 𝑒 2 𝑡 3 𝑟 2 + 3 𝑒 2 𝑡 𝑟 2 + 𝑒 2 𝑡 𝑡 𝑟 𝑡 𝑒 2 𝑡 𝑧 𝑟 𝑡 , 2 , 2 ( 𝑡 , 𝑟 ) = 4 4 𝑒 2 𝑡 + 1 5 2 4 1 5 𝑒 2 𝑡 2 4 3 𝑟 2 + 3 𝑒 2 𝑡 𝑟 2 5 2 𝑟 2 + 5 𝑒 2 𝑡 2 𝑟 2 + 𝑒 2 𝑡 𝑡 3 2 𝑡 4 + 5 𝑒 2 𝑡 2 𝑡 4 𝑟 𝑡 𝑒 2 𝑡 𝑟 𝑡 1 1 2 𝑟 𝑡 4 5 𝑒 2 𝑡 2 𝑟 𝑡 4 2 𝑡 2 2 𝑟 𝑡 2 2 , 𝑧 3 , 1 ( 𝑡 , 𝑟 ) = 3 4 3 𝑒 2 𝑡 4 + 7 𝑟 4 7 𝑒 2 𝑡 𝑟 4 𝑒 + 2 𝑡 2 𝑡 𝑡 2 𝑒 + 𝑟 𝑡 + 2 𝑡 𝑟 𝑡 2 , 𝑧 3 , 2 ( 𝑡 , 𝑟 ) = 3 4 3 𝑒 2 𝑡 4 + 2 𝑒 2 𝑡 2 + 7 𝑟 4 7 𝑒 2 𝑡 𝑟 4 + 1 1 2 𝑟 4 1 1 𝑒 2 𝑡 2 𝑟 4 𝑒 + 2 𝑡 2 𝑡 𝑡 2 + 3 2 𝑡 𝑒 2 𝑡 2 𝑒 𝑡