Next Article in Journal
Time-Truncated Group Plan under a Weibull Distribution based on Neutrosophic Statistics
Next Article in Special Issue
Blended Root Finding Algorithm Outperforms Bisection and Regula Falsi Algorithms
Previous Article in Journal
On the Non-Hypercyclicity of Normal Operators, Their Exponentials, and Symmetric Operators
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Homotopy Approach for Integrodifferential Equations

1
Institute of Mathematics, Faculty of Applied Mathematics, Silesian University of Technology, 44-100 Gliwice, Poland
2
Department of Building Structures, Faculty of Civil Engineering, Silesian University of Technology, 44-100 Gliwice, Poland
3
Department of Mechatronics, Faculty of Electrical Engineering, Silesian University of Technology, 44-100 Gliwice, Poland
*
Author to whom correspondence should be addressed.
Mathematics 2019, 7(10), 904; https://doi.org/10.3390/math7100904
Submission received: 2 September 2019 / Revised: 23 September 2019 / Accepted: 25 September 2019 / Published: 27 September 2019
(This article belongs to the Special Issue Numerical Modeling and Analysis)

Abstract

:
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.

1. 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.

2. Homotopy Analysis Method

The homotopy analysis method serves to solve the following operator equations:
N ( u ( x ) ) = 0 , x Ω ,
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:
H ( Φ , p ) ( 1 p ) L Φ ( x ; p ) u 0 ( x ) p h N Φ ( x ; p ) ,
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:
( 1 p ) L Φ ( x ; p ) u 0 ( x ) = p h N Φ ( x ; p ) .
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
Φ ( x ; p ) = m = 0 u m ( x ) p m ,
where u 0 ( x ) = Φ ( x ; 0 ) and
u m ( x ) = 1 m ! m Φ ( x ; p ) p m | p = 0 , m = 1 , 2 , 3 ,
If the above series is convergent for p = 1 , then we obtain the required solution:
u ( x ) = m = 0 u m ( x ) .
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 ):
L u m ( x ) χ m u m 1 ( x ) = h R ¯ m u ¯ m 1 , x ,
where u ¯ m 1 = { u 0 ( x ) , u 1 ( x ) , , u m 1 ( x ) } , and
χ m = 0 m 1 , 1 m > 1
and
R ¯ m u ¯ m 1 , x = 1 ( m 1 ) ! m 1 p m 1 N i = 0 u i ( x ) p i | p = 0 .
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:
u ^ n ( x ) = m = 0 n u m ( x ) .
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:
E n ( h ) = Ω N u ^ n ( x ) 2 d x .
Next, the optimal value of the convergence control parameter is obtained by minimizing the above squared residual.

3. Integrodifferential Equations

We consider the integrodifferential equations of the form
u ( x ) f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u ( t ) ) d t R 2 ( u ( x ) ) = F ( x ) ,
with condition u ( a ) = c 0 , where x [ a , b ] , c 0 R , and where, for k = 1 , 2 , R k : C [ a , b ] C [ a , b ] is the bounded linear operators f , g C [ a , b ] , a f ( x ) g ( x ) b , K 1 C ( [ a , b ] × [ a , b ] ) 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:
v 1 = sup x Ω | v ( x ) | + sup x Ω | v ( x ) | .
We assume Ω = [ a , b ] or Ω = [ a , b ] × [ a , b ] . The norm in space C 1 ( Ω ) can be also defined as given below:
v 1 = | v ( w ) | + sup x Ω | v ( x ) | ,
where w is a fixed element of set Ω .
Integrating both sides of Equation (12), we get the equivalent form of this equation:
u ( x ) = c ¯ 0 ( x ) + a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u ( t ) ) d t d x + a x R 2 ( u ( x ) ) d x ,
where we used the condition u ( a ) = c 0 and we assumed that
c ¯ 0 ( x ) = c 0 + a x F ( x ) d x .
The above equation gives the ground for defining operators L and N:
L ( v ) = v ,
N ( v ) = v ( x ) c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( v ( t ) ) d t d x a x R 2 ( v u ( x ) ) d x .
By setting u 0 C [ a , b ] and by applying the homotopy analysis method, we get the following formula for functions u m :
u m ( x ) = χ m u m 1 ( x ) + h R ¯ m ( u ¯ m 1 , x ) ,
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):
R ¯ m ( u ¯ m 1 , x ) = 1 ( m 1 ) ! m 1 p m 1 N i = 0 u i ( x ) p i p = 0 = = 1 ( m 1 ) ! m 1 p m 1 [ i = 0 u i ( x ) p i c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 i = 1 u i ( t ) p i d t d x a x R 2 i = 1 u i ( x ) p i d x ] p = 0 = = 1 ( m 1 ) ! m 1 p m 1 [ i = 0 u i ( x ) p i c ¯ 0 ( x ) i = 0 a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u i ( t ) ) p i d t d x i = 0 a x R 2 ( u i ( x ) ) p i d x ] p = 0 = = 1 ( m 1 ) ! ( ( m 1 ) ! u m 1 ( x ) ( 1 χ m ) c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) ( m 1 ) ! R 1 ( u m 1 ( t ) ) d t d x ( m 1 ) ! a x R 2 ( u m 1 ( x ) ) d x ) = = u m 1 ( x ) 1 χ m ( m 1 ) ! c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u m 1 ( t ) ) d t d x a x R 2 ( u m 1 ( x ) ) d x .
Hence, we obtain the following formulae for function u m :
u 1 ( x ) = h u 0 ( x ) c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u 0 ( t ) ) d t d x a x R 2 ( u 0 ( x ) ) d x
and for m 2 :
u m ( x ) = ( 1 + h ) u m 1 ( x ) h a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u m 1 ( t ) ) d t d x + a x R 2 ( u m 1 ( x ) ) d x .
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
lim m u m ( x ) = 0 .
Using the definition of operator L, we can write
m = 1 n L u m ( x ) χ m u m 1 ( x ) = m = 1 n u m ( x ) χ m u m 1 ( x ) = = u 1 ( x ) + ( u 2 ( x ) u 1 ( x ) ) + ( u 3 ( x ) u 2 ( x ) ) + + ( u n ( x ) u n 1 ( x ) ) = u n ( x ) .
Hence, we have
m = 1 L u m ( x ) χ m u m 1 ( x ) = lim n u n ( x ) = 0 .
From Equation (7), we get
h m = 1 R ¯ m u ¯ m 1 , x = m = 1 L u m ( x ) χ m u m 1 ( x ) .
And because h 0 , we obtain
m = 1 R ¯ m u ¯ m 1 , x = 0 .
By using the convergence of series in Equation (6) and the continuity of respective operators, we get
0 = m = 1 R ¯ m u ¯ m 1 , x = m = 1 ( u m 1 ( x ) 1 χ m ( m 1 ) ! c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u m 1 ( t ) ) d t d x a x R 2 ( u m 1 ( x ) ) d x ) = = m = 1 u m 1 ( x ) m = 1 1 χ m ( m 1 ) ! c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 m = 1 u m 1 ( t ) d t d x a x R 2 m = 1 u m 1 ( x ) d x = = u ( x ) c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u ( t ) ) d t d x a x R 2 ( u ( x ) ) d x .
Thus, we obtain
u ( x ) = c ¯ 0 ( x ) + a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u ( t ) ) d t d x + a x R 2 ( u ( x ) ) d x .
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
A : = K 1 R 1 ( b a ) 2 + R 2 ( b a ) < 1 ,
is satisfied, then with a suitable choice of the control parameter h, the series in Equation (6) is uniformly convergent in the interval [ a , b ] .
Proof. 
Let u 0 be the function of class C 1 [ a , b ] . Let us find the boundaries of functions u m in interval [ a , b ] :
| u 0 ( x ) | u 0 u 0 1 , | u 1 ( x ) | = h u 0 ( x ) c ¯ 0 ( x ) a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u 0 ( t ) ) d t d x a x R 2 ( u 0 ( x ) ) d x | h | | u 0 ( x ) | + | c ¯ 0 ( x ) | + a x f ( x ) g ( x ) | K 1 ( x , t ) | | R 1 ( u 0 ( t ) ) | d t d x + a x | R 2 ( u 0 ( x ) ) | d x | h | c ¯ 0 + u 0 + K 1 R 1 ( b a ) 2 + R 2 ( b a ) u 0 | h | c ¯ 0 + ( 1 + A ) u 0 , | u 2 ( x ) | = ( 1 + h ) u 1 ( x ) h a x f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u 1 ( t ) ) d t d x + a x R 2 ( u 1 ( x ) ) d x | 1 + h | + | h | K 1 R 1 ( b a ) 2 + R 2 ( b a ) u 1 | 1 + h | + | h | A u 1 = γ h u 1 ,
where
γ h : = | 1 + h | + | h | A .
It can be shown by induction that, for m 1 and x [ a , b ] , we have
| u m ( x ) | γ h m 1 u 1 γ h m 1 u 1 1 .
Thus, for the considered series in Equation (6), we get
m = 0 u m ( x ) m = 0 | u m ( x ) | u 0 1 + u 1 1 m = 1 γ h m 1 .
In the above estimation, the last series is the geometric series with common ratio γ h . Thus, if only γ h < 1 (of course we have γ h > 0 ), then based on the comparison test, we state that the investigated series is uniformly convergent in the interval [ a , b ] .
Let us investigate whether one can choose such value of parameter h, that γ h < 1 , that is
| 1 + h | + | h | A < 1 .
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
h > 2 1 + A a n d A < 1 ,
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.
Theorem 3.
If the inequality (19) is fulfilled and n N , then for the approximate solution u ^ n , we have the following error estimation:
| u ( x ) u ^ n ( x ) | γ h n 1 γ h u 1 γ h n 1 γ h u 1 1 ,
| u ( x ) u ^ n ( x ) | γ ¯ h n 1 γ ¯ h u 1 1
for all x [ a , b ] .
Thus, we have
u u ^ n 1 γ h n 1 γ h + γ ¯ h n 1 γ ¯ h u 1 1 ,
where γ h = | 1 + h | + | h | A and γ ¯ h = | 1 + h | + | h | A b a .
Proof. 
For all x [ a , b ] , we have
| u ( x ) u ^ n ( x ) | = | m = n + 1 u m ( x ) | m = n + 1 | u m ( x ) | u 1 m = n + 1 γ h m 1 = γ h n 1 γ h u 1 γ h n 1 γ h u 1 1 .
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:
L ( v ) = v ( x ) , N ( v ) = v ( x ) F ( x ) f ( x ) g ( x ) K 1 ( x , t ) R 1 ( v ( t ) ) d t R 2 ( v ( x ) ) .
Then, we get the following relations for the derivatives of functions u m :
u 1 ( x ) = h u 0 ( x ) F ( x ) f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u 0 ( t ) ) d t R 2 ( u 0 ( x ) )
and for m 2 :
u m ( x ) = ( 1 + h ) u m 1 ( x ) h f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u m 1 ( t ) ) d t + R 2 ( u m 1 ( x ) ) .
Now, we determine the boundaries of functions u m :
| u 0 ( x ) | u 0 1 , | u 1 ( x ) | = h u 0 ( x ) F ( x ) f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u 0 ( t ) ) d t R 2 ( u 0 ( x ) ) | h | | u 0 ( x ) | + | F ( x ) | + f ( x ) g ( x ) | K 1 ( x , t ) | | R 1 ( u 0 ( t ) ) | d t + | R 2 ( u 0 ( x ) ) | | h | F + u 0 1 + K 1 R 1 ( b a ) + R 2 u 0 1 | h | F + 1 + A b a u 0 1 ,
| u 2 ( x ) | = ( 1 + h ) u 1 ( x ) h f ( x ) g ( x ) K 1 ( x , t ) R 1 ( u 1 ( t ) ) d t + R 2 ( u 1 ( x ) ) | 1 + h | + | h | K 1 R 1 ( b a ) + R 2 u 1 1 | 1 + h | + | h | A b a u 1 1 = γ ¯ h u 1 1 .
Hence, we get in general for m 1
| u m ( x ) | γ ¯ h m 1 u 1 1 .
By using the above estimations, we have for any x [ a , b ]
| u ( x ) u ^ n ( x ) | = | m = n + 1 u m ( x ) | m = n + 1 | u m ( x ) | γ ¯ h n 1 γ ¯ h u 1 1 ,
which proves the inequality of Equation (21).
The last inequality follows from the definition of the taken norm:
u u ^ n 1 = sup x [ a , b ] | u ( x ) u ^ n ( x ) | + sup x [ a , b ] | u ( x ) u ^ n ( x ) | ,
and from the inequalities already proven. □

4. Examples

Example 1.
We start with the solution of the following Volterra–Fredholm equation:
u ( x ) + 1 2 0 x ( t x ) 2 u ( t ) d t + 1 4 0 1 ( t x ) u ( t ) d t = 3 16 + 5 3 x + x 3 6 + x 5 60
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 ) = 1 4 0 1 ( t x ) u ( t ) d t , F ( x ) = 3 16 + 5 3 x + x 3 6 + x 5 60 , and f ( x ) = 0 , g ( x ) = 1 , where x [ 0 , 1 ] , and c 0 = 1 . The solution of Equation (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
K 1 = 1 2 , R 1 = 1 , R 2 = 1 4 .
Hence, we get
A = 3 4 ,
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
u 1 ( x ) = h 362880 ( 294840 + 598905 x 350721 x 2 + + 60480 x 3 12285 x 4 + 11088 x 5 1008 x 6 + 72 x 7 + 2 x 9 ) .
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.
Example 2.
Vibrations of the freely supported reinforced concrete beam of span l = 6 m and cross section ( b / h ) 0.30 m/ 0.60 m can be described by means of the linear differential equation. Evenly distributed mass m of the beam per length unit, by assuming that the density of beam material (reinforced concrete) is of value ϱ = 2500 kg/m 3 , equals
m = b h ϱ = 450 kg / m .
The discussed beam is considered as the set of segments of length d x , and the kinetic energy of such an elementary segment, that is, the element of mass m d x (see Figure 3a) [19,20], is expressed as follows:
m d x y ˙ 2 ( x ) 2 ,
where y ˙ ( 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:
y ( x ) = 16 5 y m x l 2 x 3 l 3 + x 4 l 4 ,
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:
y ˙ ( x ) = 16 5 y ˙ m x l 2 x 3 l 3 + x 4 l 4 ,
where y ˙ m = d y m d t and the kinetic energy of the whole beam is equal to
E k = 2 0 l / 2 m y ˙ 2 ( x ) 2 d x .
By substituting Equation (26) into Equation (27) and by integrating, we get
E k = 3968 7875 m l y ˙ m 2 2 = 0.504 m l y ˙ m 2 2 .
The kinetic energy described by Equation (28) corresponds to the kinetic energy of the system with one degree of freedom, and the substitute mass equals
m = 0.504 m l = 1633 kg
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:
k = 48 E I l 3 ,
where E = 27 GPa is the modulus of elasticity of the beam material and I denotes the moment of inertia of the cross section b h 3 12 = 5.4 · 10 3 m 4 . Thus, we have k = 32.4 MN/m.
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:
m y ¨ + k y = 0 ,
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
y ¨ ( t ) + ω 0 2 y ( t ) = 0 ,
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 , y ˙ ( 0 ) = v 0 .
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, obtaining
y ˙ ( t ) + ω 0 2 0 t y ( τ ) d τ = v 0 ,
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
u 1 ( t ) = h 100 , u 2 ( t ) = 1 20000 h 200 + h ( 200 + 1985281 t 2 ) ,
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:
m y ( t ) + b 1 y ( t ) + k 1 0 t y ( t ) d t = f ( y , t ) ,
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:
y ( 0 ) = y 0 ,
y ( 0 ) = v 0 ,
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:
y ( t ) + 0 t k 1 m 0 τ y ( s ) d s + c 1 m y ( τ ) d τ + b 1 m y ( t ) = f r ( t ) ,
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
u 1 ( t ) = h t 6 7200000000 + t 5 60000 + 499997 t 4 1200000 2999 t 3 30000 + 124 t 2 125 3 t 50 + 1 10 .
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.

5. 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:
u ( x ) k = 1 n f k ( x ) g k ( x ) K k ( x , t ) R k ( u ( t ) ) d t = F ( x ) ,
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:
( b a ) 2 k = 1 n K k R k < 1 .
The homotopy analysis method can be also applied for solving the equation of the following form:
j = 1 l a j ( x ) u ( j ) ( x ) k = 1 n f k ( x ) g k ( x ) K k ( x , t ) R k ( u ( t ) ) d t = F ( x )
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
H ( Φ , p ) ( 1 p ) L Φ ( x ; p ) u 0 ( x ) p h H ( x ) N Φ ( x ; p ) ,
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 .

Author Contributions

Conceptualization, D.S. and R.W.; formal analysis, D.S.; investigation, D.S., E.H., R.W., K.G. and T.T.; methodology, D.S., E.H. and R.W.; writing—original draft, D.S. and E.H.; writing—review and editing, R.W., K.G. and T.T.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Liao, S. Homotopy analysis method: A new analytic method for nonlinear problems. Appl. Math. Mech. 1998, 19, 957–962. [Google Scholar]
  2. Liao, S. Homotopy Analysis Method in Nonlinear Differential Equations; Springer: Berlin, Germany; Higher Education Press: Beijing, China, 2012. [Google Scholar]
  3. Van Gorder, R. Control of error in the homotopy analysis of semi-linear elliptic boundary value problems. Numer. Algor. 2012, 61, 613–629. [Google Scholar] [CrossRef]
  4. Hetmaniok, E.; Słota, D.; Wituła, R.; Zielonka, A. An analytical method for solving the two-phase inverse Stefan problem. Bull. Pol. Acad. Sci. Tech. Sci. 2015, 63, 583–590. [Google Scholar] [CrossRef] [Green Version]
  5. Zurigat, M.; Momani, S.; Odibat, Z.; Alawneh, A. The homotopy analysis method for handling systems of fractional differential equations. Appl. Math. Model. 2010, 34, 24–35. [Google Scholar] [CrossRef]
  6. Zhang, X.; Tang, B.; He, Y. Homotopy analysis method for higher-order fractional integro-differential equations. Comput. Math. Appl. 2011, 62, 3194–3203. [Google Scholar] [CrossRef] [Green Version]
  7. Abbasbandy, S.; Shivanian, E. A new analytical technique to solve Fredholm’s integral equations. Numer. Algor. 2011, 56, 27–43. [Google Scholar] [CrossRef]
  8. Hetmaniok, E.; Nowak, I.; Słota, D.; Wituła, R. Convergence and error estimation of homotopy analysis method for some type of nonlinear and linear integral equations. J. Numer. Math. 2015, 23, 331–344. [Google Scholar] [CrossRef]
  9. Van Gorder, R.; Vajravelu, K. On the selection of auxiliary functions, operators, and convergence control parameters in the application of the homotopy analysis method to nonlinear differential equations: A general approach. Commun. Nonlinear Sci. Numer. Simulat. 2009, 14, 4078–4089. [Google Scholar] [CrossRef]
  10. Odibat, Z. A study on the convergence of homotopy analysis method. Appl. Math. Comput. 2010, 217, 782–789. [Google Scholar] [CrossRef]
  11. Hetmaniok, E.; Słota, D.; Trawiński, T.; Wituła, R. Usage of the homotopy analysis method for solving the nonlinear and linear integral equations of the second kind. Numer. Algor. 2014, 67, 163–185. [Google Scholar] [CrossRef]
  12. Behzadi, S.; Abbasbandy, S.; Allahviranloo, T.; Yildirim, A. Application of homotopy analysis method for solving a class of nonlinear Volterra-Fredholm integro-differential equations. J. Appl. Anal. Comput. 2012, 2, 127–136. [Google Scholar]
  13. Mohamed, M.; Gepreel, K.; Alharthi, M.; Alotabi, R. Homotopy analysis transform method for integro-differential equations. Gen. Math. Notes 2016, 32, 32–48. [Google Scholar]
  14. Ghanbari, B. The convergence study of the homotopy analysis method for solving nonlinear Volterra-Fredholm integrodifferential equations. Sci. World J. 2014, 2014. [Google Scholar] [CrossRef] [PubMed]
  15. Alipouru, M.; Vali, M.; Borzabadi, A. A direct approach for approximate optimal control of integro-differential equations based on homotopy analysis and parametrization method. IMA J. Math. Control Inf. 2017, 34, 611–630. [Google Scholar] [CrossRef]
  16. Liao, S. Notes on the homotopy analysis method: Some definitions and theorems. Commun. Nonlinear Sci. Numer. Simulat. 2009, 14, 983–997. [Google Scholar] [CrossRef]
  17. Russo, M.; Van Gorder, R. Control of error in the homotopy analysis of nonlinear Klein-Gordon initial value problems. Appl. Math. Comput. 2013, 219, 6494–6509. [Google Scholar] [CrossRef]
  18. Turkyilmazoglu, M. Convergence of the homotopy analysis method. arXiv 2010, arXiv:1006.4460v1. [Google Scholar]
  19. Gromysz, K. Investigations of Reinforced-Concrete Composite Slabs Loaded Temporarily, Cyclically And Kinematically; Publishing House of the Silesian University of Technology: Gliwice, Poland, 2013. [Google Scholar]
  20. Gromysz, K. Nonlinear analytical model of composite concrete slab free and forced vibrations. Procedia Eng. 2017, 193, 281–288. [Google Scholar] [CrossRef]
Figure 1. Logarithm of the squared residual E 9 .
Figure 1. Logarithm of the squared residual E 9 .
Mathematics 07 00904 g001
Figure 2. Distribution of error of the exact solution approximation for n = 4 (a) and n = 7 (b).
Figure 2. Distribution of error of the exact solution approximation for n = 4 (a) and n = 7 (b).
Mathematics 07 00904 g002
Figure 3. Reduction of the beam with continuous mass distribution to the system with one degree of freedom: (a) system with continuous mass distribution m and (b) substitute mass m concentrated in a point where the vibration arrow of the system with continuous mass distribution occurs.
Figure 3. Reduction of the beam with continuous mass distribution to the system with one degree of freedom: (a) system with continuous mass distribution m and (b) substitute mass m concentrated in a point where the vibration arrow of the system with continuous mass distribution occurs.
Mathematics 07 00904 g003
Figure 4. Oscillator of parameters k and m as the substitute system for the reinforced concrete beam.
Figure 4. Oscillator of parameters k and m as the substitute system for the reinforced concrete beam.
Mathematics 07 00904 g004
Figure 5. Logarithm of the squared residual E 9 .
Figure 5. Logarithm of the squared residual E 9 .
Mathematics 07 00904 g005
Figure 6. Distribution of error of the exact solution approximation for n = 10 (a) and n = 15 (b).
Figure 6. Distribution of error of the exact solution approximation for n = 10 (a) and n = 15 (b).
Mathematics 07 00904 g006
Figure 7. Distribution of error of the exact solution approximation for n = 7 (a) and n = 12 (b).
Figure 7. Distribution of error of the exact solution approximation for n = 7 (a) and n = 12 (b).
Mathematics 07 00904 g007
Table 1. Values of errors of the exact solution reconstruction ( Δ n = u e u ^ n 1 ).
Table 1. Values of errors of the exact solution reconstruction ( Δ n = u e u ^ n 1 ).
n12345
Δ n 5.101 × 10 2 1.007 × 10 3 1.963 × 10 5 3.296 × 10 7 5.946 × 10 9
Δ ( 22 ) 106.077 1 79.932 1 60.231 1 45.386 1 34.200
n 678910
Δ n 8.882 × 10 11 1.519 × 10 12 2.309 × 10 14 1.332 × 10 15 4.441 × 10 16
Δ ( 22 ) 25.770 19.419 14.633 11.026 1 8.308
Table 2. Values of errors of the exact solution reconstruction ( Δ n ).
Table 2. Values of errors of the exact solution reconstruction ( Δ n ).
n12345
Δ n 1.428 8.117 40.616 43.243 18.710
n 678910
Δ n 3.508 0.155 2.878 × 10 2 1.233 × 10 3 2.843 × 10 4
n 1112131415
Δ n 9.155 × 10 6 1.788 × 10 6 1.812 × 10 7 6.266 × 10 9 9.220 × 10 10
n 1617181920
Δ n 1.131 × 10 10 5.933 × 10 12 2.335 × 10 13 9.648 × 10 14 6.550 × 10 14
Table 3. Values of errors of the exact solution reconstruction ( Δ n = y e u ^ n 1 ).
Table 3. Values of errors of the exact solution reconstruction ( Δ n = y e u ^ n 1 ).
n12345
Δ n 3.153 2.944 1.317 0.324 4.852 × 10 2
n 678910
Δ n 4.739 × 10 5 3.132 × 10 4 1.417 × 10 5 4.290 × 10 7 7.948 × 10 9
n 1112131415
Δ n 5.889 × 10 11 7.592 × 10 13 1.704 × 10 14 2.498 × 10 16 5.551 × 10 17

Share and Cite

MDPI and ACS Style

Słota, D.; Hetmaniok, E.; Wituła, R.; Gromysz, K.; Trawiński, T. Homotopy Approach for Integrodifferential Equations. Mathematics 2019, 7, 904. https://doi.org/10.3390/math7100904

AMA Style

Słota D, Hetmaniok E, Wituła R, Gromysz K, Trawiński T. Homotopy Approach for Integrodifferential Equations. Mathematics. 2019; 7(10):904. https://doi.org/10.3390/math7100904

Chicago/Turabian Style

Słota, Damian, Edyta Hetmaniok, Roman Wituła, Krzysztof Gromysz, and Tomasz Trawiński. 2019. "Homotopy Approach for Integrodifferential Equations" Mathematics 7, no. 10: 904. https://doi.org/10.3390/math7100904

APA Style

Słota, D., Hetmaniok, E., Wituła, R., Gromysz, K., & Trawiński, T. (2019). Homotopy Approach for Integrodifferential Equations. Mathematics, 7(10), 904. https://doi.org/10.3390/math7100904

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop