On Spectral Homotopy Analysis Method for Solving Linear Volterra and Fredholm Integrodifferential Equations

and Applied Analysis 3 which can be dealt in similar way. The so-called zero-order deformation equation was constructed by Liao as follows: ( 1 − p)L[φ(τ ; p) − u0 τ ] phH τ (N[φ(τ ; p)]), 2.2 where p ∈ 0, 1 is the embedding parameter and, h is a nonzero convergence-parameter, H τ is an auxiliary function and u0 τ is called an initial guess of u τ and φ τ ; p is an unknown function. In addition, L is an auxiliary linear operator and N is nonlinear operator as follows: L ( φ ( x; p )) a1 x ∂2φ ( x; p ) ∂x2 , 2.3 with the property L c1t c2 0 where c1 and c2 are constants and N [ φ ( x; p )] a1 x ∂2φ ( x; p ) ∂x2 a2 x ∂φ ( x; p ) ∂x a3 x φ ( x; p ) − g x − ∫x −1 k x, t φ t dt 2.4 is a nonlinear operator. Obviously, when p 0 and p 1, it holds φ τ ; 0 u0 τ and φ τ ; 1 u τ . In this way, as p increase from 0 to 1, φ τ ; p alter from initial guess u0 τ to the solution u τ and φ τ ; p is expanded in Taylor series with respect to p: φ ( τ ; p ) u0 τ ∞ ∑ m 1 um τ p, 2.5 where um τ Dm [ φ ( τ ; p )] , 2.6 Dmφ 1 m! ∂ ( φ ) ∂pm ∣∣∣∣ p 0 . 2.7 The series 2.5 converges at p 1 if the auxiliary linear operator, the initial guess, the convergence parameter, and the auxiliary function are properly selected: φ τ u0 τ ∞ ∑ m 1 um τ . 2.8 The admissible and valid values of the convergence-parameter h are found from the horizontal portion of the h-curves. Liao proved that u τ is one of the solutions of original nonlinear equation. AsH τ 1 so 2.2 becomes ( 1 − p)L[φ(τ ; p) − u0 τ ] ph(N[φ(τ ; p)]). 2.9 4 Abstract and Applied Analysis Define the vector um {u0 τ , u1 τ , . . . , um τ }. Operating on both sides withDm, we have the so called mth-order deformation equation L [ um τ − χmum−1 τ ] hH τ Rm um−1 τ , 2.10 where Rm um−1 1 m − 1 ! ∂m−1N [ φ ( τ ; p )] ∂pm−1 ∣∣∣∣∣ p 0 , χm { 0, m ≤ 1 1, otherwise. 2.11 um τ for m ≥ 0 is governed by the linear equation 2.10 and can be solved by symbolic computation software such as Maple, Matlab, and so on. 3. Spectral Homotopy Analysis Solution Consider the second-order Volterra integrodifferential equation a1 x u′′ x a2 x u′ x a3 x u x g x ∫x −1 k x, t u t dt, u −1 0, u 1 0. 3.1 To obtain initial approximation we solve the following two-point boundary value problem: a1 x u′′ 0 x a2 x u ′ 0 x a3 x u0 x g x , 3.2 subject to boundary conditions u0 −1 0, u0 1 0. 3.3 We use the Chebyshev pseudospectral method to solve 3.2 . At first we approximate u0 τ by a truncated series of Chebyshev polynomial of the form u0 τ ≈ u0 ( τj ) N ∑ k 0 ûkTk ( τj ) , j 0, . . . ,N, 3.4 where Tk is the kth Chebyshev polynomials, ûk, are coefficients and τ0, τ1, . . . , τN are GaussLobatto points which are the extrema of the Nth-order Chebyshev polynomial, defined by τj cos ( πj N ) . 3.5 Abstract and Applied Analysis 5 Derivatives of the functions y0 τ at the collocation points are represented as du0 τk dτs N ∑ j 0 D kju0 ( τj ) , k 0, . . . ,N, 3.6and Applied Analysis 5 Derivatives of the functions y0 τ at the collocation points are represented as du0 τk dτs N ∑ j 0 D kju0 ( τj ) , k 0, . . . ,N, 3.6 where s is the order of differentiation andD is the Chebyshev spectral differentiation matrix. By following 24 we express the entries of the differentiation matrix D as Dkj (−1 2 ) ck cj −1 k j sin ( π ( j k ) /2N ) sin ( π ( j − k)/2N) , j / k, Dkj (−1 2 ) xk sin2 πk/N , k / 0,N, k j, D00 −DNN 2N 2 1 6 . 3.7 Substituting 3.4 – 3.6 into 3.2 results in AU0 G 3.8 subject to the boundary conditions u0 τ0 u0 τN 0, 3.9 where A a1D a2D a3, U0 u0 τ0 , u0 τ1 , . . . , u0 τN T , G [ g τ0 , g τ1 , . . . , g τN ]T , ar diag ar τ0 , ar τ1 , . . . ar τN , r 1, 2, 3. 3.10 To satisfy the boundary conditions we delete the first and the last rows and columns of A and the first and the last rows of Y0 andG. The values of u0 τi , i 0, . . . ,N are determined from the equation U0 A−1G, 3.11 which is the initial approximation for the SHAM solution of the governing equation 3.1 . The zeroth-order deformation equation is given by ( 1 − p)L[φ(τ ; p)u0 τ ] ph(N[φ(τ ; p)] − g τ ), 3.12 6 Abstract and Applied Analysis where L [ φ ( τ ; p )] a1 ∂2φ ∂τ2 a2 ∂φ ∂τ a3φ, 3.13 h is the nonzero convergence controlling auxiliary parameter, and N is a nonlinear operator given by N [ φ ( τ ; p )] a1 ∂2φ ∂τ2 a2 ∂φ ∂τ a3φ − ∫ τ −1 k τ, t φ t dt. 3.14 Differentiating 3.2 m times with respect to the embedding parameter p and then setting p 0 and finally dividing them by m!, we have the so-called mth-order deformation equation L [ um τ − χmum−1 τ ] hRm 3.15 subject to boundary conditions um −1 um 1 0, 3.16


Introduction
Functional equations such as partial differential equations, integral and integrodifferential equations and others are widely employed to model complex real-life problems.Many physical events in different fields of sciences and engineering have been formulated using integrodifferential equations 1, 2 .Integrodifferential equations are usually difficult to solve analytically so there is a need to obtain an efficient approximate solution 3, 4 .Homotopy analysis method HAM was first proposed by Liao in 5 .it is based on the homotopy concept in topology for solving nonlinear differential equations.Unlike the traditional perturbation methods like Lyapunov's artificial small parameter method 6 , Adomian decomposition method 7-10 and the δ-expansion method 11 , which are the special cases of HAM, this approach does not need a small perturbation parameter.In the HAM, original nonlinear problem is converted to infinite number of linear problems without using the perturbation techniques 12 .Homotopy analysis method is more powerful than traditional perturbation methods since it is applicable for solving problems with strong nonlinearity 13-15 , even if they do not have any small or large parameters.
We can adjust the convergence region and the rate of approximation series to give us freedom to use different base function to approximate a nonlinear problem.He in 16 proposed the so called Homotopy perturbation method HPM .Later it was pointed by Hayat et al. 17 that HPM is only a special case of the HAM h −1 .Many nonlinear problems were solved by Homotopy perturbation method 18, 19 .The HPM does not provide us with convergent series solution for strongly nonlinear problems as was indicated by Abbasbandy in 20 using a simple problem in which the physical parameters were large.In 2010, Motsa et al.21 employed Chebyshev polynomials to solve higher-order deformation equation and called his approach spectral homotopy analysis method SHAM .The spectral homotopy analysis method has been used only for solving partial and ordinary differential equations 22, 23 .Motsa et al. 21-23 found that the Spectral homotopy analysis method is more efficient than the homotopy analysis method as it does not depend on the rule of solution expression and the rule of ergodicity.This method is more flexible than homotopy analysis method since it allows for a wider range of linear and nonlinear operators and one is not restricted to use the method of higher order differential mapping for solving boundary value problems in bounded domains, unlike the homotopy analysis method.The range of admissible h values is much wider in spectral homotopy analysis method than in homotopy analysis method.The main restriction of HAM in solving integral equations is to choose the best initial guess as the series solution be convergent.In SHAM the initial approximation is taken to be the solution of the nonhomogeneous linear part of the given equation.
In this paper, we propose spectral homotopy analysis method SHAM to solve linear Volterra and Fredholm type of integrodifferential equations.Volterra integrodifferential equation is given by The paper is organized in the following way.Section 2 included a brief introduction in homotopy analysis method.Spectral homotopy analysis method for solving integral equations is presented in Section 3. We then propose the way to calculate S m−1 which is needed to obtain Y m in Section 4. Numerical examples are presented in Section 5.In Section 6, concluding remarks are given.

Homotopy Analysis Solution
In this section, we give a brief introduction to HAM.We consider the following differential equation in a general form: where N is nonlinear operator, τ denotes independent variables and u τ is an unknown function, respectively.For simplicity we disregard all initial and all boundary conditions which can be dealt in similar way.The so-called zero-order deformation equation was constructed by Liao as follows: where p ∈ 0, 1 is the embedding parameter and, h is a nonzero convergence-parameter, H τ is an auxiliary function and u 0 τ is called an initial guess of u τ and φ τ; p is an unknown function.In addition, L is an auxiliary linear operator and N is nonlinear operator as follows: with the property L c 1 t c 2 0 where c 1 and c 2 are constants and is a nonlinear operator.Obviously, when p 0 and p 1, it holds φ τ; 0 u 0 τ and φ τ; 1 u τ .In this way, as p increase from 0 to 1, φ τ; p alter from initial guess u 0 τ to the solution u τ and φ τ; p is expanded in Taylor series with respect to p: where The series 2.5 converges at p 1 if the auxiliary linear operator, the initial guess, the convergence parameter, and the auxiliary function are properly selected: The admissible and valid values of the convergence-parameter h are found from the horizontal portion of the h-curves.Liao proved that u τ is one of the solutions of original nonlinear equation.As H τ 1 so 2.2 becomes Define the vector u m {u 0 τ , u 1 τ , . . ., u m τ }.Operating on both sides with D m , we have the so called mth-order deformation equation where

2.11
u m τ for m ≥ 0 is governed by the linear equation 2.10 and can be solved by symbolic computation software such as Maple, Matlab, and so on.

Spectral Homotopy Analysis Solution
Consider the second-order Volterra integrodifferential equation

3.1
To obtain initial approximation we solve the following two-point boundary value problem: We use the Chebyshev pseudospectral method to solve 3.2 .At first we approximate u 0 τ by a truncated series of Chebyshev polynomial of the form where T k is the kth Chebyshev polynomials, u k , are coefficients and τ 0 , τ 1 , . . ., τ N are Gauss-Lobatto points which are the extrema of the Nth-order Chebyshev polynomial, defined by Derivatives of the functions y 0 τ at the collocation points are represented as where s is the order of differentiation and D is the Chebyshev spectral differentiation matrix.By following 24 we express the entries of the differentiation matrix D as

3.10
To satisfy the boundary conditions we delete the first and the last rows and columns of A and the first and the last rows of Y 0 and G.The values of u 0 τ i , i 0, . . ., N are determined from the equation where 3.17 Applying the Chebyshev pseudospectral transformation on 3.15 -3.17 results in subject to the boundary conditions where A and G are as defined in 3.10 and To implement the boundary conditions we delete the first and the last rows of S m−1 and the first and the last rows and columns of A.
Finally, this recursive formula can be written as following for m ≥ 1: so, starting from the initial approximation.We can obtain higher-order approximation U m for m ≥ 1, recursively.

Calculation of S m−1
Consider the well-known Chebyshev polynomials of first kind of degree n, defined by 25 Also they are derived from the following recursive formula: In this section we approximate the integrand in 3.15 by letting where T j t , j 0, 1, . . ., m − 1 are the basis functions.
Clearly, w mj , 0 ≤ j ≤ m − 1 can be obtained by solving the following matrix equation: where τ j t j , j 0, 1, . . ., m − 1 are Chebyshev collocation points defined in 2.6 .We can obtain w mj by solving the above system so the integral in 3.15 can be rewritten as following:  where P as in 26 , that is the N 1 × N 1 operational matrix for integration as follows:

Numerical Examples
In this section we apply the technique described in Section 3 to some illustrative examples of second-order Volterra integrodifferential equations.
Example 5.1.Consider the second-order Volterra integrodifferential equation x − 2t y t dt 5.1 subject to the boundary conditions y −1 y 1 0 with the exact solution y x x 2 − 1.
We employ HAM and SHAM to solve this example.In Table 1, there is a comparison of the numerical result against the HAM and SHAM approximation solutions at different orders.It is worth noting that the SHAM results become very highly accurate only with a few iterations and fourth-order results become very close to the exact solution.Comparison of the numerical solution with the 4th-order SHAM solution for h −0.85 is made in Figure 1.It  In Table 2, we present that the SHAM results for another value of h in h h −0.9 , are also very highly accurate, and in sixth order we obtain the exact numerical solutions.In HAM we choose y 0 x x 2 − 1 /2 as initial guess.
subject to the boundary conditions y −1 y 1 0 with the exact solution y x x 3 − x.Similar to the previous example, we employ HAM and SHAM to solve this example.In Table 3, there is a comparison of the numerical results against the HAM and SHAM approximation solutions 5.2 at different orders.It is worth noting that the rate of convergence in SHAM results is higher than the HAM results.Comparison of the numerical solution with the 16th-order SHAM solution for h −0.4 is made in Figure 4.It can be seen from the h-curves Figures 5 and 6 that if −0.4 ≤ h ≤ 0 and −0.6 ≤ h ≤ −0.3 then In this example, by a change of variables, the Fredholm-integro differential equation with nonhomogeneous conditions are transformed to the equation with homogeneous conditions and then apply the spectral homotopy analysis method for this problem.According to governing equation and the boundary condition it is reasonable to set the transformation Comparison between absolute errors in solutions by SHAM and Legendre collocation matrix method LCMM is tabulated in Table 4.It is worth noting that the SHAM results become very highly accurate only with four iterations.In this example, by a change of variables, the Fredholm-integro differential equation with nonhomogeneous conditions are transformed to the equation with homogeneous conditions and then we apply the spectral homotopy analysis method for this problem.Comparison between absolute errors in solutions by SHAM and HAM is tabulated in Table 5.It is worth noting that the SHAM results become very highly accurate only with four iterations.

Conclusion
In this paper for the first time in the literature we presented the application of spectral homotopy analysis method SHAM in solving linear integro differential equations.A comparison was made between exact analytic solutions and numerical results from the spectral homotopy analysis method and homotopy analysis method HAM solutions.The numerical results indicate that in SHAM the rate of convergence is faster than HAM.For example, in Example 5.1 we found that the fourth-order SHAM approximation sufficiently gives a match with the numerical results up to twelve decimal places.In contrast, HAM solutions have a good agreement with the numerical results in 11th-order.In this paper we choose the admissible value of h from h-curve for y 1 and y 1 .As it is shown in Figures 2 and 3 the range of admissible h values is much wider in SHAM than in HAM.In using HAM there is restriction on choosing an initial approximation that the higher-order deformation can be integrated whereas in SHAM the initial approximation is the solution of the nonhomogeneous linear part of 3.1 .The spectral homotopy analysis method is more efficient as it does not depend on the rule of coefficient ergodicity unlike the HAM.With some changes we can apply this method on Fredholm integrodifferential equations.In Example 5.3 27 , we compare the SHAM solution with Legendre collocation matrix method LCCM to solve linear Fredholm integrodifferential equation.As can be seen in Table 4 SHAM is more accurate and efficient than LCCM.In this work we obtained the numerical results up to twelve decimal places having well agreement with exact solution but reported them only up to six decimal places.
In this paper, we described the spectral homotopy analysis method to solve linear Volterra and Fredholm integrodifferential equations; however, it remains to be generalized and verified for more complicated integral equations that we consider it as future works.

1 T
j t dt PT j τ , 4.6

y 4 xFigure 1 :
Figure 1: Comparison of the numerical solution with the 4th-order SHAM solution for h −0.85.

y 16 xFigure 4 :
Figure 4: Comparison of the numerical solution with the 16th-order SHAM solution for h −0.4.

Table 1 :
Numerical result of Example 5.1 against the HAM and the SHAM solutions whit h −0.85.

Table 2 :
Numerical result of Example 5.1 against the SHAM solutions whit h −0.9.

Table 3 :
Numerical result of Example 5.2 against the HAM and the SHAM solutions whit h −0.4.

Table 4 :
Absolute errors for Example 5.3.

Table 5 :
Absolute errors for Example 5.4.