Homotopy Approach for Integrodifferential Equations

In this paper, we present the application of the homotopy analysis method for solving integrodifferential equations. In this method, a series is created, the successive elements of which are determined by calculating the appropriate integral of the previous element. In this elaboration, we prove that, if this series is convergent, then its sum is the solution of the objective equation. We formulate and prove the sufficient condition of this convergence, and we give also the estimation of error of an approximate solution obtained by taking the partial sum of the considered series. Moreover, we present in this paper the example of using the investigated method for determining the vibrations of the freely supported reinforced concrete beam as well as for solving the equation of movement of the electromagnet jumper mechanical system.


Introduction
The integrodifferential equations occur in the mathematical models of many real and scientific phenomena. The exact solutions of equations of this kind are usually difficult to find. Therefore, the various approximate methods are applied for finding their solutions. In this paper, we propose the application of the homotopy analysis method for solving these kind of equations.
The homotopy analysis method was developed in the 1990s [1,2]. This method gives the possibility to determine the solution of the wide class of problems described by means of various operator equations. In particular, it can be the ordinary and partial differential equations [3,4], the fractional differential or integrodifferential equations [5,6], and also the integral equations [7,8]. Theoretical results concerning the convergence of the homotopy analysis method in the case of differential equations can be found, among others, in References [2,9,10], whereas, in the case of integral equations, can be found in References [8,11].
The homotopy analysis method was already applied for solving the integrodifferential equations in some papers (see, for example, References [12,13]). In most of these papers, presented are only the examples of applications without any theoretical results concerning the convergence of the method. The theoretical result, the most frequently cited in the other papers, is the classical theorem stating that, if the series created in the homotopy analysis method is convergent, then its sum is the solution of the input equation. There are only two papers [14,15] that contain also some other theoretical results.
In Reference [14], the theorem defining the sufficient conditions for convergence and divergence of the created series is formulated. Next, in Reference [15], the authors used the homotopy analysis method combined with the parametrization methods to solve the optimal control problems governed by the integrodifferential equations. This paper includes also the sufficient condition for convergence of the obtained series (in addition to the classical theorem mentioned before).
In the current paper, except the classical theorem, we also formulate and prove the sufficient condition for convergence of the created series. Moreover, we give the estimation of error of the approximate solution obtained by reducing the sum of the series to the partial sum. Additionally, we present the example of using the examined method for determining the vibrations of the freely supported reinforced concrete beam as well as for solving the equation of movement of the electromagnet jumper mechanical system.

Homotopy Analysis Method
The homotopy analysis method serves to solve the following operator equations: where N denotes the operator (in particular, it may be the nonlinear operator), while u is the unknown function. In the first step, we specify the homotopy operator H as given below: where p ∈ [0, 1] is the embedding parameter, h = 0 is the convergence control parameter [2,9,16,17], u 0 denotes an initial approximation of the solution of Equation (1), and L is the auxiliary linear operator with property L(0) = 0. Considering the equation H(Φ, p) = 0, we obtain the so-called zero-order deformation equation: Substituting p = 0, we get L(Φ(x; 0) − u 0 (x)) = 0, so we have Φ(x; 0) = u 0 (x). However, if we assume p = 1, we obtain N(Φ(x; 1)) = 0. Therefore, Φ(x; 1) = u(x), where u is the desired solution of Equation (1). Thus, the change of parameter p from zero to one involves a change from the trivial problem to the original problem (and, thus, the solution from u 0 to u). If the operator N has a unique kernel, then Equation (1) has a single solution, which could be determined by using the described method. If the kernel is not unique, the method allows to determine one of the existing solutions of Equation (1).
Taking the Maclaurin series of function Φ(x; p) (with respect to parameter p), we get where u 0 (x) = Φ(x; 0) and If the above series is convergent for p = 1, then we obtain the required solution: In order to determine the function u m , we differentiate m times with respect to parameter p the left-and the right-hand sides of dependence of Equation (3). Then, we divide the result by m! and we substitute p = 0, obtaining the so-called mth-order deformation equation (m > 0): where u m−1 = {u 0 (x), u 1 (x), . . . , u m−1 (x)}, and and If we are not able to determine the sum of series in Equation (6), then we can accept the partial sum of this series as the approximate solution of considered equation: The proper selection of the convergence control parameter h affects the area of convergence of the created series as well as the rate of this convergence [2,10,18]. In order to determine the value of this parameter, we can use the so-called "optimization method" [2], in which we define the following squared residual of governing equation: Next, the optimal value of the convergence control parameter is obtained by minimizing the above squared residual.

Integrodifferential Equations
We consider the integrodifferential equations of the form with condition u(a) = c 0 , where x ∈ [a, b], c 0 ∈ R, and where, for k = 1, 2, and F ∈ C[a, b], whereas the function u is sought. In this elaboration, we will use the norm in space C 0 (Ω), where Ω is the compact set, defined as follows: v = sup x∈Ω |v(x)| and the norm in C 1 (Ω) is defined as follows: We assume Ω = [a, b] or Ω = [a, b] × [a, b]. The norm in space C 1 (Ω) can be also defined as given below: where w is a fixed element of set Ω.
Integrating both sides of Equation (12), we get the equivalent form of this equation: where we used the condition u(a) = c 0 and we assumed that The above equation gives the ground for defining operators L and N: By setting u 0 ∈ C[a, b] and by applying the homotopy analysis method, we get the following formula for functions u m : where χ m and operator R m are determined by Equations (8) and (9).
In the considered case, the operator R m takes the following form (under assumption that the series is convergent, which will be discussed later): Hence, we obtain the following formulae for function u m : and for m 2: Now, we intend to present and to prove the following three theorems describing the application of the homotopy analysis method for solving linear integrodifferential equations. Theorem 1. Let the functions u m and m ≥ 1, be defined by Equations (17) and (18). If the series in Equation (6) is uniformly convergent, then the sum of this series is the solution of Equation (12).
Proof. Let us assume that the series in Equation (6) is convergent. For any x ∈ [a, b], according to the necessary condition for convergence of the series, we have Using the definition of operator L, we can write By using the convergence of series in Equation (6) and the continuity of respective operators, we get Thus, we obtain This means that the function u(x) = ∑ ∞ m=0 u m (x) satisfies Equation (13), and therefore, it satisfies Equation (12).
In the next theorem, we prove the sufficient condition for convergence of the series created in the discussed method.

Theorem 2. If the inequality
is satisfied, then with a suitable choice of the control parameter h, the series in Equation (6) It can be shown by induction that, for m 1 and Thus, for the considered series in Equation (6), we get In the above estimation, the last series is the geometric series with common ratio γ h . Thus, if only We consider the above inequality in three cases: for h < −1, for h ∈ [−1, 0) and for h > 0. In the first case, we get the conditions In the second case, we have A < 1, whereas the third case leads to contradiction. It means that, if the inequality of Equation (19) is fulfilled, then we can select such a value of the parameter controlling convergence (it would be enough to take any h ∈ ( −2 1+A , 0)) that γ h < 1.
The estimations introduced in the previous proof will be also used in the proof of the following theorem. (19) is fulfilled and n ∈ N, then for the approximate solution u n , we have the following error estimation:

Theorem 3. If the inequality
for all x ∈ [a, b]. Thus, we have where Now, we proceed to the proof of the second inequality. For this purpose, we apply the homotopy analysis method for Equation (12) by defining the operators L and N in the following way: Then, we get the following relations for the derivatives of functions u m : and for m 2: Now, we determine the boundaries of functions u m : Hence, we get in general for m 1 By using the above estimations, we have for any which proves the inequality of Equation (21). The last inequality follows from the definition of the taken norm: and from the inequalities already proven.

Examples
Example 1. We start with the solution of the following Volterra-Fredholm equation: with condition u(0) = 1. In this equation, we have K 1 (x, t) = − 1 2 (t − x) 2 , R 1 (u) = u, R 2 (u) = −  (14) is then the function u e (x) = x 2 + 1. Let us mention that all the calculations presented here were carried out with the aid of Mathematica software.
In the considered equation, we have Hence, we get which means that the sum of the created series is the solution of the discussed equation.
Assuming the initial approximation u 0 (x) = F(x), in the first step, we get The optimal value of the convergence control parameter equals −0.9858834. Figure 1 displays the graph of logarithm of the squared residual for n = 9. Table 1 compiles the errors of reconstruction of the exact solution for the successive approximate solutions u n , n ∈ {1, 2, . . . , 10}. As revealed by the above results, together with an increase of the number of components in Equation (10), the errors quickly decrease. The errors decrease most quickly for the optimal value of parameter h. For this value, the approximate solution u 5 provides the approximation of the sought function with the error not higher than 8.363 · 10 −8 , whereas the approximate solution u 10 gives the error not higher than 4.4409 · 10 −16 . Also, for example, in the case of h = −1, we obtain, respectively, 4.7787 · 10 −7 and 3.2792 · 10 −14 . Table 1 contains also the estimation of error resulting from the inequalityof Equation (22). In the considered example, we have γ h = γ h = 0.753529, u 1 1 17.3483 for the optimal value of convergence control parameter h = −0.9858834. Equation (22) gives the estimation of error of the approximate solution (the worst possible case). In fact, the errors of approximate solutions are, in general, much smaller than the value determined in the right-hand side of this inequality. In Figure 2, the absolute error |u e (x) − u n (x)| for n = 4 and n = 7 is plotted. The obtained results show that the investigated method is very quickly convergent and that the calculation of only a few first terms of the series ensures very good approximation of the exact solution.
The discussed beam is considered as the set of segments of length dx, and the kinetic energy of such an elementary segment, that is, the element of mass m dx (see Figure 3a) [19,20], is expressed as follows: whereẏ(x) denotes the velocity of point of coordinate x (in the beam). In case of the freely supported beam wrapped by the mass m , the deflection curve y(x) is symmetrical with respect to the centre of gravity and, for 0 x 1 2 , is expressed by the following relation: where y m denotes the largest deflection of the beam (see Figure 3a). The velocity of each point of coordinate x is described by the following formula: whereẏ m = dy m dt and the kinetic energy of the whole beam is equal to By substituting Equation (26)  attached in the vibration arrow, that is, at the point where the displacement of beam is the highest (see Figure 3b).
Stiffness of the substitute system is expressed by the following formula: The equation describing the first free vibration of the reinforced concrete beam of span l = 6 m and cross section (b/h) 0.30 m/0.60 m has then the following form: where m = 1633 kg is the beam substitute mass and k = 32.4 MN/m describes the beam stiffness in the center of span (see Figure 4). Equation (31), after dividing by m, takes the form where ω 0 = k m = 140.9 rad/s denotes the circular frequency of undamped vibrations. The above equation must be completed by the initial conditions y(0) = y 0 ,ẏ(0) = v 0 . (33) Now, we will demonstrate how to apply the homotopy analysis method to solve Equation (32) with the conditions of Equation (33). For this purpose, let us integrate Equation (32) in limits from 0 to t, obtaininġ with condition y(0) = y 0 . For the initial conditions y 0 = 0.01 and v 0 = 0, the exact solution takes the form y e (t) = 0.01 cos(140.9 t).
Equation (34) can be solved by using the procedure described in this paper. If we take the function u 0 (t) = 0 as the initial approximation, then in the first steps of the method, we get . . . Figure 5 displays the logarithm of squared residual for n = 9. The optimal value of the convergence control parameter, determined numerically, is equal to −0.94423256. In Table 2, there are collected the errors (in norm · 1 ) of reconstructing the exact solution by the successive approximate solutions u n , n ∈ {1, 2, . . . , 20}. At the beginning, that is, until n = 3, the errors increase, but after that, the errors start quickly decrease. So, for n = 10, the errors are at the level of 10 −4 ; for n = 15, they are at the level of 10 −10 ; and for n = 20, the errors are at the level of 10 −14 . The main part of these errors is generated by the approximation of derivative.
For n = 4, the error in norm · 1 is equal to 43.243, whereas in norm · , it is only equal to 0.351. A similar situation can be observed in the other cases. As revealed by the above results, starting from some moment, we can observe that an increase of the number of components in Equation (10) causes the errors to quickly decrease (see also Figure 6).
The investigated equation shows that the homotopy analysis method can be successfully applied also for solving the equations not satisfying the condition in Equation (19). In the above considered case, we have A = 39.478 because the condition of Equation (19) is only the sufficient condition and it is not the necessary condition for convergence of the method.

Example 3.
The equation of movement of the electromagnet jumper mechanical system can be presented in the following way: where m is the mass, b 1 denotes the damping, y describes the displacement of the electromagnet jumper, and f denotes the force occurring in consequence to the electrical system action. In this example, we take f (y, t) = f 1 − c 1 y(t). The third term in the left-hand side of the equation is the result of taking into account in the mechanical system the nonlinear spring, which generates the force component depending on the displacement. The above differential equation is completed by the following initial conditions: Now: we present the application of the homotopy analysis method for solving Equation (32) with the conditions of Equation (33). Equation (35), after integration in limits from 0 to t and some transformations, can be written in the following form: where f r (t) is the function depending on the initial conditions. In the considered example, we take m = 1, b 1 = 2/10, k 1 = 1/1000, f 1 = 10, c 1 = 10, v 0 = 1/10, y 0 = 0. In order to test the usage of investigated method, we select the function f r (t) = 1/10 + t/50 + t 2 /2 + t 3 /60000 so that the exact solution of above equation is known (y(t) = t/10).
By applying the homotopy analysis method for solving Equation (38) and by taking u 0 (t) = f r (t), we obtain in the first step The optimal value of the convergence control parameter, determined numerically, is equal to −0.99. In Table 3, there are compiled the errors (in norm · 1 ) of reconstructing the exact solution by the successive approximate solutions u n , n ∈ {1, 2, . . . , 20}. As revealed by the above results, starting from some moment, the errors quickly decrease together with an increase of the number of components in Equation (10). Next, in Figure 7, the absolute error y e − u n is plotted for n = 7 and n = 12. Again, we may conclude, on the basis of the obtained results, that the method is very quickly convergent to the exact solution.

Conclusions
This paper presents the application of the homotopy analysis method for solving linear integrodifferential equations. In this method, we create a series, the successive elements of which are calculated by integrating properly the previous elements. In this paper, it is proven that, if the created series is convergent, then its sum is the solution of the input equation. The sufficient condition of this convergence is also formulated and proven. Moreover, the estimation of error of approximate solution, obtained by reducing the considered series to its partial sum, is also given.
The series created with the discussed method is, in general, very quickly convergent. After selecting the optimal value of the convergence control parameter, calculation of only a few or several first terms of the series ensures very good approximation of the sought solution. The executed calculations show that the homotopy analysis method is effective in solving the problems under consideration.
The results obtained in this paper can be easily generalized to the case of equations having the following form: where all the functions are defined analogously to Equation (12). In this case, the condition of Equation (19) should be replaced by the condition in the following form: The homotopy analysis method can be also applied for solving the equation of the following form: in a similar way as described above, with the appropriate initial conditions. One has to only integrate l times the given equation and to repeat the procedure described in this paper.
It is also possible to find another formulation of the homotopy operator [2], namely where H is an auxiliary function. Also in this case, all the results obtained and presented in this work are true. The only change concerns the formula for constant γ h . Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.