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

Abstract

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In this paper, we consider the following initial boundary value problem for a viscoelastic wave equation with nonlinear damping:$$\begin{array}{}\text{(1)}& u_{tt}-u_{xx}+\lambda {u_t}^{q-2}{u}_t+{\displaystyle {\int }_0^{t}gt-s{u}_{xx}x,s\text{d}s}=fx,t,u,\hspace{1em}0<x<1,0<t<T,\text{(2)}& u_{x}0,t-h_{0}u0,t=u_{x}1,t+h_{1}u1,t=0,\text{(3)}& ux,0={\tilde{u}}_{0}x,u_{t}x,0={\tilde{u}}_{1}x,\end{array}$$where $\lambda >0,q\ge 2,h_0\ge 0,\text{\hspace{0.17em}and\hspace{0.17em}}{h}_1\ge 0$ are constants, with $h_0+h_1>0$, and $f,g$, ${\tilde{u}}_0$, and ${\tilde{u}}_1$ are given functions.

Equation (1) arises naturally within frameworks of mathematical models in engineering and physical sciences. The 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–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$$\begin{array}{}\text{(4)}& \begin{array}{c}{u}_0\equiv 0,\hfill \\ u_m^{″}-\Delta u_m+\lambda {u_m^{\prime }}^{q-2}{u}_m^{\prime }+{\displaystyle {\int }_0^{t}gt-s\Delta u_{m}x,s\text{d}s}\hfill \\ ={\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}{D}_3^{i}fx,t,u_{m-1}{u_m-u_{m-1}}^i},\hfill \end{array}\end{array}$$$0<x<\mathrm{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\le C{u_{m-1}-u}_X^N$, for some $C>0$, all integer numbers $N\ge 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. The 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–17]. In [17], Truong et al. considered the nonlinear wave equation of Kirchhoff-Carrier type as follows:$$\begin{array}{}\text{(5)}& \begin{array}{c}{u}_{tt}-\mu t,{ut}^2,{u_{x}t}^2{u}_{xx}=fu,\hspace{1em}0<x<\mathrm{1,0}<t<T,\hfill \\ u_{x}0,t-hu0,t=gt,u1,t=0,\hfill \\ ux,0={\tilde{u}}_{0}x,u_{t}x,0={\tilde{u}}_{1}x,\hfill \end{array}\end{array}$$where $\mu ,f,{\tilde{u}}_0$, and ${\tilde{u}}_1$ are given functions, $h\ge 0$ is a given constant, and $\mu t,{ut}^2,{u_{x}t}^2$ depends on the integrals ${ut}^2={\displaystyle {\int }_0^1{u}^{2}x,t\text{d}x}$ and ${u_{x}t}^2={\displaystyle {\int }_0^1{u}_x^{2}x,t\text{d}x}$. The authors associated the first equation in (5) with a recurrent sequence $u_m$ defined by$$\begin{array}{}\text{(6)}& \frac{{\partial }^2{u}_m}{\partial t^2}-\mu t,{u_{m}t}^2,{u_{mx}t}^2{u}_{mxx}={\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}\frac{{\partial }^{i}f}{\partial u^i}{u}_{m-1}{u_m-u_{m-1}}^i},\hspace{1em}0<x<1,0<t<T,\end{array}$$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:$$\begin{array}{}\text{(7)}& \begin{array}{c}{u}_{tt}-\frac{\partial }{\partial x}\mu x,t{u}_x+\lambda u_t=fx,t,u+{\displaystyle {\int }_0^{t}gx,t,s,ux,s\text{d}s},\hspace{1em}0<x<\mathrm{1,0}<t<T,\hfill \\ u0,t=u1,t=0,\hfill \\ ux,0={\tilde{u}}_{0}x,u_{t}x,0={\tilde{u}}_{1}x,\hfill \end{array}\end{array}$$where $\mu$, $f$, $g$, ${\tilde{u}}_0$, and ${\tilde{u}}_1$ are given functions and $\lambda \ne 0$ is a given constant. If $\mu \in C^1\mathrm{0,1}\times {ℝ}_{+}$, $f\in C^N\mathrm{0,1}\times {ℝ}_{+}\times ℝ$, and $g\in C^{N-1}\mathrm{0,1}\times \mathrm{\Delta }\times ℝ$, with $\mathrm{\Delta }=t,s\in {ℝ}_{+}^2:s\le 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$$\begin{array}{}\text{(8)}& \begin{array}{c}{u}_0\equiv 0,\hfill \\ \begin{array}{c}\frac{{\partial }^2{u}_m}{\partial t^2}-\frac{\partial }{\partial x}\mu x,t\frac{\partial u_m}{\partial x}+\lambda \frac{\partial u_m}{\partial t}={\displaystyle \sum _{i=0}^{N-1}\frac{1}{k!}\frac{{\partial }^{k}f}{\partial u^k}x,t,u_{m-1}{u_m-u_{m-1}}^k}\\ +{\displaystyle \sum _{k=0}^{N-1}\frac{1}{k!}{\displaystyle {\int }_0^t\frac{{\partial }^{k}g}{\partial u^k}x,t,s,u_{m-1}x,s}}{u_{m}x,s-u_{m-1}x,s}^k\text{d}s,\end{array}\hfill \end{array}\end{array}$$$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\text{\hspace{0.17em}and\hspace{0.17em}}f\equiv b{u}^{p-2}u$, equation (1) is reduced to the following nonlinear wave equation:$$\begin{array}{}\text{(9)}& u_{tt}-\Delta u+a{u_t}^{q-2}{u}_t=b{u}^{p-2}u,\end{array}$$where $a,b>0$ and $p,q\ge 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–23]. If $b=0$, the damped term $a{u_t}^{q-2}{u}_t$ assures the global existence for suitably initial data; see [24–29]. The 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\ge 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\ne 0$ and $a=b=1$ in (9), Messaoudi [7] considered a Dirichlet problem for a nonlinear wave equation as follows:$$\begin{array}{}\text{(10)}& \begin{array}{c}{u}_{tt}-\Delta u+{{\displaystyle \int }}_0^{t}gt-s\Delta ux,s\text{d}s+{u_t}^{m-2}{u}_t=u^{p-2}u,\text{in\hspace{0.17em}}\Omega \times 0,\infty ,\hfill \\ ux,t=0,x\in \partial \Omega ,t\ge 0,\hfill \\ ux,0=u_{0}x,u_{t}x,0=u_{1}x,x\in \Omega ,\hfill \end{array}\end{array}$$where $\mathrm{\Omega }$ is a bounded domain of ${ℝ}^n$ ($n\ge 1$) with a smooth boundary $\partial \mathrm{\Omega }$, $p>2$, $m\ge 1,$ and $g:{ℝ}^{+}⟶{ℝ}^{+}$ 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\le m$ satisfied the condition$$\begin{array}{}\text{(11)}& \max m,p\le \frac{2n-1}{n-2}\text{,}\hspace{1em}\text{with\hspace{0.17em}}n\ge 3.\end{array}$$

Later, in case of $m=2$ and $x\in {ℝ}^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 $-\mathrm{\Delta }{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$$\begin{array}{}\text{(12)}& u_{x}0,t=u0,t,u_{x}1,t+\eta u1,t=ht,\end{array}$$and initial condition (3). With $q=4$ and $f=u^2+Fx,t$, the derivatives were first approximated by finite-difference schemes. Then, 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:$$\begin{array}{}\text{(13)}& \begin{array}{c}{u}_{tt}-u_{xx}+u+\lambda u_t=a{u}^{p-2}u+fx,t,\hspace{1em}0<x<1,t>0,\hfill \\ u_{x}0,t=g_{0}t-{\displaystyle {\int }_0^t{H}_{0}t-su0,s\text{d}s}+{\displaystyle {\int }_0^t{k}_{0}x,tux,t\text{d}x},\hfill \\ -u_{x}1,t=g_{1}t-{\displaystyle {\int }_0^t{H}_{1}t-su1,s\text{d}s}+{\displaystyle {\int }_0^1{k}_{1}x,tux,t\text{d}x},\hfill \\ ux,0={\tilde{u}}_{0}x,u_{t}x,0={\tilde{u}}_{1}x,\hspace{1em}0<x<1,\hfill \end{array}\end{array}$$where $a=0$, $\lambda >0$, and $p>2$ are given constants and ${\tilde{u}}_0$, ${\tilde{u}}_1$, $f$, $g_i$, $k_i$, and $H_{i}i=\mathrm{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:$$\begin{array}{}\text{(14)}& u_{tt}+2\alpha u_t+{\beta }^{2}u=u_{xx}+fx,t,\hspace{1em}0<x<1,t>0,\end{array}$$where $\alpha >0\text{\hspace{0.17em}and\hspace{0.17em}}\beta \ge 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. The 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. The 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.

The 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$, $\lambda =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:$$\begin{array}{}\text{(15)}& \begin{array}{c}{u}_0\equiv 0,\hfill \\ u_m^{″}t+u_m^{\prime }t-\Delta u_{m}t+{{\displaystyle \int }}_0^{t}gt-s\Delta u_{m}sds-b_{m}x,t{u}_{m}x,t=F_{m}x,t,0<x<\mathrm{1,0}<t<T,\hfill \\ u_{mx}0,t-u_{m}0,t=u_{mx}1,t+u_{m}1,t=0,\hfill \\ u_{m}x,0={\tilde{u}}_{0}x,u_m^{\prime }x,0={\tilde{u}}_{1}x,\hfill \end{array}\end{array}$$where$$\begin{array}{c}\text{(16)}& F_{m}x,t=fx,t,u_{m-1}x,t-b_{m}x,t{u}_{m-1}x,t,b_{m}x,t=D_{3}fx,t,u_{m-1}x,t.\end{array}$$

Furthermore, in a specific case of (15) with $b_{m}x,t=0$, the corresponding scheme called the single-iterative scheme is considered. Then, 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=\mathrm{0,1},\dots$, $N$, and the forward difference formulas (see [38], pages 36 and 43) to approximate the $k^{\text{th}}$ derivatives. Then, problems (15) and (16) will be changed into a system of second-order integrodifferential equations (in the time variable) of the unknown functions $u_i^{m}t=u_m{x}_i,t$, $i=\mathrm{0,1},\dots$, $N$. Normally, this system will be converted into a system of $2N$ first-order integrodifferential equations by using the auxiliary functions ${\dot{v}}_i^{m}t=u_i^{m}t$; 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\mathrm{\Delta }t$, $\mathrm{\Delta }t=T/M$, $j=\mathrm{0,1},\dots$, $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):$$\begin{array}{}\text{(17)}& {\stackrel{⟶}{u}}_{j+1}^m={\stackrel{⟶}{u}}_j^m+\Delta t{\mathrm{\Psi }}_j^{\ast ,m},\end{array}$$where ${\mathrm{\Psi }}_j^{\ast ,m}$ 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. The 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 finite-difference 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.

2. Preliminaries

First, we put $\mathrm{\Omega }=\mathrm{0,1}$ and denote the usual function spaces used in this paper by the notations $L^p=L^p\mathrm{\Omega }$ and $H^m=H^m\mathrm{\Omega }$. Let $\cdot ,\cdot$ be either the scalar product in $L^2$ or the dual pairing of a continuous linear functional and an element of a function space. The notation $\cdot$ stands for the norm in $L^2$, ${\cdot }_X$ is the norm in the Banach space $X$, and $X^{\prime }$ is the dual space of $X$.

We denote by $L^{p}0,T;X$, $1\le p\le \infty$, for the Banach space of real functions $u:0,T⟶X$ measurable such that$$\begin{array}{c}\text{(18)}& u_{L^{p}0,T;X}={{{\displaystyle \int }}_0^T{ut}_X^p\text{d}t}^{1/p}<\infty \hspace{1em}\text{for\hspace{0.17em}}1\le p<\infty ,u_{L^{\infty }0,T;X}=\underset{0<t<T}{ess\sup }{ut}_X\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\text{for\hspace{0.17em}}p=\infty .\end{array}$$

With $f\in C^k\mathrm{0,1}\times {ℝ}_{+}\times ℝ$, $f=fx,t,u$, we put $D_{1}f=\partial f/\partial x$, $D_{2}f=\partial f/\partial t$, $D_{3}f=\partial f/\partial u$, and $D^{\alpha }f=D_1^{{\alpha }_1}{D}_2^{{\alpha }_2}{D}_3^{{\alpha }_3}f$; $\alpha ={\alpha }_1,{\alpha }_2,{\alpha }_3\in {\mathbb{Z}}_{+}^3$, $\alpha ={\alpha }_1+{\alpha }_2+{\alpha }_3\le k$, and $D^{\mathrm{0,0,0}}f=D^{0}f=f$.

Next, we put$$\begin{array}{c}\text{(19)}& au,v={\displaystyle {\int }_0^1{u}_{x}x{v}_{x}x\text{d}x+h_{0}u0v0+h_{1}u1v1}\text{,}\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}u,v\in H^1,v_a=\sqrt{av,v},\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}v\in H^1,\\ \text{(20)}& v_i={v^{2}i+{\displaystyle {\int }_0^1{v}_x^2}x\text{d}x}^{1/2},\hspace{1em}i=\mathrm{0,1.}\end{array}$$

On $H^1$, three norms $v_{H^1},v_a$, and $v_i$ are equivalent norms.

The weak solution of problems (1)–(3) can be defined as follows. A function $u=ux,t$ is a weak solution of problems (1)–(3) if$$\begin{array}{}\text{(21)}& u\in L^{\infty }0,T;H^2,u^{\prime }\in L^{\infty }0,T;H^1,u^{″}\in L^{\infty }0,T;L^2,\end{array}$$and $u$ satisfies the following variational equation:$$\begin{array}{}\text{(22)}& u^{″}t,w+aut,w+\lambda {u^{\prime }t}^{q-2}{u}^{\prime }t,w={\displaystyle {\int }_0^{t}gt-saus,w\text{d}s}+fx,t,u,w,\end{array}$$for all $w\in H^1$, and a.e., $t\in 0,T$, together with the initial conditions$$\begin{array}{}\text{(23)}& u0={\tilde{u}}_0,u^{\prime }0={\tilde{u}}_1,\end{array}$$where $a\cdot ,\cdot$ 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.

Lemma 1.

The imbedding $H^1\hookrightarrow C^0\overline{\mathrm{\Omega }}$ is compact, and$$\begin{array}{c}\text{(24)}& \text{i}{v}_{C^0\overline{\Omega }}\le \sqrt{2}{v}_{H^1},\text{ii}{v}_{C^0\overline{\Omega }}\le \sqrt{2}{v}_i,\\ \text{iii}\frac{1}{\sqrt{3}}{v}_{H^1}\le v_i\le \sqrt{3}{v}_{H^1},\hfill \end{array}$$for all $v\in H^1$, $i=\mathrm{0,1}$.

Lemma 2.

Let $h_0\ge 0$ and $h_1\ge 0$ with $h_0+h_1>0$. Then, the symmetric bilinear form $a\cdot ,\cdot$ defined by (19) is continuous on $H^1\times H^1$ and coercive on $H^1$, i.e.,$$\begin{array}{c}\text{(25)}& \text{i}au,v\le a_1{u}_{H^1}{v}_{H^1},\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}u,v\in H^1,\text{ii}av,v\ge a_0{v}_{H^1}^2,\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}v\in H^1,\hfill \end{array}$$where $a_0=1/3\min 1,\max h_0,h_1$ and $a_1=1+2h_0+h_1$.

Lemma 3.

There exists the Hilbert orthonormal base $w_j$ of $L^2$ consisting of the eigenfunctions $w_j$ corresponding to the eigenvalue ${\lambda }_j$ such that$$\begin{array}{}\text{(26)}& \begin{array}{c}0<{\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_j\le \cdots ,\hspace{1em}\underset{j⟶+\infty }{\text{lim}}{\lambda }_j=+\infty ,\hfill \\ a{w}_j,v={\lambda }_j{w}_j,v\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}v\in H^1,j=\mathrm{1,2},\dots .\hfill \end{array}\end{array}$$

Furthermore, the sequence $w_j/\sqrt{{\lambda }_j}$ is also a Hilbert orthonormal base of $H^1$ with respect to the scalar product $a\cdot ,\cdot$.

On the contrary, we have $w_j$ satisfying the following boundary value problem:$$\begin{array}{}\text{(27)}& \begin{array}{c}-\Delta w_j={\lambda }_j{w}_j,\text{in\hspace{0.17em}}\mathrm{0,1},\hfill \\ w_{jx}0-h_0{w}_{j}0=w_{jx}1+h_1{w}_{j}1=0,w_j\in C^{\infty }\overline{\Omega }.\hfill \end{array}\end{array}$$

The proof of Lemma 3 can be found in Theorem 7.7 of [39], p.87, with $H=L^2$ and $a\cdot ,\cdot$ as defined by (19).

3. The High-Order Iterative Method

First, we make the following assumptions:$$\begin{array}{c}\text{(28)}& H_1{\tilde{u}}_0,{\tilde{u}}_1\in H^2\times H^1,H_{2}g\in H^{1}0,T^{\ast },\hfill \\ H_{3}f\in C^0\mathrm{0,1}\times {ℝ}_{+}\times ℝ\text{such\hspace{0.17em}that}\hfill \\ \text{i}{D}_3^{i}f\in C^0\mathrm{0,1}\times {ℝ}_{+}\times ℝ,\hspace{1em}1\le i\le N,\hfill \\ \text{ii}{D}_1{D}_3^{i}f\in C^0\mathrm{0,1}\times {ℝ}_{+}\times ℝ,\hspace{1em}0\le i\le N-1.\hfill \end{array}$$

Fix$T^{\ast }>0$. For each $M>0$ given, we set the constant $K_{M}f$ as follows:$$\begin{array}{}\text{(29)}& \begin{array}{c}{K}_{M}f={\displaystyle \sum _{i=0}^N{D_3^{i}f}_{C^0{\Omega }_M}}+{\displaystyle \sum _{i=1}^{N-1}{D_1{D}_3^{i}f}_{C^0{\Omega }_M}},\hfill \\ f_{C^0{\Omega }_M}=\text{sup}fx,t,u:x,t,u\in {\Omega }_M,\hfill \\ {\Omega }_M=\mathrm{0,1}\times 0,T^{\ast }\times -\sqrt{2}M,\sqrt{2}M.\hfill \end{array}\end{array}$$

For every $T\in 0,T^{\ast }$ and$M>0$, we put$$\begin{array}{}\text{(30)}& \begin{array}{c}WM,T=\begin{array}{c}v\in L^{\infty }0,T;H^2:v^{\prime }\in L^{\infty }0,T;H^1,v^{″}\in L^2{Q}_T,\\ \text{with\hspace{0.17em}}{v}_{L^{\infty }0,T;H^2},{v^{\prime }}_{L^{\infty }0,T;H^1},{v^{″}}_{L^2{Q}_T}\le M\end{array},\\ W_{1}M,T=v\in WM,T:v^{″}\in L^{\infty }0,T;L^2.\end{array}\end{array}$$

Now, we establish a recurrent sequence $u_m$ in which the first term is chosen as $u_0\equiv 0$, and suppose that$$\begin{array}{}\text{(31)}& u_{m-1}\in W_{1}M,T.\end{array}$$

The $m^{\text{th}}$ iteration $u_m$ can be defined by the following problem: find $u_m\in W_{1}M,Tm\ge 1$ satisfying the nonlinear variational problem$$\begin{array}{}\text{(32)}& \begin{array}{c}{u}_m^{″}t,w+a{u}_{m}t,w+\lambda {u_m^{\prime }t}^{q-2}{u}_m^{\prime }t,w={\displaystyle {\int }_0^{t}gt-sa{u}_{m}s,w\text{d}s}+F_{m}t,w,\forall w\in H^1,\hfill \\ u_{m}0={\tilde{u}}_0,u_m^{\prime }0={\tilde{u}}_1,\hfill \end{array}\end{array}$$where$$\begin{array}{}\text{(33)}& F_{m}x,t={\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}{D}_3^{i}fx,t,u_{m-1}{u_m-u_{m-1}}^i}.\end{array}$$

The existence of the above sequence $u_m$ is claimed by the following theorem.

Theorem 1.

Let $H_1-H_3$ hold. Then, there exist a constant $M>0$ depending on ${\tilde{u}}_0$, ${\tilde{u}}_1$, $h_0$, and $h_1$ and a constant $T>0$ depending on ${\tilde{u}}_0$, ${\tilde{u}}_1$, $g$, $f$, $h_0$, $h_1$, $q$, and $\lambda$ such that, for $u_0\equiv 0$, there exists a recurrent sequence $u_m\subset W_{1}M,T$ defined by (32) and (33).

Proof.

The proof of Theorem 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$$\begin{array}{}\text{(34)}& u_m^{k}t={\displaystyle {\sum }_{j=1}^k{c}_{mj}^k}t{w}_j,\end{array}$$where the coefficients $c_{mj}^k$ satisfy the following system of nonlinear differential equations:$$\begin{array}{}\text{(35)}& \begin{array}{c}{\ddot{u}}_m^{k}t,w_j+a{u}_m^{k}t,w_j+\lambda {{\dot{u}}_m^{k}t}^{q-2}{\dot{u}}_m^{k}t,w_j\hfill \\ ={{\displaystyle \int }}_0^{t}gt-sa{u}_m^{k}s,w_j\text{d}s+F_m^{k}t,w_j,\hspace{1em}1\le j\le k,\hfill \\ u_m^{k}0={\tilde{u}}_{0k},{\dot{u}}_m^{k}0={\tilde{u}}_{1k},\hfill \end{array}\end{array}$$in which$$\begin{array}{}\text{(36)}& \begin{array}{c}{\tilde{u}}_{0k}={\displaystyle \sum _{j=1}^k{\alpha }_j^k}{w}_j⟶{\tilde{u}}_0\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{H}^2,\hfill \\ {\tilde{u}}_{1k}={\displaystyle \sum _{j=1}^k{\beta }_j^k}{w}_j⟶{\tilde{u}}_1\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{H}^1,\hfill \end{array}\end{array}$$and$$\begin{array}{}\text{(37)}& F_m^{k}x,t={\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}{D}_3^{i}fx,t,u_{m-1}{u_m^k-u_{m-1}}^i}={\displaystyle \sum _{j=0}^{N-1}{A}_{j}x,t,u_{m-1}{u_m^k}^j},\end{array}$$with$$\begin{array}{}\text{(38)}& A_{j}x,t,u_{m-1}={\displaystyle \sum _{i=j}^{N-1}\frac{{-1}^{i-j}}{j!i-j!}}{D}_3^{i}fx,t,u_{m-1}{u}_{m-1}^{i-j}.\end{array}$$

Note that, by using (31) and standard methods in ordinary differential equations, we can prove that system (35) admits a unique solution $c_{mj}^{k}t$, $1\le j\le k$, on interval $0,T_m^k\subset 0,T$. The following estimates allow one to take $T_m^k=T$ independent of $m$ and $k$.

Step 2.

(a priori estimates). First, for all $j=1,\dots ,k$, multiplying the first equation in (35) by ${\dot{c}}_{mj}^{k}t$, summing on $j$, and integrating with respect to the time variable from 0 to $t$, we have$$\begin{array}{}\text{(39)}& X_m^{k}t=X_m^{k}0+2{\displaystyle {\int }_0^t\text{d}s{\displaystyle {\int }_0^{s}gs-\tau a{u}_m^k\tau ,{\dot{u}}_m^{k}s\text{d}\tau }}+2{\displaystyle {\int }_0^t{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s},\end{array}$$where$$\begin{array}{}\text{(40)}& X_m^{k}t={{\dot{u}}_m^{k}t}^2+{u_m^{k}t}_a^2+2\lambda {{\displaystyle \int }}_0^t{{\dot{u}}_m^{k}s}_{L^q}^q\text{d}s.\end{array}$$

Next, by replacing $w_j$ in the first equation in (35) by $-\mathrm{\Delta }{w}_j$, we obtain that$$\begin{array}{}\text{(41)}& a{\ddot{u}}_m^{k}t,w_j+\Delta u_m^{k}t,\Delta w_j+\lambda a{{\dot{u}}_m^{k}t}^{q-2}{\dot{u}}_m^{k}t,w_j={{\displaystyle \int }}_0^{t}gt-\tau \Delta u_m^k\tau ,\Delta w_j\text{d}\tau +a{F}_m^{k}t,w_j,\hspace{1em}1\le j\le k,\end{array}$$similar to the first equation in (35), and it yields$$\begin{array}{c}\text{(42)}& Y_m^{k}t=Y_m^{k}0+2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^{s}gs-\tau \Delta u_m^k\tau ,\Delta {\dot{u}}_m^{k}s\text{d}\tau }+2{\displaystyle {\int }_0^{t}a{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s},=Y_m^{k}0+2{\displaystyle {\int }_0^{t}gt-\tau \Delta u_m^k\tau ,\Delta u_m^{k}t\text{d}\tau }-2g0{\displaystyle {\int }_0^t{\Delta u_m^{k}s}^2\text{d}s}-2{\displaystyle {\int }_0^t\text{d}s}{{\displaystyle \int }}_0^s{g}^{\prime }s-\tau \Delta u_m^k\tau ,\Delta u_m^{k}s\text{d}\tau +2{\displaystyle {\int }_0^{t}a{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s},\end{array}$$with$$\begin{array}{}\text{(43)}& Y_m^{k}t={{\dot{u}}_m^{k}t}_a^2+{\Delta u_m^{k}t}^2+\frac{8\lambda }{q^2}q-1{\displaystyle {\int }_0^t{\frac{\partial }{\partial x}{{\dot{u}}_m^{k}s}^{q-2/2}{\dot{u}}_m^{k}s}^2\text{d}s}+2\lambda h_0{\displaystyle {\int }_0^t{{\dot{u}}_m^{k}0,s}^q\text{d}s}+2\lambda h_1{\displaystyle {\int }_0^t{{\dot{u}}_m^{k}1,s}^q\text{d}s}.\end{array}$$

Putting$$\begin{array}{c}\text{(44)}& S_m^{k}t=X_m^{k}t+Y_m^{k}t+{\displaystyle {\int }_0^t{{\ddot{u}}_m^{k}s}^2\text{d}s}={{\dot{u}}_m^{k}t}^2+{{\dot{u}}_m^{k}t}_a^2+{u_m^{k}t}_a^2+{\Delta u_m^{k}t}^2\\ +2\lambda {\displaystyle {\int }_0^t{{\dot{u}}_m^{k}s}_{L^q}^q\text{d}s}+\frac{8\lambda }{q^2}q-1{\displaystyle {\int }_0^t{\frac{\partial }{\partial x}{{\dot{u}}_m^{k}s}^{q-2/2}{\dot{u}}_m^{k}s}^2\text{d}s}\\ +2\lambda h_0{\displaystyle {\int }_0^t{{\dot{u}}_m^{k}0,s}^q\text{d}s}+2\lambda h_1{\displaystyle {\int }_0^t{{\dot{u}}_m^{k}1,s}^q\text{d}s}+{\displaystyle {\int }_0^t{{\ddot{u}}_m^{k}s}^2\text{d}s},\end{array}$$it follows from (39), (42), and (44) that$$\begin{array}{c}\text{(45)}& S_m^{k}t=S_m^{k}0-2g0{\displaystyle {\int }_0^t{\Delta u_m^{k}s}^2\text{d}s}+2{\displaystyle {\int }_0^{t}gt-\tau \Delta u_m^k\tau ,\Delta u_m^{k}t\text{d}\tau }+2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^{s}gs-\tau a{u}_m^k\tau ,{\dot{u}}_m^{k}s\text{d}\tau }-2{\displaystyle {\int }_0^t\text{d}s}{{\displaystyle \int }}_0^s{g}^{\prime }s-\tau \Delta u_m^k\tau ,\Delta u_m^{k}s\text{d}\tau \\ +{\displaystyle {\int }_0^t{{\ddot{u}}_m^{k}s}^2\text{d}s}+2{\displaystyle {\int }_0^t{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s}+2{\displaystyle {\int }_0^{t}a{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s}=S_m^{k}0+{\displaystyle \sum _{j=1}^7{I}_j}.\end{array}$$

We estimate $S_m^{k}0$ and the integrals on the right-hand side of (45) as follows.

First integral $I_1$: it is clear that$$\begin{array}{}\text{(46)}& I_1=-2g0{\displaystyle {\int }_0^t{\Delta u_m^{k}s}^2\text{d}s}\le 2g0{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}.\end{array}$$

Second integral $I_2$: by the inequality $2ab\le 1/2a^2+2b^2$, $\forall a$, $b\in ℝ$, we have$$\begin{array}{}\text{(47)}& I_2=2{\displaystyle {\int }_0^{t}gt-\tau \Delta u_m^k\tau ,\Delta u_m^{k}t\text{d}\tau }\le 1/2S_m^{k}t+2g_{L^{2}0,T^{\ast }}^2{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}.\end{array}$$

Third integral $I_3$: using the following inequality $au,v\le \sqrt{au,u}\sqrt{av,v}=u_a{v}_a$ for all $u$, $v\in H^1$, we obtain$$\begin{array}{c}\text{(48)}& I_3=2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^{s}gs-\tau a{u}_m^k\tau ,{\dot{u}}_m^{k}s\text{d}\tau }\le 2{\displaystyle {\int }_0^t\text{d}s}{{\displaystyle \int }}_0^{s}gs-\tau {u_m^k\tau }_a{{\dot{u}}_m^{k}s}_a\text{d}\tau \\ \le 2{{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}}^{1/2}{{\displaystyle {\int }_0^t\text{d}s}{{\displaystyle {\int }_0^{s}gs-\tau \sqrt{S_m^k\tau }\text{d}\tau }}^2}^{1/2}\\ \le 2\sqrt{T^{\ast }}{g}_{L^{2}0,T^{\ast }}{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}.\end{array}$$

Fourth integral $I_4$:$$\begin{array}{c}\text{(49)}& I_4=-2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^s{g}^{\prime }s-\tau \Delta u_m^k\tau ,\Delta u_m^{k}s\text{d}\tau }\le 2\sqrt{T^{\ast }}{g^{\prime }}_{L^{2}0,T^{\ast }}{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}.\end{array}$$

Fifth integral $I_5$: note that the first equation in (35) can be written as follows:$$\begin{array}{c}\text{(50)}& {\ddot{u}}_m^{k}t,w_j+\lambda {{\dot{u}}_m^{k}t}^{q-2}{\dot{u}}_m^{k}t,w_j=\Delta u_m^{k}t,w_j-{\displaystyle {\int }_0^{t}gt-s\Delta u_m^{k}s,w_j\text{d}s}+F_m^{k}t,w_j,\hspace{1em}1\le j\le k.\end{array}$$

Hence, replacing $w_j$ with ${\ddot{u}}_m^{k}t$ and using the inequality ${a+b+c}^2\le 3a^2+b^2+c^2$ for all $a$, $b$, $c\in ℝ$, we obtain$$\begin{array}{c}\text{(51)}& {{\ddot{u}}_m^{k}t}^2+\frac{\lambda }{q}\frac{d}{\text{d}t}{{\dot{u}}_m^{k}t}_{L^q}^q=\Delta u_m^{k}t,{\ddot{u}}_m^{k}t-{{\displaystyle \int }}_0^{t}gt-\tau \Delta u_m^k\tau ,{\ddot{u}}_m^{k}t\text{d}\tau +F_m^{k}t,{\ddot{u}}_m^{k}t\le \Delta u_m^{k}t+{{\displaystyle \int }}_0^{t}gt-\tau \Delta u_m^k\tau \text{d}\tau +F_m^{k}t{\ddot{u}}_m^{k}t\\ \le \frac{3}{2}{\Delta u_m^{k}t}^2+g_{L^{2}0,T^{\ast }}^2{{\displaystyle \int }}_0^t{\Delta u_m^k\tau }^2\text{d}\tau +{F_m^{k}t}^2+\frac{1}{2}{{\ddot{u}}_m^{k}t}^2\\ \le \frac{3}{2}{S}_m^{k}t+g_{L^{2}0,T^{\ast }}^2{{\displaystyle \int }}_0^t{S}_m^k\tau \text{d}\tau +{F_m^{k}t}^2+\frac{1}{2}{{\ddot{u}}_m^{k}t}^2.\end{array}$$

This implies that$$\begin{array}{}\text{(52)}& {{\ddot{u}}_m^{k}t}^2+\frac{2\lambda }{q}\frac{d}{dt}{{\dot{u}}_m^{k}t}_{L^q}^q\le 3S_m^{k}t+g_{L^{2}0,T^{\ast }}^2{\displaystyle {\int }_0^t{S}_m^k\tau \text{d}\tau }+{F_m^{k}t}^2.\end{array}$$

Integrating in $t$, we have$$\begin{array}{c}\text{(53)}& I_5={\displaystyle {\int }_0^t{{\ddot{u}}_m^{k}s}^2\text{d}s}\le {\displaystyle {\int }_0^t{{\ddot{u}}_m^{k}s}^2\text{d}s}+\frac{2\lambda }{q}{{\dot{u}}_m^{k}t}_{L^q}^q\le \frac{2\lambda }{q}{{\tilde{u}}_{1k}}_{L^q}^q+3{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}+3T^{\ast }{g}_{L^{2}0,T^{\ast }}^2{\displaystyle {\int }_0^t{S}_m^k\tau \text{d}\tau }+3{\displaystyle {\int }_0^t{F_m^{k}s}^2\text{d}s}.\end{array}$$

Sixth integral $I_6$:$$\begin{array}{}\text{(54)}& I_6=2{\displaystyle {\int }_0^t{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s}\le 2{\displaystyle {\int }_0^t{F}_m^{k}s{\dot{u}}_m^{k}s\text{d}s}\le {\displaystyle {\int }_0^t{F_m^{k}s}^2\text{d}s}+{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}.\end{array}$$

Seventh integral $I_7$:$$\begin{array}{}\text{(55)}& I_7=2{\displaystyle {\int }_0^{t}a{F}_m^{k}s,{\dot{u}}_m^{k}s\text{d}s}\le 2{\displaystyle {\int }_0^t{F_m^{k}s}_a{{\dot{u}}_m^{k}s}_a\text{d}s}\le {\displaystyle {\int }_0^t{F_m^{k}s}_a^2\text{d}s}+{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}.\end{array}$$

Combining (45), (46)–(49), and (53)–(55), it leads to$$\begin{array}{}\text{(56)}& S_m^{k}t\le 2S_m^{k}0+\frac{4\lambda }{q}{{\tilde{u}}_{1k}}_{L^q}^q+D_T{\displaystyle {\int }_0^t{S}_m^{k}s\text{d}s}+8{\displaystyle {\int }_0^t{F_m^{k}s}^2\text{d}s}+2{\displaystyle {\int }_0^t{F_m^{k}s}_a^2\text{d}s},\end{array}$$where$$\begin{array}{}\text{(57)}& D_T=10+4g0+22+3T^{\ast }{g}_{L^{2}0,T^{\ast }}^2+4\sqrt{T^{\ast }}{g}_{L^{2}0,T^{\ast }}+{g^{\prime }}_{L^{2}0,T^{\ast }}.\end{array}$$

To estimate integrals ${\displaystyle {\int }_0^t{F_m^{k}s}^2\text{d}s}$ and ${\displaystyle {\int }_0^t{F_m^{k}s}_a^2\text{d}s}$, we use the following lemma.

Lemma 4.

The following inequalities are valid:$$\begin{array}{c}\text{(58)}& \text{i}{F_m^{k}t}_{L^{\infty }}\le {\tilde{A}}_{M}1+{\sqrt{S_m^{k}t}}^{N-1},\text{ii}{F_m^{k}t}_a\le \sqrt{a_1}{\tilde{A}}_M+{\tilde{B}}_{M}1+{\sqrt{S_m^{k}t}}^{N-1},\hfill \end{array}$$where ${\tilde{A}}_M$ and ${\tilde{B}}_M$ are defined as follows:$$\begin{array}{}\text{(59)}& \begin{array}{c}{\tilde{A}}_M=K_{M}f{\displaystyle \sum _{i=0}^{N-1}{\tilde{a}}_i},{\tilde{B}}_M=K_{M}f{\displaystyle \sum _{i=0}^{N-1}{\tilde{b}}_i},\hfill \\ {\tilde{a}}_i=\begin{array}{cc}1+\frac{1}{2}{\displaystyle \sum _{i=1}^{N-1}\frac{{2\sqrt{2}}^i}{i!}{M}^i},\hfill & i=0,\hfill \\ \frac{1}{2}\frac{1}{i!}{\frac{2\sqrt{2}}{\sqrt{a_0}}}^i,\hfill & i\ge 1,\hfill \end{array}\hfill \\ {\tilde{b}}_i=\frac{2^{i-1}}{i!}1+M{\sqrt{\frac{2}{a_0}}}^i+i{\sqrt{\frac{2}{a_0}}}^{i-1},\hspace{1em}i\ge 1,\hfill \\ {\tilde{b}}_0=1+M+{\displaystyle \sum _{i=1}^{N-1}{\tilde{b}}_i{M}^i}.\hfill \end{array}\end{array}$$

Proof of Lemma 4.

 

proof.

of (i). Using the inequalities ${a+b}^i\le 2^{i-1}{a}^i+b^i$, for all $a$, $b\ge 0$, $i\ge 1$, and $s^i\le 1+s^q$, $\forall s\ge 0$, $\forall i,q$, $0\le i\le q$, we have$$\begin{array}{c}\text{(60)}& {\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}{u_m^k+u_{m-1}}^i}\le 1+{\displaystyle \sum _{i=1}^{N-1}\frac{\sqrt{2^i}}{i!}{{u_m^{k}t}_{H^1}+M}^i}\le 1+{\displaystyle \sum _{i=1}^{N-1}\frac{{\sqrt{2}}^i}{i!}{2}^{i-1}{u_m^{k}t}_{H^1}^i+M^i}\\ \le 1+\frac{1}{2}{\displaystyle \sum _{i=1}^{N-1}\frac{{2\sqrt{2}}^i}{i!}{M}^i}+\frac{1}{2}{\displaystyle \sum _{i=1}^{N-1}\frac{{2\sqrt{2}}^i}{i!}\frac{1}{\sqrt{a_0^i}}{\sqrt{S_m^{k}t}}^i}\\ ={\displaystyle \sum _{i=0}^{N-1}{\tilde{a}}_i{\sqrt{S_m^{k}t}}^i\le {\displaystyle \sum _{i=0}^{N-1}{\tilde{a}}_{i}1+{\sqrt{S_m^{k}t}}^{N-1}}},\end{array}$$with ${\tilde{a}}_i$, $0\le i\le N-1$, defined by (59). Hence,$$\begin{array}{}\text{(61)}& F_m^{k}x,t\le {\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}{D}_3^{i}fx,t,u_{m-1}{u_m^k-u_{m-1}}^i}\le K_{M}f{\displaystyle \sum _{i=0}^{N-1}\frac{1}{i!}{u_m^k+u_{m-1}}^i}\le K_{M}f{\displaystyle \sum _{i=0}^{N-1}{\tilde{a}}_{i}1+{\sqrt{S_m^{k}t}}^{N-1}}={\tilde{A}}_{M}1+{\sqrt{S_m^{k}t}}^{N-1},\end{array}$$and it implies that equation (i) in (58) holds.

proof.

of (ii). We also have$$\begin{array}{}\text{(62)}& F_{mx}^{k}x,t=D_{1}fx,t,u_{m-1}+D_{3}fx,t,u_{m-1}\nabla u_{m-1}+{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}{D}_1{D}_3^{i}fx,t,u_{m-1}+D_3^{i+1}fx,t,u_{m-1}\nabla u_{m-1}{u_m^k-u_{m-1}}^i+D_3^{i}fx,t,u_{m-1}i{u_m^k-u_{m-1}}^{i-1}{u}_{mx}^k-\nabla u_{m-1}}.\end{array}$$

Hence,$$\begin{array}{}\text{(63)}& F_{mx}^{k}x,t\le K_{M}f1+\nabla u_{m-1}+K_{M}f{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}1+\nabla u_{m-1}{u_m^k-u_{m-1}}^i+i{u_m^k-u_{m-1}}^{i-1}{u}_{mx}^k-\nabla u_{m-1}}.\end{array}$$

It follows that$$\begin{array}{c}\text{(64)}& F_{mx}^{k}t\le K_{M}f1+\nabla u_{m-1}+K_{M}f{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}1+\nabla u_{m-1}{u_m^k-u_{m-1}}_{C^0\mathrm{0,1}}^i+i{u_m^k-u_{m-1}}_{C^0\mathrm{0,1}}^{i-1}{u}_{mx}^k-\nabla u_{m-1}}\le K_{M}f1+M+{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}1+M{\sqrt{\frac{2}{a_0}}}^i+i{\sqrt{\frac{2}{a_0}}}^{i-1}{u_m^k-u_{m-1}}_a^i}\\ \le K_{M}f1+M+{\displaystyle \sum _{i=1}^{N-1}\frac{2^{i-1}}{i!}1+M{\sqrt{\frac{2}{a_0}}}^i+i{\sqrt{\frac{2}{a_0}}}^{i-1}{u_m^k}_a^i+{u_{m-1}}_a^i}\\ \le K_{M}f1+M+{\displaystyle \sum _{i=1}^{N-1}{\tilde{b}}_i{\sqrt{S_m^{k}t}}^i+M^i}\\ =K_{M}f1+M+{\displaystyle \sum _{i=1}^{N-1}{\tilde{b}}_i{M}^i+{\displaystyle \sum _{i=1}^{N-1}{\tilde{b}}_i{\sqrt{S_m^{k}t}}^i}}\\ =K_{M}f{\displaystyle \sum _{i=0}^{N-1}{\tilde{b}}_i{\sqrt{S_m^{k}t}}^i}\le K_{M}f{\displaystyle \sum _{i=0}^{N-1}{\tilde{b}}_{i}1+{\sqrt{S_m^{k}t}}^{N-1}}\\ ={\tilde{B}}_{M}1+{\sqrt{S_m^{k}t}}^{N-1},\end{array}$$with ${\tilde{b}}_i$, $0\le i\le N-1$, and ${\tilde{B}}_M$ defined by (59).

Hence,$$\begin{array}{c}\text{(65)}& {F_m^{k}t}_a\le \sqrt{a_1}{F_m^{k}t}_{H^1}=\sqrt{a_1}{{F_m^{k}t}^2+{F_{mx}^{k}t}^2}^{1/2}\le \sqrt{a_1}{F}_m^{k}t+F_{mx}^{k}t\le \sqrt{a_1}{\tilde{A}}_M+{\tilde{B}}_{M}1+{\sqrt{S_m^{k}t}}^{N-1}.\end{array}$$

Therefore, equation (ii) in (58) follows. Lemma 4 is proved.

Now, the integrals ${\displaystyle {\int }_0^t{F_m^{k}s}^2\text{d}s}$ and ${\displaystyle {\int }_0^t{F_m^{k}s}_a^2\text{d}s}$ are estimated as follows:$$\begin{array}{c}\text{(66)}& {\displaystyle {\int }_0^t{F_m^{k}s}^2\text{d}s\le {\tilde{A}}_M^2{\displaystyle {\int }_0^t{1+{\sqrt{S_m^{k}s}}^{N-1}}^2\text{d}s}\le 2{\tilde{A}}_M^2{\displaystyle {\int }_0^{t}1+{S_m^{k}s}^{N-1}\text{d}s},}{\displaystyle {\int }_0^t{F_m^{k}s}_a^2\text{d}s}\le 2a_1{{\tilde{A}}_M+{\tilde{B}}_M}^2{\displaystyle {\int }_0^{t}1+{S_m^{k}s}^{N-1}\text{d}s}.\end{array}$$

Combining (60) and (66), we have$$\begin{array}{}\text{(67)}& S_m^{k}t\le 2S_m^{k}0+\frac{4\lambda }{q}{{\tilde{u}}_{1k}}_{L^q}^q+{\text{TC}}_{1}M,T+C_{1}M,T{\displaystyle {\int }_0^t{S_m^{k}s}^N\text{d}s},\end{array}$$where$$\begin{array}{}\text{(68)}& C_{1}M,T=D_T+16{\tilde{A}}_M^2+4a_1{{\tilde{A}}_M+{\tilde{B}}_M}^2.\end{array}$$

By means of convergences (36), we can deduce the existence of a constant $M>0$ independent of $k$ and $m$ such that$$\begin{array}{}\text{(69)}& 2S_m^{k}0+\frac{4\lambda }{q}{{\tilde{u}}_{1k}}_{L^q}^q\le \frac{M^2}{4},\hspace{1em}\forall m,k\in \mathbb{N}.\end{array}$$

Finally, it follows from (67) and (69) that$$\begin{array}{}\text{(70)}& S_m^{k}t\le \frac{M^2}{4}+{\text{TC}}_{1}M,T+C_{1}M,T{\displaystyle {\int }_0^t{S_m^{k}s}^N\text{d}s},\hspace{1em}0\le t\le T_m^k\le T.\end{array}$$

Then, by solving nonlinear Volterra integral inequality (70) (based on the methods in [40]), the following lemma is proved.

Lemma 5.

There exists a constant $T>0$ independent of $k$ and $m$ such that$$\begin{array}{}\text{(71)}& S_m^{k}t\le M^2,\hspace{1em}\forall t\in 0,T,\forall m,k\in \mathbb{N}.\end{array}$$

By Lemma 5, we can take constant $T_m^k=T$ for all $k$ and $m\in \mathbb{N}$. Thus, we have$$\begin{array}{}\text{(72)}& u_m^k\in W_{1}M,T,\hspace{1em}\forall m,k\in \mathbb{N}.\end{array}$$

Step 3.

(limiting process). Thanks to (72), there exists a subsequence $u_m^{k_j}$ of $u_m^k$, still denoted by $u_m^k$, such that$$\begin{array}{}\text{(73)}& \begin{array}{cc}{u}_m^k⟶u_m\hfill & \text{in\hspace{0.17em}}{L}^{\infty }0,T;H^2{\mathrm{weakly}}^{\ast },\hfill \\ {\dot{u}}_m^k⟶u_m^{\prime }\hfill & \text{in\hspace{0.17em}}{L}^{\infty }0,T;H^1{\mathrm{weakly}}^{\ast },\hfill \\ {\ddot{u}}_m^k⟶u_m^{″}\hfill & \text{in\hspace{0.17em}}{L}^2{Q}_T\text{\hspace{0.17em}weakly,}\hfill \\ u_m\in WM,T.\hfill & \hfill \end{array}\end{array}$$

Using the compactness lemma of Lions ([41], p.57) and applying the Fischer–Riesz theorem, from (73), there exists a subsequence of $u_m^k$, denoted by the same symbol satisfying$$\begin{array}{}\text{(74)}& \begin{array}{cc}{u}_m^k⟶u_m\hfill & \text{strongly\hspace{0.17em}in\hspace{0.17em}}{L}^{2}0,T;H^1\text{\hspace{0.17em}and\hspace{0.17em}a}.\text{e}.\text{\hspace{0.17em}in\hspace{0.17em}}{Q}_T\text{,}\hfill \\ {\dot{u}}_m^k⟶u_m^{\prime }\hfill & \text{strongly\hspace{0.17em}in\hspace{0.17em}}{L}^2{Q}_T\text{\hspace{0.17em}and\hspace{0.17em}a}.\text{e}.\text{\hspace{0.17em}in\hspace{0.17em}}{Q}_T.\hfill \end{array}\underset{x⟶\infty }{\lim }\end{array}$$

By using the following inequality,$$\begin{array}{}\text{(75)}& {{\dot{u}}_m^{k}t}^{q-2}{\dot{u}}_m^{k}t-{u_m^{\prime }t}^{q-2}{u}_m^{\prime }t\le q-1{\sqrt{2}M}^{q-2}{\dot{u}}_m^{k}t-u_m^{\prime }t,\hspace{1em}q\ge 2,\end{array}$$it follows from (44), (71), and (74)${}_2$ that$$\begin{array}{}\text{(76)}& {{\dot{u}}_m^{k}t}^{q-2}{\dot{u}}_m^{k}t\hookrightarrow {u_m^{\prime }t}^{q-2}{u}_m^{\prime }t\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{L}^2{Q}_T.\end{array}$$

On the contrary, by $L^{\infty }0,T;H^1\hookrightarrow L^{\infty }{Q}_T$ and using the inequality$$\begin{array}{}\text{(77)}& a^j-b^j\le j{M}^{j-1}a-b,\hspace{1em}\forall a,b\in -M,M,\forall M>0,\forall j\in \mathbb{N},\end{array}$$we deduce from (71) that$$\begin{array}{}\text{(78)}& {u_m^k}^j-u_m^j\le j{M}^{j-1}{u}_m^k-u_m,\hspace{1em}0\le j\le N-1.\end{array}$$

Therefore, (74)${}_1$ and (78) lead to$$\begin{array}{}\text{(79)}& {u_m^k}^j⟶u_m^j\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{L}^2{Q}_T.\end{array}$$

We note that$$\begin{array}{}\text{(80)}& A_{j}x,t,u_{m-1}t\le K_{M}f{\displaystyle \sum _{i=j}^{N-1}\frac{1}{j!i-j!}{\sqrt{2}M}^{i-j}}\equiv {\overline{D}}_{j}M.\end{array}$$

By (33), (37), and (71), it gives$$\begin{array}{}\text{(81)}& F_m^{k}t-F_{m}t\le {\displaystyle \sum _{j=0}^{N-1}{\overline{D}}_{j}M}{u_m^{k}t}^j-u_m^{j}t.\end{array}$$

Hence, we have$$\begin{array}{}\text{(82)}& {F_m^k-F_m}_{L^2{Q}_T}^2\le N{\displaystyle \sum _{j=0}^{N-1}{\overline{D}}_j^{2}M{{u_m^k}^j-u_m^j}_{L^2{Q}_T}^2},\end{array}$$so$$\begin{array}{}\text{(83)}& F_m^k⟶F_m\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{L}^2{Q}_T.\end{array}$$

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$$\begin{array}{}\text{(84)}& u_m^{\prime \prime }=\Delta u_m-\lambda {u_m^{\prime }}^{q-2}{u}_m^{\prime }-{\displaystyle {\int }_0^{t}gt-s\Delta u_{m}s\text{d}s}+F_m\in L^{\infty }0,T;L^2.\end{array}$$

Hence, $u_m\in W_{1}M,T$, and Theorem 1 is proved.

Next, the main result is also given by the following theorem. We consider the space $W_{1}T$, defined by$$\begin{array}{}\text{(85)}& W_{1}T=C0,T;H^1\cap C^{1}0,T;L^2,\end{array}$$and then $W_{1}T$ is a Banach space with respect to the norm$$\begin{array}{}\text{(86)}& v_{W_{1}T}=v_{C0,T;H^1}+{v^{\prime }}_{C^{1}0,T;L^2}.\end{array}$$

Theorem 2.

Let $H_1-H_3$ hold. Then, there exist constants $M>0$ and $T>0$ such that problems (1)–(3) have a unique weak solution $u\in W_{1}M,T$ and the recurrent sequence $u_m$ defined by (32) and (33) converges at the $N$-order rate to the solution $u$ strongly in $W_{1}T$ in the sense$$\begin{array}{}\text{(87)}& {u_m-u}_{W_{1}T}\le C{u_{m-1}-u}_{W_{1}T}^N,\end{array}$$for all $m\ge 1$, where $C$ is a suitable constant. Furthermore, the following estimate is fulfilled:$$\begin{array}{}\text{(88)}& {u_m-u}_{W_{1}T}\le C_T{{\gamma }_T}^{N^m},\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}m\in \mathbb{N},\end{array}$$where $C_T$ and $0<{\gamma }_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$. Then, $v_m$ satisfies the variational problem$$\begin{array}{}\text{(89)}& \begin{array}{c}\begin{array}{c}{v}_m^{″}t,w+a{v}_{m}t,w+\lambda {u_{m+1}^{\prime }t}^{q-2}{u}_{m+1}^{\prime }t-{u_m^{\prime }t}^{q-2}{u}_m^{\prime }t,w\\ ={\displaystyle {\int }_0^{t}gt-sa{v}_{m}s,w\text{d}s}+F_{m+1}t-F_{m}t,w,\hspace{1em}\forall w\in H^1,\end{array}\hfill \\ v_{m}0=v_m^{\prime }0=0.\hfill \end{array}\end{array}$$

Taking $w=v_m^{\prime }$ in (89), after integrating in $t$, we get$$\begin{array}{c}\text{(90)}& {\rho }_{m}t=-2\lambda {\displaystyle {\int }_0^t{u_{m+1}^{\prime }s}^{q-2}{u}_{m+1}^{\prime }s-{u_m^{\prime }s}^{q-2}{u}_m^{\prime }s,v_m^{\prime }s\text{d}s+2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^{s}gs-\tau a{v}_m\tau ,v_m^{\prime }s\text{d}\tau }+2{\displaystyle {\int }_0^t{F}_{m+1}s-F_{m}s,v_m^{\prime }s\text{d}s}}=-2\lambda {\displaystyle {\int }_0^t{u_{m+1}^{\prime }s}^{q-2}{u}_{m+1}^{\prime }s-{u_m^{\prime }s}^{q-2}{u}_m^{\prime }s,v_m^{\prime }s\text{d}s}+2{\displaystyle {\int }_0^{t}gt-\tau a{v}_m\tau ,v_{m}t\text{d}\tau }-2{\displaystyle {\int }_0^{t}g0a{v}_{m}s,v_{m}s\text{d}s}-2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^s{g}^{\prime }s-\tau a{v}_m\tau ,v_{m}s\text{d}\tau }+2{\displaystyle {\int }_0^t{F}_{m+1}s-F_{m}s,v_m^{\prime }s\text{d}s}\\ \equiv {\displaystyle \sum _{k=1}^5{J}_k},\end{array}$$with$$\begin{array}{}\text{(91)}& {\rho }_{m}t={v_m^{\prime }t}^2+{v_{m}t}_a^2.\end{array}$$

Next, we need to estimate the integrals on the right side of (90).

By the inequality$$\begin{array}{}\text{(92)}& {u_{m+1}^{\prime }x,s}^{q-2}{u}_{m+1}^{\prime }x,s-{u_m^{\prime }x,s}^{q-2}{u}_m^{\prime }x,s\le q-1{\sqrt{2}M}^{q-2}{v}_m^{\prime }x,s,\hspace{1em}q\ge 2,\end{array}$$we have$$\begin{array}{c}\text{(93)}& J_1=-2\lambda {\displaystyle {\int }_0^t{u_{m+1}^{\prime }s}^{q-2}{u}_{m+1}^{\prime }s-{u_m^{\prime }s}^{q-2}{u}_m^{\prime }s,v_m^{\prime }s\text{d}s}\le 2\lambda {\displaystyle {\int }_0^t{u_{m+1}^{\prime }s}^{q-2}{u}_{m+1}^{\prime }s-{u_m^{\prime }s}^{q-2}{u}_m^{\prime }s{v}_m^{\prime }s\text{d}s}\\ \le 2\lambda q-1{\sqrt{2}M}^{q-2}{\displaystyle {\int }_0^t{v_m^{\prime }s}^2\text{d}s}\\ \le 2\lambda q-1{\sqrt{2}M}^{q-2}{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s}.\end{array}$$

It is not difficult to estimate $J_2$, $J_3$, and $J_4$ as follows:$$\begin{array}{}\text{(94)}& J_2=2{\displaystyle {\int }_0^{t}gt-\tau a{v}_m\tau ,v_{m}t\text{d}\tau }\le \frac{1}{2}{\rho }_{m}t+2g_{L^{2}0,T^{\ast }}^2{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s},\text{(95)}& J_3=-2{\displaystyle {\int }_0^{t}g0a{v}_{m}s,v_{m}s\text{d}s}\le 2g0{\displaystyle {\int }_0^t{v_{m}s}_a^2\text{d}s}\le 2g0{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s},\text{(96)}& J_4=-2{\displaystyle {\int }_0^t\text{d}s}{\displaystyle {\int }_0^s{g}^{\prime }s-\tau a{v}_m\tau ,v_{m}s\text{d}\tau }\le 2\sqrt{T^{\ast }}{g^{\prime }}_{L^{2}0,T^{\ast }}{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s}.\end{array}$$

Using Taylor’s expansion of the function $fx,t,u_m=fx,t,u_{m-1}+v_{m-1}$ around the point $u_{m-1}$ up to order $N$, we obtain$$\begin{array}{}\text{(97)}& fx,t,u_m-fx,t,u_{m-1}={\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}{D}_3^{i}fx,t,u_{m-1}{v}_{m-1}^i+\frac{1}{N!}{D}_3^{N}fx,t,{\tilde{\theta }}_m{v}_{m-1}^N},\end{array}$$where ${\tilde{\theta }}_m={\tilde{\theta }}_{m}x,t=u_m+{\theta }_1{v}_{m-1}$, $0<{\theta }_1<1$.

Hence, it follows from (33) and (97) that$$\begin{array}{}\text{(98)}& F_{m+1}x,t-F_{m}x,t={\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}{D}_3^{i}fx,t,u_m{v}_m^i+\frac{1}{N!}{D}_3^{N}fx,t,{\tilde{\theta }}_m{v}_{m-1}^N},\end{array}$$and then$$\begin{array}{c}\text{(99)}& F_{m+1}t-F_{m}t\le K_{M}f{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}{\sqrt{2}{v_{m}t}_{H^1}}^i+\frac{1}{N!}{K}_{M}f{\sqrt{2}{v_{m-1}t}_{H^1}}^N}\le {\eta }_T^1\sqrt{{\rho }_{m}t}+{\eta }_T^2{v_{m-1}}_{W_{1}T}^N,\end{array}$$in which ${\eta }_T^1=K_{M}f/\sqrt{a_0}{\displaystyle {\sum }_{i=1}^{N-1}\sqrt{2}{\sqrt{2}M}^{i-1}/i!}$ and ${\eta }_T^2={\sqrt{2}}^N/N!K_{M}f$.

So,$$\begin{array}{c}\text{(100)}& J_5=2{\displaystyle {\int }_0^t{F}_{m+1}s-F_{m}s,v_m^{\prime }s\text{d}s}\le 2{\displaystyle {\int }_0^t{F}_{m+1}s-F_{m}s{v}_m^{\prime }s\text{d}s}\\ \le 2{\eta }_T^1{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s}+T{\eta }_T^2{v_{m-1}}_{W_{1}T}^{2N}+{\eta }_T^2{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s}.\end{array}$$

Then, it implies from (90), (93)–(96), and (100) that$$\begin{array}{}\text{(101)}& {\rho }_{m}t\le 2T{\eta }_T^2{v_{m-1}}_{W_{1}T}^{2N}+{\eta }_T^3{\displaystyle {\int }_0^t{\rho }_{m}s\text{d}s},\end{array}$$where ${\eta }_T^3=4\lambda q-1{\sqrt{2}M}^{q-2}+g0+g_{L^{2}0,T^{\ast }}^2+\sqrt{T^{\ast }}{g^{\prime }}_{L^{2}0,T^{\ast }}+{\eta }_T^1+1/2{\eta }_T^2$.

By using Gronwall’s lemma, (101) gives$$\begin{array}{}\text{(102)}& {v_m}_{W_{1}T}\le {\mu }_T{v_{m-1}}_{W_{1}T}^N,\end{array}$$where ${\mu }_T=1+1/\sqrt{a_0}\sqrt{2T{\eta }_T^2\exp T{\eta }_T^3}$.

Choosing $T>0$ small enough such that ${\gamma }_T=M{\mu }_T^{1/N-1}<1$, it follows from (102) that$$\begin{array}{}\text{(103)}& {u_m-u_{m+p}}_{W_{1}T}\le {1-{\gamma }_T}^{-1}{{\mu }_T}^{-1/N-1}{{\gamma }_T}^{N^m},\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}m\text{\hspace{0.17em}and\hspace{0.17em}}p\in \mathbb{N}.\end{array}$$

Hence, $u_m$ is a Cauchy sequence in $W_{1}T$. Then, there exists $u\in W_{1}T$ such that$$\begin{array}{}\text{(104)}& u_m⟶u\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{W}_{1}T.\end{array}$$

Note that $u_m\in W_{1}M,T$; then, there exists a subsequence $u_{m_j}$ of $u_m$ such that$$\begin{array}{}\text{(105)}& \begin{array}{cc}{u}_{m_j}⟶u\hfill & \text{in\hspace{0.17em}}{L}^{\infty }0,T;H^2{\mathrm{weakly}}^{\ast },\hfill \\ u_{m_j}^{\prime }⟶u^{\prime }\hfill & \text{in\hspace{0.17em}}{L}^{\infty }0,T;H^1{\mathrm{weakly}}^{\ast },\hfill \\ u_{m_j}^{\prime \prime }⟶u^{\prime \prime }\hfill & \text{in\hspace{0.17em}}{L}^2{Q}_T\text{\hspace{0.17em}weakly,}\hfill \\ u\in WM,T.\hfill & \hfill \end{array}\end{array}$$

Moreover, by (104) and the inequality$$\begin{array}{c}\text{(106)}& {u_m^{\prime }t}^{q-2}{u}_m^{\prime }t-{u^{\prime }t}^{q-2}{u}^{\prime }t\le q-1{\sqrt{2}M}^{q-2}{u}_m^{\prime }t-u^{\prime }t\le q-1{\sqrt{2}M}^{q-2}{u_m^{\prime }-u^{\prime }}_{C^{1}0,T;L^2}\\ \le q-1{\sqrt{2}M}^{q-2}{u_m-u}_{W_{1}T},\end{array}$$we have$$\begin{array}{}\text{(107)}& {u_m^{\prime }t}^{q-2}{u}_m^{\prime }t⟶{u^{\prime }t}^{q-2}{u}^{\prime }t\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{L}^{\infty }0,T;L^2.\end{array}$$

On the contrary,$$\begin{array}{c}\text{(108)}& F_m\cdot ,t-f\cdot ,t,ut\le f\cdot ,t,u_{m-1}t-f\cdot ,t,ut+{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}{D}_3^{i}f\cdot ,t,u_{m-1}{u_m-u_{m-1}}^i}\\ \le K_{M}f{u_{m-1}-u}_{W_{1}T}+{\displaystyle \sum _{i=1}^{N-1}\frac{1}{i!}{u_m-u_{m-1}}_{W_{1}T}^i}.\end{array}$$

Therefore, it implies from (104) and (38) that$$\begin{array}{}\text{(109)}& F_{m}t⟶f\cdot ,t,ut\text{\hspace{0.17em}strongly\hspace{0.17em}in\hspace{0.17em}}{L}^{\infty }0,T;L^2.\end{array}$$

Finally, passing to the limit in (32) and (33) as $m=m_j⟶\infty$, there exists $u\in WM,T$ satisfying the equation$$\begin{array}{}\text{(110)}& u^{\prime \prime }t,w+aut,w+\lambda {u^{\prime }t}^{q-2}{u}^{\prime }t,w+{\displaystyle {\int }_0^{t}gt-saus,w\text{d}s}=f\cdot ,t,ut,w,\end{array}$$for all $w\in H^1$, and the initial condition$$\begin{array}{}\text{(111)}& u0={\tilde{u}}_0,u^{\prime }0={\tilde{u}}_1.\end{array}$$

On the contrary, it follows from the fourth equation in (105) and (110) that$$\begin{array}{}\text{(112)}& u^{″}=\Delta u-\lambda {u^{\prime }}^{q-2}{u}^{\prime }+{\displaystyle {\int }_0^{t}gt-s\Delta us\text{d}s+fx,t,u\in L^{\infty }0,T;L^2}.\end{array}$$

Hence, $u\in W_{1}M,T$.

Proof.

(ii) Uniqueness: applying the similar estimations used in the proof of Theorem 1, it is easy to prove that $u\in W_{1}M,T$ is a unique local weak solution of problems (11)–(13).

Passing to the limit in (103) as $p⟶\infty$ for fixed $m$, we get (88). Also, with a similar argument, (87) follows. Theorem 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:$$\begin{array}{}\text{(113)}& f\in C^1\mathrm{0,1}\times {ℝ}_{+}\times ℝ,\end{array}$$see [42].

4. Numerical Results

In this section, we shall construct a numerical algorithm for the 2-order iterative scheme obtained by (14) and (15). The finite-difference method and some techniques for approximating double integrals will be used to find the finite-difference approximate solutions of $u^{m}x,t$ of this scheme. Moreover, a numerical algorithm for the single-iterative 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 single-iterative scheme and the 2-order iterative scheme.

Consider problems (11)–(13) with $\lambda =h_0=h_1=1$ and $q=2$; then, they are transformed into the following problem:$$\begin{array}{}\text{(114)}& \begin{array}{c}{u}_{tt}-\Delta u+u_t+{\displaystyle {\int }_0^{t}gt-s\Delta ux,s\text{d}s}=fx,t,u,\hspace{1em}0<x<\mathrm{1,0}<t<T,\hfill \\ u_{x}0,t-u0,t=u_{x}1,t+u1,t=0,\hfill \\ ux,0={\tilde{u}}_{0}x,u_{t}x,0={\tilde{u}}_{1}x,\hfill \end{array}\end{array}$$where the functions $g$, $f$, ${\tilde{u}}_0$, and ${\tilde{u}}_1$ are defined by$$\begin{array}{}\text{(115)}& \begin{array}{c}gt=e^{-2t},\hfill \\ fx,t,u=u^3+4e^{-2t}-8{1+x-x^2}^3{e}^{-3t},\hfill \\ {\tilde{u}}_{0}x=21+x-x^2,{\tilde{u}}_{1}x=-21+x-x^2.\hfill \end{array}\end{array}$$

The exact solution of problem (114), with $g$, $f$, ${\tilde{u}}_0$, and ${\tilde{u}}_1$ defined by (115), is a function $u_{ex}$ given by$$\begin{array}{}\text{(116)}& u_{ex}x,t=21+x-x^2{e}^{-t}.\end{array}$$

With datum (115) and $0\le x\le 1$ and $0\le t\le 0.5$, Figure 1 describes the surface of the exact solution $u_{ex}x,t$ of problem (114) below.

[figure omitted; refer to PDF]

We consider the 2-order iterative scheme of problem (114) as follows:$$\begin{array}{}\text{(117)}& \begin{array}{c}{u}_0\equiv 0,\hfill \\ {\ddot{u}}^{m}t+{\dot{u}}^{m}t-\Delta u^{m}t+{\displaystyle {\int }_0^{t}gt-s\Delta u^{m}s\text{d}s}-b^{m}x,t{u}^{m}x,t=F^{m}x,t,\hspace{1em}0<x<\mathrm{1,0}<t<T,\hfill \\ u_x^{m}0,t-u^{m}0,t=u_x^{m}1,t+u^{m}1,t=0,\hfill \\ u^{m}x,0={\tilde{u}}_{0}x,{\dot{u}}^{m}x,0={\tilde{u}}_{1}x,\hfill \end{array}\end{array}$$where$$\begin{array}{c}\text{(118)}& F^{m}x,t=fx,t,u^{m-1}x,t-b^{m}x,t{u}^{m-1}x,t,b^{m}x,t=D_{3}fx,t,u^{m-1}x,t.\end{array}$$

In order to solve problems (117) and (118) numerically, we first use the spatial difference grid $u_i^{m}t\equiv u^m{x}_i,t$, with $x_i=ih$, $h=1/N$, $i=\overline{0,N}$.

Replacing the derivatives in spatial variable $x$ of (117) by the following approximations,$$\begin{array}{c}\text{(119)}& \frac{\partial u^m}{\partial x}{x}_i,t\approx \frac{u_i^m-u_{i-1}^{mm}}{h},\Delta u^m{x}_i,t=\frac{{\partial }^2{u}^m}{\partial x^2}{x}_i,t\approx \frac{u_{i-1}^m-2u_i^m+u_{i+1}^m}{h^2},\end{array}$$problems (117) and (118) are turned into the following system of intergodifferential equations:$$\begin{array}{}\text{(120)}& \begin{array}{c}{\ddot{u}}_i^{m}t+{\dot{u}}_i^{m}t-\frac{u_{i-1}^m-2u_i^m+u_{i+1}^m}{h^2}\hfill \\ +{\displaystyle {\int }_0^{t}gt-s\frac{u_{i-1}^{m}s-2u_i^{m}s+u_{i+1}^{m}s}{h^2}\text{d}s}-b_i^{m}t{u}_i^{m}t=F_i^{m}t,\hspace{1em}i=\overline{1,N},\hfill \\ \frac{u_1^m-u_0^m}{h}-u_0^m=\frac{u_{N+1}^m-u_N^m}{h}+u_N^m=0,\hfill \\ u_i^{m}0={\tilde{u}}_0{x}_i\equiv {\tilde{u}}_{0i},{\dot{u}}_i^{m}0={\tilde{u}}_1{x}_i\equiv {\tilde{u}}_{1i},i=\overline{1,N},\hfill \end{array}\end{array}$$where$$\begin{array}{c}\text{(121)}& b_i^{m}t=b^m{x}_i,t=D_{3}f{x}_i,t,u_i^{m-1}t,F_i^{m}t=F^m{x}_i,t=f{x}_i,t,u_i^{m-1}t-b_i^{m}t{u}_i^{m-1}t,\hspace{1em}i=\overline{0,N}.\end{array}$$

Boundary conditions in the second equation in (120) lead to$$\begin{array}{}\text{(122)}& u_0^m=\frac{u_1^m}{1+h},u_{N+1}^m=1-h{u}_N^m.\end{array}$$

Eliminating the unknown functions $u_0^m$ and $u_{N+1}^m$ in the first equation ($i=1$) and the last equation ($i=N$), respectively, system (120) is equivalent to$$\begin{array}{}\text{(123)}& \begin{array}{c}{\ddot{u}}_1^{m}t+{\dot{u}}_1^{m}t+{\tilde{a}}_1{u}_1^m-\frac{u_2^m}{h^2}+{{\displaystyle \int }}_0^{t}gt-s-{\tilde{a}}_1{u}_1^{m}s+\frac{u_2^m}{h^2}s\text{d}s-b_1^{m}t{u}_1^{m}t=F_1^{m}t,\hspace{1em}i=1,\\ {\ddot{u}}_i^{m}t+{\dot{u}}_i^{m}t-\frac{u_{i-1}^m-2u_i^m+u_{i+1}^m}{h^2}+{{\displaystyle \int }}_0^{t}gt-s\frac{u_{i-1}^{m}s-2u_i^{m}s+u_{i+1}^{m}s}{h^2}\text{d}s-b_i^{m}t{u}_i^{m}t=F_i^{m}t,\hspace{1em}i=\overline{2,N-1},\\ {\ddot{u}}_N^{m}t+{\dot{u}}_N^{m}t-\frac{u_{N-1}^m-1+h{u}_N^m}{h^2}+{{\displaystyle \int }}_0^{t}gt-s\frac{u_{N-1}^{m}s-1+h{u}_N^{m}s}{h^2}\text{d}s-b_N^{m}t{u}_N^{m}t=F_N^{m}t,\hspace{1em}i=N,\\ u_i^{m}0={\tilde{u}}_0{x}_i\equiv {\tilde{u}}_{0i},{\dot{u}}_i^{m}0={\tilde{u}}_1{x}_i\equiv {\tilde{u}}_{1i},i=\overline{1,N},\end{array}\end{array}$$where ${\tilde{a}}_1=1+2h/h^{2}1+h$.

Rewrite (123) into a vector equation:$$\begin{array}{}\text{(124)}& \begin{array}{c}E\frac{{\text{d}}^2{\stackrel{⟶}{u}}^m}{\text{d}{t}^2}t+E\frac{\text{d}{\stackrel{⟶}{u}}^m}{\text{d}t}t+B^{m}t{\stackrel{⟶}{u}}^{m}t-{\displaystyle {\int }_0^{t}gt-sC{\stackrel{⟶}{u}}^{m}s\text{d}s}={\stackrel{⟶}{F}}^{m}t,\\ {\stackrel{⟶}{u}}^{m}0={u_1^{m}0,\dots ,u_N^{m}0}^T={{\tilde{u}}_0{x}_1,\dots ,{\tilde{u}}_0{x}_N}^T\equiv {\tilde{X}}_0\in {ℝ}^N,\\ \frac{\text{d}{\stackrel{⟶}{u}}^m}{\text{d}t}0={{\dot{u}}_1^{m}0,\dots ,{\dot{u}}_N^{m}0}^T={{\tilde{u}}_1{x}_1,\dots ,{\tilde{u}}_1{x}_N}^T\equiv {\tilde{X}}_1\in {ℝ}^N,\end{array}\end{array}$$where$$\begin{array}{}\text{(125)}& \begin{array}{c}{\stackrel{⟶}{u}}^{m}t={u_1^{m}t,\dots ,u_N^{m}t}^T={u^m{x}_1,t,\dots ,u^m{x}_N,t}^T,\\ {\stackrel{⟶}{F}}^{m}t={F_1^{m}t,\dots ,F_N^{m}t}^T,\end{array}\end{array}$$and $E$, $B^{m}t$, and $C\in {\mathfrak{M}}_N$ (${\mathfrak{M}}_N$ is the set of real $N$-order matrices) are defined by$$\begin{array}{c}\text{(126)}& E=\begin{array}{ccccccc}1& & & & & & \\ & 1& & & & & \\ & & 1& & & & \\ & & & \ddots & & & \\ & & & & \ddots & & \\ & & & & & 1& \\ & & & & & & 1\end{array},C=\begin{array}{ccccccc}{\tilde{a}}_1& \gamma & & & & & \\ \gamma & \tilde{a}& \gamma & & & & \\ & \gamma & \tilde{a}& \gamma & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & \ddots & \ddots & \ddots & \\ & & & & \gamma & \tilde{a}& \gamma \\ & & & & & \gamma & {\tilde{a}}_N\end{array},\\ {\tilde{B}}_{m}t=\begin{array}{ccccc}{b}_1^{m}t& & & & \\ & b_2^{m}t& & & \\ & & \ddots & & \\ & & & \ddots & \\ & & & & b_N^{m}t\end{array},\\ B^{m}t=C-{\tilde{B}}_{m}t,\\ {\tilde{a}}_1=\frac{1+2h}{h^{2}1+h},\\ {\tilde{a}}_N=\frac{1+h}{h^2},\\ \tilde{a}=\frac{2}{h^2},\\ \gamma =\frac{-1}{h^2}.\end{array}$$

Integrating in time variable $t$, we have$$\begin{array}{}\text{(127)}& \begin{array}{c}\frac{\text{d}{\stackrel{⟶}{u}}^m}{\text{d}t}t={\stackrel{⟶}{H}}^{m}t-e^{-t}{\displaystyle {\int }_0^t{e}^s{B}^{m}s{\stackrel{⟶}{u}}^{m}s\text{d}s}+e^{-t}{\displaystyle {\int }_0^t{e}^s\text{d}s}{\displaystyle {\int }_0^{t}gs-rC{\stackrel{⟶}{u}}^{m}r\text{d}r},\\ {\stackrel{⟶}{u}}^{m}0={\tilde{X}}_0,\end{array}\end{array}$$where$$\begin{array}{c}\text{(128)}& {\stackrel{⟶}{H}}^{m}t={\tilde{X}}_1{e}^{-t}+e^{-t}{\displaystyle {\int }_0^t{e}^s{\stackrel{⟶}{F}}^{m}s\text{d}s},{\tilde{X}}_0={{\tilde{u}}_0{x}_1,\dots ,{\tilde{u}}_0{x}_N}^T,\\ {\tilde{X}}_1={{\tilde{u}}_1{x}_1,\dots ,{\tilde{u}}_1{x}_N}^T.\end{array}$$

Approximating the derivatives $\text{d}{\stackrel{⟶}{u}}^m/\text{d}t{t}_j$ by the following time differences,$$\begin{array}{c}\text{(129)}& \frac{\text{d}{\stackrel{⟶}{u}}^m}{\text{d}t{t}_j}\simeq \frac{{\stackrel{⟶}{u}}_{j+1}^m-{\stackrel{⟶}{u}}_j^m}{\Delta t},{\stackrel{⟶}{u}}_j^m={\stackrel{⟶}{u}}^m{t}_j,t_j=j\Delta t,j=0,\dots ,M,\Delta t=\frac{T}{M},\\ {\stackrel{⟶}{H}}_j^m={\stackrel{⟶}{H}}^m{t}_j,\\ B_j^m=B^m{t}_j,\end{array}$$equation (127)${}_1$ is turned into$$\begin{array}{}\text{(130)}& \frac{{\stackrel{⟶}{u}}_{j+1}^m-{\stackrel{⟶}{u}}_j^m}{\Delta t}={\stackrel{⟶}{H}}_j^m-e^{-t_j}{\displaystyle {\int }_0^{t_j}{e}^s{B}^{m}s{\stackrel{⟶}{u}}^{m}s\text{d}s}+e^{-t_j}{\displaystyle {\int }_0^{t_j}{e}^s\text{d}s}{\displaystyle {\int }_0^{s}gs-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}.\end{array}$$

Note that the integral ${\displaystyle {\int }_0^{t_j}\mathrm{\Phi }s\text{d}s}$ can be approximated by the trapezoidal formula:$$\begin{array}{c}\text{(131)}& {\displaystyle {\int }_0^{t_j}\Phi s\text{d}s}\simeq \begin{array}{cc}\Delta t\frac{{\Phi }_0+{\Phi }_1}{2},\hfill & j=1,\hfill \\ \Delta t\frac{{\Phi }_0+{\Phi }_j}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{\Phi }_{\nu }},\hfill & 2\le j\le M,\hfill \end{array}{\Phi }_j=\Phi t_j.\end{array}$$

Applying (131) with $\mathrm{\Phi }s=e^s{B}^{m}s{\stackrel{⟶}{u}}^{m}s$ and $\mathrm{\Phi }s=g{t}_j-rC{\stackrel{⟶}{u}}^{m}r\text{d}r$, respectively, we have the following approximations:$$\begin{array}{c}\text{(132)}& {\displaystyle {\int }_0^{t_j}{e}^s{B}^{m}s{\stackrel{⟶}{u}}^{m}s\text{d}s\simeq {\Psi }_j^{1,m}}\equiv \begin{array}{c}\Delta t\frac{e^{t_0}{B}_0^m{\stackrel{⟶}{u}}_0^m+e^{t_1}{B}_1^m{\stackrel{⟶}{u}}_1^m}{2},\hspace{1em}j=1,\\ \Delta t\frac{e^{t_0}{B}_0^m{\stackrel{⟶}{u}}_0^m+e^{t_j}{B}_j^m{\stackrel{⟶}{u}}_j^m}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}{B}_{\nu }^m{\stackrel{⟶}{u}}_{\nu }^m},\hspace{1em}2\le j\le M,\end{array}{\displaystyle {\int }_0^{t_j}g{t}_j-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}\simeq \begin{array}{c}\Delta t\frac{g_{1}C{\stackrel{⟶}{u}}_0^m+g_{0}C{\stackrel{⟶}{u}}_1^m}{2},\hspace{1em}j=1,\\ \Delta t\frac{g_{j}C{\stackrel{⟶}{u}}_0^m+g_{0}C{\stackrel{⟶}{u}}_j^m}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{g}_{j-\nu }C{\stackrel{⟶}{u}}_{\nu }^m},\hspace{1em}2\le j\le M,\end{array}\\ g_j=g{t}_j,\hspace{1em}0\le j\le M.\end{array}$$

Again using (131) with $\mathrm{\Phi }s=e^s{\displaystyle {\int }_0^{s}gs-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}$, the double integral ${\displaystyle {\int }_0^{t_j}{e}^s\text{d}s}{\displaystyle {\int }_0^{s}gs-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}$ can be approximated as follows:$$\begin{array}{c}\text{(133)}& {\displaystyle {\int }_0^{t_j}{e}^s\text{d}s}{\displaystyle {\int }_0^{s}gs-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}={\displaystyle {\int }_0^{t_j}{e}^s{\displaystyle {\int }_0^{s}gs-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}\text{d}s}\simeq {\Psi }_j^{2,m}\\ \equiv \begin{array}{c}\frac{\Delta t}{2}{e}^{t_1}{\displaystyle {\int }_0^{t_1}g{t}_1-rC{\stackrel{⟶}{u}}^{m}r\text{d}r},\hspace{1em}j=1,\\ \Delta t\frac{1}{2}{e}^{t_j}{\displaystyle {\int }_0^{t_j}g{t}_j-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}}{\displaystyle {\int }_0^{t_{\nu }}g{t}_{\nu }-rC{\stackrel{⟶}{u}}^{m}r\text{d}r},\hspace{1em}2\le j\le M.\end{array}\end{array}$$

Similarly, we also obtain the following approximations:$$\begin{array}{c}\text{(134)}& {\stackrel{⟶}{H}}_j^m={\tilde{X}}_1{e}^{-t_j}+e^{-t_j}{\displaystyle {\int }_0^{t_j}{e}^s{\stackrel{⟶}{F}}^{m}s\text{d}s}\simeq {\Psi }_j^{3,m}\equiv \begin{array}{c}{\tilde{X}}_1{e}^{-t_1}+e^{-t_1}\Delta t\frac{e^{t_0}{\stackrel{⟶}{F}}_0^m+e^{t_1}{\stackrel{⟶}{F}}_1^m}{2},\hspace{1em}j=1,\\ {\tilde{X}}_1{e}^{-t_j}+e^{-t_j}\Delta t\frac{e^{t_0}{\stackrel{⟶}{F}}_0^m+e^{t_j}{\stackrel{⟶}{F}}_j^m}{2}+{\displaystyle \sum }_{\nu =1}^{j-1}{e}^{t_{\nu }}{\stackrel{⟶}{F}}_{\nu }^m,\hspace{1em}2\le j\le M,\end{array}\\ {\stackrel{⟶}{F}}_j^m={\stackrel{⟶}{F}}^m{t}_j.\end{array}$$

So, equation (130) can be rewritten by$$\begin{array}{}\text{(135)}& {\stackrel{⟶}{u}}_{j+1}^m={\stackrel{⟶}{u}}_j^m+\Delta t{\Psi }_j^{\ast ,m},\end{array}$$where$$\begin{array}{}\text{(136)}& {\Psi }_j^{\ast ,m}=-e^{-t_j}{\Psi }_j^{1,m}+e^{-t_j}{\Psi }_j^{2,m}+{\Psi }_j^{3,m}.\end{array}$$

Note that ${\mathrm{\Psi }}_j^{1,m}$ and ${\mathrm{\Psi }}_j^{2,m}$ depend on ${\stackrel{⟶}{u}}_0^m,{\stackrel{⟶}{u}}_1^m,\dots ,{\stackrel{⟶}{u}}_j^m$ and ${\mathrm{\Psi }}_j^{3,m}$ on ${\stackrel{⟶}{u}}_0^{m-1},{\stackrel{⟶}{u}}_1^{m-1},\dots ,{\stackrel{⟶}{u}}_j^{m-1}$, so$$\begin{array}{}\text{(137)}& {\Psi }_j^{\ast ,m}=-e^{-t_j}{\Psi }_j^{1,m}+e^{-t_j}{\Psi }_j^{2,m}+{\Psi }_j^{3,m},\end{array}$$where$$\begin{array}{c}\text{(138)}& {\Psi }_j^{3,m}={\tilde{X}}_1{e}^{-t_j}+e^{-t_j}\Delta t\frac{e^{t_0}{\stackrel{⟶}{F}}_0^m+e^{t_j}{\stackrel{⟶}{F}}_j^m}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}{\stackrel{⟶}{F}}_{\nu }^m},{\Psi }_j^{1,m}=\Delta t\frac{e^{t_0}{B}_0^m{\stackrel{⟶}{u}}_0^m+e^{t_j}{B}_j^m{\stackrel{⟶}{u}}_j^m}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}{B}_{\nu }^m{\stackrel{⟶}{u}}_{\nu }^m},\\ {\Psi }_j^{2,m}=\Delta t\frac{1}{2}{e}^{t_j}{\displaystyle {\int }_0^{t_j}g{t}_j-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}}{\displaystyle {\int }_0^{t_{\nu }}g{t}_{\nu }-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}.\end{array}$$

5. Description of Finite-Difference Algorithm (135) and (136)

(i) Let $M\text{\hspace{0.17em}and\hspace{0.17em}}N$ be fixed constants. At the first iteration with $m=0$, we set up the given vector

$$\begin{array}{}\text{(139)}& {\stackrel{⟶}{u}}_j^0={\stackrel{⟶}{u}}^0{t}_j={u_1^0{t}_j,\dots ,u_N^0{t}_j}^T={u^0{x}_1,t_j,\dots ,u^0{x}_N,t_j}^T\equiv 0,\hspace{1em}j=1,\dots ,M.\end{array}$$

(ii) At the ${m-1}^{\text{th}}$ iteration, suppose that we can compute

$$\begin{array}{}\text{(140)}& {\stackrel{⟶}{u}}_j^{m-1}={\stackrel{⟶}{u}}^{m-1}{t}_j={u_1^{m-1}{t}_j,\dots ,u_N^{m-1}{t}_j}^T\text{,}\hspace{1em}j=1,\dots ,M.\end{array}$$

(iii) The computation of the vectors ${\stackrel{⟶}{u}}_j^m={\stackrel{⟶}{u}}^m{t}_j={u_1^m{t}_j,\dots ,u_N^m{t}_j}^T$, $j=1,\dots ,M$, can be done consecutively by the following steps:

  C1. The computation of ${\stackrel{⟶}{u}}_1^m={u_1^m{t}_1,\dots ,u_N^m{t}_1}^T$.

(i) With the first given vector ${\stackrel{⟶}{u}}_0^m={\stackrel{⟶}{u}}^m{t}_0={{\tilde{u}}_0{x}_1,\dots ,{\tilde{u}}_0{x}_N}^T={\tilde{X}}_0$, we calculate ${\mathrm{\Psi }}_0^{\ast ,m}$ as follows:

$$\begin{array}{c}\text{(141)}& {\Psi }_0^{3,m}={\tilde{X}}_1,{\Psi }_0^{1,m}={\Psi }_0^{2,m}\equiv 0.\end{array}$$

  Hence,

$$\begin{array}{}\text{(142)}& {\Psi }_0^{\ast ,m}=-e^{-t_0}{\Psi }_0^{1,m}+e^{-t_0}{\Psi }_0^{2,m}+{\Psi }_0^{3,m}={\tilde{X}}_1.\end{array}$$

(ii) Finding ${\stackrel{⟶}{u}}_1^m={u_1^m{t}_1,\dots ,u_N^m{t}_1}^T$ by

$$\begin{array}{}\text{(143)}& {\stackrel{⟶}{u}}_1^m={\stackrel{⟶}{u}}_0^m+\Delta t{\Psi }_0^{\ast ,m}={\tilde{X}}_0+\Delta t{\tilde{X}}_1.\end{array}$$

  C2. The computation of ${\stackrel{⟶}{u}}_2^m={u_1^m{t}_2,\dots ,u_N^m{t}_2}^T$.

(i) Calculating ${\mathrm{\Psi }}_1^{\ast ,m}$:

$$\begin{array}{c}\text{(144)}& {\Psi }_1^{3,m}={\tilde{X}}_1{e}^{-t_1}+e^{-t_1}\Delta t\frac{e^{t_0}{\stackrel{⟶}{F}}_0^m+e^{t_1}{\stackrel{⟶}{F}}_1^m}{2},{\Psi }_1^{1,m}=\Delta t\frac{e^{t_0}{B}_0^m{\stackrel{⟶}{u}}_0^m+e^{t_1}{B}_1^m{\stackrel{⟶}{u}}_1^m}{2},\\ {\Psi }_1^{2,m}=\frac{\Delta t}{2}{e}^{t_1}{\displaystyle {\int }_0^{t_1}g{t}_1-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}.\end{array}$$

  Hence,$$\begin{array}{}\text{(145)}& {\Psi }_1^{\ast ,m}=-e^{-t_1}{\Psi }_1^{1,m}+e^{-t_1}{\Psi }_1^{2,m}+{\Psi }_1^{3,m}.\end{array}$$

(ii) Finding ${\stackrel{⟶}{u}}_2^m={u_1^m{t}_2,\dots ,u_N^m{t}_2}^T$ by

$$\begin{array}{}\text{(146)}& {\stackrel{⟶}{u}}_2^m={\stackrel{⟶}{u}}_1^m+\Delta t{\Psi }_1^{\ast ,m}.\end{array}$$

  C3. The computation of ${\stackrel{⟶}{u}}_{j+1}^m={u_1^m{t}_{j+1},\dots ,u_N^m{t}_{j+1}}^T$.

  Suppose that ${\stackrel{⟶}{u}}_1^m$, ${\stackrel{⟶}{u}}_2^m,\dots ,{\stackrel{⟶}{u}}_j^m$ have been calculated; then, we define ${\stackrel{⟶}{u}}_{j+1}^m={u_1^m{t}_{j+1},\dots ,u_N^m{t}_{j+1}}^T$ determined by recurrence as follows:

(i) Calculating ${\mathrm{\Psi }}_j^{\ast ,m}$:

$$\begin{array}{c}\text{(147)}& {\Psi }_j^{3,m}={\tilde{X}}_1{e}^{-t_j}+e^{-t_j}\Delta t\frac{e^{t_0}{\stackrel{⟶}{F}}_0^m+e^{t_j}{\stackrel{⟶}{F}}_j^m}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}{\stackrel{⟶}{F}}_{\nu }^m},{\Psi }_j^{1,m}=\Delta t\frac{e^{t_0}{B}_0^m{\stackrel{⟶}{u}}_0^m+e^{t_j}{B}_j^m{\stackrel{⟶}{u}}_j^m}{2}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}{B}_{\nu }^m{\stackrel{⟶}{u}}_{\nu }^m},\\ {\Psi }_j^{2,m}=\Delta t\frac{1}{2}{e}^{t_j}{\displaystyle {\int }_0^{t_j}g{t}_j-rC{\stackrel{⟶}{u}}^{m}r\text{d}r}+{\displaystyle \sum _{\nu =1}^{j-1}{e}^{t_{\nu }}}{\displaystyle {\int }_0^{t_{\nu }}g{t}_{\nu }-rC{\stackrel{⟶}{u}}^{m}r\text{d}r},\end{array}$$

  and

$$\begin{array}{}\text{(148)}& {\Psi }_j^{\ast ,m}=-e^{-t_j}{\Psi }_j^{1,m}+e^{-t_j}{\Psi }_j^{2,m}+{\Psi }_j^{3,m}.\end{array}$$

(ii) Finding ${\stackrel{⟶}{u}}_{j+1}^m={u_1^m{t}_{j+1},\dots ,u_N^m{t}_{j+1}}^T$ by

$$\begin{array}{}\text{(149)}& {\stackrel{⟶}{u}}_{j+1}^m={\stackrel{⟶}{u}}_j^m+\Delta t{\Psi }_j^{\ast ,m}.\end{array}$$

  When the process of the computations is reached to $j=M-1$, we obtain

$$\begin{array}{}\text{(150)}& {\stackrel{⟶}{u}}_j^m={u^m{x}_1,t_j,\dots ,u^m{x}_N,t_j}^T,\hspace{1em}1\le j\le M.\end{array}$$

  C4. The error of two successively approximate iterations.

$$\begin{array}{}\text{(151)}& {u^m-u^{m-1}}_{M,N}=\underset{1\le j\le M}{\max }\underset{1\le i\le N}{\max }{u}^m{x}_i,t_j-u^{m-1}{x}_i,t_j.\end{array}$$

  The process of iteration will be stopped at the $m^{\text{th}}$ step when the following estimate is satisfied:

$$\begin{array}{}\text{(152)}& {u^m-u^{m-1}}_{M,N}<{10}^{-4}.\end{array}$$

  C5. The error between the approximate solution (at the $m^{\text{th}}$ step) and the exact solution is defined by

$$\begin{array}{}\text{(153)}& E_{M,N}={u^m-u_{ex}}_{M,N}=\underset{1\le j\le M}{\max }\underset{1\le i\le N}{\max }{u}^m{x}_i,t_j-u_{ex}{x}_i,t_j,\end{array}$$

  in which $u_{ex}x,t=21+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).

[figure omitted; refer to PDF]

Table 1

Errors of the approximate solution and the exact solution.

		Single-iterative scheme	2-order iterative scheme
$N$	$M$	$E_{M,N}={u^m-u_{ex}}_{M,N}$	$E_{M,N}={u^m-u_{ex}}_{M,N}$
10	20	0.0657021647807872	0.0657021388668704
20	40	0.0343282999318164	0.0343282410061467
30	60	0.0231960539051621	0.0231959763238274
40	80	0.0174021796754980	0.0174020827329464
50	100	0.0145961553028378	0.0145960139057242

Table 2

Errors of the approximate solution (at the $m^{\text{th}}$ step) and the exact solution, with $N=10$ and $M=20$.

Number of iterations	Single-iterative scheme	2-order iterative scheme
$m$	$E_{M,N}={u^m-u_{ex}}_{M,N}$	$E_{M,N}={u^m-u_{ex}}_{M,N}$
1	0.9418510828896527	0.9418510828896527
2	0.1916822555160533	0.0846870773982589
3	0.0723290237392826	0.0657041906876774
4	0.0658909334218025	0.0657021388668704
5	0.0657050043944243	0.0657021388668677
6	0.0657021647807872	0.0657021388668675
7	0.0657021390166173	0.0657021388668675
8	0.0657021388674475	0.0657021388668675
9	0.0657021388668695	0.0657021388668675
10	0.0657021388668677	0.0657021388668675
11	0.0657021388668677	0.0657021388668675
12	0.0657021388668677	0.0657021388668675
13	0.0657021388668677	0.0657021388668675
14	0.0657021388668677	0.0657021388668675
15	0.0657021388668677	0.0657021388668675

Table 3

Errors of two consecutive iterations, with $N=10$ and $M=20$.

Number of iterations	Single-iterative scheme	2-order iterative scheme
$m$	$D_{M,N}^m={u^m-u^{m-1}}_{M,N}$	$D_{M,N}^m={u^m-u^{m-1}}_{M,N}$
1	2.5	2.5
2	0.7578162458413513	0.8858616655396145
3	0.1501823972253367	0.0356453241520862
4	0.0129970934628956	9.425500365978223e(−06)
5	5.099302560935826e(−04)	8.992806499463768e(−14)
6	1.061574497773776e(−05)	2.220446049250313e(−16)
7	1.316120679106803e(−07)	0
8	1.043858777194373e(−09)	0
9	5.554223747594733e(−12)	0
10	2.042810365310288e(−14)	0
11	0	0
12	0	0
13	0	0
14	0	0
15	0	0

Table 4

Errors of the approximate solution (at the $m^{\text{th}}$ step) and the exact solution, with $N=50$ and $M=100$.

Number of iterations	Single-iterative scheme	2-order iterative scheme
$m$	$E_{M,N}^m={u^m-u_{ex}}_{M,N}$	$E_{M,N}^m={u^m-u_{ex}}_{M,N}$
1	0.9473484068123711	0.9473484068123711
2	0.1901034101701440	0.0478701831481729
3	0.0258575358259656	0.0146003885629640
4	0.0149815977268044	0.0145960139057242
5	0.0146051199031552	0.0145960139056984
6	0.0145961553028378	0.0145960139056975
7	0.0145960154543019	0.0145960139056975
8	0.0145960139182475	0.0145960139056978
9	0.0145960139057733	0.0145960139056978
10	0.0145960139056933	0.0145960139056978
11	0.0145960139056931	0.0145960139056978
12	0.0145960139056931	0.0145960139056978
13	0.0145960139056931	0.0145960139056978
14	0.0145960139056931	0.0145960139056978
15	0.0145960139056931	0.0145960139056978

Table 5

Errors of two consecutive iterations, with $N=50$ and $M=100$.

Number of iterations	Single-iterative scheme	2-order iterative scheme
$m$	$D_{M,N}^m={u^m-u^{m-1}}_{M,N}$	$D_{M,N}^m={u^m-u^{m-1}}_{M,N}$
1	2.5	2.5
2	0.7572449966422271	0.9008606827792921
3	0.1675736824247880	0.0428472126493304
4	0.0179199191769293	1.897245848470064e(−05)
5	9.579076481733839e(−04)	7.116529587847253e(−13)
6	2.975046644859702e(−05)	1.088018564132653e(−14)
7	6.028582577588537e(−07)	3.996802888650564e(−15)
8	8.580162136340164e(−09)	8.881784197001252e(−16)
9	9.013212398656378e(−11)	2.220446049250313e(−16)
10	7.265299473147024e(−13)	0
11	4.884981308350689e(−15)	0
12	2.220446049250313e(−16)	0
13	0	0
14	0	0
15	0	0

6. 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 2-order 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 single-order iterative scheme.

Authors’ Contributions

All authors contributed equally to this article. They read and approved the final manuscript.

References

[1] S. Ahmad, K. Ijaz, A. Waleed, T. Khan, T. Hayat, A. Alsaedi, "Impact of arrhenius activation energy in viscoelastic nanomaterial flow subject to binary chemical reaction and non-linear mixed convection," Thermal Science, vol. 24 no. 2, pp. 1143-1155, DOI: 10.2298/tsci180524212a, 2020.

[2] N. El-Dabe, M. Gabr, A.-E. Elshekhipy, S. Zaher, "The motion of a non-Newtonian nanofluid over a semi-infinite moving vertical plate through porous medium with heat and mass transfer," Thermal Science, vol. 24 no. 2, pp. 1311-1321, DOI: 10.2298/tsci180604021e, 2020.

[3] N. Ijaz, M. Bhatti, A. Zeeshan, "Heat transfer analysis in magnetohydrodynamic flow of solid particles in non-Newtonian Ree-Eyring fluid due to peristaltic wave in a channel," Thermal Science, vol. 23 no. 2, pp. 1017-1026, DOI: 10.2298/tsci170220155i, 2019.

[4] M. M. Rashidi, M. Ali, B. Rostami, P. Rostami, G. Xie, "Heat and mass transfer for MHD viscoelastic fluid flow over a vertical stretching sheet with considering Soret and Dufour effects," Mathematical Problems in Engineering, vol. 2015, 2015.

[5] A. Zeeshan, M. Bhatti, N. Ijaz, O. Bég, A. Kadir, "Biologically inspired transport of solid spherical nanoparticles in an electrically-conducting viscoelastic fluid with heat transfer," Thermal Science, vol. 24 no. 2, pp. 1251-1260, DOI: 10.2298/tsci180706324z, 2020.

[6] Y. Zuo, H. Liu, "A fractal rheological model for SiC paste using a fractal derivative," Journal of Applied and Computational Mechanics, vol. 7 no. 1, pp. 13-18, 2021.

[7] S. A. Messaoudi, "Blow up and global existence in a nonlinear viscoelastic wave equation," Mathematische Nachrichten, vol. 260 no. 1, pp. 58-66, DOI: 10.1002/mana.200310104, 2003.

[8] Q. Li, L. He, "General decay and blow-up of solutions for a nonlinear viscoelastic wave equation with strong damping," Bound Value Problems, vol. 2018,DOI: 10.1186/s13661-018-1072-1, 2018.

[9] N. Mezouar, S. Boulaaras, "Global existence and decay of solutions for a class of viscoelastic Kirchhoff equation," Bulletin of the Malaysian Mathematical Sciences Society, vol. 43 no. 1, pp. 725-755, DOI: 10.1007/s40840-018-00708-2, 2020.

[10] N. Mezouar, S. Boulaaras, "Global existence and decay of solutions of a singular nonlocal viscoelastic system with damping terms," Topol. Methods Nonlinear Anal, vol. 56 no. 1, pp. 283-312, 2020.

[11] P. K. Parida, D. K. Gupta, "Recurrence relations for a Newton-like method in Banach spaces," Journal of Computational and Applied Mathematics, vol. 206 no. 2, pp. 873-887, DOI: 10.1016/j.cam.2006.08.027, 2007.

[12] L. T. P. Ngoc, N. T. Long, "Existence, blow-up, and exponential decay estimates for a system of nonlinear wave equations with nonlinear boundary conditions," Mathematical Methods in the Applied Sciences, vol. 37 no. 4, pp. 464-487, DOI: 10.1002/mma.2803, 2014.

[13] L. T. P. Ngoc, B. M. Tri, N. T. Long, "An n-order iterative scheme for a nonlinear wave equation containing a nonlocal term," Filomat, vol. 32 no. 6, pp. 1755-1767, DOI: 10.2298/fil1706755n, 2017.

[14] L. T. P. Ngoc, N. A. Triet, A. P. Ngoc Dinh, N. T. Long, "Existence and exponential decay of solutions for a wave equation with integral nonlocal boundary conditions of memory type," Numerical Functional Analysis and Optimization, vol. 38 no. 9, pp. 1173-1207, DOI: 10.1080/01630563.2017.1320672, 2017.

[15] N. H. Nhan, L. T. P. Ngoc, T. M. Thuyet, N. T. Long, "On a high order iterative scheme for a nonlinear wave equation with the source term containing a nonlinear integral," Nonlinear Functional Analysis and Applications, vol. 21 no. 1, pp. 65-84, 2016.

[16] L. X. Truong, L. T. P. Ngoc, N. T. Long, "High-order iterative schemes for a nonlinear Kirchhoff-Carrier wave equation associated with the mixed homogeneous conditions," Nonlinear Analysis TMA, vol. 71 no. 1-2, pp. 467-484, DOI: 10.1016/j.na.2008.10.086, 2009.

[17] L. X. Truong, L. T. P. Ngoc, N. T. Long, "The n-order iterative schemes for a nonlinear Kirchhoff-Carrier wave equation associated with the mixed inhomogeneous conditions," Applied Mathematics and Computation, vol. 215 no. 5, pp. 1908-1925, DOI: 10.1016/j.amc.2009.07.056, 2009.

[18] X. X. Li, C. H. He, "Homotopy perturbation method coupled with the enhanced perturbation method," Journal of Low Frequation Noise V. A, vol. 38 no. 3-4, pp. 1399-1403, DOI: 10.1177/1461348418800554, 2019.

[19] Q.-P. Ji, J. Wang, L.-X. Lu, C.-F. Ge, "Li-He’s modified homotopy perturbation method coupled with the energy method for the dropping shock response of a tangent nonlinear packaging system," Journal of Low Frequency Noise, Vibration and Active Control, vol. 146,DOI: 10.1177/1461348420914457, 2020.

[20] H. A. Levine, "Instability and nonexistence of global solutions to nonlinear wave equations of the form $Pu\sb{tt} = −Au + \mathcal{F}(u)$ P u t t = A u + F u," Transactions of the American Mathematical Society, vol. 192,DOI: 10.1090/s0002-9947-1974-0344697-2, 1974.

[21] H. A. Levine, "Some additional remarks on the nonexistence of global solutions to nonlinear wave equations," SIAM Journal on Mathematical Analysis, vol. 5 no. 1, pp. 138-146, DOI: 10.1137/0505015, 1974.

[22] L. T. P. Ngoc, L. H. K. Son, T. M. Thuyet, N. T. Long, "An N -order iterative scheme for a nonlinear Carrier wave equation in the annular with Robin-Dirichlet conditions," Nonlinear Functional Analysis and Applications, vol. 22 no. 1, pp. 147-169, 2017.

[23] L. X. Truong, L. T. P. Ngoc, A. P. N. Dinh, N. T. Long, "Existence, blow-up and exponential decay estimates for a nonlinear wave equation with boundary conditions of two-point type," Nonlinear Analysis: Theory, Methods & Applications, vol. 74 no. 18, pp. 6933-6949, DOI: 10.1016/j.na.2011.07.015, 2011.

[24] V. Georgiev, G. Todorova, "Existence of a solution of the wave equation with nonlinear damping and source terms," Journal of Differential Equations, vol. 109 no. 2, pp. 295-308, DOI: 10.1006/jdeq.1994.1051, 1994.

[25] O. M. Jokhadze, "Global Cauchy problem for wave equations with a nonlinear damping term," Differential Equations, vol. 50 no. 1, pp. 57-65, DOI: 10.1134/s001226611401008x, 2014.

[26] N. T. Long, L. T. P. Ngoc, "On a nonlinear wave equation with boundary conditions of two-point type," Journal of Mathematical Analysis and Applications, vol. 385 no. 2, pp. 1070-1093, DOI: 10.1016/j.jmaa.2011.07.034, 2012.

[27] L. T. P. Ngoc, L. X. Truong, N. T. Long, "An N -order iterative scheme for a nonlinear Kirchhoff-Carrier wave equation associated with mixed homogeneous conditions," Acta Mathematica Vietnamica, vol. 35 no. 2, pp. 207-227, 2010.

[28] G. Todorova, "Cauchy problem for a nonlinear wave equation with nonlinear damping and source terms," Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, vol. 326 no. 2, pp. 191-196, DOI: 10.1016/s0764-4442(97)89469-4, 1998.

[29] N. A. Triet, L. T. P. Ngoc, A. P. N. Dinh, N. T. Long, "Existence and exponential decay for a nonlinear wave equation with nonlocal boundary conditions of 2N-point type," Mathematical Methods in the Applied Sciences, vol. 44 no. 1, pp. 668-692, 2021.

[30] M. Kafini, S. A. Messaoudi, "A blow-up result in a Cauchy viscoelastic problem," Applied Mathematics Letters, vol. 21 no. 6, pp. 549-553, DOI: 10.1016/j.aml.2007.07.004, 2008.

[31] N. T. Long, A. P. N. Dinh, L. X. Truong, "Existence and decay of solutions of a nonlinear viscoelastic problem with a mixed nonhomogeneous condition," Numerical Functional Analysis and Optimization, vol. 29 no. 11-12, pp. 1363-1393, DOI: 10.1080/01630560802605955, 2008.

[32] S. Boulaaras, "Some new properties of asynchronous algorithms of theta scheme combined with finite elements methods for an evolutionary implicit 2-sided obstacle problem," Mathematical Methods in the Applied Sciences, vol. 40 no. 18, pp. 7231-7239, DOI: 10.1002/mma.4525, 2017.

[33] S. Boulaaras, M. Haiour, M. A. Bencheick Le hocine, "L∞-error estimates of discontinuous Galerkin methods with theta time discretization scheme for an evolutionary HJB equations," Mathematical Methods in the Applied Sciences, vol. 40 no. 12, pp. 4310-4319, DOI: 10.1002/mma.4306, 2017.

[34] R. K. Mohanty, V. Gopal, "High accuracy cubic spline finite difference approximation for the solution of one-space dimensional non-linear wave equations," Applied Mathematics and Computation, vol. 218 no. 8, pp. 4234-4244, DOI: 10.1016/j.amc.2011.09.054, 2011.

[35] M. M. Alsuyuti, E. H. Doha, S. S. Ezz‐Eldien, B. I. Bayoumi, D. Baleanu, "Modified Galerkin algorithm for solving multitype fractional differential equations," Mathematical Methods in the Applied Sciences, vol. 42 no. 5, pp. 1389-1412, DOI: 10.1002/mma.5431, 2019.

[36] Ö. Oruç, "Two meshless methods based on local radial basis function and barycentric rational interpolation for solving 2D viscoelastic wave equation," Computers & Mathematics with Applications, vol. 79 no. 12, pp. 3272-3288, DOI: 10.1016/j.camwa.2020.01.025, 2020.

[37] V. Singh, S. Islam, R. K. Mohanty, "Local meshless method for convection dominated steady and unsteady partial differential equations," Engineering with Computers, vol. 35 no. 6, pp. 803-812, DOI: 10.1007/s00366-018-0632-4, 2019.

[38] G. F. Pinder, Numerical Methods for Solving Partial Differential Equations, 2018.

[39] R. E. Showalter, "Hilbert space methods for partial differential equations," Electronic Journal of Differential Equations Monograph, vol. 1, 1994.

[40] V. Lakshmikantham, S. Leela, Differential and Integral Inequalities, 1969.

[41] J. L. Lions, Quelques méthodes de résolution des problèmes aux limites nonlinéaires, Dunod, 1969.

[42] N. T. Long, T. N. Diem, "On the nonlinear wave equation u t t − u x x = f x , t , u x , u t , associated with the mixed homogeneous conditions," Nonlinear Anal. TMA, vol. 29 no. 11, pp. 1217-1230, 1997.