High-Order Iterative Scheme for a Viscoelastic Wave Equation and Numerical Results

In this paper, we consider a Robin problem for a viscoelastic wave equation. First, by the high-order iterative method coupled with the Galerkin method, the existence of a recurrent sequence via an N-order iterative scheme is established, and then the N-order convergent rate of the obtained sequence to the unique weak solution of the proposed model is also proved. Next, with N � 2, a numerical algorithm given by the finite-difference method is constructed to approximate the solution via the 2-order iterative scheme. Moreover, the same algorithm for the single-iterative scheme generated by the 2-order iterative scheme is also considered. Finally, comparison with errors of the numerical solutions obtained by the single-iterative scheme and the 2-order iterative scheme shows that the convergent rate of the 2-order iterative scheme is faster than that of the single-iterative scheme.


Introduction
In this paper, we consider the following initial boundary value problem for a viscoelastic wave equation with nonlinear damping: u(x, 0) � u 0 (x), u t (x, 0) � u 1 (x), where λ > 0, q ≥ 2, h 0 ≥ 0, and h 1 ≥ 0 are constants, with h 0 + h 1 > 0, and f, g, u 0 , and u 1 are given functions. Equation (1) arises naturally within frameworks of mathematical models in engineering and physical sciences. e left-hand integral of equation (1) stands for the characters of viscoelastic materials. Many researchers have paid attention to viscoelastic materials for a quite long time, especially in the last two decades, and have made a lot of progress, taking into account viscoelastic fluid, which achieved major attention due to its application in different physiological and industrial processes. In the same content, nanofluid has become an interesting objective which describes various phenomena such as electrical conductivity, especially in the bubble electrospinning [1][2][3], heat transfer on solid particle motion [4], biologically inspired peristaltic transport [5], and rheology controlled by the concentration of the added particles (such as SiC) [6]. In addition to studying the specific properties of viscoelastic materials, numerous researchers have considered the extensions of the mathematical model for viscoelastic problems and have obtained many interesting properties of solutions such as global existence, decay, and blow-up result. One of the problems similar to problems (1)-(3) was considered by Messaoudi [7], in which the blow-up result of solutions with negative initial energy was established. After that, Li and He [8] proved, under suitable conditions, the global existence and the general decay of solutions for the same model. Recently, the results obtained in [8] have also been investigated by Mezouar and Boulaaras [9,10] for the proposed nonlocal viscoelastic problems. Consider a recurrent sequence u m associated with equation (1) and defined by 0 < x < 1, 0 < t < T, with u m satisfying (2) and (3). If the sequence u m converges to the weak solution u of problems (1)- (3) and satisfies an estimate of N order in the form of ‖u m − u‖ X ≤ C‖u m− 1 − u‖ N X , for some C > 0, all integer numbers N ≥ 2, and X is a certain Banach space, such a method for finding the solution of problems (1)-(3) is called high-order iterative method. e original idea of this method is based on investigating the recurrent relations of a Newton-like method in Banach spaces [11]. After that, this approach was also applied successfully to [12][13][14][15][16][17]. In [17], Truong et al. considered the nonlinear wave equation of Kirchhoff-Carrier type as follows: where μ, f, u 0 , and u 1 are given functions, h ≥ 0 is a given constant, and μ(t, ‖u(t)‖ 2 , ‖u x (t)‖ 2 ) depends on the integrals ‖u(t)‖ 2 � 1 0 u 2 (x, t)dx and ‖u x (t)‖ 2 � 1 0 u 2 x (x, t)dx. e authors associated the first equation in (5) with a recurrent sequence u m defined by where u m satisfies second and third equations in (5), and proved that u m converges to the unique weak solution of problem (5). Moreover, with N � 3, the 3-order iterative scheme was established, and some numerical results of finite-difference approximate solutions were presented. In [15], the authors investigated the Dirichlet problem for a wave equation with linear damping and nonlinear integral as follows: g(x, t, s, u(x, s))ds, 0 < x < 1, 0 < t < T, where μ, f, g, u 0 , and u 1 are given functions and λ ≠ 0 is a given constant. If μ ∈ C 1 ([0, 1] × R + ), f ∈ C N ([0, 1] × R + × R), and g ∈ C N− 1 ([0, 1] × Δ × R), with Δ � (t, s) { ∈ R 2 + : s ≤ t}, by the high-order iterative method coupled with the Galerkin method, the existence of a recurrent sequence u m associated with equation (7) and defined by 0 < x < 1, 0 < t < T, with u m satisfying second and third equations in (7), was proved, and u m convergence with the N-order rate to the unique weak solution of problem (7) was also confirmed. Some other iteration methods can be widely used to solve various nonlinear problems, for example, the variational iteration method and the homotopy perturbation method. Based on the basic idea of the enhanced perturbation method and the homotopy technology, Li and He [18] constructed a modified homotopy perturbation method, making a very high accuracy of the obtained approximate solution. Recently, Ji et al. [19] successfully applied Li and He's modified homotopy perturbation method coupled with the energy method for the dropping shock response of a tangent nonlinear packaging system. When g � 0 and f ≡ b|u| p− 2 u, equation (1) is reduced to the following nonlinear wave equation: where a, b > 0 and p, q ≥ 2, which has been extensively studied and obtained many interesting results such as global existence, exponential decay, and finite-time blow-up result.
If a � 0, it is well known that the source term b|u| p− 2 u causes a finite-time blow-up phenomenon of the solution with negative initial energy; see [20][21][22][23]. If b � 0, the damped term a|u t | q− 2 u t assures the global existence for suitably initial data; see [24][25][26][27][28][29]. e interactions between the damping and the source terms give the results that the solutions with any initial data continue to exist globally in time if q ≥ p and blow up in finite time if q < p and the initial energy is sufficiently negative.
In the presence of the viscoelastic term (g ≠ 0) and a � b � 1 in (9), Messaoudi [7] considered a Dirichlet problem for a nonlinear wave equation as follows: where Ω is a bounded domain of R n (n ≥ 1) with a smooth boundary zΩ, p > 2, m ≥ 1, and g: R + ⟶ R + is a positive nonincreasing function. With suitable conditions on g, he proved that the solutions with initial negative energy blow up in finite time if p > m and continue to exist if p ≤ m satisfied the condition Later, in case of m � 2 and x ∈ R n , Kafini and Messaoudi [30] also established the finite-time blow-up result with suitable conditions on the initial data and the relaxation function g. In the presence of the strong damping − Δu t and the linear damping u t (m � 2), Li and He [8] proved the global existence of solutions and established a general decay rate estimate for problem (10). Moreover, they showed the finite-time blow-up results of solutions with both negative initial energy and positive initial energy.
Although there are many studies of solution properties of viscoelastic problems, however, it seems that few works related to numerical algorithms for this type were published. In [31], Long et al. proved the global existence and exponential decay of equation (1) associated with a mixed nonhomogeneous condition and initial condition (3). With q � 4 and f � u 2 + F(x, t), the derivatives were first approximated by finite-difference schemes. en, a linear recursive scheme generated by the nonlinear difference equation was constructed. Finally, the exact solution and the approximate solution were illustrated numerically. In [14], Ngoc et al. also obtained the same results given in [31] for the following nonlinear wave equation associated with nonlocal boundary conditions:

Mathematical Problems in Engineering
where a � 0, λ > 0, and p > 2 are given constants and u 0 , u 1 , f, g i , k i , and H i (i � 0, 1) are given functions satisfying conditions specified later.
In some recent literature studies, various difference methods have been applied to studying the consistency, stability, efficiency, and convergence of the proposed schemes such as Boulaaras [32] used finite element methods to prove the existence and uniqueness of the discrete solution for an evolutionary implicit 2-sided obstacle problem. Boulaaras and Haiour [33] analysed the convergence and regularity of the proposed algorithm via the finite element methods coupled with a theta time discretization scheme for evolutionary Hamilton-Jacobi-Bellman equations. Mohanty and Gopal [34] used a cubic spline finite difference method of Numerov type with an accuracy of order two in time and four in space directions for the solution of the telegraphic equation with the forcing function: (14) where α > 0 and β ≥ 0 are real parameters. In [35], Alsuyuti et al. applied a new modified Galerkin algorithm based on the shifted Jacobi polynomials to solve fractional differential equations and a system of fractional differential equations governed by homogeneous and nonhomogeneous initial and boundary conditions. e new algorithm was also applied for fractional partial differential equations with Robin boundary conditions and the time-fractional telegraph equation. In addition, some illustrative examples were presented to confirm the validity and efficiency of the proposed algorithm. Recently, capability of the meshless methods has been proved by using them for many different problems; see [36,37] and the references therein. In [36], Oruç used two meshless methods based on the local radial basis function and barycentric rational interpolation for solving a two-dimensional viscoelastic wave equation. e stability of the methods, the accuracy, and the computational cost were considered. From the comparisons, it was shown that the performance of the barycentric rational interpolation in the sense of accuracy is slightly better than the performance of the local radial basis function; however, the computational cost of the local radial basis function is less than the computational cost of the barycentric rational interpolation. Before, a local meshless method was proposed in [37] for convection-dominated steady and unsteady partial differential equations in which numerical results have confirmed that the new approach is accurate and efficient for solving a wide class of one-and two-dimensional convection-dominated problems having sharp corners and jump discontinuities.
e first goal of our present paper is devoted to studying the existence and the N-order convergence of the high-order iterative scheme defined by (4). In case N � 2, λ � 1, and q � 2, the second goal is mentioned to building a numerical algorithm in order to approximate the successive solutions of the 2-order iterative scheme as follows: Furthermore, in a specific case of (15) with b m (x, t) � 0, the corresponding scheme called the single-iterative scheme is considered. en, the numerical algorithms of both schemes are established, and the results of errors are presented to compare their convergent rates also.
For the first purpose, by using the high-order iterative method coupled with the Galerkin method, we shall prove the existence of a recurrent sequence u m associated with equation (1) and defined by (4), and then u m convergence with the N-order rate to the unique weak solution of problems (1)-(3) will also be claimed. For the second purpose, first, we shall use the uniform spatial partition x i � ih, h � 1/N, i � 0, 1, . . ., N, and the forward difference formulas (see [38], pages 36 and 43) to approximate the k th derivatives. en, problems (15) and (16) will be changed into a system of second-order integrodifferential equations (in the time variable) of the unknown functions . ., N. Normally, this system will be converted into a system of 2N first-order integrodifferential equations by using the auxiliary functions ; see the same transformations in [14,17,23,26,29,31]; however, these transformations will increase computations. To reduce the computations, the above second-order integrodifferential system will be transformed into a system of N first-order integrodifferential equations by integrating in time on interval (0, t). After that, to approximate double integrals, the trapezoidal formula will be successively used twice. It can be said with much confidence that this technique has never been used before. Next, by using uniform partition t j � jΔt, Δt � T/M, j � 0, 1, . . ., M, for discretization in time variable t, we will obtain an algorithm to determine the finite-difference approximate solutions of u (m) via 2-order iterative schemes (15) and (16) given by the following difference equation (formula (149) below): where Ψ [ * ,m] j is defined by (147) and (148). Finally, we will construct the algorithm to find the finite-difference approximate solutions of u (m) given by the single-iterative scheme (formulas (157)-(159) below) and present a numerical example to compare convergent rates of two schemes. e errors of computations of the numerical solutions given by two schemes show that the convergent rate of the 2-order iterative scheme is faster than that of the single-iterative scheme.
Our paper is organized as follows. In Section 2, we introduce some notations and modified lemmas. In Section 3, by using the Galerkin method and standard arguments of compactness, we prove the existence and convergence of the high-order sequence defined by (4). In Section 4, by using the finite-difference method and some new techniques to reduce computations and to approximate a double integral, we construct a numerical algorithm to determine the finitedifference approximate solutions of u (m) via 2-order iterative schemes (15) and (16). Moreover, a concrete example is numerically illustrated to compare the convergent rate of the single-iterative scheme with that of the 2-order iterative scheme. In Section 5, we summarize the main outcomes of our paper.

Preliminaries
First, we put Ω � (0, 1) and denote the usual function spaces used in this paper by the notations L p � L p (Ω) and H m � H m (Ω). Let 〈·, ·〉 be either the scalar product in L 2 or the dual pairing of a continuous linear functional and an element of a function space. e notation ‖·‖ stands for the norm in L 2 , ‖·‖ X is the norm in the Banach space X, and X ′ is the dual space of X.
We denote by L p (0, T; X), 1 ≤ p ≤ ∞, for the Banach space of real functions u: (0, T) ⟶ X measurable such that On H 1 , three norms ‖v‖ H 1 , ‖v‖ a , and ‖v‖ i are equivalent norms. e weak solution of problems (1)-(3) can be defined as and u satisfies the following variational equation: for all w ∈ H 1 , and a.e., t ∈ (0, T), together with the initial conditions where a(·, ·) is the symmetric bilinear form on H 1 defined by (19).
We now have the following lemmas, the proofs of which are straightforward, so we omit the details.

Mathematical Problems in Engineering
ere exists the Hilbert orthonormal base w j of L 2 consisting of the eigenfunctions w j corresponding to the eigenvalue λ j such that Furthermore, the sequence w j / �� λ j is also a Hilbert orthonormal base of H 1 with respect to the scalar product a(·, ·).
On the contrary, we have w j satisfying the following boundary value problem: e proof of Lemma 3 can be found in eorem 7.7 of [39], p.87, with H � L 2 and a(·, ·) as defined by (19).

The High-Order Iterative Method
First, we make the following assumptions: Fix T * > 0. For each M > 0 given, we set the constant K M (f) as follows: For every T ∈ (0, T * ] and M > 0, we put Now, we establish a recurrent sequence u m in which the first term is chosen as u 0 ≡ 0, and suppose that e m th iteration u m can be defined by the following problem: find u m ∈ W 1 (M, T)(m ≥ 1) satisfying the nonlinear variational problem e existence of the above sequence u m is claimed by the following theorem.
Proof. e proof of eorem 1 consists of three steps. □ Step 1 (Faedo-Galerkin approximation). Let w j be a basis of H 1 as in Lemma 3; we find an approximate solution of problems (32) and (33) in the form of where the coefficients c (k) mj satisfy the following system of nonlinear differential equations: with Note that, by using (31) and standard methods in ordinary differential equations, we can prove that system (35) admits a unique solution c (k) . e following estimates allow one to take T (k) m � T independent of m and k.
Step 2 (a priori estimates). First, for all j � 1, . . . , k, multiplying the first equation in (35) by _ c (k) mj (t), summing on j, and integrating with respect to the time variable from 0 to t, we have where Next, by replacing w j in the first equation in (35) similar to the first equation in (35), and it yields with (43)

Mathematical Problems in Engineering
Putting it follows from (39), (42), and (44) that We estimate S (k) m (0) and the integrals on the right-hand side of (45) as follows.
First integral I 1 : it is clear that Second integral I 2 : by the inequality 2ab ≤ 1/2a 2 + 2b 2 , ∀a, b ∈ R, we have ird integral I 3 : using the following inequality Fourth integral I 4 : Fifth integral I 5 : note that the first equation in (35) can be written as follows: 8

Mathematical Problems in Engineering
Hence, replacing w j with € u (k) m (t) and using the inequality Integrating in t, we have Sixth integral I 6 : Seventh integral I 7 : Combining (45), (46)-(49), and (53)-(55), it leads to where To estimate integrals e following inequalities are valid: where A M and B M are defined as follows: Proof of Lemma 4 □ proof of (i). Using the inequalities (a + b) i ≤ 2 i− 1 (a i + b i ), for all a, b ≥ 0, i ≥ 1, and s i ≤ 1 + s q , ∀s ≥ 0, ∀i, q, 0 ≤ i ≤ q, we have with a i , 0 ≤ i ≤ N − 1, defined by (59). Hence, and it implies that equation (i) in (58) holds.
□ proof of (ii). We also have Hence, It follows that with b i , 0 ≤ i ≤ N − 1, and B M defined by (59). Hence,

Mathematical Problems in Engineering 11
Combining (60) and (66), we have where By means of convergences (36), we can deduce the existence of a constant M > 0 independent of k and m such that Finally, it follows from (67) and (69) that (70) en, by solving nonlinear Volterra integral inequality (70) (based on the methods in [40]), the following lemma is proved.

Lemma 5.
ere exists a constant T > 0 independent of k and m such that By Lemma 5, we can take constant T (k) m � T for all k and m ∈ N. us, we have Step 3 (limiting process). anks to (72), there exists a sub- By using the following inequality, it follows from (44), (71), and (74) 2 that On the contrary, by L ∞ (0, T; H 1 )↪L ∞ (Q T ) and using the inequality we deduce from (71) that erefore, (74) 1 and (78) lead to We note that By (33), (37), and (71), it gives Hence, we have so Passing to the limit in (35) and (36), we have u m satisfying (32) and (33) in L 2 (0, T).
On the contrary, it follows from the first equation in (32) and the fourth equation in (73) that Hence, u m ∈ W 1 (M, T), and eorem 1 is proved. Next, the main result is also given by the following theorem. We consider the space W 1 (T), defined by and then W 1 (T) is a Banach space with respect to the norm  (32) and (33) converges at the N-order rate to the solution u strongly in W 1 (T) in the sense for all m ≥ 1, where C is a suitable constant. Furthermore, the following estimate is fulfilled: where C T and 0 < c T < 1 are the constants only depending on T.
Proof. (i) Existence of solutions: we shall prove that u m is a Cauchy sequence in W 1 (T). Indeed, we put v m � u m+1 − u m . en, v m satisfies the variational problem Taking w � v m ′ in (89), after integrating in t, we get with Next, we need to estimate the integrals on the right side of (90).
By the inequality we have It is not difficult to estimate J 2 , J 3 , and J 4 as follows: Mathematical Problems in Engineering 13 Using Taylor's expansion of the function f( Hence, it follows from (33) and (97) that and then (2) T � ( (100) en, it implies from (90), (93)-(96), and (100) that where for all m and p ∈ N.

(103)
Hence, u m is a Cauchy sequence in W 1 (T). en, there exists u ∈ W 1 (T) such that u m ⟶ u strongly in W 1 (T). (104) Note that u m ∈ W 1 (M, T); then, there exists a subsequence u m j of u m such that Moreover, by (104) and the inequality we have On the contrary, (108) erefore, it implies from (104) and (38) that Finally, passing to the limit in (32) and (33) as m � m j ⟶ ∞, there exists u ∈ W(M, T) satisfying the equation for all w ∈ H 1 , and the initial condition On the contrary, it follows from the fourth equation in (105) and (110) that Hence, u ∈ W 1 (M, T). □ Proof. (ii) Uniqueness: applying the similar estimations used in the proof of eorem 1, it is easy to prove that u ∈ W 1 (M, T) is a unique local weak solution of problems (11)-(13).
Passing to the limit in (103) as p ⟶ ∞ for fixed m, we get (88). Also, with a similar argument, (87) follows. eorem 2 is proved completely. □ Remark 1. In order to construct the N− order iterative scheme defined by (14), (H 3 ) is a necessary assumption. If we only consider the existence of solutions of problems (11)- (13), this condition can be weakened as follows: see [42].

Numerical Results
In this section, we shall construct a numerical algorithm for the 2-order iterative scheme obtained by (14) and (15). e finite-difference method and some techniques for approximating double integrals will be used to find the finitedifference approximate solutions of u (m) (x, t) of this scheme. Moreover, a numerical algorithm for the singleiterative scheme obtained by (154) and (155) below will also be constructed, and a concrete example will be numerically illustrated to compare the convergent rates of the singleiterative scheme and the 2-order iterative scheme. Consider problems (11)-(13) with λ � h 0 � h 1 � 1 and q � 2; then, they are transformed into the following problem: where the functions g, f, u 0 , and u 1 are defined by e exact solution of problem (114), with g, f, u 0 , and u 1 defined by (115), is a function u ex given by With datum (115) and 0 ≤ x ≤ 1 and 0 ≤ t ≤ 0.5, Figure 1 describes the surface of the exact solution u ex (x, t) of problem (114) below.
We consider the 2-order iterative scheme of problem (114) as follows: In order to solve problems (117) and (118) numerically, we first use the spatial difference grid u (m) Replacing the derivatives in spatial variable x of (117) by the following approximations, problems (117) and (118) are turned into the following system of intergodifferential equations:  16 Mathematical Problems in Engineering Boundary conditions in the second equation in (120) lead to Eliminating the unknown functions u (m) 0 and u (m) N+1 in the first equation (i � 1) and the last equation (i � N), respectively, system (120) is equivalent to where a 1 � 1 + 2h/h 2 (1 + h). Rewrite (123) into a vector equation: Mathematical Problems in Engineering and E, B (m) (t), and C ∈ M N (M N is the set of real N-order matrices) are defined by Integrating in time variable t, we have (128) Approximating the derivatives d u →(m) /dt(t j ) by the following time differences, (130) Note that the integral t j 0 Φ(s)ds can be approximated by the trapezoidal formula: (131) Applying (131) Similarly, we also obtain the following approximations: So, equation (130) can be rewritten by where Note that Ψ [1,m] j and Ψ [2,m] where

Description of Finite-Difference Algorithm (135) and (136)
(i) Let M and N be fixed constants. At the first iteration with m � 0, we set up the given vector (ii) At the (m − 1) th iteration, suppose that we can compute as follows: Hence, Hence, When the process of the computations is reached to C4. e error of two successively approximate iterations.
e process of iteration will be stopped at the m th step when the following estimate is satisfied: C5. e error between the approximate solution (at the m th step) and the exact solution is defined by in which u ex (x, t) � 2(1 + x − x 2 )e − t is the exact solution of (114) satisfying datum (115).
With datum (115) and the grid of N � 50 and M � 100, Figure 2 describes the surface of the finite-difference approximate solution of u (m) (x, t) defined by 2-order iterative scheme (117) and (118), with respect to algorithm (135) and (136).
With T � 0.5, the computational results presented in Tables 1-5 get along with algorithms (135) and (136) and (157)-(159) of iterative schemes (117) and (118) and (154) and (155), respectively. Table 1 (Table 3) are also decreased when the number of iterations is increased. Comparing columns 2 and 3 of Table 3, the errors of the 2-order iterative scheme are still smaller than those of the single-iterative scheme on the same number of iterations.
Note that, with the grid of N � 10 and M � 20, the above remarks on Tables 2 and 3 are also valid on Tables 4 and 5.
Finally, with the datum as in (115) and the grid of N � 50 and M � 100, we plot the curved surfaces of the finitedifference approximate solutions defined by algorithms (135) and (136) and (157)-(159), respectively, and the exact solution u ex of problem (115).

Conclusion
In this article, an initial boundary value problem for a viscoelastic wave equation with nonlinear damping is investigated, and its main outcomes are summarized in two parts. In part 1, theoretically, the existence of a recurrent sequence via a high-order iterative scheme is proved, and the high-order convergence of this sequence to the unique weak solution of the proposed model is also claimed. In part 2, two specific cases of the high-order iterative scheme called the 2order iterative scheme and the single-iterative scheme are considered, in which two numerical algorithms for finding the approximate solutions corresponding to these schemes are constructed by the finite-difference method and the techniques approximating double integrals. To close this part, a concrete example is numerically considered. And the studied results of the errors show that the convergent rate of the 2-order iterative scheme is faster than that of the singleorder iterative scheme.

Data Availability
e data used to support the findings of this study are included within the article.