Solution of Second-Order IVP and BVP of Matrix Differential Models Using Matrix DTM

Abstract

We introduce a matrix form of differential transformation method (DTM) and apply for nonlinear second-order initial value problems (IVPs) and boundary value problems (BVPs) of matrix models which are given by 𝐮^′′(𝑡)=𝑓(𝑡,𝐮(𝑡),𝐮^′(𝑡)) and subject to initial conditions 𝐮(𝑎)=𝐮₀,𝐮^′(𝑎)=𝐮₁ and boundary conditions 𝐮(𝑎)=𝐮₀,𝐮(𝑏)=𝐮₁, where 𝐮₀,𝐮₁∈𝑅^𝑟×𝑞. Also the convergence of present method is established. Several illustrative examples are given to demonstrate the effectiveness of the present method.

1. Introduction

In this work, we study two important cases of nonlinear second-order matrix models given by the following.(1)Matrix initial-value problems are of the following form: 𝐮^′′(𝑡)=𝑓𝑡,𝐮(𝑡),𝐮^′(𝑡),𝐮(𝑎)=𝐮₀,𝐮^′(𝑎)=𝐮₁,𝑎≤𝑡≤𝑏,[𝑎,𝑏]⊂ℝ. (1.1)(2)Matrix boundary-value problems are of the following form: 𝐮^′′(𝑡)=𝑓𝑡,𝐮(𝑡),𝐮^′(𝑡),𝐮(𝑎)=𝐮₀,𝐮(𝑏)=𝐮₁,𝑎≤𝑡≤𝑏,[𝑎,𝑏]⊂ℝ. (1.2)

Matrices 𝐮₀,𝐮₁,𝐮(𝑡)∈ℝ^𝑟×𝑞 and matrix function 𝑓∶[𝑎,𝑏]×ℝ^𝑟×𝑞×ℝ^𝑟×𝑞↦ℝ^𝑟×𝑞 are frequent in different fields in physics and engineering. Note that (1.1) could be the statement of Newton second law of motion for a coupled mechanical system. These problems generally arise frequently in many areas of science and engineering, for example, fluid mechanics, quantum mechanics, optimal control, chemical reactor theory, aerodynamics, reaction-diffusion process, geophysics, and so forth, where one solves scalar or vectorial problems with boundary-value conditions [1–3].

Many authors studied (1.1), (1.2), and other similar forms of these equation, by different numerical methods, such as Jodar et al. [4–6], Al-Said [7], Al-Said and Noor [8], and Hargreaves and Higham [9]. In [10], (1.1) studied was subject to initial conditions using Cubic Matrix Splines method.

The differential transform method (DTM) is a seminumerical/analytic technique that formalizes the Taylor series in a totally different manner. It was first introduced by Zhou in a study about electrical circuits [11]. The DTM obtains an analytical solution in the form of a polynomial. It is different from the traditional high-order Taylor’s series method, which requires symbolic competition of the necessary derivatives of the data functions. The Taylor series method is computationally taken long time for large orders. With this method, it is possible to obtain highly accurate results or exact solutions for differential equations. The differential transform is an iterative procedure for obtaining analytic Taylor series solutions of differential equations. DTM has been successfully applied to solve many nonlinear problems arising in engineering, physics, mechanics, biology, and so forth. Ayaz [12] applied DTM for solution of system of ordinary differential equations. Arikoglu and Ozkol [13] employed DTM on differential-difference equations. Furthermore, the method may be employed for the solution of partial differential equations. Abazari et al. employed DTM on some partial differential equations (PDEs) and their coupled version in [14–18] and extended DTM to solve the first and second kinds of the Riccati matrix differential equations in [19]. Recently, Abazari and Ganji [20] extended to the reduced differential transform method (RDTM) to study the partial differential equation with-proportional delay, and [21] applied RDTM on simulation of generalized Hirota-Satsuma coupled KdV equation. We also note that the differential transformation method was applied for solution of fractional order nonlinear boundary value problems in [22] and the nonlinear higher-order boundary value problems; see [23].

The purpose of this paper is extended DTM, mentioned in [19], to apply for (1.1) subject to initial conditions and (1.2) subject to boundary conditions.

The layout of the paper is as follows. In Section 2, the matrix differential transform method will be introduced. In Section 3, the convergence of present method is illustrated. In Section 4, some numerical results are given to clarify the method and a comparison is made with the existing results. Section 5 is the brief conclusion of this paper. Finally some references are listed in the end. Note that we have computed the numerical results by Maple programming.

2. Basic Definitions

With reference to [5–9, 13, 14], we introduce in this section the basic definition of the matrix form of one-dimensional DTM.

Definition 2.1. If 𝐮(𝑡)∈ℝ^𝑛×𝑛 is matrix analytical function in the domain 𝑇, then it will be differentiated continuously with respect to time 𝑡∈𝑇, 𝑑^𝑘𝐮(𝑡)𝑑𝑡^𝑘=𝜙(𝑡,𝑘), (2.1) where 𝑘 belongs to the set of nonnegative integer, denoted by the 𝐾 domain. Therefore, for 𝑡=𝑡_𝑖, (2.1) can be written as 𝐔_𝑖(𝑘)=𝜙𝑡_𝑖=𝑑,𝑘^𝑘𝐮(𝑡)𝑑𝑡^𝑘_{𝑡=𝑡𝑖},∀𝑘∈𝐾, (2.2) where 𝐔_𝑖(𝑘)∈ℝ^𝑛×𝑛 is called the spectrum of 𝐮(𝑡) at 𝑡=𝑡_𝑖, in the 𝐾 domain.

Definition 2.2. If 𝐮(𝑡)∈ℝ^𝑛×𝑛 can be expressed by Taylor’s series about fixed point 𝑡_𝑖, then 𝐮(𝑡) can be represented as 𝐮(𝑡)=^∞_𝑘=0𝐮^(𝑘)𝑡_𝑖𝑘!𝑡−𝑡_𝑖^𝑘. (2.3) If 𝐮_𝑛(𝑡)=∑^𝑛_𝑘=0(𝐮^(𝑘)(𝑡_𝑖)/𝑘!)(𝑡−𝑡_𝑖)^𝑘 is the 𝑛-partial sums of a Taylor’s series (2.3), then 𝐮(𝑡)=^𝑛_𝑘=0𝐮^(𝑘)𝑡_𝑖𝑘!𝑡−𝑡_𝑖^𝑘+𝑅_𝑛(𝑡), (2.4) where 𝐮_𝑛(𝑡) is called the 𝑛th Taylor polynomial for 𝐮(𝑡) about 𝑡_𝑖 and 𝑅_𝑛(𝑡) is remainder term. If 𝐔(𝑘) is defined as 𝐔(𝑘)=𝐮^(𝑘)𝑡_𝑖𝑘!,𝑘=0,1,…, (2.5) then (2.3) reduces to 𝐮(𝑡)=^∞_𝑘=0𝐔(𝑘)𝑡−𝑡_𝑖^𝑘, (2.6) and the 𝑛-partial sums of a Taylor’s series (2.6) reduce to 𝐮_𝑛(𝑡)=^𝑛_𝑘=0𝐔(𝑘)𝑡−𝑡_𝑖^𝑘. (2.7) The 𝐔(𝑘) defined in (2.5) is called the matrix differential transform of matrix function 𝑢(𝑡). For a special case, when 𝑡₀=0, then solution (2.6) reduces to 𝐮(𝑡)=^𝑛_𝑘=0𝐔(𝑘)𝑡^𝑘+𝑅_𝑛(𝑡). (2.8) From the aforementioned definitions, it can be found that the concept of the one-dimensional matrix differential transform is derived from the Taylor series expansion. With (2.5) and (2.6), the fundamental mathematical operations performed by one-dimensional matrix differential transform can readily be obtained and listed in Table 1.

Table 1 

The fundamental operations of one-dimensional matrix differential transform method.

3. Convergence Analysis

In this section, we show that the presented matrix differential transformation method is convergence.

Theorem 3.1. Let the matrix function 𝐮∶ℝ^𝑛×𝑛→ℝ^𝑛×𝑛 and its first 𝑛 derived functions be relatively continuous and finite on an interval 𝑇 and differentiable on 𝑇−𝑄 (𝑄 countable). Let 𝑡₀,𝑡∈𝑇. Then (2.3) and (2.4) hold, with 𝑅_𝑛(𝑥)=1𝑛!^𝑥_𝑡0𝐮^(𝑛+1)(𝑡)(𝑥−𝑡)^𝑛𝑅𝑑𝑡,_𝑛(𝑥)_∞≤𝑥−𝑡₀^𝑛+1()!𝐮𝑛+1^(𝑛+1)(𝑡)_∞, (3.1) where ‖⋅‖_∞ is infinity matrix norm.

Proof. By Definition 2.2, we get 𝑅_𝑛(𝑥)=𝐮(𝑥)−𝐮𝑡₀−^𝑛_𝑘=1𝐮^(𝑘)𝑡₀(𝑥−𝑡₀)^𝑘.𝑘! (3.2) We use the right-hand side as a “pattern” to define a function 𝐡∶ℝ^𝑛×𝑛→ℝ^𝑛×𝑛. This time, we keep 𝑥 fixed (say, 𝑥=𝑎∈𝑇) and replace 𝑡₀ by a variable 𝑡. Thus we set 𝐡(𝑡)(𝑎)(𝑡)−𝐮=𝐮−𝐮^′(𝑡)()𝐮1!𝑎−𝑡−⋯−^(𝑛)(𝑡)()𝑛!𝑎−𝑡^𝑛,𝑡∈ℝ. (3.3) Then 𝐡(𝑡₀)=𝑅_𝑛(𝑎) and 𝐡(𝑎)=0. Our assumptions imply that 𝐡 is relatively continuous and finite on 𝑇 and differentiable on 𝑇−𝑄. Differentiating (3.3), we see that all cancels out except for one term: 𝐡^′(𝑡)𝐮=−^(𝑛+1)(𝑡)()𝑛!𝑎−𝑡^𝑛,𝑡∈𝑇−𝑄. (3.4) Then for 𝑡∈𝑇 we get −𝐡(𝑡)=∫^𝑎_𝑡(𝐮^(𝑛+1)(𝑡)/𝑛!)(𝑎−𝑡)^𝑛𝑑𝑡, and ^𝑎_𝑡0𝐮^(𝑛+1)(𝑡)()𝑛!𝑎−𝑡^𝑛𝑑𝑡=−𝐡(𝑎)+𝐡𝑡₀=𝑅_𝑛(𝑎). (3.5) As 𝑥=𝑎, (3.1) is proved. Next, let 𝑀=‖𝐮^(𝑛+1)(𝑡)‖_∞. If 𝑀=+∞, (3.1) is valid. If 𝑀<+∞, we define 𝑔(𝑡)=𝑀((𝑡−𝑎)^𝑛+1/(𝑛+1)!) for 𝑡≥𝑎, and 𝑔(𝑡)=−𝑀((𝑎−𝑡)^𝑛+1/(𝑛+1)!) for 𝑡≤𝑎. In both cases, for 𝑡∈𝑇−𝑄 we have 𝑔^′(𝑡)=𝑀(|𝑡−𝑎|^𝑛/𝑛!)≥‖𝐡^′(𝑡)‖_∞; then we get ‖𝐡(𝑡₀)−𝐡(𝑎)‖_∞≤‖𝑔(𝑡₀)−𝑔(𝑎)‖_∞, or ‖𝑅_𝑛(𝑎)‖_∞≤𝑀(|𝑎−𝑡₀|^𝑛+1/(𝑛+1)!). Thus (3.1) follows, because 𝑎 is arbitrary value.

From Theorem 3.1, we get that if 𝑛→∞, then ‖𝑅_𝑛(𝑥)‖_∞→0_𝑛×𝑛. Then by this theorem, the convergence of this method is investigated.

4. Applications and Numerical Results

This section is devoted to computational results. We applied the method presented in this paper and solved six examples. The first two examples are allotted to matrix initial-value problem and the next examples are about matrix boundary-value form. These examples are chosen such that there exist exact solutions for them. All calculations are implemented by Maple.

Example 4.1. In the first example, we consider the following non linear second-order matrix initial-value problem: 𝐮^′′(𝑡)+𝐮^′(𝑡)𝐮(𝑡)(𝑡),𝐮(0)==𝐜0−101,𝐮^′(0)=,2001 (4.1) where 𝐜(𝑡)=^{2𝑡3−6𝑡2+4𝑡−22𝑡−2}_{0𝑒𝑡+𝑒2𝑡}. By applying matrix differential transform operator Table 1 on nonlinear system (4.1), for 𝑘=0,1,2,…,𝑛, we get (𝑘+2)!𝐔()+𝑘!𝑘+2^𝑘_𝑙=0()𝐔()𝐔(𝑙)(𝑘),𝐔(0)=(1)=,𝑘−𝑙+1𝑘−𝑙+1=𝐂0−101,𝐔2001 (4.2) where 𝐔(𝑘) and 𝐂(𝑘) are the differential transform of 𝐮(𝑡) and 𝐜(𝑡), respectively. From (4.2), and for 𝑘=0,1,2,…, we get 𝐔(𝑘+2)=−𝑘!()!,𝑘+2^𝑘_𝑙=0(𝑘−𝑙+1)𝐔(𝑘−𝑙+1)𝐔(𝑙)−𝐂(𝑘), (4.3) we substitute the initial condition (4.2), in recursive equation (4.3), for 𝑘=0,1,2, we get 𝐔(2)=⁻¹⁰_01/2,𝐔(3)=⁰⁰_01/6, and 𝐔(4)=⁰⁰_01/24, and then from the inverse differential transform operator (2.7), we get 𝐔₄(𝑡)=+010−1012001𝑡+−102𝑡²+01006𝑡³+0100𝑡24⁴. (4.4) The closed form of above solution is 𝐮(𝑡)=^{2𝑡−𝑡2−1}_0𝑒𝑡, which is exactly the same as the exact solution.

Example 4.2. In this example, we consider the following nonlinear second-order matrix differential equations: 𝐮^′′(𝑡)+𝐮(𝑡)𝐮^′(𝑡)−𝐮²(𝑡)(𝑡),𝐮(0)==𝑐1−101,𝐮^′(0)=,1−111 (4.5) where 𝑐(𝑡)=^{𝑒𝑡−𝑒2𝑡−𝑒𝑡}_{2𝑒𝑡+𝑡𝑒𝑡+𝑒2𝑡𝑒𝑡}. The differential transform version of nonlinear system (4.5), for 𝑘=0,1,2,…,𝑛, is (𝑘+2)!𝐔()+𝑘!𝑘+2^𝑘_𝑙=0(𝑙+1)𝐔(𝑘−𝑙)𝐔(𝑙+1)−^𝑘_𝑙=0𝐔()𝐔(𝑙)(𝑘),𝐔(0)=(1)=,𝑘−𝑙=𝐂1−101,𝐔1−111 (4.6) where 𝐔(𝑘) and 𝐂(𝑘) are the matrix differential transform of 𝐮(𝑡) and 𝐜(𝑡), respectively. From (4.6), and for 𝑘=0,1,2,…, we get following iteration equation: 𝐔(𝑘+2)=−𝑘!(𝑘+2)!^𝑘_𝑙=0(𝑙+1)𝐔(𝑘−𝑙)𝐔(𝑙+1)−^𝑘_𝑙=0𝐔(𝑘−𝑙)𝐔(𝑙)−𝐂(𝑘). (4.7) By utilizing the initial values 𝐔(0) and 𝐔(1) in recursive equations (4.7), for 𝑘=0,1,2, the first four terms of 𝐔(𝑘) are obtained as follows: 𝐔(2)=^1/2−1/2_11/2,𝐔(3)=^1/6−1/6_1/21/6, and 𝐔(4)=^1/24−1/2_1/61/24, and then from (2.7), we get 𝐔₄(𝑡)=+11−1011−111𝑡+2−12112𝑡²+16−161216𝑡³+1−1242161𝑡24⁴. (4.8) The closed form of the previous solution is 𝐮(𝑡)=^{𝑒𝑡−𝑒𝑡}_{𝑡𝑒𝑡𝑒𝑡}, which is exactly the same as the exact solution.

Example 4.3. In the third example, consider the following linear second-order boundary-value problem: 𝐮^′′(𝑡)+𝐚(𝑡)𝐮^′(𝑡)(𝑡)𝐮(𝑡)(𝑡),𝐮(0)(1)=,+𝐛=𝐜=𝐮0000 (4.9) where 𝐚(𝑡)=^−3𝑡−3₋₃₀,𝐛(𝑡)=⁹⁰_−𝑡1, and 𝐜(𝑡)^{9𝑡2−2−6𝑡2−6𝑡+5}_{−𝑡4−8𝑡2+2𝑡2−𝑡+3}. By applying matrix differential transform operator listed in Table 1, on linear system (4.9), for 𝑘=0,1,2,…,𝑛, we obtain (𝑘+2)!𝐔()+𝑘!𝑘+2^𝑘_𝑙=0(𝑙+1)𝐀(𝑘−𝑙)𝐔(𝑙+1)+^𝑘_𝑙=0𝐁()𝐔(𝑙)(𝑘),𝐔(0)=,𝑘−𝑙=𝐂0000^𝑛𝑘=0𝐔(𝑘)=,0000 (4.10) where 𝐔(𝑘),𝐀(𝑘),𝐁(𝑘), and 𝐂(𝑘) are the matrix differential transform of 𝐮(𝑡),𝐚(𝑡),𝐛(𝑡), and 𝐜(𝑡), respectively. For 𝑘=0,1,2,…, (4.10) can be rewriten as in the following iterative equation: 𝐔(𝑘+2)=−𝑘!(𝑘+2)!^𝑘_𝑙∑₌₀(𝑙+1)𝐀(𝑘−𝑙)𝐔(𝑙+1)−∑^𝑘_𝑙=0𝐁(𝑘−𝑙)𝐔(𝑙)−𝐂(𝑘). (4.11) Then from (4.11), for 𝑘=0,1,2, we get 𝐔(2)1=−2(𝐀(0)𝐔(1)(0)𝐔(0)(0,𝐔(3)1+𝐁−𝐂))=−6((0)𝐔(2)(1)𝐔(1)(0)𝐔(1)(1)𝐔(0)(1,𝐔(4)12𝐀+𝐀+𝐁+𝐁−𝐂))=−(𝐀(0)𝐔(3)(1)𝐔(2)(2)𝐔(1)(0)𝐔(2)(1)𝐔(1)(2)𝐔(0)(2.12+2𝐀+𝐀+𝐁+𝐁+𝐁−𝐂)) (4.12) In the same manner, the rest of components were obtained using the MAPLE Package. From boundary condition (4.10) we get 𝐔(0)=,0000^𝑛𝑘=0𝐔(𝑘)=0000. (4.13) To obtain the remain coefficient component of 𝐔(𝑘), for 𝑘≥1, it is enough to find 𝐔(1) from differential transform version of boundary conditions (4.13). Therefore, assume that 𝐘(1)=^{𝑎1,1𝑎1,2}_{𝑎2,1𝑎2,2}; then for 𝑛=3, from ∑³_𝑘=0𝐔(𝑘)=⁰⁰₀₀, we get 32𝑎_1,1+32𝑎_2,132𝑎_1,2+32𝑎_2,27+33𝑎_2,1+32𝑎_1,173𝑎_2,2+32𝑎_1,2+236=0000. (4.14) By solving (4.14), we get 𝐔(1)=⁰⁻¹₀₋₁. Note that for 𝑛≥3, by solving ∑^𝑛_𝑘=0𝐔(𝑘)=⁰⁰₀₀, we always get 𝐔(1)=⁰⁻¹₀₋₁. By substituting 𝐔(0) and 𝐔(1) in list (4.12), we get 𝐔(2)=⁻¹¹₁₀,𝐔(3)=¹⁰₋₁₁, and for 𝑛≥4, we have 𝐔(𝑛)=⁰⁰₀₀. Therefore from inverse differential transform operator (2.7), the four-term approximation solutions is obtained as follows: 𝐔₄(𝑡)=+𝑡00000−10−1𝑡+−1110²+𝑡10−11³, (4.15) which is exactly the same as the exact solution 𝐮(𝑡)=^{𝑡3−𝑡2𝑡2−𝑡}_{−𝑡3+𝑡2𝑡3−𝑡}.

Example 4.4. Finally, we consider the following nonlinear second-order matrix boundary-value problem: 𝐮^′′(𝑡)−𝐮^′(𝑡)𝐮(𝑡)(𝑡)(𝑡),𝐮(0)(1)=,−𝐮=𝐜=𝐮00000≤𝑡≤1, (4.16) where 𝐜(𝑡)=^{2+4𝑡2−7𝑡3+3𝑡4−2+6𝑡−2𝑡2+5𝑡3−3𝑡5}_{−2+𝑡−3𝑡2−𝑡3+3𝑡46𝑡+𝑡2+2𝑡4−3𝑡5}. Then for 𝑘=0,1,2,…,𝑛, we have (𝑘+2)!𝐔()−𝑘!𝑘+2^𝑘_𝑙=0()𝐔()𝐔(𝑙)(𝑘)(𝑘),𝐔(0)=,𝑘−𝑙+1𝑘−𝑙+1−𝐔=𝐂0000^𝑛𝑘=0𝐔(𝑘)=.0000 (4.17) Therefore from (4.17), some of the first quantities are obtained as follows: 𝐔(2)=12𝐔(1)𝐔(0)+12𝐔(0)+12𝐂(0),𝐔(3)=13𝐔(2)𝐔(0)+16𝐔(1)²+16𝐔(1)+16𝐂(1),𝐔(4)=14𝐔(3)𝐔(0)+14𝐔(2)𝐔(1)+1𝐔(2)+112𝐂(2),𝐔(5)=1125𝐔(4)𝐔(0)+15𝐔(3)𝐔(1)+1𝐔(2)10²+1𝐔(3)+120𝐂(3).20 (4.18) In the same manner, the rest of components were obtained using the MAPLE Package. Similar to the previous examples, assume that 𝐔(1)=^{𝑎1,1𝑎1,2}_{𝑎2,1𝑎2,2}; then by assuming the various 𝑛=3, in the boundary condition (4.17), we get 71+6𝑎_1,1+16𝑎²_1,1+16𝑎_1,2𝑎_2,176𝑎_1,2+16𝑎_1,1𝑎_1,2+16𝑎_1,2𝑎_2,2−56+76𝑎_2,1+16𝑎_1,1𝑎_2,1+16𝑎_2,1𝑎_2,271+6𝑎_2,2+16𝑎_1,2𝑎_2,1+16𝑎²_2,2=0000. (4.19) From the common solution of (4.19), we obtain 𝐔(1)=⁻¹⁰₁₋₁, and by setting 𝑛=4 we get 17𝑎12^1,1+17+1126𝑎²_1,143𝑎_1,2+16𝑎_1,1𝑎_1,2+16𝑎_1,2𝑎_2,2+16𝑎_1,2𝑎_2,1−16𝑎_2,1−1𝑎12^1,2−16𝑎_2,2−1𝑎12^1,1−1454𝑎_2,1−76+16𝑎_2,1𝑎_1,176𝑎_2,2+16𝑎_1,2𝑎_2,1+16𝑎²_2,2+16𝑎_2,2𝑎_2,1−16𝑎_1,1−1𝑎12^2,2+13−1126𝑎_1,2−1𝑎12^2,1=[0]_2×2. (4.20) In this case also from the common solution of (4.20), we get 𝐔(1)=⁻¹⁰₁₋₁. By a sample computation for 𝑛=5,6,…, it is easy to obtain that the common solution of ∑^𝑛_𝑘=0𝐔(𝑘)=⁰⁰₀₀ is 𝐔(1)=⁻¹⁰₁₋₁. By substituting the 𝐔(0) and 𝐔(1) in (4.18) and utilizing the obtained quantities coefficients 𝐔(∶) in (2.7), the series form of solution of equation (4.16) is obtained as follows: 𝐔₄(𝑡)=+𝑡0000−101−1𝑡+1−1−10²+𝑡0101³. (4.21) The closed form of previous solution is 𝐮(𝑡)=^{−𝑡+𝑡2−𝑡2+𝑡3}_{𝑡−𝑡2−𝑡+𝑡3}, which is exactly the same as the exact solution.

5. Conclusions

In this paper, we have shown that DTM can successfully be used for solving the linear and nonlinear second-order Matrix IVPs and Matrix BVPs. This method is simple and easy to use and solves the problem without any need for discretizing the variables. Also this method is useful for finding an accurate approximation of the exact solution. A symbolic calculation software package, Maple, is used in the derivations. The method gives rapidly converging series solutions. The accuracy of the obtained solution can be improved by taking more terms in the solution. In many cases, the series solutions obtained with DTM can be written in exact closed form. Further we note that DTM can also be applied to solve the nonlinear matrix differential Riccati equations (first and second kinds of Riccati matrix differential equations) where DTM technique provides a sequence of matrix functions which converges to the exact solution of the problem; see [19]. In fact the present work is the extension of [19].