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

Authors: Z. Pashazadeh Atabakan; A. Kılıçman; A. Kazemi Nasab
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

A modification of homotopy analysis method (HAM) known as spectral homotopy analysis method (SHAM) is proposed to solve linear Volterra integrodifferential equations. Some examples are given in order to test the efficiency and the accuracy of the proposed method. The SHAM results show that the proposed approach is quite reasonable when compared to SHAM results and exact solutions.


---

## Body

## 1. 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(1.1)a1(x)u′′(x)+a2(x)u′(x)+a3(x)u(x)=g(x)+∫-1x‍k(x,t)u(t)dt,u(-1)=0,u(1)=0.The paper is organized in the following way. Section2 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 Sm-1 which is needed to obtain Ym in Section 4. Numerical examples are presented in Section 5. In Section 6, concluding remarks are given.


## 2. Homotopy Analysis Solution

In this section, we give a brief introduction to HAM. We consider the following differential equation in a general form:(2.1)N[u(τ)]=0,
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:
(2.2)(1-p)L[ϕ(τ;p)-u0(τ)]=phH(τ)(N[ϕ(τ;p)]),
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:
(2.3)L(ϕ(x;p))=a1(x)∂2ϕ(x;p)∂x2,
with the property L(c1t+c2)=0 where c1 and c2 are constants and
(2.4)N[ϕ(x;p)]=a1(x)∂2ϕ(x;p)∂x2+a2(x)∂ϕ(x;p)∂x+a3(x)ϕ(x;p)-g(x)-∫-1x‍k(x,t)ϕ(t)dt
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:
(2.5)ϕ(τ;p)=u0(τ)+∑m=1+∞‍um(τ)pm,
where
(2.6)um(τ)=Dm[ϕ(τ;p)],(2.7)Dmϕ=1m!∂m(ϕ)∂pm|p=0.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:
(2.8)ϕ(τ)=u0(τ)+∑m=1+∞‍um(τ).The admissible and valid values of the convergence-parameterh 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
(2.9)(1-p)L[ϕ(τ;p)-u0(τ)]=ph(N[ϕ(τ;p)]).Define the vectorum={u0(τ),u1(τ),…,um(τ)}. Operating on both sides with Dm, we have the so called mth-order deformation equation
(2.10)L[um(τ)-χmum-1(τ)]=hH(τ)Rm(um-1(τ)),
where
(2.11)Rm(um-1)=1(m-1)!∂m-1N[ϕ(τ;p)]∂pm-1|p=0,χm={0,m≤11,otherwise.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(3.1)a1(x)u′′(x)+a2(x)u′(x)+a3(x)u(x)=g(x)+∫-1x‍k(x,t)u(t)dt,u(-1)=0,u(1)=0.To obtain initial approximation we solve the following two-point boundary value problem:(3.2)a1(x)u0′′(x)+a2(x)u0′(x)+a3(x)u0(x)=g(x),
subject to boundary conditions
(3.3)u0(-1)=0,u0(1)=0.
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
(3.4)u0(τ)≈u0N(τj)=∑k=0N‍u^kTk(τj),j=0,…,N,
where Tk 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
(3.5)τj=cos(πjN).Derivatives of the functionsy0(τ) at the collocation points are represented as
(3.6)dsu0(τk)dτs=∑j=0N‍Dkjsu0(τj),k=0,…,N,
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.7)Dkj=(-12)ckcj(-1)k+jsin(π(j+k)/2N)sin(π(j-k)/2N),j≠k,Dkj=(-12)xksin2(πk/N),k≠0,N,k=j,D00=-DNN=2N2+16.
Substituting (3.4)–(3.6) into (3.2) results in
(3.8)AU0=G
subject to the boundary conditions
(3.9)u0(τ0)=u0(τN)=0,
where
(3.10)A=a1D2+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.To satisfy the boundary conditions we delete the first and the last rows and columns ofA and the first and the last rows of Y0 and G. The values of u0(τi), i=0,…,N are determined from the equation
(3.11)U0=A-1G,
which is the initial approximation for the SHAM solution of the governing equation (3.1). The zeroth-order deformation equation is given by
(3.12)(1-p)L[ϕ(τ;p)u0(τ)]=ph(N[ϕ(τ;p)]-g(τ)),
where
(3.13)L[ϕ(τ;p)]=a1∂2ϕ∂τ2+a2∂ϕ∂τ+a3ϕ,h is the nonzero convergence controlling auxiliary parameter, and N is a nonlinear operator given by
(3.14)N[ϕ(τ;p)]=a1∂2ϕ∂τ2+a2∂ϕ∂τ+a3ϕ-∫-1τ‍k(τ,t)ϕ(t)dt.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
(3.15)L[um(τ)-χmum-1(τ)]=hRm
subject to boundary conditions
(3.16)um(-1)=um(1)=0,
where
(3.17)Rm(τ)=a1um-1′′+a2u′m-1+a3um-1-g(τ)(1-χm)-∫-1τ‍k(τ,t)um-1(t)dt.Applying the Chebyshev pseudospectral transformation on (3.15)–(3.17) results in
(3.18)AUm=(χm+h)AUm-1-h[Sm-1+(1-χm)G],
subject to the boundary conditions
(3.19)um(τ0)=um(τN)=0,
where A and G are as defined in (3.10) and
(3.20)Um=[um(τ0),um(τ1),…,um(τN)]T,Sm-1=∫-1τ‍k(τ,t)Um-1dt.
To implement the boundary conditions we delete the first and the last rows of Sm-1 and the first and the last rows and columns of A.Finally, this recursive formula can be written as following form≥1:
(3.21)Um=(χm+h)A-1AUm-1-hA-1[Sm-1+(1-χm)G],
so, starting from the initial approximation. We can obtain higher-order approximation Um for m≥1, recursively.


## 4. Calculation ofSm-1

Consider the well-known Chebyshev polynomials of first kind of degreen, defined by [25]
(4.1)Tn(x)=cos(ncos-1(x)),n≥0.Also they are derived from the following recursive formula:(4.2)T0(x)=1,T1(x)=x,Tn+1(x)=2xTn(x)-Tn-1(x),n=1,2,3,…In this section we approximate the integrand in (3.15) by letting
(4.3)k(τm,t,ym-1(t))=∑j=0m-1‍Tj(t)wmj,t=τ0,τ1,…,τm-1,m=1,2,…,N,
where Tj(t),  j=0,1,…,m-1 are the basis functions.Clearly,wmj,  0≤j≤m-1 can be obtained by solving the following matrix equation:
(4.4)(T0(τ0)T1(τ0)…Tm-1(τ0)T0(τ1)T1(τ1)…Tm-1(τ1)⋮⋮⋮T0(τm-1)T1(τm-1)…Tm-1(τm-1))(wm,0wm,1⋮wm,m-1)=(k(tm,τ0,y0)k(tm,τ1,y1)⋮k(tm,τm-1,ym-1)),
where τj=tj, j=0,1,…,m-1 are Chebyshev collocation points defined in (2.6). We can obtain wmj by solving the above system so the integral in (3.15) can be rewritten as following:
(4.5)∫-1τ‍k(τ,t)ym-1(t)dt=∑j=0m-1‍am,jwm,j,
where
(4.6)am,j=∫-1τ‍Tj(t)dt=PTj(τ),
where P as in [26], that is the (N+1)×(N+1) operational matrix for integration as follows:
(4.7)P=(1100…00-140-140…00-13-12016…00180-140…00⋮⋮⋮⋮⋱⋮⋮(-1)N(N-1)2-1000…012N(-1)N+1N2-1000…-12N-20).


## 5. Numerical Examples

In this section we apply the technique described in Section3 to some illustrative examples of second-order Volterra integrodifferential equations.Example 5.1.
Consider the second-order Volterra integrodifferential equation(5.1)y′′(x)-xy′(x)+xy(x)=x46+x3-2x2-x3+52+∫-1x‍(x-2t)y(t)dt
subject to the boundary conditions y(-1)=y(1)=0 with the exact solution y(x)=x2-1.We employ HAM and SHAM to solve this example. In Table1, 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 can be seen from the h-curves (Figures 2 and 3) that if -1≤h≤-0.2 and -1≤h≤-0.4 then the SHAM an HAM approximate solutions are in well agreement with the exact solution. 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 y0(x)=(x2-1)/2 as initial guess.Table 1
Numerical result of Example5.1 against the HAM and the SHAM solutions whit h=-0.85.
x
SHAM
HAM
Numerical
2nd order
4th order
5th order
2nd order
5th order
11th order
−0.8
−0.360376
−0.360001
−0.36000000
−0.370525
−0.360086
−0.36000000
−0.36000000
−0.6
−0.640699
−0.640003
−0.64000000
−0.658320
−0.640123
−0.64000000
−0.64000000
−0.4
−0.840976
−0.840005
−0.84000000
−0.860438
−0.840105
−0.84000000
−0.84000000
−0.2
−0.961198
−0.960005
−0.96000000
−0.979684
−0.960078
−0.96000000
−0.96000000
0.0
−1.001354
−1.000005
−1.00000000
−1.018938
−1.000066
−1.00000000
−1.00000000
0.2
−0.961421
−0.960005
−0.96000000
−0.978828
−0.960067
−0.96000000
−0.96000000
0.4
−0.841370
−0.840005
−0.84000000
−0.857858
−0.840073
−0.84000000
−0.84000000
0.6
−0.641159
−0.640004
−0.64000000
−0.654041
−0.64048
−0.64000000
−0.64000000
0.8
−0.360727
−0.360002
−0.36000000
−0.367054
−0.360011
−0.36000000
−0.36000000Table 2
Numerical result of Example5.1 against the SHAM solutions whit h=-0.9.
x
5th order
6th order
7th order
Numerical
-
0.8
-
0.35999595
-
0.36000000
-
0.36000000
-
0.36000000
-
0.6
-
0.63999432
-
0.64000000
-
0.64000000
-
0.64000000
-
0.4
-
0.83999202
-
0.84000000
-
0.84000000
-
0.84000000
-
0.2
-
0.95999009
-
0.96000000
-
0.96000000
-
0.96000000
0.0
-
0.99998868
-
1.00000000
-
1.00000000
-
1.00000000
0.2
-
0.95998890
-
0.96000000
-
0.96000000
-
0.96000000
0.4
-
0.83998833
-
0.84000000
-
0.84000000
-
0.84000000
0.6
-
0.63999070
-
0.64000000
-
0.64000000
-
0.64000000
0.8
-
0.35999373
-
0.36000000
-
0.36000000
-
0.36000000Figure 1
Comparison of the numerical solution with the 4th-order SHAM solution forh=-0.85.Figure 2
4th-order SHAMh-curve.Figure 3
11th-order HAMh-curve.Example 5.2.
Consider the second-order Volterra integrodifferential equation(5.2)y′′(x)+y(x)=-4x6+6x4+x3+5x-2+∫-1x‍24t2y(t)dt
subject to the boundary conditions y(-1)=y(1)=0 with the exact solution y(x)=x3-x.Similar to the previous example, we employ HAM and SHAM to solve this example. In Table3, 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 the approximate solution is in well agreement with the exact solution. In HAM, we choose y0(x)=-x2+1 as initial guess.Table 3
Numerical result of Example5.2 against the HAM and the SHAM solutions whit h=-0.4.
x
SHAM
Numerical
HAM
8th order
14th order
16th order
8th order
23th order
-
0.8
0.28818730
0.28800272
0.28800000
0.28800000
0.28786022
0.28800097
-
0.6
0.38426358
0.38400174
0.38400000
0.38400000
0.38548670
0.3840003
-
0.4
0.33618044
0.33599640
0.33600000
0.33600000
0.34009630
0.33600088
-
0.2
0.19199056
0.19198946
0.19200000
0.19200000
0.19905046
0.1920014
0.0
-
0.00024441
-
0.00000165
0.00000000
0.00000000
-
0.00041662
0.00000166
0.2
-
0.19249771
-
0.19202092
-
0.19200000
-
0.19200000
-
0.18206491
-
0.19209884
0.4
-
0.33673440
-
0.33602266
-
0.33600000
-
0.33600000
-
0.32825115
-
0.33600119
0.6
-
0.36072793
-
0.38401942
-
0.38400000
-
0.38400000
-
0.38049994
-
0.38400153
0.8
-
0.28863110
-
0.28801009
-
0.28800000
-
0.28800000
-
0.28821582
-
0.28800159Figure 4
Comparison of the numerical solution with the 16th-order SHAM solution forh=-0.4.Figure 5
16th-order SHAMh-curve3.Figure 6
23th-order HAMh-curve.Example 5.3.
Consider the first-order Fredholm integrodifferential equation [27]
(5.3)y′(x)-y(x)=sin(4πx)-cos(2πx)-2πsin(2πx)-∫-11‍sin(4πx+2πt)y(t)dt
subject to the boundary conditions y(-1)=y(1)=1.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(5.4)y(x)=cos(2πx)+Y(x)
to make the governing boundary conditions homogeneous. Substituting (5.4) in (5.3) and boundary conditions gives
(5.5)Y′(x)-Y(x)=∫-11‍sin(4πx+2πt)Y(t)dt
subject to the boundary condition
(5.6)Y(-1)=Y(1)=0.Comparison between absolute errors in solutions by SHAM and Legendre collocation matrix method (LCMM) is tabulated in Table4. It is worth noting that the SHAM results become very highly accurate only with four iterations.Table 4
Absolute errors for Example5.3.
x
SHAM
LCMM
-
1.0
3.30
×
10−19
2.49
×
100
-
0.8
1.73
×
10−9
3.56
×
10−1
-
0.6
6.69
×
10−9
2.66
×
10−1
-
0.4
2.91
×
10−10
5.85
×
10−2
-
0.2
6.93
×
10−9
1.74
×
10−2
0.0
9.16
×
10−9
0.00
0.2
6.93
×
10−9
1.97
×
10−2
0.4
2.91
×
10−10
7.50
×
10−2
0.6
6.69
×
10−9
2.59
×
10−2
0.8
1.73
×
10−9
4.16
×
10−1
1.0
6.38
×
10−17
2.58
×
100Example 5.4.
Consider the following second-order Fredholm integrodifferential equation [28]:
(5.7)y′′(x)+xy′(x)-xy(x)=ex-2sin(x)+∫-11‍sin(x)e-ty(t)dt,
where y(-1)=e-1,  y(1)=e.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. According to governing equation and the boundary condition it is reasonable to set the transformation(5.8)y(x)=ex+Y(x)
to make the governing boundary conditions homogeneous. Substituting (5.8) in (5.7) and boundary conditions gives
(5.9)Y′′(x)+xY′(x)-xY(x)=∫-11‍sin(x)e-tY(t)dt
subject to boundary condition
(5.10)Y(-1)=Y(1)=0.Comparison between absolute errors in solutions by SHAM and HAM is tabulated in Table5. It is worth noting that the SHAM results become very highly accurate only with four iterations.Table 5
Absolute errors for Example5.4.
x
SHAM
HAM
4th order
20th order
0.0
0
0
0.2
0
2.2760
×
10−14
0.4
0
3.5971
×
10−13
0.6
1
×
10−19
1.8209
×
10−12
0.8
1
×
10−19
5.7547
×
10−12
1.0
0
1.4550
×
10−11
1.2
1.2241
×
10−15
2.9134
×
10−11
1.4
1.1239
×
10−13
5.3974
×
10−11


## 6. 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 Example5.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.



---

*Source: 960289-2012-11-05.xml*