1 Introduction↩︎

1.1 Formulation of the Problem↩︎

Motivated by the theory of beams, plates, and shells, this work presents a natural abstract generalization of the nonlinear dynamic Timoshenko beam system. As a second-generation beam theory, the Timoshenko model [1] extends beyond the classical Euler-Bernoulli (Kirchhoff-Love) theory by accounting for shear deformation and rotational bending. This approach is particularly suitable for analyzing thick beams, sandwich composite beams, and beams subjected to high-frequency excitation where the wavelength is comparable to the beam thickness. We paid attention to this model, since it has wide-ranging applications in structural and mechanical engineering, stemming from the nonlinear theory of elasticity.

Let us consider the Cauchy problem associated with a nonlinear abstract system of coupled hyperbolic equations defined in the real Hilbert space \(\mathcal{H}\): \[\begin{gather} \frac{\mathrm{d}^2 u\left( t \right)}{\mathrm{d}t^2} + \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u \right\rVert}^2 \right) A u\left( t \right) + a_1 B v\left( t \right) = f_1\left( t \right)\,,\quad t \in \left[ 0,T \right]\,,\tag{1}\\ \frac{\mathrm{d}^2 v\left( t \right)}{\mathrm{d}t^2} + \gamma A v\left( t \right) + \delta C v\left( t \right) + a_2 B u\left( t \right) = f_2\left( t \right)\,,\tag{2}\\ u\left( 0 \right) = \varphi_0\,,\quad u^{\prime}\left( 0 \right) = \varphi_1\,,\quad v\left( 0 \right) = \psi_0\,,\quad v^{\prime}\left( 0 \right) = \psi_1\,.\tag{3} \end{gather}\]

Let \(\alpha\), \(\beta\), \(\gamma\), and \(\delta\) be positive constants, and let \(a_1\) and \(a_2\) be constants unrestricted in sign. Consider \(A\) to be a self-adjoint, positive-definite operator (independent of \(t\) and generally unbounded) with domain \(D\left( A \right)\), which is everywhere dense in \(\mathcal{H}\), i.e. \(\overline{D\left( A \right)} = \mathcal{H}\), \(A = A^{\ast} \geq \nu I\) with \(\nu > 0\) (where \(I\) denotes the identity operator). Let \(C\) be a symmetric, bounded operator such that \(\left( C \varphi,\varphi \right) \geq 0\) for all \(\varphi \in \mathcal{H}\). The operator \(B\) is assumed to be a closed linear operator satisfying the following condition: \[\label{eq:B95cond} {\left\lVert B \varphi \right\rVert}^2 \leq \cstBoper^2 \left( A \varphi,\varphi \right)\,,\quad \forall \varphi \in D\left( A \right) \subset D\left( B \right)\,,\quad \cstBoper > 0\,,\tag{4}\] here, \(\left( \cdot,\cdot \right)\) denotes the inner product defined on the Hilbert space \(\mathcal{H}\), and \(\left\lVert \cdot \right\rVert\) denotes the norm associated with this inner product. Moreover, let \(\varphi_0\), \(\varphi_1\), \(\psi_0\), and \(\psi_1\) be specified vectors in the Hilbert space \(\mathcal{H}\). The unknown functions \(u\left( t \right)\) and \(v\left( t \right)\) are continuous and twice continuously differentiable, mapping into \(\mathcal{H}\). Additionally, \(f_1\left( t \right)\) and \(f_2\left( t \right)\) are specified continuous functions taking values in \(\mathcal{H}\).

We now comment on the appearance of the square root in equation 1 .

In accordance with the linear case (cf. Kreı̆n Chapter III

[3] Chapter III ), the vector functions \(u\left( t \right)\) and \(v\left( t \right)\), taking values in \(\mathcal{H}\) and defined on the interval \(\left[ 0,T \right]\), are called solutions to the problem 1 3 if they fulfill the following conditions:

In this context, continuity and differentiability are considered in relation to a metric defined on the Hilbert space \(\mathcal{H}\).

Moving from theory to application, we focus on a specific case of the problem 1 3 , corresponding to the spatially one-dimensional initial-boundary value problem 65 68 . This model describes the geometrically nonlinear vibrations of a beam and is formulated as a coupled system of nonlinear integro-partial differential equations with homogeneous Dirichlet boundary conditions. The homogeneous version of this nonlinear system was first proposed by Sapir and Reiss in the appendix of [4] as a mathematical model for the free vibrations of an elastic beam with fixed endpoints. The unknown functions \(u\) and \(v\) represent the transverse displacement of the beam centerline and the rotational displacement of the cross section, respectively. The constants in the model characterize the physical and mechanical properties of the beam and are determined by its material parameters (the significance of these constants is discussed in [4] and [5]).

The same mechanical problem has been analyzed within the framework of the nonlinear Euler-Bernoulli beam theory (cf. the articles by Dickey [6], Nayfeh and Mook [7], and Woinowsky-Krieger [8]). However, unlike the Timoshenko model, this theory is limited because it does not account for shear deformation and rotary inertia effects. For the resulting initial and boundary value problems, global existence and uniqueness results in appropriate Sobolev spaces were established by Ball [9] and Dickey [6] using the Galerkin method. The equations examined in these studies can be formulated as a semilinear evolution equation in a Hilbert space. This formulation enables the application of semigroup theory to reduce the problem to an integral equation (cf. Ball [10] and Fitzgibbon [11]). Although this approach is not applicable to the study of 1 3 because the corresponding evolution equation is nonlinear. Issues of existence and uniqueness of solutions to the abstract nonlinear wave equation have been addressed, for instance, in papers of Hochstein-Kind [12] and Pohožaev [13].

As far as we know, the issue of solvability in the abstract setting of the dynamic nonlinear Timoshenko beam problem 1 3 in a Hilbert space has not been studied. Regarding the Timoshenko beam system 65 68 , Tucsnak first raised the question of existence and uniqueness of solutions in [5] and proved local-in-time solvability under suitable regularity assumptions on the initial data. Subsequently, Ammari [14] investigated the global existence and large-time behavior of the system governing the nonlinear vibrations of a Timoshenko beam and, under small initial data, established the global existence of strong solutions together with exponential decay of the associated energy. More recently, Narciso and Cousin [15] obtained existence and uniqueness results for local solutions in noncylindrical domains by employing the Faedo–Galerkin method. Furthermore, in [16], Aouragh, Segaoui, and Soufyane analyzed the stability of a nonlinear shear beam system and established the well-posedness of the considered model through the Faedo–Galerkin method. In the subsequent work [17], Aouragh, El Baz, and Soufyane investigated a thermoelastic nonlinear shear beam model with thermal dissipation. Using the Faedo–Galerkin method, they established the well-posedness of the system and proved its exponential stability via the multiplier method. For numerical purposes, they employed a finite element discretization in space, combined with the Euler and Crank–Nicolson schemes for temporal discretization.

We would like to state that we are not aware of the works that directly address the design of numerical algorithms and the construction of approximate solutions for the abstract analogue of the Timoshenko system 1 3 . Concerning the construction and analysis of numerical schemes for the concrete Timoshenko system 65 68 , several articles are devoted to these issues. In particular, Bernardi and Copetti [18], [19] proposed a numerical approach for a contact problem arising in the nonlinear dynamic thermoviscoelastic Timoshenko beam model. After establishing the well-posedness of the corresponding system of three equations, they applied finite element discretization in space combined with Euler and Crank–Nicolson time-stepping schemes, derived a priori error estimates for the discrete problem, and presented results from numerical experiments. Peradze [20] investigated the homogeneous version of problem 65 68 under the assumption that the associated Cauchy data are analytic. The numerical solution was obtained by combining the Galerkin method with a Crank–Nicolson-type difference scheme and a Picard iteration procedure, and convergence and accuracy properties of the resulting algorithm were analyzed. In a subsequent work [21], Peradze considers problem 65 68 in a homogeneous setting and proposes a numerical solution based on the finite element method, a modified Crank–Nicolson-type difference scheme, and a Picard-type iteration process. The total error of the proposed method is estimated. In the paper by Hauck, Målqvist, and Rupp [22], two-level domain decomposition methods for spatial network models are developed, and an application of this method to the Timoshenko beam network is presented.

In our opinion, viewing the system 1 2 as an abstract analogue of the nonlinear Timoshenko model has certain advantages. In this context, we can employ a unified approach for a class of related problems. In particular, the formulation 1 3 includes, within the framework of the Timoshenko model, the Dirichlet, Neumann, Robin, and mixed boundary value problems. Furthermore, by using the abstract setting, we can also address the spatially multidimensional case on a Lipschitz domain.

The approach used in the following discussion also covers the case where the coupling coefficients \(a_1\) and \(a_2\) vanish. When \(a_1 = a_2 = 0\), system 1 2 decouples into two independent equations. The first equation is an abstract counterpart of the classical Kirchhoff equation and includes both one-dimensional and multidimensional settings; see Lions [23]. Moreover, our approach can be extended to the case where, in equation 1 , the term \({\left\lVert A^{\frac{1}{2}} u \right\rVert}^2\) is replaced by \({\left\lVert u \right\rVert}^2\). This formulation includes, as a particular case, the Carrier equation [24].

The second equation provides an abstract formulation of the linear model of wave propagation. Clearly, it covers both one-dimensional and multidimensional cases. In addition, this equation also includes dynamic plane and spatial problems arising in the theory of elasticity. In this setting, \(A\) represents the Friedrichs extension of the symmetric positive definite elasticity operator associated with the displacement formulation of the elasticity system under classical boundary conditions; see Mikhlin Section 29

[25, pp. Section 29] .

1.2 Outline of the Problem↩︎

We develop a time-stepping scheme for the nonlinear problem 1 3 . To achieve this objective, a uniform temporal grid is introduced, and a symmetric three-layer semi-discrete scheme is proposed, in which the nonlinear term is evaluated at a temporal midpoint. Such a treatment permits an approximate solution of the original nonlinear problem by solving a linear problem in parallel at each temporal layer.

The investigation of the stability and convergence of the proposed scheme is based on the following fact: the sequences of vectors \(A^{\frac{s}{2}} \left( u_k - u_{k - 1} \right) / \tau\) and \(A^{\frac{\left( s + 1 \right)}{2}} u_k\), with \(s = 0, 1\), as well as \(\left( v_k - v_{k - 1} \right) / \tau\) and \(L^{\frac{1}{2}} v_k\), where \(L = \delta A + \gamma C\), are uniformly bounded. Here, \(u_k\) and \(v_k\) denote approximate solutions, and \(\tau > 0\) is the time step. These facts allow us to estimate the error of the approximate solution and to show that the order of convergence is \(\bigO \left( \tau^2 \right)\) in the class of smooth solutions. Moreover, the approximation error of the first derivative, obtained by applying the central finite difference formula to the approximate solution, is also of order \(\bigO \left( \tau^2 \right)\).

We subsequently consider a nonlinear, spatially one-dimensional Timoshenko beam system 65 68 as a special case of the problem 1 3 , and the results obtained within the abstract framework are extended to this setting. In this case, the proposed scheme yields second-order linear ordinary differential equations for the unknowns \(u_k \left( x \right)\) and \(v_k \left( x \right)\) at each temporal layer. Note that these differential equations are not coupled and can therefore be solved in parallel. To compute numerical solutions of these equations, we employ the Legendre–Galerkin spectral method. The chosen basis functions ensure that the resulting Galerkin linear system is sparse, which can be decoupled into two linear subsystems. More precisely, the coefficient matrix of the corresponding Galerkin linear system is a symmetric positive-definite tridiagonal matrix with a gap, in the sense that it has only nonzero entries on the main diagonal and the second sub- and superdiagonals, while the first sub- and superdiagonals vanish. By exploiting this structure, the associated linear system can be decomposed into two independent tridiagonal subsystems: one corresponding to the odd-indexed unknowns and the other to the even-indexed unknowns. This property is particularly important for numerical implementation.

An important step in the numerical computation of the Timoshenko model using the proposed scheme is the estimation of the error of the Legendre–Galerkin spectral method. The approximation error of this method is estimated using both the \(L^2\)-norm and the uniform norm for the ordinary differential equations resulting from time discretization.

The crowning stage of each algorithm’s development is numerical experiments. It is well known that the construction of algorithms for the numerical computation of mathematical models describing various processes, together with their computer implementation, enables the use of numerical experiments that replace many expensive real-world experiments spanning various fields of natural science. To assess the practical value of an algorithm, it is important to conduct numerical experiments on model problems that fully demonstrate the role of each component of the algorithm.

The paper discusses four benchmark problems whose solutions exhibit oscillatory behavior. Numerical results obtained using the proposed algorithm indicate that the combined algorithm is stable and achieves high practical accuracy for these problems. The proposed algorithm is implemented in Python, and the corresponding source code is archived in an open-access Zenodo repository [26].

2 Symmetric Three-Layer Semi-Discrete Scheme↩︎

We aim to find a solution to the problem 1 3 using the following semi-discrete scheme: \[\begin{gather} \frac{\Delta^2 u_{k - 1}}{\tau^2} + \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^2 \right) \frac{A u_{k + 1} + A u_{k - 1}}{2} = f_{1,k} - a_1 B v_k\,,\tag{5} \\ \frac{\Delta^2 v_{k - 1}}{\tau^2} + \frac{L v_{k + 1} + L v_{k - 1}}{2} = f_{2,k} - a_2 B u_k\,,\tag{6} \end{gather}\] where \(k = 1,2,\ldots,n - 1\) and \(\tau = T/n\) (with \(n > 1\)), \(\Delta u_{k - 1} = u_k - u_{k - 1}\) and \(\Delta^2 u_{k - 1} = \Delta\left( \Delta u_{k - 1} \right)\). Additionally, let \(f_{1,k} = f_1\left( t_k \right)\) and \(f_{2,k} = f_2\left( t_k \right)\), where \(t_k = k\tau\). The initial conditions are given by \(u_0 = \varphi_0\) and \(v_0 = \psi_0\). The operator \(L\) is defined as \(L = \gamma A + \delta C\).

Let \(u\left( t_k \right)\) and \(v\left( t_k \right)\) be the values of the solutions to problem 1 3 at the discrete time points \(t = t_k\). Their corresponding numerical approximations at these points are denoted by \(u_k\) and \(v_k\), respectively; thus, \(u\left( t_k \right) \approx u_k\) and \(v\left( t_k \right) \approx v_k\).

To perform computations using schemes 5 and 6 , it is essential to ascertain the initial vectors \(u_0\), \(v_0\), \(u_1\), and \(v_1\). The vectors \(u_0\) and \(v_0\) are prescribed by \(u_0 = \varphi_0\) and \(v_0 = \psi_0\), respectively. The vectors \(u_1\) and \(v_1\) need to be approximated. It is well-established that to approximate \(u_1\) and \(v_1\), one should expand the exact solutions \(u\left( t \right)\) and \(v\left( t \right)\) in a Taylor series around \(t = 0\), retaining at least the first two terms. This yields: \(u_1 = \varphi_0 + \tau \varphi_1\) and \(v_1 = \psi_0 + \tau \psi_1\). To achieve second-order accuracy, it is necessary to include the first three terms in the Taylor expansions of \(u\left( \tau \right)\) and \(v\left( \tau \right)\), which involves the second-order derivatives. These second-order derivatives of \(u\left( t \right)\) and \(v\left( t \right)\) at \(t = 0\) can be determined from equations 1 and 2 , considering the initial conditions specified in 3 . Thus, we obtain: \[\begin{gather} u_1 = \varphi_0 + \tau \varphi_1 + \frac{\tau^2}{2} \varphi_2\,,\quad \varphi_2 = f_{1,0} - a_1 B \psi_0 - \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} \varphi_0 \right\rVert}^2 \right) A \varphi_0\,,\tag{7} \\ v_1 = \psi_0 + \tau \psi_1 + \frac{\tau^2}{2} \psi_2\,,\quad \psi_2 = f_{2,0} - a_2 B \varphi_0 - L \psi_0\,.\tag{8} \end{gather}\] By substituting the values of the starting vectors \(u_0\), \(v_0\), \(u_1\), and \(v_1\) into equations 5 and 6 , we derive the subsequent linear system of equations that defines the vectors \(u_2\) and \(v_2\): \[\begin{align} \left( I + \frac{\tau^2}{2} q_1 A \right) u_2 &= g_1\,,\tag{9} \\ \left( I + \frac{\tau^2}{2} L \right) v_2 &= g_2\,,\tag{10} \end{align}\] where the scalar \(q_1 = \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_1 \right\rVert}^2\) and the right-hand sides \(g_1\) and \(g_2\) are given.

Given that \(A\) is a self-adjoint, positively-definite operator and \(q_1 > 0\), it follows that the operator \(I + 0.5 \tau^2 q_1 A\) is also self-adjoint and positively definite. Consequently, the operator \(I + 0.5 \tau^2 q_1 A\) is continuously invertible, ensuring that equation 9 has a unique solution for each \(g_1\) from \(\mathcal{H}\), which depends continuously on the right-hand side. The same reasoning applies to equation 10 . Moreover, we consider the following fact: the sum of a self-adjoint operator and a symmetric, bounded operator is also self-adjoint (cf. Kato Chapter IV

[2] Chapter IV ). Therefore, \(u_2 = {\left( I + 0.5 \tau^2 q_1 A \right)}^{-1} g_1\) and \(v_2 = {\left( I + 0.5 \tau^2 L \right)}^{-1} g_2\), where \({\left( I + 0.5 \tau^2 q_1 A \right)}^{-1}\) and \({\left( I + 0.5 \tau^2 L \right)}^{-1}\) are bounded, self-adjoint operators defined over the entire Hilbert space \(\mathcal{H}\). In a way analogous to the determination of \(u_2\) and \(v_2\), the vectors \(u_k\) and \(v_k\) for \(k > 2\) are computed using \(u_{k - 1}\), \(v_{k - 1}\), \(u_{k - 2}\), and \(v_{k - 2}\). Consequently, the application of schemes 5 and 6 is reduced to solving linear problems at each temporal layer. The proposed algorithm enables the parallel computation of the vectors \(u_k\) and \(v_k\).

3 Preliminaries↩︎

In the following, we present two lemmas, stated without proof, which we require to establish the convergence of the scheme 5 6 .

The subsequent lemma is stated with an adjustment to the proposed scheme 6 . The lemma is formulated as follows.

For our purposes, in the estimate of ?? , we need to replace the operator \(L\) with the operator \(A\). This requires dual assessments of the operator \(L\) in terms of the operator \(A\). We shall state these estimates in the form of remarks.

Let us now demonstrate the left-hand inequality of ?? . By applying a straightforward transformation, it is obtained \[\label{eq:remark195lower95bound} {\left\lVert L\varphi \right\rVert}^{2} = \gamma^2 {\left\lVert A\varphi \right\rVert}^2 + 2 \gamma \delta \left( A\varphi,C\varphi \right) + \delta^2 {\left\lVert C\varphi \right\rVert}^2 \geq \gamma^2 {\left\lVert A\varphi \right\rVert}^2 - 2 \gamma \delta \left\lVert A\varphi \right\rVert\left\lVert C\varphi \right\rVert + \delta^2 {\left\lVert C\varphi \right\rVert}^2\,.\tag{11}\] In the derivation, we have applied the Cauchy-Schwarz inequality.

Applying Young’s inequality with \(\varepsilon\) (valid for every \(\varepsilon > 0\)) to the term \(2\left\lVert A\varphi \right\rVert\left\lVert C\varphi \right\rVert\) yields the following result for assessment 11 \[\label{eq:remark195Lu95ineqt} {\left\lVert L\varphi \right\rVert}^{2} + \frac{\delta}{\varepsilon} \left( \gamma - \varepsilon \delta \right) {\left\lVert C\varphi \right\rVert}^2 \geq \gamma \left( \gamma - \varepsilon \delta \right) {\left\lVert A\varphi \right\rVert}^2\,,\quad 0 < \varepsilon < \frac{\gamma}{\delta}\,.\tag{12}\] Given the positive definiteness of the operator \(A\) and using the Cauchy-Schwarz inequality, we obtain the following result: \[\left\lVert L\varphi \right\rVert \left\lVert \varphi \right\rVert \geq \left( L\varphi,\varphi \right) \geq \gamma \left( A\varphi,\varphi \right) \geq \gamma \nu {\left\lVert \varphi \right\rVert}^2\,,\] which leads to the conclusion that \[\left\lVert \varphi \right\rVert \leq \frac{1}{\gamma \nu} \left\lVert L\varphi \right\rVert\,.\] Furthermore, under the assumption that \(\left\lVert C\varphi \right\rVert \leq c \left\lVert \varphi \right\rVert\), with \(c = \left\lVert C \right\rVert\), we obtain the resulting inequality from 12 \[\left\lVert L\varphi \right\rVert \geq \hat{\gamma} \left\lVert A\varphi \right\rVert\,,\quad \hat{\gamma} = \sqrt{\frac{\gamma \left( \gamma - \varepsilon \delta \right)}{1 + \frac{c^2 \delta}{\varepsilon \gamma^2 \nu^2}\left( \gamma - \varepsilon \delta \right)}}\,.\]

Proceeding to demonstrate the right-hand inequality of ?? , we account for the fact that \(\left\lVert C\varphi \right\rVert \leq c \left\lVert \varphi \right\rVert\) and thus arrive at the following conclusion: \[\label{eq:remark195norm95L} \left\lVert L\varphi \right\rVert \leq \gamma \left\lVert A\varphi \right\rVert + c \delta \left\lVert \varphi \right\rVert\,.\tag{13}\] By taking an additional step, which involves applying the Cauchy-Schwarz inequality and employing the positive definiteness of the operator \(A\), we obtain: \[\left\lVert A\varphi \right\rVert\left\lVert \varphi \right\rVert \geq \left( A\varphi,\varphi \right) \geq \nu {\left\lVert \varphi \right\rVert}^2\,.\] Based on the aforementioned inequality in 13 , the desired result follows \[\left\lVert L\varphi \right\rVert \leq \nu_0 \left\lVert A\varphi \right\rVert\,,\quad \nu_0 = \gamma + \frac{c\delta}{\nu}\,.\]

Due to the positive definiteness of the operator \(A\), it follows directly that \[\label{eq:remark295result95post95def} {\left\lVert A^{\frac{1}{2}} \varphi \right\rVert}^2 \geq \nu {\left\lVert \varphi \right\rVert}^2\,,\quad \varphi \in D\left( A \right)\,.\tag{14}\] On the other hand, considering inequality 14 along with the fact that \(\left\lVert C\varphi \right\rVert \leq c \left\lVert \varphi \right\rVert\), we obtain \[\begin{align} \left( L\varphi,\varphi \right) &= \gamma \left( A\varphi,\varphi \right) + \delta\left( C\varphi,\varphi \right) \leq \gamma{\left\lVert A^{\frac{1}{2}} \varphi \right\rVert}^2 + c\delta{\left\lVert \varphi \right\rVert}^2 \\ &\leq \nu_0 {\left\lVert A^{\frac{1}{2}} \varphi \right\rVert}^2\,,\quad \varphi \in D\left( A \right)\,,\quad \nu_0 = \gamma + \frac{c\delta}{\nu}\,. \end{align}\] Consequently, we derive that \[\label{eq:remark295final95res} \left\lVert L^{\frac{1}{2}}\varphi \right\rVert \leq \nu_0 \left\lVert A^{\frac{1}{2}} \varphi \right\rVert\,,\quad \varphi \in D\left( A \right)\,.\tag{15}\] Given that the operator \(A^{\frac{1}{2}}\) maps \(D\left( A \right)\) onto \(D\left( A^{\frac{1}{2}} \right)\) and that the range of \(A^{\frac{1}{2}}\), \(R\left( A^{\frac{1}{2}} \right) = \mathcal{H}\). Under these conditions, it follows that inequality 15 can be extended to the entire domain \(D\left( A^{\frac{1}{2}} \right)\). Consequently, the right-hand inequality of ?? is satisfied.

The left-hand inequality of ?? is derived through analogous reasoning.

It is evident that from 4 , the following implication arises: \[\label{eq:remark695relat95norm95B95norm95halfA} \left\lVert B \varphi \right\rVert \leq \cstBoper \left\lVert A^{\frac{1}{2}} \varphi \right\rVert\,,\quad \forall \varphi \in D\left( A \right) \subset D\left( B \right)\,.\tag{16}\] It is known that \(D\left( A \right)\) is a core of \(A^{\frac{1}{2}}\) (see, Kato Lemma 3.38

[2, pp. Lemma 3.38] ). This implies that for every \(\varphi \in D\left( A^{\frac{1}{2}} \right)\), there exists a sequence \(\varphi_n \in D\left( A \right)\) such that \(\varphi_n \to \varphi\) and \(A^{\frac{1}{2}} \varphi_n \to A^{\frac{1}{2}} \varphi\). Consequently, by 16 , it follows that \(B \varphi_n\) forms a Cauchy sequence. Given that \(\mathcal{H}\) is a complete space, it is evident that this sequence converges. Furthermore, since the operator \(B\) is closed, we conclude that \(\varphi \in D\left( B \right)\) and \(B \varphi_n \to B \varphi\). As a result of this reasoning and inequality 16 , the relation ?? is established.

In the continuous problem, the coupling of equations 1 and 2 gives rise to the terms \(a_1 B v\left( t \right)\) and \(a_2 B u\left( t \right)\). Similarly, in the discrete problem, the terms \(a_1 B v_k\) and \(a_2 B u_k\) are responsible for the coupling of schemes 5 and 6 . Hence, to naturally connect the estimates of schemes 5 and 6 , it is necessary to fulfil a certain condition for the vector \(B \varphi\). The following remark addresses this matter.

From the second inequality in Remark [4](#prop:remark3){reference-type=“ref” reference=“prop:remark3”}, it follows immediately that \[\begin{align} \left\lVert A v_{k + 1} \right\rVert &\leq \frac{\nu_0}{\hat{\gamma}}\left[ \sqrt{2} \left( \left\lVert A v_0 \right\rVert + \left\lVert A^{\frac{1}{2}} \frac{\Delta v_0}{\tau} \right\rVert \right) + \tau\left\lVert A\frac{\Delta v_0}{\tau} \right\rVert\right. \\ &+ \left.\tau\sum_{i = 1}^{k} \left\lVert A^{\frac{1}{2}}f_{2,i} \right\rVert + \left\lvert a_2 \right\rvert\tau\sum_{i = 1}^{k} \left\lVert A^{\frac{1}{2}}B u_i \right\rVert \right]\,. \end{align}\] By considering inequality ?? from Remark [6](#prop:remark4){reference-type=“ref” reference=“prop:remark4”}, we obtain the desired estimate.

4 Uniform Boundedness of the Solution to a Discrete Problem and Its Corresponding Difference Analogues of the Derivative↩︎

It should be noted that, throughout this text and in all sections, the letters \(c\) and \(M\), indexed with lower subscripts, denote positive constants.

In this section, we demonstrate that the vectors \(\Delta u_{k - 1} / \tau\), \(\Delta v_{k - 1} / \tau\), \(A^{\frac{1}{2}} u_k\), and \(L^{\frac{1}{2}} v_k\) are uniformly bounded, which is essential for proving the convergence of the approximate solution. Additionally, establishing convergence requires us to prove the uniform boundedness of the vectors \(A u_k\) and \(A^{\frac{1}{2}} \Delta u_{k - 1} / {\tau}\). It is important to note that this issue is nontrivial, as the energy method does not yield a recurrence inequality that would allow the application of the telescoping series cancellation technique. The subsequent section addresses this by demonstrating that the vectors \(A u_k\) and \(A^{\frac{1}{2}} \Delta u_{k - 1} / {\tau}\) are locally uniformly bounded.

By evaluating the inner product of both sides of equation 5 with \(u_{k + 1} - u_{k - 1} = \Delta u_k + \Delta u_{k - 1}\) and considering the properties of the operator \(A\), which is self-adjoint and positive-definite, we arrive at \[\begin{align} \label{eq:lemma195main95equality} {\left\lVert \frac{\Delta u_k}{\tau} \right\rVert}^{2} + \frac{1}{2} \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2} \right) {\left\lVert A^{\frac{1}{2}} u_{k + 1} \right\rVert}^{2} &= {\left\lVert \frac{\Delta u_{k - 1}}{\tau} \right\rVert}^{2} + \frac{1}{2} \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2} \right) {\left\lVert A^{\frac{1}{2}} u_{k - 1} \right\rVert}^{2}\nonumber \\ &+ \left( f_{1,k} - a_1 B v_k,\Delta u_k \right) + \left( f_{1,k} - a_1 B v_k,\Delta u_{k - 1} \right)\,, \end{align}\tag{17}\] Let us denote \[\alpha_{1,k} = {\left\lVert \frac{\Delta u_{k - 1}}{\tau} \right\rVert}^{2}\,,\quad \gamma_{1,k} = {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2}\,.\] Employing these notations, equality 17 can be expressed as follows: \[\begin{align} \alpha_{1,k + 1} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \gamma_{1,k + 1} &= \alpha_{1,k} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \gamma_{1,k - 1} \\ &+ \left( f_{1,k} - a_1 B v_k,\Delta u_k \right) + \left( f_{1,k} - a_1 B v_k,\Delta u_{k - 1} \right)\,. \end{align}\] By applying the Cauchy-Schwarz inequality to the right-hand side of the given equality, one can conclude that \[\begin{align} \alpha_{1,k + 1} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \gamma_{1,k + 1} &\leq \alpha_{1,k} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \gamma_{1,k - 1} \\ &+ \tau \left( \sqrt{\alpha_{1,k}} + \sqrt{\alpha_{1,k + 1}} \right) \left\lVert f_{1,k} - a_1 B v_k \right\rVert\,. \end{align}\] Consequently, we derive the following result \[\label{eq:lemma195ineq95lambda195eps1} \lambda_{1,k + 1} \leq \lambda_{1,k} + \varepsilon_{1,k}\,,\tag{18}\] where \[\lambda_{1,k} = \alpha_{1,k} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k - 1} \right) \gamma_{1,k}\,, \quad \varepsilon_{1,k} = \frac{1}{2} \alpha \left( \gamma_{1,k - 1} - \gamma_{1,k} \right) + \tau \left( \sqrt{\alpha_{1,k}} + \sqrt{\alpha_{1,k + 1}} \right) \left\lVert f_{1,k} - a_1 B v_k \right\rVert\,.\] By involving the telescoping series cancellation technique to inequality 18 , the following result is obtained \[\begin{align} \lambda_{1,k + 1} &\leq \lambda_{1,1} + \sum_{i = 1}^{k} \varepsilon_{1,i} \\ &= \lambda_{1,1} + \frac{1}{2} \alpha \sum_{i = 1}^{k} \left( \gamma_{1,i - 1} - \gamma_{1,i} \right) + \tau \sum_{i = 1}^{k} \left( \sqrt{\alpha_{1,i}} + \sqrt{\alpha_{1,i + 1}} \right) \left\lVert f_{1,i} - a_1 B v_i \right\rVert \\ &= \lambda_{1,1} + \frac{1}{2} \alpha \left( \gamma_{1,0} - \gamma_{1,k} \right) + \tau \sum_{i = 1}^{k} \left( \sqrt{\alpha_{1,i}} + \sqrt{\alpha_{1,i + 1}} \right) \left\lVert f_{1,i} - a_1 B v_i \right\rVert\,. \end{align}\] Subsequently, by rearranging the terms and taking into account that \(\alpha_{1,i} \leq \lambda_{1,i}\), we derive the following result \[\begin{align} \lambda_{1,k + 1} + \frac{1}{2} \alpha \gamma_{1,k} \leq \lambda_{1,1} + \frac{1}{2} \alpha \gamma_{1,0} + \tau \sum_{i = 1}^{k} \left( \sqrt{\lambda_{1,i}} + \sqrt{\lambda_{1,i + 1}} \right) \left\lVert f_{1,i} - a_1 B v_i \right\rVert\,. \end{align}\] Upon introducing the notation \[\delta_{1,k} = \sqrt{\lambda_{1,k} + \frac{1}{2} \alpha \gamma_{1,k - 1}}\,,\] it follows that \[\label{eq:lemma195ineq95delta1} \delta_{1,k + 1}^{2} \leq \delta_{1,1}^{2} + \tau \sum_{i = 1}^{k} \left( \delta_{1,i} + \delta_{1,i + 1} \right) \left\lVert f_{1,i} - a_1 B v_i \right\rVert\,.\tag{19}\] By employing the technique for inequality 19 as described in Theorem 2.1

[28, pp. Theorem 2.1] , we arrive at the following conclusion \[\label{eq:lemma195final95ineq95dlt1} \delta_{1,k + 1} \leq \delta_{1,1} + 2 \tau \sum_{i = 1}^{k} \left\lVert f_{1,i} - a_1 B v_i \right\rVert\,.\tag{20}\]

By performing the inner product of both sides of equation 6 with \(v_{k + 1} - v_{k - 1} = \Delta v_k + \Delta v_{k - 1}\), and considering the self-adjoint and positive-definite nature of the operator \(A\), one can immediately deduce that \[\alpha_{2,k + 1} + \frac{1}{2} \gamma_{2,k + 1} = \alpha_{2,k} + \frac{1}{2} \gamma_{2,k - 1} + \left( f_{2,k} - a_2 B u_k,\Delta v_k \right) + \left( f_{2,k} - a_2 B u_k,\Delta v_{k - 1} \right)\,,\] where \[\alpha_{2,k} = {\left\lVert \frac{\Delta v_{k - 1}}{\tau} \right\rVert}^{2}\,,\quad \gamma_{2,k} = {\left\lVert L^{\frac{1}{2}} v_k \right\rVert}^{2}\,,\quad L = \gamma A + \delta C\,.\] If we adopt the approach employed to establish inequality 20 , the ensuing result follows \[\label{eq:lemma195final95ineq95dlt2} \delta_{2,k + 1} \leq \delta_{2,1} + 2 \tau \sum_{i = 1}^{k} \left\lVert f_{2,i} - a_2 B u_i \right\rVert\,.\tag{21}\] Here, \[\delta_{2,k} = \sqrt{\lambda_{2,k} + \frac{1}{2} \gamma_{2,k - 1}}\,,\quad \lambda_{2,k} = \alpha_{2,k} + \frac{1}{2} \gamma_{2,k}\,.\] By summing inequalities 20 and 21 , taking into account condition 4 , and introducing the notation \(\delta_k = \delta_{1,k} + \delta_{2,k}\), it consequently follows that \[\begin{align} \label{eq:lemma195ineq95delta} \delta_{k + 1} &\leq \delta_1 + 2 \tau \sum_{i = 1}^{k} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right) + 2 \tau \sum_{i = 1}^{k}\left( \left\lvert a_1 \right\rvert \left\lVert B v_i \right\rVert + \left\lvert a_2 \right\rvert \left\lVert B u_i \right\rVert \right)\nonumber \\ &\leq \delta_1 + 2 \tau \sum_{i = 1}^{k} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right) + \cstDeltaIneq \tau \sum_{i = 1}^{k} \left( \sqrt{\gamma_{1,i}} + \sqrt{\gamma_{2,i}} \right)\,,\quad \cstDeltaIneq = 2 \cstBoper \max\left( \left\lvert a_1 \right\rvert,\left\lvert a_2 \right\rvert \right)\,. \end{align}\tag{22}\] Observe that the following straightforward inequality holds \[\begin{align} \delta_k = \delta_{1,k} + \delta_{2,k} \geq \frac{1}{\sqrt{2}} \left( \sqrt{\alpha}\sqrt{\gamma_{1,k}} + \sqrt{\gamma_{2,k}} \right) \geq \frac{1}{\sqrt{2}} \min\left( 1,\sqrt{\alpha} \right)\left( \sqrt{\gamma_{1,k}} + \sqrt{\gamma_{2,k}} \right)\,, \end{align}\] Hence, it follows that \[\label{eq:lemma195sqrt95gamma} \sqrt{\gamma_{1,k}} + \sqrt{\gamma_{2,k}} \leq \cstGammaIneq \delta_k\,,\quad \cstGammaIneq = \sqrt{2} \max\left( 1,\frac{1}{\sqrt{\alpha}} \right)\,.\tag{23}\] Considering inequality 23 in the estimation of 22 yields the finding that \[\delta_{k + 1} \leq \delta_1 + \cstFinalGammaIneq \tau \sum_{i = 1}^{k} \delta_i + 2 \tau \sum_{i = 1}^{k} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right)\,.\] Hence, through employing the discrete Grönwall-type inequality (cf., e.g., Lemma 3.1

[29, pp. Lemma 3.1] ), it can be established that \[\label{eq:lemma195gronwall95ineqt} \delta_{k + 1} \leq e^{\cstFinalGammaIneq t_k} \left( \delta_1 + 2 \tau \sum_{i = 1}^{k} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right) \right)\,.\tag{24}\] Upon consideration of the estimate, \[\sum_{i = 1}^{k} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right) \leq k \max_{1 \leq i \leq k} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right)\,,\] it is evident that inequality 24 can be formulated as \[\delta_{k + 1} \leq e^{\cstFinalGammaIneq T} \left( \delta_1 + 2 T \max_{1 \leq i \leq n} \left( \left\lVert f_{1,i} \right\rVert + \left\lVert f_{2,i} \right\rVert \right) \right)\,.\] Thus, it can be inferred that \(\alpha_{1,k}\), \(\alpha_{2,k}\), \(\gamma_{1,k}\), and \(\gamma_{2,k}\) are uniformly bounded.

5 Uniform Boundedness of High-Order Terms Corresponding to the Solution of the Discrete Problem↩︎

It should be emphasized that the proof of the convergence of the approximate solution depends, among other considerations, on the uniform boundedness of the vectors \(A^{\frac{1}{2}} \Delta u_{k - 1} / {\tau}\) and \(A u_k\). This is quite natural, as the behavior of the solution to the discrete problem is predominantly governed by the first equation 5 of the system, which represents a difference analogue of the nonlinear Kirchhoff equation (see [29]), excluding the term \(a_1 B v_k\).

In this section, we establish that the vectors \(A^{\frac{1}{2}} \Delta u_{k - 1} / {\tau}\) and \(A u_k\) are locally uniformly bounded. The proof of this result relies on the nonlinear inequality outlined in Lemma [1](#lemma:rogava-tsiklauri1){reference-type=“ref” reference=“lemma:rogava-tsiklauri1”}.

Consider taking the inner product of both sides of equation 5 with \(A\left( u_{k + 1} - u_{k - 1} \right) = A\left( \Delta u_k \right) + A\left( \Delta u_{k - 1} \right)\). By leveraging the properties of the operator \(A\), which is self-adjoint and positive-definite, we obtain: \[\begin{align} \label{eq:lemma295first95norm95eq} {\left\lVert \frac{1}{\tau} A^{\frac{1}{2}} \left( \Delta u_k \right) \right\rVert}^2 &+ \frac{1}{2} \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2} \right) {\left\lVert A u_{k + 1} \right\rVert}^2\nonumber \\ &= {\left\lVert \frac{1}{\tau} A^{\frac{1}{2}} \left( \Delta u_{k - 1} \right) \right\rVert}^2 + \frac{1}{2} \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2} \right) {\left\lVert A u_{k - 1} \right\rVert}^2\nonumber \\ &+ \left( A^{\frac{1}{2}} \left( f_{1,k} - a_1 B v_k \right),A^{\frac{1}{2}}\left( \Delta u_k \right) \right) + \left( A^{\frac{1}{2}} \left( f_{1,k} - a_1 B v_k \right),A^{\frac{1}{2}}\left( \Delta u_{k - 1} \right) \right)\,. \end{align}\tag{25}\] In this context, we assume that \(f_{1,k} - a_1 B v_k\) belongs to \(D\left( A^{\frac{1}{2}} \right)\).

By employing the Cauchy-Schwarz inequality along with inequality ?? from Remark [6](#prop:remark4){reference-type=“ref” reference=“prop:remark4”}, the following result can be derived \[\begin{align} \label{eq:lemma295f95k95cauchy95schwarz} &\left\lvert \left( A^{\frac{1}{2}} \left( f_{1,k} - a_1 B v_k \right),A^{\frac{1}{2}}\left( \Delta u_k \right) \right) + \left( A^{\frac{1}{2}} \left( f_{1,k} - a_1 B v_k \right),A^{\frac{1}{2}}\left( \Delta u_{k - 1} \right) \right) \right\rvert\nonumber \\ \leq& \left( \left\lVert A^{\frac{1}{2}} f_{1,k} \right\rVert + \left\lvert a_1 \right\rvert \cstrmkthree \left\lVert A v_k \right\rVert \right) \tau \left( \left\lVert \frac{1}{\tau} A^{\frac{1}{2}} \left( \Delta u_k \right) \right\rVert + \left\lVert \frac{1}{\tau} A^{\frac{1}{2}} \left( \Delta u_{k - 1} \right) \right\rVert \right)\,. \end{align}\tag{26}\] Let us introduce the following denotations: \[\widetilde{\alpha}_{1,k} = {\left\lVert \frac{1}{\tau} A^{\frac{1}{2}} \left( \Delta u_{k - 1} \right) \right\rVert}^2\,,\quad \beta_{1,k} = {\left\lVert A u_k \right\rVert}^2\,,\quad \gamma_{1,k} = {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2}\,,\quad \sigma_{1,k} = \max_{1 \leq i \leq k} \left\lVert A^{\frac{1}{2}} f_{1,i} \right\rVert\,.\] By substituting the notations introduced in the preceding step and employing inequality 26 , we can rewrite equality 25 as follows \[\begin{align} \label{eq:lemma295notat95ineq95main} \widetilde{\alpha}_{1,k + 1} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \beta_{1,k + 1} &\leq \widetilde{\alpha}_{1,k} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \beta_{1,k - 1}\nonumber \\ &+ \left( \sigma_{1,k} + \left\lvert a_1 \right\rvert \cstrmkthree \left\lVert A v_k \right\rVert \right) \tau \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right)\,. \end{align}\tag{27}\] By incorporating Remark [7](#prop:remark5){reference-type=“ref” reference=“prop:remark5”} into inequality 27 , we obtain the subsequent result \[\begin{align} \widetilde{\alpha}_{1,k + 1} &+ \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \beta_{1,k + 1} \leq \widetilde{\alpha}_{1,k} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \beta_{1,k - 1} \\ &+ \widetilde{M}_{k - 1} \tau \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right) + \cstmainalph \tau^2 \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right) \sum_{i = 1}^{k - 1} \sqrt{\beta_{1,i}}\,, \end{align}\] where \(\widetilde{M}_{k - 1} = \sigma_{1,k} + \left\lvert a_1 \right\rvert \cstrmkthree \hat{M}_{k - 1}\) and \(\cstmainalph = \left\lvert a_1 \right\rvert \cstrmkthree \cstrmkfora = \left\lvert a_1 a_2 \right\rvert \cstrmkthree^2 \cstrmkforb\).

Continuing from the previous step, by adding the term \(\frac{1}{2} \left( \alpha + \beta \gamma_{1,k} \right) \beta_{1,k}\) to both sides of the previous inequality, and introducing the notation \(\widetilde{\beta}_{1,k} = \frac{1}{2} \left( \beta_{1,k} + \beta_{1,k - 1} \right)\), we obtain \[\begin{align} \label{eq:lemma295tilde95M} \widetilde{\alpha}_{1,k + 1} &+ \left( \alpha + \beta \gamma_{1,k} \right) \widetilde{\beta}_{1,k + 1} \leq \widetilde{\alpha}_{1,k} + \left( \alpha + \beta \gamma_{1,k - 1} \right) \widetilde{\beta}_{1,k} + \beta\left( \gamma_{1,k} - \gamma_{1,k - 1} \right)\widetilde{\beta}_{1,k}\nonumber \\ &+ \widetilde{M}_{k - 1} \tau \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right) + \cstmainalph \tau^2 \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right) \sum_{i = 1}^{k - 1} \sqrt{\beta_{1,i}}\,. \end{align}\tag{28}\] To evaluate the absolute value of the difference \(\left\lvert \gamma_{1,k} - \gamma_{1,k - 1} \right\rvert\), we shall employ Lemma [3](#prop:lemma1){reference-type=“ref” reference=“prop:lemma1”} \[\label{eq:lemma295diff95gamma} \begin{align} \left\lvert \gamma_{1,k} - \gamma_{1,k - 1} \right\rvert &\leq 2 \cstmthree \left\lVert A^{\frac{1}{2}} \left( \Delta u_{k - 1} \right) \right\rVert = 2 \cstmthree \tau \sqrt{\widetilde{\alpha}_{1,k}}\,. \end{align}\tag{29}\] Through the use of inequalities 29 and \(\beta_{1,i} \leq 2 \widetilde{\beta}_{1,i}\), we can reformulate inequality 28 as follows \[\begin{align} \label{eq:lemma295tilde95lam95ineq} \widetilde{\lambda}_{1,k + 1} &\leq \widetilde{\lambda}_{1,k} + 2 \cstmthree \beta \tau \sqrt{\widetilde{\alpha}_{1,k}} \widetilde{\beta}_{1,k} + \widetilde{M}_{k - 1} \tau \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right)\nonumber \\ &+ \sqrt{2} \cstmainalph \left( \sqrt{\widetilde{\alpha}_{1,k + 1}} + \sqrt{\widetilde{\alpha}_{1,k}} \right) \tau^2 \sum_{i = 1}^{k - 1} \sqrt{\widetilde{\beta}_{1,i}}\,, \end{align}\tag{30}\] where \(\widetilde{\lambda}_{1,k} = \widetilde{\alpha}_{1,k} + \left( \alpha + \beta \gamma_{1,k - 1} \right) \widetilde{\beta}_{1,k}\).

Given the inequalities \(\widetilde{\alpha}_{1,k} \leq \widetilde{\lambda}_{1,k}\) and \(\widetilde{\beta}_{1,k} \leq \frac{1}{\alpha} \widetilde{\lambda}_{1,k}\), inequality 30 should be rewritten as \[\begin{align} \widetilde{\lambda}_{1,k + 1} &\leq \widetilde{\lambda}_{1,k} + \csttildelama \tau \sqrt{\widetilde{\lambda}_{1,k}} \widetilde{\lambda}_{1,k} + \widetilde{M}_{k} \tau \left( \sqrt{\widetilde{\lambda}_{1,k + 1}} + \sqrt{\widetilde{\lambda}_{1,k}} \right) \\ &+ \csttildelamb \tau^2 \left( \sqrt{\widetilde{\lambda}_{1,k + 1}} + \sqrt{\widetilde{\lambda}_{1,k}} \right) \sum_{i = 1}^{k} \sqrt{\widetilde{\lambda}_{1,i}}\,,\quad \csttildelama = \frac{2 \cstmthree \beta}{\alpha}\,,\quad \csttildelamb = \frac{2 \cstmainalph}{\sqrt{2 \alpha}}\,. \end{align}\]

Let us denote the maximum of \(\widetilde{\lambda}_{1,i}\) for \(1 \leq i \leq k\) by \(\hat{\lambda}_{1,k}\), i.e., \(\hat{\lambda}_{1,k} = \max_{1 \leq i \leq k} \widetilde{\lambda}_{1,i}\). Let \(\widetilde{\lambda}_{1,i + 1}\) (for \(i \leq k\)) attain its maximum at \(i = j\). Thus, we have \(\hat{\lambda}_{1,k + 1} = \widetilde{\lambda}_{1,j + 1}\). Clearly, from the aforementioned inequality, we find that \[\begin{align} \widetilde{\lambda}_{1,j + 1} &\leq \hat{\lambda}_{1,j} + \csttildelama \tau \sqrt{\hat{\lambda}_{1,j}} \hat{\lambda}_{1,j} + \widetilde{M}_{j} \tau \left( \sqrt{\hat{\lambda}_{1,j + 1}} + \sqrt{\hat{\lambda}_{1,j}} \right) \\ &+ \csttildelamb \tau^2 \left( \sqrt{\hat{\lambda}_{1,j + 1}} + \sqrt{\hat{\lambda}_{1,j}} \right) \left( k \sqrt{\hat{\lambda}_{1,j}} \right)\,. \end{align}\] Hence, it follows that \[\begin{align} \hat{\lambda}_{1,k + 1} &\leq \hat{\lambda}_{1,k} + \csttildelama \tau \sqrt{\hat{\lambda}_{1,k}} \hat{\lambda}_{1,k} + \widetilde{M}_{k} \tau \left( \sqrt{\hat{\lambda}_{1,k + 1}} + \sqrt{\hat{\lambda}_{1,k}} \right) \\ &+ \csthatlam \tau \left( \sqrt{\hat{\lambda}_{1,k + 1}} + \sqrt{\hat{\lambda}_{1,k}} \right) \sqrt{\hat{\lambda}_{1,k}}\,,\quad \csthatlam = T \csttildelamb\,. \end{align}\] If we rearrange the terms on the right-hand side of the previously stated inequality, it takes the form \[\label{eq:lemma295hat95lam95rearranged} \hat{\lambda}_{1,k + 1} \leq \left( 1 + \csthatlam \tau + \csttildelama \tau \sqrt{\hat{\lambda}_{1,k}} \right) \hat{\lambda}_{1,k} + \widetilde{M}_{k} \tau \left( \sqrt{\hat{\lambda}_{1,k + 1}} + \sqrt{\hat{\lambda}_{1,k}} \right) + \csthatlam \tau \sqrt{\hat{\lambda}_{1,k} \hat{\lambda}_{1,k + 1}}\,.\tag{31}\] By employing the well-known Young’s inequality for products, one can deduce that: \[\begin{align} \sqrt{\hat{\lambda}_{1,k} \hat{\lambda}_{1,k + 1}} &\leq \frac{1}{3} \left( \hat{\lambda}_{1,k}^{\frac{3}{2}} + 2 \hat{\lambda}_{1,k + 1}^{\frac{3}{4}} \right)\,,\tag{32} \\ \sqrt{\hat{\lambda}_{1,k + 1}} + \sqrt{\hat{\lambda}_{1,k}} &\leq \frac{5}{6} + \frac{1}{2} \hat{\lambda}_{1,k} + \frac{2}{3} \hat{\lambda}_{1,k + 1}^{\frac{3}{4}}\,.\tag{33} \end{align}\] By incorporating inequalities 32 and 33 into bound 31 , one can immediately deduce the following result \[\label{eq:lemma295hat95lam95final} \hat{\lambda}_{1,k + 1} \leq \left( 1 + \cstlemtwoa \tau + \cstlemtwob \tau \sqrt{\hat{\lambda}_{1,k}} \right) \hat{\lambda}_{1,k} + \cstlemtwoc \tau \hat{\lambda}_{1,k + 1}^{\frac{3}{4}} + \cstlemtwod \tau\,.\tag{34}\] Here, \[\cstmaxMk = \max\limits_{1 \leq k \leq n - 1} \widetilde{M}_{k}\,,\quad \cstlemtwoa = \csthatlam + \frac{\cstmaxMk}{2}\,,\quad \cstlemtwob = \csttildelama + \frac{\csthatlam}{3}\,,\quad \cstlemtwoc = \frac{2}{3}\left( \csthatlam + \cstmaxMk \right)\,,\quad \cstlemtwod = \frac{5 \cstmaxMk}{6}\,.\] By introducing the notations \[w_k = \left( 1 + \cstlemtwoa \tau + \cstlemtwob \tau \sqrt{\hat{\lambda}_{1,k}} \right) \hat{\lambda}_{1,k} + \cstlemtwod \tau\] and \(y_{k + 1} = \sqrt[4]{\hat{\lambda}_{1,k + 1}}\), the estimate 34 should be rewritten as follows \[y_{k + 1}^4 - \cstlemtwoc \tau y_{k + 1}^3 - w_k \leq 0\,.\] From here, it follows \[\label{eq:lemma295fourth95ord95ineq} \left( \frac{y_{k + 1}}{\cstlemtwoc \tau} \right)^4 - \left( \frac{y_{k + 1}}{\cstlemtwoc \tau} \right)^3 - \frac{w_k}{\left( \cstlemtwoc \tau \right)^4} \leq 0\,.\tag{35}\]

Consider the polynomial \(P\left( \xi \right) = \xi^4 - \xi^3 - b\), where \(b = w_k / \left( \cstlemtwoc \tau \right)^4\). This polynomial corresponds to the left-hand side of inequality 35 . It is evident that \(P\left( \xi \right)\) has exactly two real roots, one positive and one negative. Since \(P\left( 0 \right) = -b < 0\) and \(P\left( 1 + \sqrt[4]{b} \right) > 0\), it follows that the positive root of \(P\left( \xi \right)\) must be less than \(1 + \sqrt[4]{b}\). Considering these facts, from inequality 35 , we have \[\frac{y_{k + 1}}{\cstlemtwoc \tau} \leq 1 + \sqrt[4]{b}\,,\] or which is the same \[y_{k + 1} \leq \cstlemtwoc \tau + \sqrt[4]{w_k}\,.\] By squaring both sides of the given inequality and substituting the product term with the sum of squares, we obtain the following \[y_{k + 1}^2 \leq \left( 1 + \cstlemtwoc \tau \right) \left( \cstlemtwoc \tau + \sqrt{w_k} \right)\,.\] If we apply the same transformation as in the previous case, we have \[\label{eq:lemma295fourth95y} y_{k + 1}^4 \leq \left( 1 + \cstlemtwoc \tau \right)^3 \left( \cstlemtwoc \tau + w_k \right)\,.\tag{36}\] We may now revert to the previous notation. In this case, inequality 36 takes the form \[\label{eq:lemma295hat95lam95bef95div} \hat{\lambda}_{1,k + 1} \leq \left( 1 + \cstmaxab \tau \right)^4 \left( 1 + \frac{\cstlemtwob \tau}{1 + \cstlemtwoa \tau} \sqrt{\hat{\lambda}_{1,k}} \right) \hat{\lambda}_{1,k} + \left( 1 + \cstmaxab \tau \right)^3 \cstykplusonefin \tau\,,\tag{37}\] where \(\cstykplusonefin = \cstlemtwoc + \cstlemtwod\) and \(\cstmaxab = \max\left( \cstlemtwoa,\cstlemtwoc \right)\).

If we replace \(\left( 1 + \cstmaxab \tau \right)^4\) with \(1 + \cstfourthpoly \tau\) in inequality 37 , and subsequently divide both sides of the resulting inequality by \(\left( 1 + \cstfourthpoly \tau \right)^{k + 1}\), performing straightforward transformations, we obtain the following result \[\xi_{k + 1} \leq \xi_k \left( 1 + \cstsqrtfracinxi \tau \sqrt{\xi_k} \right) + \cstxikfin \tau\,,\] where \[\xi_k = \frac{\hat{\lambda}_{1,k}}{\left( 1 + \cstfourthpoly \tau \right)^k}\,.\] Consider the transformation where \(\overline{\tau} = \cstoverlinetau \tau\) with \(\cstoverlinetau = \max\left( \cstsqrtfracinxi,\cstxikfin \right)\). Then we have \[\xi_{k + 1} \leq \xi_k \left( 1 + \overline{\tau} \sqrt{\xi_k} \right) + \overline{\tau}\,.\] From Lemma [1](#lemma:rogava-tsiklauri1){reference-type=“ref” reference=“lemma:rogava-tsiklauri1”}, it follows that \[\label{eq:lemma295xi95and95hat95gamma} \xi_k \leq \frac{\xi}{\left( 1 - \overline{t}_k \sqrt{\xi} \right)^2} \leq \frac{\hat{\lambda}}{\left( 1 - \cstoverlinetau \sqrt{\hat{\lambda}} t_k \right)^2}\,,\quad k = 1,2,\ldots,m\,,\tag{38}\] where \[\xi = \max\left( 1,\xi_1 \right) \leq \max\left( 1,\hat{\lambda}_{1,1} \right) = \hat{\lambda}\,,\quad \overline{t}_k = k \overline{\tau} = \cstoverlinetau t_k < \frac{1}{\sqrt{\hat{\lambda}}} \leq \frac{1}{\sqrt{\xi}}\,.\] Following the established notation, we have \[\label{eq:lemma295xi95with95exp95denom} \xi_k = \frac{\hat{\lambda}_{1,k}}{\left( 1 + \cstfourthpoly \tau \right)^k} \geq \frac{\hat{\lambda}_{1,k}}{e^{\cstfourthpoly t_k}}\,.\tag{39}\] Using inequality 38 along with 39 , we obtain the following estimate \[\label{eq:lemma295final95ineq} \hat{\lambda}_{1,k} \leq \frac{\hat{\lambda}}{\left( 1 - \cstoverlinetau \sqrt{\hat{\lambda}} t_k \right)^2} e^{\cstfourthpoly t_k}\,,\quad k = 1,2,\ldots,m\,.\tag{40}\]

It should be noted that the value of \(m\) in the inequality 40 is influenced both by the coefficient of \(t_k\) (appearing in the denominator of the fraction) and by the number of subdivisions of the time interval, denoted by \(n\). The coefficient of \(t_k\) can be explicitly determined based on the data provided in problem 1 3 , while also taking into account the value of \(T\). Furthermore, the inequality \({\hat{\lambda}} \leq {\widetilde{M}}\) is satisfied, where \({\widetilde{M}}\) represents a positive constant that is dependent on the original data from the stated problem 1 3 , as well as the value of \(T\).

The following inequality is derived from 40 \[\label{eq:lemma295final95result95loc95bound} \hat{\lambda}_{1,k} \leq \frac{\widetilde{M}}{\left( 1 - \overline{M}{\,}\overline{T} \right)^2} e^{\cstfourthpoly \overline{T}}\,,\quad k = 1, 2, \ldots, \left[ \frac{\overline{T}}{\tau} \right]\,,\tag{41}\] where \(\overline{M} = \cstoverlinetau \sqrt{\widetilde{M}}\) and \(\displaystyle \overline{T} = \frac{q}{\overline{M}}\), with \(0 < q < 1\).

The inequality given in 41 implies that the vectors \(A u_k\) and \(A^{\frac{1}{2}} \Delta u_{k - 1} / {\tau}\) are uniformly bounded over the local interval \(\left[ 0,\overline{T} \right]\). One should observe that the local uniform boundedness of the vectors \(A v_k\) follows directly from Remark [7](#prop:remark5){reference-type=“ref” reference=“prop:remark5”}.

6 Error Estimates for Approximate Solutions of the System↩︎

Before the theorem concerning the convergence of the scheme 5 6 is presented, a remark regarding the smoothness of the solutions to the problem 1 3 is made to clarify the order of convergence of the proposed symmetric three-layer semi-discrete scheme 5 6 . A minimum degree of smoothness in the solutions is required to ensure the well-posedness of the problem. It should be noted that this condition guarantees convergence but is insufficient for determining the order of convergence. When the smoothness of the solutions is increased by one degree, an order of convergence equal to one is achieved. Nevertheless, in this case, as well as in the previous case, the following initial conditions should be taken: \(u_1 = \varphi_0 + \tau \varphi_1\) and \(v_1 = \psi_0 + \tau \psi_1\). Furthermore, if the smoothness is increased by two degrees and the initial functions are specified according to formulas 7 and 8 , an additional degree of convergence is attained, resulting in a total order of two. However, any further increase in smoothness would be regarded as superfluous, since the approximation order of the scheme 5 6 does not exceed two.

The following theorem is formulated to address the convergence of the scheme 5 6 .

The initial step is a standard procedure. We shall derive the system corresponding to the scheme 5 6 for the errors associated with the approximate solutions \(z_{1,k}\) and \(z_{2,k}\), and proceed to evaluate the remainder terms. Following this, equations 1 and 2 are reformulated at the discrete time points \(t = t_k\), where \(k = 1,2,\ldots,n - 1\), as follows: \[\begin{align}\tag{42} \frac{\Delta^2 u \left( t_{k - 1} \right)}{\tau^2} &+ \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u \left( t_k \right) \right\rVert}^2 \right) \frac{A u \left( t_{k + 1} \right) + A u \left( t_{k - 1} \right)}{2}\nonumber \\ &= \left( f_1\left( t_k \right) - a_1 B v\left( t_k \right) \right) + R_{1,k}\left( \tau \right) + R_{2,k}\left( \tau \right)\,, \end{align} \begin{equation}\tag{43} \frac{\Delta^2 v \left( t_{k - 1} \right)}{\tau^2} + \frac{L v\left( t_{k + 1} \right) + L v\left( t_{k - 1} \right)}{2} = \left( f_2\left( t_k \right) - a_2 B u\left( t_k \right) \right) + R_{3,k}\left( \tau \right) + R_{4,k}\left( \tau \right)\,. \end{equation}\] Here, \[\begin{align} R_{1,k}\left( \tau \right) &= \frac{\Delta^2 u \left( t_{k - 1} \right)}{\tau^2} - \frac{\mathrm{d}^2 u\left( t_k \right)}{\mathrm{d}t^2}\,,\quad R_{2,k}\left( \tau \right) = \frac{1}{2} \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u \left( t_k \right) \right\rVert}^2 \right) A\left( \Delta^2 u \left( t_{k - 1} \right) \right)\,, \\ R_{3,k}\left( \tau \right) &= \frac{\Delta^2 v \left( t_{k - 1} \right)}{\tau^2} - \frac{\mathrm{d}^2 v\left( t_k \right)}{\mathrm{d}t^2}\,,\quad R_{4,k}\left( \tau \right) = \frac{1}{2} L\left( \Delta^2 v \left( t_{k - 1} \right) \right)\,. \end{align}\]

Based on conditions [itm95theorem95b] and [itm95theorem95c] of Theorem [8](#prop:main95thm95convrg){reference-type=“ref” reference=“prop:main95thm95convrg”}, the following estimates can be derived for the remainder terms: \[\label{eq:main95thm95convrg95est95rems} \left\lVert R_{ cboundremaindersvariable\endcsname ,k}\left( \tau \right) \right\rVert \leq \boundremainders \tau^2\,,\quad cboundremaindersvariable\endcsname = 1,2,3,4\,.\tag{44}\]

By subtracting equations 42 and 43 from equations 5 and 6 , respectively, we obtain the following system of equations: \[\begin{gather} \frac{\Delta^2 z_{1,k - 1}}{\tau^2} + \left( \alpha + \beta {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^2 \right) \frac{A z_{1,k + 1} + A z_{1,k - 1}}{2} = \widetilde{g}_k - a_1 B z_{2,k}\,,\tag{45} \\ \frac{\Delta^2 z_{2,k - 1}}{\tau^2} + \frac{L z_{2,k + 1} + L z_{2,k - 1}}{2} = R_k\left( \tau \right) - a_2 B z_{1,k}\,,\tag{46} \end{gather}\] where \[\begin{gather} \widetilde{g}_k = \beta \left( {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^2 - {\left\lVert A^{\frac{1}{2}} u \left( t_k \right) \right\rVert}^2 \right) \frac{A u \left( t_{k + 1} \right) + A u \left( t_{k - 1} \right)}{2} + R_{1,k}\left( \tau \right) + R_{2,k}\left( \tau \right)\,, \\ R_k\left( \tau \right) = R_{3,k}\left( \tau \right) + R_{4,k}\left( \tau \right)\,. \end{gather}\]

Taking the inner product of both sides of equality 45 with \(z_{1,k + 1} - z_{1,k - 1} = \Delta z_{1,k} + \Delta z_{1,k - 1}\), and considering the self-adjointness and positive definiteness of the operator \(A\), we arrive at the following result \[\label{eq:main95thm95convrg95inner95prod95eqt} \overline{\lambda}_{1,k + 1} = \overline{\lambda}_{1,k} + \left( \overline{\varepsilon}_{1,k} + \eta_{1,k} + \eta_{2,k} \right)\,,\tag{47}\] where: \[\begin{gather} \overline{\lambda}_{1,k} = \overline{\alpha}_{1,k}^{2} + \frac{1}{2}\left( \alpha + \beta \gamma_{1,k - 1} \right) \overline{\gamma}_{1,k}^{2}\,,\quad \overline{\alpha}_{1,k} = \left\lVert \frac{\Delta z_{1,k - 1}}{\tau} \right\rVert\,,\quad \overline{\gamma}_{1,k} = \left\lVert A^{\frac{1}{2}} z_{1,k} \right\rVert\,, \\ \gamma_{1,k} = {\left\lVert A^{\frac{1}{2}} u_k \right\rVert}^{2}\,,\quad \overline{\varepsilon}_{1,k} = \frac{1}{2} \alpha \left( \overline{\gamma}_{1,k - 1}^{2} - \overline{\gamma}_{1,k}^{2} \right) + \frac{1}{2} \beta \left( \gamma_{1,k} \overline{\gamma}_{1,k - 1}^{2} - \gamma_{1,k - 1} \overline{\gamma}_{1,k}^{2} \right)\,, \\ \eta_{1,k} = \left( \widetilde{g}_k,\Delta z_{1,k} + \Delta z_{1,k - 1} \right)\,,\quad \eta_{2,k} = -a_1 \left( B z_{2,k},\Delta z_{1,k} + \Delta z_{1,k - 1} \right)\,. \end{gather}\]

The resulting recurrence relation 47 provides a “descent” from \(k\) to \(1\). Therefore, it follows from 47 that \[\label{eq:main95thm95convrg95sum95inner95prod} \begin{align} \overline{\lambda}_{1,k + 1} &= \overline{\lambda}_{1,1} + \sum_{i = 1}^{k} \left( \overline{\varepsilon}_{1,i} + \eta_{1,i} + \eta_{2,i} \right) \\ &= \overline{\lambda}_{1,1} + \frac{1}{2} \alpha \left( \overline{\gamma}_{1,0}^{2} - \overline{\gamma}_{1,k}^{2} \right) + \frac{1}{2} \beta \sum_{i = 1}^{k} \left( \gamma_{1,i} \overline{\gamma}_{1,i - 1}^{2} - \gamma_{1,i - 1} \overline{\gamma}_{1,i}^{2} \right) + \sum_{i = 1}^{k} \left( \eta_{1,i} + \eta_{2,i} \right)\,. \end{align}\tag{48}\] We observe that for the following sum, the subsequent representation is satisfied \[\label{eq:main95thm95convrg95sum95term95repr} \sum_{i = 1}^{k} \left( \gamma_{1,i} \overline{\gamma}_{1,i - 1}^{2} - \gamma_{1,i - 1} \overline{\gamma}_{1,i}^{2} \right) = \gamma_{1,1} \overline{\gamma}_{1,0}^{2} + \sum_{i = 1}^{k - 1} \overline{\gamma}_{1,i}^{2}\left( \gamma_{1,i + 1} - \gamma_{1,i - 1} \right) - \gamma_{1,k - 1} \overline{\gamma}_{1,k}^{2}\,.\tag{49}\] By substituting the corresponding term in the equality 48 with the sum representation given in expression 49 , the following relation is established \[\label{eq:main95thm95convrg95overl95lamb95eq} \begin{align} \overline{\lambda}_{1,k + 1} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k - 1} \right) \overline{\gamma}_{1,k}^{2} &= \overline{\lambda}_{1,1} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,1} \right) \overline{\gamma}_{1,0}^{2} \\ &+ \frac{1}{2} \beta \sum_{i = 1}^{k - 1} \overline{\gamma}_{1,i}^{2}\left( \gamma_{1,i + 1} - \gamma_{1,i - 1} \right) + \sum_{i = 1}^{k} \left( \eta_{1,i} + \eta_{2,i} \right)\,. \end{align}\tag{50}\] We shall introduce the following notation \[\overline{\delta}_{1,k} = \sqrt{\overline{\lambda}_{1,k} + \frac{1}{2} \left( \alpha + \beta \gamma_{1,k - 2} \right) \overline{\gamma}_{1,k - 1}^{2}}\,.\] Under this denotation, assuming \(\gamma_{1,-1} = \gamma_{1,1}\), equality 50 can be expressed as follows \[\label{eq:main95thm95convrg95tild95delta95one} \overline{\delta}_{1,k + 1}^{2} = \overline{\delta}_{1,1}^{2} + \frac{1}{2} \beta \sum_{i = 1}^{k - 1} \overline{\gamma}_{1,i}^{2}\left( \gamma_{1,i + 1} - \gamma_{1,i - 1} \right) + \sum_{i = 1}^{k} \left( \eta_{1,i} + \eta_{2,i} \right)\,.\tag{51}\]

Following Lemma [3](#prop:lemma1){reference-type=“ref” reference=“prop:lemma1”} and Lemma [4](#prop:lemma2){reference-type=“ref” reference=“prop:lemma2”}, an estimate for the difference \(\gamma_{1,i + 1} - \gamma_{1,i - 1}\) can be derived \[\begin{align} \label{eq:main95thm95convrg95abs95diff95gamma} \left\lvert \gamma_{1,i + 1} - \gamma_{1,i - 1} \right\rvert &\leq \left( \left\lvert \sqrt{\gamma_{1,i + 1}} - \sqrt{\gamma_{1,i}} \right\rvert + \left\lvert \sqrt{\gamma_{1,i}} - \sqrt{\gamma_{1,i - 1}} \right\rvert \right) \left( \sqrt{\gamma_{1,i + 1}} + \sqrt{\gamma_{1,i - 1}} \right)\nonumber \\ &\leq \tau \left( \left\lVert A^{\frac{1}{2}} \frac{\Delta u_{i}}{\tau} \right\rVert + \left\lVert A^{\frac{1}{2}} \frac{\Delta u_{i - 1}}{\tau} \right\rVert \right) \left( \left\lVert A^{\frac{1}{2}} u_{i + 1} \right\rVert + \left\lVert A^{\frac{1}{2}} u_{i - 1} \right\rVert \right) \leq 4 \cstmthree \cstmfive \tau\,. \end{align}\tag{52}\] On the other hand, it is easy to see that \[\label{eq:main95thm95convrg95overl95gamma95est} \overline{\gamma}_{1,i}^{2} \leq \frac{2}{\alpha} \overline{\delta}_{1,i}^{2}\,.\tag{53}\]

By employing the Cauchy-Schwarz inequality, we can ascertain that \[\label{eq:main95thm95convrg95etaone} \left\lvert \eta_{1,i} \right\rvert \leq \left\lVert \widetilde{g}_i \right\rVert \left( \left\lVert \Delta z_{1,i} \right\rVert + \left\lVert \Delta z_{1,i - 1} \right\rVert \right) = \tau \left\lVert \widetilde{g}_i \right\rVert \left( \overline{\alpha}_{1,i + 1} + \overline{\alpha}_{1,i} \right)\,.\tag{54}\] We shall proceed to estimate \(\left\lVert \widetilde{g}_i \right\rVert\). To achieve this, it is necessary to assess certain quantities. By virtue of inequality 44 , we conclude that \(\left\lVert R_{1,k}\left( \tau \right) \right\rVert + \left\lVert R_{2,k}\left( \tau \right) \right\rVert \leq \inthmsumonetwo \tau^2\). Taking into account this fact, together with estimate 53 , Lemma [3](#prop:lemma1){reference-type=“ref” reference=“prop:lemma1”}, and condition [itm95theorem95c] of Theorem [8](#prop:main95thm95convrg){reference-type=“ref” reference=“prop:main95thm95convrg”}, we establish that the following inequality is valid \[\begin{align} \label{eq:main95thm95convrg95norm95tild95g} \left\lVert \widetilde{g}_i \right\rVert &\leq \frac{\beta}{2} \left\lvert {\left\lVert A^{\frac{1}{2}} u_i \right\rVert}^2 - {\left\lVert A^{\frac{1}{2}} u \left( t_i \right) \right\rVert}^2 \right\rvert \left( \left\lVert A u \left( t_{i + 1} \right) \right\rVert + \left\lVert A u \left( t_{i - 1} \right) \right\rVert \right) + \left\lVert R_{1,i}\left( \tau \right) \right\rVert + \left\lVert R_{2,i}\left( \tau \right) \right\rVert\nonumber \\ &\leq \beta \overline{\gamma}_{1,i} \left( \cstmthree + \inthmmaxhalfa \right) \inthmmaxa + \inthmsumonetwo \tau^2 \leq \beta \sqrt{2 {\alpha}^{-1}} \overline{\delta}_{1,i} \left( \cstmthree + \inthmmaxhalfa \right) \inthmmaxa + \inthmsumonetwo \tau^2 = \inthmtilgfin \overline{\delta}_{1,i} + \inthmsumonetwo \tau^2\,, \end{align}\tag{55}\] where \(\inthmmaxhalfa = \nu^{-\frac{1}{2}} \inthmmaxa\) and \(\inthmmaxa = \max_{0 \leq t \leq T} \left\lVert A u \left( t \right) \right\rVert\).

Given that \(\overline{\alpha}_{1,i}^{2} \leq \overline{\lambda}_{1,i} \leq \overline{\delta}_{1,i}^{2}\), and in view of the derived inequality 55 , the following estimate for 54 can be established \[\label{eq:main95thm95convrg95abs95eta95one} \left\lvert \eta_{1,i} \right\rvert \leq \tau \left( \inthmtilgfin \overline{\delta}_{1,i} + \inthmsumonetwo \tau^2 \right) \left( \overline{\delta}_{1,i + 1} + \overline{\delta}_{1,i} \right)\,.\tag{56}\] By applying the Cauchy-Schwarz inequality along with the estimate ?? , we easily find that \[\begin{align} \left\lvert \eta_{2,i} \right\rvert &\leq \left\lvert a_1 \right\rvert \left\lVert B z_{2,i} \right\rVert \left( \left\lVert \Delta z_{1,i} \right\rVert + \left\lVert \Delta z_{1,i - 1} \right\rVert \right) = \tau \left\lvert a_1 \right\rvert \left\lVert B z_{2,i} \right\rVert \left( \overline{\alpha}_{1,i + 1} + \overline{\alpha}_{1,i} \right) \\ &\leq \tau \left\lvert a_1 \right\rvert \left\lVert B z_{2,i} \right\rVert \left( \overline{\delta}_{1,i + 1} + \overline{\delta}_{1,i} \right) \leq \tau \left\lvert a_1 \right\rvert \cstBoper \left\lVert A^{\frac{1}{2}} z_{2,i} \right\rVert \left( \overline{\delta}_{1,i + 1} + \overline{\delta}_{1,i} \right)\,. \end{align}\] Hence, from Remark [3](#prop:remark2){reference-type=“ref” reference=“prop:remark2”}, the following inequality follows: \[\label{eq:main95thm95convrg95abs95eta95two} \left\lvert \eta_{2,i} \right\rvert \leq \inthmnutwo \tau \overline{\gamma}_{2,i} \left( \overline{\delta}_{1,i + 1} + \overline{\delta}_{1,i} \right)\,,\quad \overline{\gamma}_{2,i} = \left\lVert L^{\frac{1}{2}} z_{2,i} \right\rVert\,,\quad \inthmnutwo = \frac{\left\lvert a_1 \right\rvert \cstBoper}{\sqrt{\gamma}}\,.\tag{57}\] By incorporating 52 , 53 , 56 , and 57 into expression 51 , we derive the following conclusion \[\begin{align} \overline{\delta}_{1,k + 1}^{2} \leq \overline{\delta}_{1,1}^{2} &+ \inthmsumovrldlt \tau \sum_{i = 1}^{k - 1} \overline{\delta}_{1,i}^{2} \\ &+ \tau \sum_{i = 1}^{k} \left( \inthmtilgfin \overline{\delta}_{1,i} +\inthmnutwo \overline{\gamma}_{2,i} + \inthmsumonetwo \tau^2 \right) \left( \overline{\delta}_{1,i + 1} + \overline{\delta}_{1,i} \right)\,,\quad \inthmsumovrldlt = \frac{4 \cstmthree \cstmfive \beta}{\alpha}\,. \end{align}\] By implementing the approach for the previously established inequality, as indicated in Theorem 2.1

[28, pp. Theorem 2.1] , we arrive at the following result: \[\begin{align} \label{eq:main95thm95convrg95overl95deltone95final} \overline{\delta}_{1,k + 1} &\leq \overline{\delta}_{1,1} + \inthmsumovrldlt \tau \sum_{i = 1}^{k} \overline{\delta}_{1,i} + 2 \tau \sum_{i = 1}^{k} \left( \inthmtilgfin \overline{\delta}_{1,i} +\inthmnutwo \overline{\gamma}_{2,i} + \inthmsumonetwo \tau^2 \right) \nonumber \\ &= \overline{\delta}_{1,1} + \left( 2 \inthmtilgfin + \inthmsumovrldlt \right) \tau \sum_{i = 1}^{k} \overline{\delta}_{1,i} + 2 \inthmnutwo \tau \sum_{i = 1}^{k} \overline{\gamma}_{2,i} + 2 \inthmsumonetwo t_k \tau^2\,. \end{align}\tag{58}\]

Let us now move forward to derive an estimate of the type 58 for the equation 46 . To do this, we adopt a similar approach as was used for the previous equation 45 . Specifically, we apply the inner product to both sides of equality 46 with \(z_{2,k + 1} - z_{2,k - 1} = \Delta z_{2,k} + \Delta z_{2,k - 1}\), employing the properties of the operator \(L = \gamma A + \delta C\), which is self-adjoint and positive definite. Furthermore, reasoning analogous to that used to establish the bound 21 in Lemma [3](#prop:lemma1){reference-type=“ref” reference=“prop:lemma1”} leads us to the following relation: \[\label{eq:main95thm95convrg95overl95delt95z95two} \overline{\delta}_{2,k + 1} \leq \overline{\delta}_{2,1} + 2 \tau \sum_{i = 1}^{k} \left\lVert R_i\left( \tau \right) \right\rVert + 2 \left\lvert a_2 \right\rvert \tau \sum_{i = 1}^{k} \left\lVert B z_{1,i} \right\rVert\,,\tag{59}\] for which \[\overline{\delta}_{2,k} = \sqrt{\overline{\lambda}_{2,k} + \frac{1}{2} \overline{\gamma}_{2,k - 1}^{2}}\,,\quad \overline{\lambda}_{2,k} = \overline{\alpha}_{2,k}^{2} + \frac{1}{2} \overline{\gamma}_{2,k}^{2}\,,\quad \overline{\alpha}_{2,k} = \left\lVert \frac{\Delta z_{2,k - 1}}{\tau} \right\rVert\,,\quad \overline{\gamma}_{2,k} = \left\lVert L^{\frac{1}{2}} z_{2,k} \right\rVert\,.\]

It can be observed that by employing the estimates 44 and 53 alongside Remark [5](#prop:remark6){reference-type=“ref” reference=“prop:remark6”}, from inequality 59 we can immediately deduce the following result: \[\label{eq:main95thm95convrg95overl95delttwo95final} \overline{\delta}_{2,k + 1} \leq \overline{\delta}_{2,1} + \inthmfrsumoverldlt \tau \sum_{i = 1}^{k} \overline{\delta}_{1,i} + 2 \inthmsumthreefour t_k \tau^2\,,\quad \inthmfrsumoverldlt = \frac{2 \sqrt{2} \left\lvert a_2 \right\rvert \cstBoper}{\sqrt{\alpha}}\,.\tag{60}\]

It is evident that a similar inequality holds for \(\overline{\gamma}_{2,i}\), as given in 53 , specifically \(\overline{\gamma}_{2,i}^{2} \leq 2 \overline{\delta}_{2,i}^{2}\). Taking this inequality into consideration in the bound 58 , we arrive at the following result: \[\label{eq:main95thm95convrg95overl95deltone95latest} \overline{\delta}_{1,k + 1} \leq \overline{\delta}_{1,1} + \left( 2 \inthmtilgfin + \inthmsumovrldlt \right) \tau \sum_{i = 1}^{k} \overline{\delta}_{1,i} + 2 \sqrt{2} \inthmnutwo \tau \sum_{i = 1}^{k} \overline{\delta}_{2,i} + 2 \inthmsumonetwo t_k \tau^2\,.\tag{61}\] By summing inequalities 60 and 61 and introducing the notation \(\overline{\delta}_{k} = \overline{\delta}_{1,k} + \overline{\delta}_{2,k}\), we can infer that: \[\overline{\delta}_{k + 1} \leq \overline{\delta}_{1} + \inthmsumovrlndlt \tau \sum_{i = 1}^{k} \overline{\delta}_{i} + \inthmassoctau t_k \tau^2\,,\] in which \[\inthmsumovrlndlt = \max\left( 2 \inthmtilgfin + \inthmsumovrldlt + \inthmfrsumoverldlt, 2 \sqrt{2} \inthmnutwo \right)\, \text{and}\,\, \inthmassoctau = 2\left( \inthmsumonetwo + \inthmsumthreefour \right)\,.\] Thus, by employing the discrete Grönwall-type inequality (cf. Lemma 3.1

[29, pp. Lemma 3.1] ), it can be concluded that: \[\label{eq:main95thm95convrg95reslt95gronwl} \overline{\delta}_{k + 1} \leq {e}^{\inthmsumovrlndlt t_k} \left( \overline{\delta}_{1} + \inthmassoctau t_k \tau^2 \right) \leq {e}^{\inthmsumovrlndlt T} \left( \overline{\delta}_{1} + \inthmassoctau T \tau^2 \right)\,.\tag{62}\]

In order to derive the desired estimates, it becomes essential to take an additional step. Specifically, this entails estimating the term \(\overline{\delta}_{1}\) in inequality 62 . On the other hand, the term \(\overline{\delta}_{1}\) is expressed as the sum of \(\overline{\delta}_{1,1}\) and \(\overline{\delta}_{2,1}\). The first summand involves the quantities \(\overline{\alpha}_{1,1}\), \(\overline{\gamma}_{1,0}\), and \(\overline{\gamma}_{1,1}\), whereas the second summand contains the quantities \(\overline{\alpha}_{2,1}\), \(\overline{\gamma}_{2,0}\), and \(\overline{\gamma}_{2,1}\).

In view of assumptions [itm95theorem95a]–[itm95theorem95c] of Theorem [8](#prop:main95thm95convrg){reference-type=“ref” reference=“prop:main95thm95convrg”}, the following estimates hold: \[\label{eq:main95thm95convrg95overln95alphbet} \overline{\alpha}_{1,1} \leq \inthmalphoneone \tau^2\,,\quad \overline{\gamma}_{1,1} \leq \inthmagammoneone \tau^2\,,\quad \overline{\alpha}_{2,1} \leq \inthmalphtwoone \tau^2\,,\quad \overline{\gamma}_{2,1} \leq \inthmagammtwoone \tau^2\,.\tag{63}\]

By employing these estimates given in 63 , the following bound is established: \[\label{eq:main95thm95convrg95overl95delta95one} \overline{\delta}_{1} = \overline{\delta}_{1,1} + \overline{\delta}_{2,1} \leq \inthmoverldltone \tau^2\,.\tag{64}\]

From 62 , and considering the bound provided in 64 , the estimates for Theorem [8](#prop:main95thm95convrg){reference-type=“ref” reference=“prop:main95thm95convrg”} are derived.

7 Spatial One-Dimensional Nonlinear Dynamic Timoshenko Model↩︎

Observe that the system 1 2 is an abstract analogue of the spatial one-dimensional nonlinear dynamic Timoshenko model, which takes the following form: \[\begin{gather} \frac{\partial^2 u\left( x,t \right)}{\partial t^2} - \left( \alpha + \beta \int\limits_{0}^{\ell} \left[ \frac{\partial u\left( x,t \right)}{\partial x} \right]^2 \mathrm{d}x \right) \frac{\partial^2 u\left( x,t \right)}{\partial x^2} + a_1 \frac{\partial v\left( x,t \right)}{\partial x} = f_1\left( x,t \right)\,,\tag{65} \\ \frac{\partial^2 v\left( x,t \right)}{\partial t^2} - \gamma \frac{\partial^2 v\left( x,t \right)}{\partial x^2} + \delta v\left( x,t \right) - a_2 \frac{\partial u\left( x,t \right)}{\partial x} = f_2\left( x,t \right)\,.\tag{66} \end{gather}\] In this setting, \(\left( x,t \right) \in \left( 0,\ell \right) \times \left( 0,T \right]\); the constants \(\alpha\), \(\beta\), \(\gamma\), \(\delta\), \(a_1\), and \(a_2\) are positive; the functions \(f_1\left( x,t \right)\) and \(f_2\left( x,t \right)\) are continuous over the prescribed domain.

For the system 65 66 , we pose the following initial-boundary value problem: \[\begin{gather} u\left( x,0 \right) = \varphi_0\left( x \right)\,,\quad u_{t}^{\prime}\left( x,0 \right) = \varphi_1\left( x \right)\,,\quad v\left( x,0 \right) = \psi_0\left( x \right)\,,\quad v_{t}^{\prime}\left( x,0 \right) = \psi_1\left( x \right)\,,\tag{67} \\ u\left( 0,t \right) = 0\,,\quad u\left( \ell,t \right) = 0\,,\quad v\left( 0,t \right) = 0\,,\quad v\left( \ell,t \right) = 0\,.\tag{68} \end{gather}\] Here, the functions \(\varphi_0\left( x \right)\), \(\varphi_1\left( x \right)\), \(\psi_0\left( x \right)\), and \(\psi_1\left( x \right)\) are continuous. Moreover, compatibility conditions are assumed to hold for \(\varphi_0\left( x \right)\) and \(\psi_0\left( x \right)\), specifically that \(\varphi_0\left( 0 \right) = \varphi_0\left( \ell \right) = 0\) and \(\psi_0\left( 0 \right) = \psi_0\left( \ell \right) = 0\); the functions \(u\left( x,t \right)\) and \(v\left( x,t \right)\) are unknown.

7.1 Semi-Discrete Scheme↩︎

We shall establish a uniform grid for the time domain \(\left[ 0,T \right]\) with a step size of \(\tau\), defined as follows: \[0 < t_0 < t_1 < \cdots < t_{n - 1} < t_n = T\,,\quad t_k = k \tau\,,\quad k = 0,1,\ldots,n\,,\quad \tau = \frac{T}{n}\,.\] Let us write a semi-discrete scheme for the problem 65 68 based on the proposed scheme 5 6 for the abstract coupled system. It has the following form: \[\begin{align} \tag{69} \frac{{\Delta}^2 u_{k - 1}\left( x \right)}{\tau^2} - \frac{1}{2} q_k \left( \frac{\mathrm{d}^2 u_{k + 1}\left( x \right)}{\mathrm{d}x^2} + \frac{\mathrm{d}^2 u_{k - 1}\left( x \right)}{\mathrm{d}x^2} \right) = f_{1,k}\left( x \right) - a_1 \frac{\mathrm{d}v_k\left( x \right)}{\mathrm{d}x}\,,\end{align} \begin{align}\tag{70} \frac{{\Delta}^2 v_{k - 1}\left( x \right)}{\tau^2} - \frac{1}{2} \gamma \left( \frac{\mathrm{d}^2 v_{k + 1}\left( x \right)}{\mathrm{d}x^2} + \frac{\mathrm{d}^2 v_{k - 1}\left( x \right)}{\mathrm{d}x^2} \right) + \frac{1}{2} \delta \left( v_{k + 1}\left( x \right) + v_{k - 1}\left( x \right) \right)\nonumber \\ = f_{2,k}\left( x \right) + a_2 \frac{\mathrm{d}u_k\left( x \right)}{\mathrm{d}x}\,, \end{align}\] where \(k = 1,2,\ldots,n - 1\), \(f_1\left( x,t_k \right) = f_{1,k}\left( x \right)\), and \(f_2\left( x,t_k \right) = f_{2,k}\left( x \right)\).

The nonlinear term in equation 69 is evaluated at the middle node point and denoted by \(q_k\), specifically: \[\label{eq:tim95spec95q95k95nonlin95term} q_k = \alpha + \beta \int\limits_{0}^{\ell} \left( \frac{\mathrm{d}u_k \left( x \right)}{\mathrm{d}x} \right)^2 \mathrm{d}x\,.\tag{71}\]

For each discrete time layer, the boundary conditions prescribed in 68 are rewritten as follows: \[\label{eq:tim95spec95bound95conds} u_{k + 1}\left( 0 \right) = 0\,,\quad u_{k + 1}\left( \ell \right) = 0\,,\quad v_{k + 1}\left( 0 \right) = 0\,,\quad v_{k + 1}\left( \ell \right) = 0\,.\tag{72}\]

The values of the unknown functions at the zeroth and first temporal layers are determined by the initial conditions specified in 67 and the system given by 65 66 as follows: \[\begin{gather} u_0\left( x \right) = \varphi_0\left( x \right)\,, \\ u_1\left( x \right) = \varphi_0\left( x \right) + \tau \varphi_1\left( x \right) + \frac{\tau^2}{2} \varphi_2\left( x \right)\,,\quad \varphi_2\left( x \right) = f_{1,0}\left( x \right) - a_1 \frac{\mathrm{d}\psi_0\left( x \right)}{\mathrm{d}x} + q_0 \frac{\mathrm{d}^2 \varphi_0\left( x \right)}{\mathrm{d}x^2}\,, \end{gather}\] and \[\begin{gather} v_0\left( x \right) = \psi_0\left( x \right)\,, \\ v_1\left( x \right) = \psi_0\left( x \right) + \tau \psi_1\left( x \right) + \frac{\tau^2}{2} \psi_2\left( x \right)\,,\quad \psi_2\left( x \right) = f_{2,0}\left( x \right) + a_2 \frac{\mathrm{d}\varphi_0\left( x \right)}{\mathrm{d}x} + \gamma \frac{\mathrm{d}^2 \psi_0\left( x \right)}{\mathrm{d}x^2} - \delta \psi_0\left( x \right)\,. \end{gather}\]

In the subsequent discussion, let \(u_k\left( x \right)\) and \(v_k\left( x \right)\) denote the solutions to the system of differential-difference equations 69 70 , subject to the boundary conditions 72 . Consequently, these solutions are declared as approximate values of the exact solution \(u\left( x,t \right)\) and \(v\left( x,t \right)\) of the problem 65 68 at the discrete time instants \(t = t_k\), respectively, so that \(u\left( x,t_k \right) \approx u_k\left( x \right)\) and \(v\left( x,t_k \right) \approx v_k\left( x \right)\).

Let us now define the notation for the first-order and second-order differential operators as follows: \[\begin{gather} {\mathcal{B}}_0 = \frac{\mathrm{d}}{\mathrm{d}x}\,,\quad D\left( {\mathcal{B}}_0 \right) = C^1\left( \left[ 0,\ell \right] \right)\,,\tag{73} \\ {\mathcal{A}}_0 = - \frac{\mathrm{d}^2}{\mathrm{d}x^2}\,,\quad D\left( {\mathcal{A}}_0 \right) = \left\{ u\left( x \right) \in C^2\left( \left[ 0,\ell \right] \right) \mid u\left( 0 \right) = u\left( \ell \right) = 0 \right\}\,.\tag{74} \end{gather}\] Employing the notations defined in 73 and 74 , the system 69 70 can be reformulated as follows: \[\begin{align} \tag{75} \frac{{\Delta}^2 u_{k - 1}\left( x \right)}{\tau^2} + q_k \frac{{\mathcal{A}}_0 u_{k + 1} + {\mathcal{A}}_0 u_{k - 1}}{2} = f_{1,k}\left( x \right) - a_1 {\mathcal{B}}_0 v_k\left( x \right)\,,\end{align} \begin{align} \tag{76} \frac{{\Delta}^2 v_{k - 1}\left( x \right)}{\tau^2} + \frac{{\mathcal{L}}_0 v_{k + 1} + {\mathcal{L}}_0 v_{k - 1}}{2} = f_{2,k}\left( x \right) + a_2 {\mathcal{B}}_0 u_k\left( x \right)\,,\end{align}\] for \(k = 1,2,\ldots,n - 1\), and \({\mathcal{L}}_0 = \gamma {\mathcal{A}}_0 + \delta {\mathcal{I}}\), where \({\mathcal{I}}\) represents the identity operator.

It follows easily from the application of integration by parts that the quantity \(q_k\) can be represented in the subsequent form, using the notation introduced in 74 , specifically: \[q_k = \alpha + \beta \left( {\mathcal{A}}_0 u_k,u_k \right)\,.\] In the following discussions, the notation \(\left( \cdot,\cdot \right)\) is used to denote the inner product in the space \(L^2\left( 0,\ell \right)\), whereas the corresponding norm is represented by \(\left\lVert \cdot \right\rVert\).

Since \({\mathcal{A}}_{0}\) is a symmetric and positive definite operator (see, e.g., Chap. 8 in Rektorys [30]), it admits an extension to a self-adjoint and positive definite operator (cf. Chap. V, Sec. 4 in Maurin [31]). We denote this operator by \(\mathcal{A}\) (\({\mathcal{A}}_{0} \subset {\mathcal{A}}\)). It is evident that the operator \({\mathcal{B}}_{0}\) is closable and can be extended to a closed operator. The closure of \({\mathcal{B}}_{0}\) is denoted by \(\mathcal{B}\) (\({\mathcal{B}}_{0} \subset {\mathcal{B}}\)). We shall rewrite equations 75 and 76 in terms of the operators \(\mathcal{A}\) and \(\mathcal{B}\) as follows: \[\begin{align} \tag{77} \frac{{\Delta}^2 u_{k - 1}\left( x \right)}{\tau^2} + q_k \frac{{\mathcal{A}} u_{k + 1} + {\mathcal{A}} u_{k - 1}}{2} = f_{1,k}\left( x \right) - a_1 {\mathcal{B}} v_k\left( x \right)\,,\end{align} \begin{align} \tag{78} \frac{{\Delta}^2 v_{k - 1}\left( x \right)}{\tau^2} + \frac{{\mathcal{L}} v_{k + 1} + {\mathcal{L}} v_{k - 1}}{2} = f_{2,k}\left( x \right) + a_2 {\mathcal{B}} u_k\left( x \right)\,,\quad {\mathcal{L}} = \gamma {\mathcal{A}} + \delta {\mathcal{I}}\,,\end{align}\] where \[q_k = \alpha + \beta {\left\lVert {\mathcal{A}}^{\frac{1}{2}} u_k \right\rVert}^2\,.\]

Our objective is to extend a theorem analogous to Theorem [8](#prop:main95thm95convrg){reference-type=“ref” reference=“prop:main95thm95convrg”} to the scheme 77 78 . To achieve this, it requires demonstrating that inequalities analogous to 4 and ?? are retained for the operators \(\mathcal{A}\) and \(\mathcal{B}\).

It is easy to obtain \[{\left\lVert {\mathcal{B}}_{0} u \right\rVert}^2 = \int\limits_{0}^{\ell} {\left[ u^{\prime}\left( x \right) \right]}^2 \mathrm{d}x = - \int\limits_{0}^{\ell} u^{\prime\prime}\left( x \right) u\left( x \right) \mathrm{d}x = \left( {\mathcal{A}}_{0} u,u \right)\,,\quad \forall u \in D\left( {\mathcal{A}}_{0} \right) \subset D\left( {\mathcal{B}}_{0} \right)\,.\] Hence, it follows that \[\label{eq:rmk95spec95semidiscrete95square95norm} {\left\lVert {\mathcal{B}} u \right\rVert}^2 = \left( {\mathcal{A}} u,u \right)\,,\quad \forall u \in D\left( {\mathcal{A}}_{0} \right) \subset D\left( {\mathcal{B}}_{0} \right)\,.\tag{79}\] Define \({\mathcal{A}} u = v\). In this notation, 79 takes the following form: \[\label{eq:rmk95spec95semidiscrete95inv95oper} {\left\lVert {\mathcal{B}} {\mathcal{A}}^{-1} v \right\rVert}^2 = \left( v,{\mathcal{A}}^{-1} v \right)\,,\quad \forall v \in {R}_{0} = {\mathcal{A}} D\left( {\mathcal{A}}_{0} \right)\,.\tag{80}\]

It is evident that \(R_0\) is dense in \(L^2\). Indeed, \(R_0\) represents the set of functions \(f\) for which the differential equation \(-u^{\prime \prime}\left( x \right) = f\left( x \right)\), subject to homogeneous boundary conditions, admits a unique classical solution. It is well established that, for such a simple differential equation, the solution \(u\) lies in \(D\left( {\mathcal{A}}_0 \right)\) whenever \(f \in C\left( \left[ 0,\ell \right] \right)\) (cf. Chap. II, Sec. 10 in book by Fučı́k and Kufner [32]). Moreover, the density of \(C\left( \left[ 0,\ell \right] \right)\) in \(L^2\left( 0,\ell \right)\) is well-known.

The relation 80 may now be extended over the entire space \(L^2\). Given that \(R_0\) is dense in \(L^2\), for each \(v \in L^2\), there exists a sequence \(v_n \in R_0\) such that \(v_n \to v\). Accordingly, by 80 , for the sequence \(v_n\) we have: \[\label{eq:rmk95spec95semidiscrete95inv95oper95vn} {\left\lVert {\mathcal{B}} {\mathcal{A}}^{-1} v_n \right\rVert}^2 = \left( v_n,{\mathcal{A}}^{-1} v_n \right)\,.\tag{81}\] Applying the Cauchy–Schwarz inequality to 81 , we obtain: \[{\left\lVert {\mathcal{B}} {\mathcal{A}}^{-1} v_n \right\rVert}^2 \leq \left\lVert v_n \right\rVert \left\lVert {\mathcal{A}}^{-1} v_n \right\rVert \leq \left\lVert {\mathcal{A}}^{-1} \right\rVert {\left\lVert v_n \right\rVert}^2\,.\] From this, it follows that \({\mathcal{B}} {\mathcal{A}}^{-1} v_n\) forms the Cauchy sequence. Taking into account that \({\mathcal{B}}\) is a closed operator and that \({\mathcal{A}}^{-1} v_n \to {\mathcal{A}}^{-1} v\), we conclude that \({\mathcal{A}}^{-1} v \in D\left( \mathcal{B} \right)\) and that \({\mathcal{B}} {\mathcal{A}}^{-1} v_n \to {\mathcal{B}} {\mathcal{A}}^{-1} v\). Considering this fact, it follows from 81 that: \[{\left\lVert {\mathcal{B}} {\mathcal{A}}^{-1} v \right\rVert}^2 = \left( v,{\mathcal{A}}^{-1} v \right)\,,\quad \forall v \in L^2\,,\] or, which is the same \[{\left\lVert {\mathcal{B}} u \right\rVert}^2 = \left( {\mathcal{A}} u,u \right)\,,\quad \forall u \in D\left( {\mathcal{A}} \right)\,.\]

Indeed, by an application of integration by parts, we obtain: \[\left( {\mathcal{A}}_0 {\mathcal{B}}_0 u, {\mathcal{B}}_0 u \right) = -\int_{0}^{\ell} u^{\prime \prime \prime} \left( x \right) u^{\prime} \left( x \right) \mathrm{d}x = \int_{0}^{\ell} {\left[ u^{\prime \prime} \left( x \right) \right]}^2 \mathrm{d}x = {\left\lVert {\mathcal{A}}_0 u \right\rVert}^2\,,\quad \forall u \in D_0\,.\] Thus, it is evident that ?? follows since the operators \(\mathcal{A}\) and \(\mathcal{B}\) are extensions of \({\mathcal{A}}_0\) and \({\mathcal{B}}_0\), respectively.

Due to Remark [6](#prop:remark4){reference-type=“ref” reference=“prop:remark4”}, the following relation follows from Remark [12](#prop:rmk95operat95ab){reference-type=“ref” reference=“prop:rmk95operat95ab”}: \[\left\lVert {\mathcal{A}}^{\frac{1}{2}} {\mathcal{B}} u \right\rVert = \left\lVert {\mathcal{A}} u \right\rVert\,,\quad \forall u \in D\left( {\mathcal{A}} \right)\,.\]

Based on Theorem [8](#prop:main95thm95convrg){reference-type=“ref” reference=“prop:main95thm95convrg”}, we now formulate a theorem concerning the convergence of the scheme 77 78 .

7.2 Legendre–Galerkin Spectral Method↩︎

It should be observed that solving the system 69 70 reduces to solving the following intermediate system: \[\begin{align} w_{1,k} \left( x \right) - \frac{\tau^2}{2} q_k \frac{\mathrm{d}^2 w_{1,k}}{\mathrm{d}x^2} &= \tau^2 f_{1,k}\left( x \right) + 2 u_k\left( x \right) - a_1 \tau^2 \frac{\mathrm{d}v_k}{\mathrm{d}x}\,,\tag{82} \\ \left( 1 + \frac{\tau^2}{2} \delta \right) w_{2,k} \left( x \right) - \frac{\tau^2}{2} \gamma \frac{\mathrm{d}^2 w_{2,k}}{\mathrm{d}x^2} &= \tau^2 f_{2,k}\left( x \right) + 2 v_k\left( x \right) + a_2 \tau^2 \frac{\mathrm{d}u_k}{\mathrm{d}x}\,,\tag{83} \\ w_{1,k} \left( 0 \right) = w_{1,k} \left( \ell \right) &= 0\,,\quad w_{2,k} \left( 0 \right) = w_{2,k} \left( \ell \right) = 0\,,\tag{84} \end{align}\] for which \(k = 1,2,\ldots,n - 1\).

By solving the system 82 83 , the unknowns \(u_k\left( x \right)\) and \(v_k\left( x \right)\) are determined by the following formulas: \[\label{eq:tim95spec95inv95operat95ref95uv} u_{k + 1}\left( x \right) = w_{1,k} \left( x \right) - u_{k - 1}\left( x \right)\,,\quad v_{k + 1}\left( x \right) = w_{2,k} \left( x \right) - v_{k - 1}\left( x \right)\,.\tag{85}\]

Before presenting the method in detail, we shall first introduce the shifted Legendre polynomials. These polynomials are obtained by an affine transformation that involves “shifting” and “scaling” the argument of the standard Legendre polynomials, mapping the interval \(\left[ 0,\ell \right]\) onto \(\left[ -1,1 \right]\). Specifically, the shifted Legendre polynomial of degree \(m\) is defined as \[\widetilde{P}_m \left( x \right) = P_m \left( \frac{2}{\ell} x - 1 \right)\,,\quad x \in \left[ 0,\ell \right]\,,\] where \(P_m \left( x \right)\) represents the \(m\)-th standard Legendre polynomial. In the following discussion, the shifted Legendre polynomial of degree \(m\) is denoted by \(\widetilde{P}_m \left( x \right)\).

Recall that the shifted Legendre polynomials retain an orthogonality property, similar to their unshifted counterparts, but adjusted to the interval \(\left[ 0,\ell \right]\), yielding a dependence on the interval length \(\ell\), namely: \[\label{eq:galerkin95orthog95prop} \left( \widetilde{P}_i,\widetilde{P}_m \right) = \ell A_i A_m \delta_{im}\,,\quad A_m = \frac{1}{\sqrt{2m + 1}}\,,\tag{86}\] where \(\delta_{im}\) denotes the Kronecker delta, taking the value \(1\) when \(i = m\) and \(0\) otherwise.

In order to approximate the solutions \(w_{1,k} \left( x \right)\) and \(w_{2,k} \left( x \right)\) of the system 82 84 at each temporal layer, we seek them in the form of the following linear combinations: \[\label{eq:galerkin95intermed95syst95ansatzes} {w}_{k,N}^{\left( j \right)} \left( x \right) = \sum_{m = 1}^{N} {w}_{j,m}^{k} \phi_m \left( x \right)\,,\quad j = 1,2\,.\tag{87}\] Here, we choose the differences of the shifted Legendre polynomials as ansatz (trial) functions: \[\phi_m \left( x \right) = \frac{1}{A_m \sqrt{\ell}} \int\limits_{0}^{x} \widetilde{P}_m \left( s \right) \mathrm{d}s = \frac{\sqrt{\ell}}{2} A_m \left( \widetilde{P}_{m + 1} \left( x \right) - \widetilde{P}_{m - 1} \left( x \right) \right)\,.\]

Observe that the functions \(\phi_m \left( x \right)\) inherently vanish at the endpoints of the interval \(\left[ 0,\ell \right]\). Furthermore, it is obvious that: \[\phi_{m}^{\prime} \left( x \right) = \widehat{P}_m \left( x \right)\,,\quad \widehat{P}_m \left( x \right) = \frac{1}{A_m \sqrt{\ell}} \widetilde{P}_m \left( x \right)\,.\] Evidently, \(\widehat{P}_m \left( x \right)\) denotes the orthonormal shifted Legendre polynomial on the interval \(\left[ 0,\ell \right]\).

To approximate the solutions \(u_k\left( x \right)\) and \(v_k\left( x \right)\) of the system 69 70 at each temporal layer, we employ the same linear combination as in 87 , that is, \[\label{eq:galerkin95ansatzes} \widetilde{u}_{k,N} \left( x \right) = \sum_{m = 1}^{N} \tilde{u}_{m}^{k} \phi_m \left( x \right)\,,\quad \widetilde{v}_{k,N} \left( x \right) = \sum_{m = 1}^{N} \tilde{v}_{m}^{k} \phi_m \left( x \right)\,,\quad k = 2,3,\ldots,n\,.\tag{88}\]

Substituting the Galerkin approximations defined by 87 and 88 into the relations provided in 85 , we derive the following expressions for the expansion coefficients: \[\tilde{u}_{m}^{k + 1} = {w}_{1,m}^{k} - \tilde{u}_{m}^{k - 1} \quad \text{and} \quad \tilde{v}_{m}^{k + 1} = {w}_{2,m}^{k} - \tilde{v}_{m}^{k - 1}\,.\]

Let the test functions be chosen to coincide with the trial functions \(\phi_m \left( x \right)\). The derivation of the Galerkin system involves substituting the Galerkin approximations from 87 into equations 82 and 83 , followed by taking the inner product of both sides of the resulting equations with the test functions \(\phi_m \left( x \right)\) for \(m = 1,2,\ldots,N\). Finally, the resulting Galerkin subsystems can be reformulated in the following matrix-vector form: \[\begin{align} \boldsymbol{\T}_{N}^{k} \boldsymbol{w}_{1}^{k} &= \frac{4 \tau^2}{\ell^2} \boldsymbol{I}_{k}^{\left( 1 \right)} + 2 \boldsymbol{\mathcal{H}}_N \boldsymbol{\tilde{u}}^k - \frac{2 a_1 \tau^2}{\ell} \boldsymbol{\mathcal{B}}_N \boldsymbol{\tilde{v}}^k\,,\tag{89} \\ \boldsymbol{\T}_{N} \boldsymbol{w}_{2}^{k} &= \frac{2 a_0 \tau^2}{\ell^2} \boldsymbol{I}_{k}^{\left( 2 \right)} + a_0 \boldsymbol{\mathcal{H}}_N \boldsymbol{\tilde{v}}^k + \frac{a_0 a_2 \tau^2}{\ell} \boldsymbol{\mathcal{B}}_N \boldsymbol{\tilde{u}}^k\,,\quad a_0 = \frac{4}{2 + \delta \tau^2}\,.\tag{90} \end{align}\]

The coefficient matrices corresponding to the subsystems 89 and 90 are defined by \[\boldsymbol{\T}_{N}^{k} = \boldsymbol{\mathcal{H}}_N + \frac{2 \tau^2}{\ell^2} q_{k,N} \boldsymbol{\mathcal{I}}_N \quad \text{and} \quad \boldsymbol{\T}_{N} = \boldsymbol{\mathcal{H}}_N + \frac{a_0 \tau^2}{\ell^2} \gamma \boldsymbol{\mathcal{I}}_N\,,\] for which \[q_{k,N} = \alpha + \beta \sum_{m = 1}^{N} \left( \tilde{u}_{m}^{k} \right)^2\,.\]

The matrix \(\boldsymbol{\mathcal{H}}_N = \left( h_{i,j} \right)_{1 \leq i,j \leq N}\) is symmetric and sparse. Its entries are defined as follows: the main diagonal entries are given by \(h_{i,i} = 2 A_{i - 1}^{2} A_{i + 1}^{2}\), for \(i = 1,2,\ldots,N\). The nonzero entries on the second sub- and super-diagonals are \(h_{i,i + 2} = h_{i + 2,i} = - A_{i} A_{i + 1}^{2} A_{i + 2}\), for \(i = 1,2,\ldots,N - 2\). All other entries of this matrix are zero. The matrix \(\boldsymbol{\mathcal{B}}_N = \left( b_{i,j} \right)_{1 \leq i,j \leq N}\) is sparse and skew-symmetric, with nonzero entries confined to the first sub- and super-diagonals. Its off-diagonal entries satisfy \(b_{i,i + 1} = A_i A_{i + 1}\) and \(b_{i + 1,i} = -A_i A_{i + 1}\), for \(i = 1,2,\ldots,N - 1\), while all other entries, including those on the main diagonal, are equal to zero. Moreover, \(\boldsymbol{\mathcal{I}}_N\) denotes the identity matrix of order \(N\).

The column vectors \(\boldsymbol{\tilde{u}}^k\), \(\boldsymbol{\tilde{v}}^k\), \(\boldsymbol{w}_{1}^{k}\), and \(\boldsymbol{w}_{2}^{k}\), appearing in the subsystems of the Galerkin linear equations 89 and 90 , consist of the expansion coefficients given in 88 and 87 , respectively.

For each \(j = 1,2\), the vector \[\boldsymbol{I}_{k}^{\left( j \right)} = {\left( I_{k,1}^{\left( j \right)},I_{k,2}^{\left( j \right)},\ldots,I_{k,N}^{\left( j \right)} \right)}^{\top}\] is defined by \(I_{k,m}^{\left( j \right)} = \left( f_{j,k},\phi_m \right)\). Thus, the components \(I_{k,m}^{\left( j \right)}\) are given by the inner products of the functions \(f_{j,k}\left( x \right)\) with the basis functions \(\phi_m\left( x \right)\).

This Theorem follows from the subsequent Lemma.

Let us introduce the vector: \[\boldsymbol{v} = {\left( c_1,c_2,\ldots,c_N \right)}^{\top}\,.\] Here, \(c_1,c_2,\ldots,c_N\) are the coefficients of the Galerkin approximation, corresponding to the basis functions \(\phi_1,\phi_2,\ldots,\phi_N\). This vector constitutes the solution to the Galerkin system of linear equations within the finite-dimensional subspace spanned by the chosen \(N\) basis functions.

It follows directly that: \[\boldsymbol{\mathcal{S}}_N \boldsymbol{v} = {\left( \left( A u_N,\phi_1 \right)_{\mathcal{H}},\left( A u_N,\phi_2 \right)_{\mathcal{H}},\ldots,\left( A u_N,\phi_N \right)_{\mathcal{H}} \right)}^{\top}\,,\] where \[\label{eq:galerkin95lemma95gal95approx} u_N = \sum_{i = 1}^{N} c_i \phi_i\,.\tag{91}\] In fact, we have: \[\label{eq:galerkin95lemma95inner95prod95oper} \left( A u_N,\phi_m \right)_{\mathcal{H}} = \sum_{i = 1}^{N} c_i \left( A \phi_i,\phi_m \right)_{\mathcal{H}}\,,\quad m = 1,2,\ldots,N\,.\tag{92}\] As a result of 92 , we obtain: \[\label{eq:galerkin95lemma95equality95inn95prod} \left( \boldsymbol{\mathcal{S}}_N \boldsymbol{v},\boldsymbol{v} \right) = \sum_{i = 1}^{N} c_i \left( A u_N,\phi_i \right)_{\mathcal{H}} = \left( A u_N,\sum_{i = 1}^{N} c_i \phi_i \right)_{\mathcal{H}} = \left( A u_N,u_N \right)_{\mathcal{H}}\,.\tag{93}\]

By combining relations ?? and 93 , it follows that: \[\label{eq:galerkin95lemma95finitedim95condition} \left( \boldsymbol{\mathcal{S}}_N \boldsymbol{v},\boldsymbol{v} \right) \geq \alpha \left( B u_N,u_N \right)_{\mathcal{H}} + \nu \left\lVert u_N \right\rVert_{\mathcal{H}}^2\,.\tag{94}\] Substituting 91 into inequality 94 and leveraging the \(B\)-orthogonality property ?? , we arrive at: \[\begin{align} \left( \boldsymbol{\mathcal{S}}_N \boldsymbol{v},\boldsymbol{v} \right) &\geq \alpha \left( \sum_{i = 1}^{N} c_i B \phi_i,\sum_{m = 1}^{N} c_m \phi_m \right)_{\mathcal{H}} + \nu \left\lVert u_N \right\rVert_{\mathcal{H}}^2 \\ &\geq \alpha \sum_{i = 1}^{N} \sum_{m = 1}^{N} c_i c_m \left( B \phi_i,\phi_m \right)_{\mathcal{H}} = \alpha \sum_{m = 1}^{N} b_m c_m^2 \geq \widehat{\alpha} \left\lVert \boldsymbol{v} \right\rVert_{2}^2\,,\quad \widehat{\alpha} = \alpha \min_{1 \leq m \leq N} b_m\,. \end{align}\] Finally, from the inequality stated above, it follows that the matrix \(\boldsymbol{\mathcal{S}}_N\) is positive definite.

7.3 Some Notes on Solving the Derived Subsystems of the Galerkin Linear Equations↩︎

In this section, our objective is to establish an efficient approach for addressing the subsystems of the linear equations 89 and 90 that originate from the Legendre–Galerkin spectral approximation. Specifically, we focus on solving the subsystems of Galerkin linear equations at each temporal step. For simplicity, the temporal layer index \(k\) is suppressed. Each subsystem is subsequently represented in matrix-vector form as follows: \[\label{eq:galerkin95eff95alg95unif95system} \boldsymbol{\mathcal{A}}_{N} \boldsymbol{w} = \boldsymbol{f}\,,\tag{95}\] where the coefficient matrix of the system is expressed as: \[\boldsymbol{\mathcal{A}}_{N} = \boldsymbol{\mathcal{H}}_N + \frac{4}{a \ell^2} b \boldsymbol{\mathcal{I}}_N\,,\] and the unknown vector \(\boldsymbol{w}\), which has \(N\) components, is defined as: \[\boldsymbol{w} = {\left( w_1,w_2,\ldots,w_N \right)}^{\top}\,.\] Furthermore, the vector \(\boldsymbol{f} = {\left( f_1,f_2,\ldots,f_N \right)}^{\top}\) assumes values as specified in 89 and 90 , depending on the subsystem being solved. For subsystem 89 , the coefficients \(\left( a,b \right)\) are assigned the values \(\left( 1,\tau^2 q_{k,N} / 2 \right)\), whereas for subsystem 90 , the coefficients are \(\left( 1 + \tau^2 \delta / 2,\tau^2 \gamma / 2 \right)\).

The matrix \(\boldsymbol{\mathcal{A}}_{N} \in \mathbb{R}^{N \times N}\) is distinguished by a sparsity pattern in which nonzero entries appear exclusively on the main diagonal, the second sub-diagonal, and the second super-diagonal (cf. Subsec. 7.2). Moreover, \(\boldsymbol{\mathcal{A}}_{N}\) is symmetric, and as shown in Theorem [14](#prop:thm95galerkin95positive95def95syst){reference-type=“ref” reference=“prop:thm95galerkin95positive95def95syst”}, it is positive-definite. We shall occasionally refer to this matrix as a tridiagonal matrix with a gap, where the “gap” refers to the two skipped diagonals: one below and one above the main diagonal.

We propose the following method to solve the system of Galerkin linear equations 95 . Leveraging the structure of the matrix \(\boldsymbol{\mathcal{A}}_{N}\), the system can be decomposed into two independent subsystems, as described below.

If the number of basis functions is even, i.e., \(N = 2s\), where \(s \in \mathbb{N}\), the unknowns are enumerated in the following manner: \[w_{2j - 1} = \widetilde{w}_j \quad \text{and}\quad w_{2j} = \widehat{w}_j\,,\quad \text{for}\,\, j = 1,2,\ldots,s\,.\]

Similarly, when the number of basis functions is odd, that is, \(N = 2s - 1\) with \(s \in \mathbb{N}\), the unknowns are reordered as follows: \[\begin{align} w_{2j - 1} &= \widetilde{w}_j\,,\quad \text{for}\,\, j = 1,2,\ldots,s\,, \\ w_{2j} &= \widehat{w}_j\,,\quad \text{for}\,\, j = 1,2,\ldots,s - 1\,. \end{align}\]

The matrices associated with the two independent subsystems, resulting from the considered decomposition, are tridiagonal. These systems can be solved in parallel for the unknowns \(\widetilde{w}_j\) and \(\widehat{w}_j\) using the Thomas algorithm (tridiagonal matrix algorithm). Although this algorithm is generally not stable, stability is assured when the matrix is diagonally dominant (either row- or column-wise) or symmetric positive definite (see Theorem 9.12 in Higham [33]). It is straightforward to establish that, in this case, the matrices corresponding to both decomposed tridiagonal subsystems are positive definite. In order to derive closed-form expressions for the unknowns \(\widetilde{w}_j\) and \(\widehat{w}_j\), one may apply the Cholesky decomposition to the coefficient matrix associated with the linear system 95 .

Consider a variant of the classical Cholesky decomposition, referred to as the square-root-free Cholesky decomposition, applied to the real symmetric positive-definite matrix \(\boldsymbol{\mathcal{A}}_{N}\), which corresponds to the coefficient matrix of the system of linear equations 95 , and is formulated as follows: \[\boldsymbol{\mathcal{A}}_{N} = \boldsymbol{L} \boldsymbol{D} {\boldsymbol{L}}^{\top}\,.\] In the case of the square-root-free Cholesky decomposition of a tridiagonal matrix with a gap, the matrix \(\boldsymbol{L}\) is defined such that all entries along its main diagonal are equal to unity, with nonzero off-diagonal entries appearing only on the second sub-diagonal. The matrix \(\boldsymbol{D}\) is diagonal, and \({\boldsymbol{L}}^{\top}\) denotes the transpose of \(\boldsymbol{L}\).

It should be observed that this factorization not only provides closed-form solutions to the system 95 , but also permits the parallel computation of solutions corresponding to the unknowns with odd and even indices.

Under this decomposition, the system 95 takes the following form: \[\boldsymbol{L} \boldsymbol{D} {\boldsymbol{L}}^{\top} \boldsymbol{w} = \boldsymbol{f}\,.\] This system evidently allows for decomposition into the ensuing subsystems by introducing auxiliary unknowns \(\boldsymbol{y}\) and \(\boldsymbol{z}\), as follows: \[\label{eq:galerkin95cholesky95decom} \begin{cases} \boldsymbol{L} \boldsymbol{z} &= \boldsymbol{f}\,, \\ \boldsymbol{D} \boldsymbol{y} &= \boldsymbol{z}\,, \\ {\boldsymbol{L}}^{\top} \boldsymbol{w} &= \boldsymbol{y}\,. \end{cases}\tag{96}\] As in the preceding case, two distinct scenarios should be considered, depending on whether \(N\), the number of basis functions, is even or odd, which in turn determines the dimension of the associated matrix. Thus, the problem of solving system 95 is reduced to that of solving system 96 , for which the corresponding closed-form solutions are subsequently discussed.

Consider the case where \(N\) is even, i.e., \(N = 2s\) with \(s \in \mathbb{N}\). In this case, the solution to the system 96 can be expressed as follows: \[\begin{align} z_1 = f_1\,, \qquad z_2 &= f_2\,, \\ d_1 = C_1 + \frac{4}{a \ell^2} b\,,\qquad d_2 &= C_2 + \frac{4}{a \ell^2} b\,. \end{align}\] For \(j = 2,3,\ldots,s\), the following holds: \[\begin{align} z_{2j - 1} = f_{2j - 1} + \frac{B_{2 \left( j - 1 \right)}}{d_{2j - 3}} z_{2j - 3}\,,\qquad z_{2j} &= f_{2j} + \frac{B_{2j - 1}}{d_{2 \left( j - 1 \right)}} z_{2 \left( j - 1 \right)}\,, \\ d_{2j - 1} = \left( C_{2j - 1} + \frac{4}{a \ell^2} b \right) - \frac{B_{2 \left( j - 1 \right)}^2}{d_{2j - 3}}\,,\qquad d_{2j} &= \left( C_{2j} + \frac{4}{a \ell^2} b \right) - \frac{B_{2j - 1}^2}{d_{2 \left( j - 1 \right)}}\,. \end{align}\] The vector \(\boldsymbol{y}\) is computed component-wise by multiplying each entry of the vector \(\boldsymbol{z}\) by the reciprocal of the corresponding diagonal element of the matrix \(\boldsymbol{D}\). Therefore, for \(j = 1,2,\ldots,s\), we have: \[y_{2j - 1} = \frac{z_{2j - 1}}{d_{2j - 1}}\,,\qquad y_{2j} = \frac{z_{2j}}{d_{2j}}\,.\] Observing that \(w_{2s} = y_{2s}\) and \(w_{2s - 1} = y_{2s - 1}\), the solution to the system 95 is subsequently determined for \(j = s - 1,s - 2,\ldots,1\) and is represented by: \[w_{2j - 1} = y_{2j - 1} + \frac{B_{2j}}{d_{2j - 1}} w_{2j + 1}\,,\qquad w_{2j} = y_{2j} + \frac{B_{2j + 1}}{d_{2j}} w_{2 \left( j + 1 \right)}\,.\]

Let \(N\) be an odd number, so that \(N = 2s - 1\), where \(s \in \mathbb{N}\). Even in this setting, the same formulas for determining the components of the unknown vectors \(\boldsymbol{z}\), \(\boldsymbol{y}\) and \(\boldsymbol{w}\), as well as the diagonal matrix \(\boldsymbol{D} = \operatorname{diag}\!\left( d_1,\ldots,d_N \right)\), continue to be valid. The only distinction occurs in the computation of even-indexed components within a for-loop. Specifically, for \(j = 2,3,\ldots,s - 1\), the components \(z_{2j}\) and \(d_{2j}\) are determined using the closed-form expressions derived previously. Likewise, for \(j = 1,2,\ldots,s - 1\), the components \(y_{2j}\) are defined as in the prior case. Finally, noting that \(w_{2 \left( s - 1 \right)} = d_{2 \left( s - 1 \right)}\), the remaining components \(w_{2j}\) for \(j = s - 2, s - 3,\ldots,1\) are computed recursively by applying the formula established in the preceding case.

Thus, as shown above, the Galerkin system is solved explicitly, and the corresponding computational cost is optimal, namely of order \(\bigO \left( N \right)\).

If sine functions are chosen as both trial and test functions, then the matrix of the resulting Galerkin system is diagonal. However, in order to form the right-hand side of this system, one must multiply the matrix associated with the coupling terms by the relevant coefficient vectors. Since this matrix is almost dense, this requires \(\bigO \left( N^2 \right)\) operations. Hence, the computation of the solution at each time layer requires \(\bigO \left( N^2 \right)\) arithmetic operations.

Consequently, the Galerkin system can be solved at each time layer with \(\bigO \left( N \right)\) operations for the Legendre–Galerkin spectral method, whereas the corresponding sine-based formulation requires \(\bigO \left( N^2 \right)\) operations. Since we tackle a dynamic problem, the total number of arithmetic operations is \(\bigO \left( n \times N \right)\) and \(\bigO \left( n \times N^2 \right)\), respectively, where \(n\) denotes the number of subintervals in the partition of the interval \(\left[ 0,T \right]\). Therefore, for relatively large values of \(n\), the difference between the computational costs becomes noticeable. Thus, in this sense, the Legendre–Galerkin spectral method has a computational advantage in the present setting.

In addition to the above, in order to implement the three-layer semi-discrete scheme 69 –70 efficiently, it is necessary to obtain a simple recurrence relation between the Galerkin expansion coefficients of the corresponding approximate solutions. This is achieved by the proposed method (see Subsection 7.2). If, instead of employing the Legendre–Galerkin spectral method, one represents the solution of problem 69 –70 , subject to the homogeneous boundary conditions 72 , by means of Green’s formula, then the derivation of a recurrence relation similar to the one mentioned above becomes substantially more involved. Consequently, efficient computation is no longer feasible.

At the end of this subsection, we would like to draw attention to one nuance that illustrates the advantages of the proposed method. Suppose that the data of the discrete problem are continuous functions. Clearly, these functions can be approximated by Bernstein polynomials. In this case, the solution of the discrete problem will be a polynomial, since our scheme is locally linear. Moreover, this solution can be constructed exactly by means of the Legendre–Galerkin spectral method. Since the discrete problem is well posed, the polynomial solution provides an approximation to the exact solution of the same order of accuracy as the approximation of the prescribed data.

7.4 Estimate of the Error in the Legendre–Galerkin Spectral Method↩︎

In this subsection, we shall estimate the error incurred by the Legendre–Galerkin spectral method when applied to equations 82 and 83 . For this purpose, we reformulate these equations in the following manner: \[\label{eq:ritz95gal95main} a w \left( x \right) - b \frac{\mathrm{d}^2 w \left( x \right)}{\mathrm{d}x^2} = f \left( x \right)\,,\quad x \in \left[ 0,\ell \right]\,.\tag{97}\] Here, the pair of coefficients \(\left( a,b \right)\) is defined as \(\left( a,b \right) = \left( 1,\tau^2 q_{k,N} / 2 \right)\) for equation 82 and \(\left( a,b \right) = \left( 1 + \tau^2 \delta / 2,\tau^2 \gamma / 2 \right)\) for equation 83 . In equation 97 , the right-hand side \(f \left( x \right)\) is a continuous function, coinciding with the right-hand side of either equation 82 or equation 83 , depending on the context under consideration.

Equation 97 is considered subject to homogeneous boundary conditions, which are specified as follows: \[\label{eq:ritz95gal95bound95cond} w \left( 0 \right) = w \left( \ell \right) = 0\,.\tag{98}\]

The variational formulation of the problem 97 98 is well-established in the literature (cf., e.g., Kantorovich and Krylov [34]). The problem 97 98 is thus equivalent to the problem of minimizing the following functional: \[\label{eq:ritz95gal95variat95form} I \left( w \right) = \int\limits_{0}^{\ell} \left( b w^{\prime 2} + a w^2 - 2 f w \right) \mathrm{d}x\,,\quad w \left( 0 \right) = w \left( \ell \right) = 0\,.\tag{99}\]

The subsequent discussion follows the proof technique outlined in the aforementioned book (see Chap. IV, Sec. 4 in [34]), which concerns error estimation in the variational method for ordinary differential equations, with sines taken as the trial (basis) functions.

Let \(\widetilde{w} \left( x \right)\) be any function that satisfies the homogeneous boundary conditions \(\widetilde{w} \left( 0 \right) = \widetilde{w} \left( \ell \right) = 0\). Define the function \(\eta \left( x \right) = \widetilde{w} \left( x \right) - w \left( x \right)\). It follows immediately that \(\eta \left( 0 \right) = \eta \left( \ell \right) = 0\).

From 99 , a straightforward transformation yields: \[\begin{align} I \left( \widetilde{w} \right) - I \left( w \right) &= I \left( w + \eta \right) - I \left( w \right) \\ &= 2 \int\limits_{0}^{\ell} \left( b w^{\prime} \eta^{\prime} + a w \eta - f \eta \right) \mathrm{d}x + \int\limits_{0}^{\ell} \left( b \eta^{\prime 2} + a {\eta}^2 \right) \mathrm{d}x\,. \end{align}\] In the resulting expression, the first summand is simply the variation of the integral \(I\), which vanishes, i.e., \(\delta I = 0\). Consequently, we obtain: \[\label{eq:ritz95gal95diff95var} I \left( \widetilde{w} \right) - I \left( w \right) = \int\limits_{0}^{\ell} \left( b \left( \widetilde{w}^{\prime} - w^{\prime} \right)^{2} + a \left( \widetilde{w} - w \right)^2 \right) \mathrm{d}x\,.\tag{100}\]

Recall that we have chosen the following system of functions as the basis functions: \[\phi_m \left( x \right) = \int\limits_{0}^{x} \widehat{P}_m \left( s \right) \mathrm{d}s\,,\quad m = 1,2,\ldots\,,\] where \(\widehat{P}_m \left( x \right)\) represents the orthonormal shifted Legendre polynomial defined on the interval \(\left[ 0,\ell \right]\).

It is known that if \(u \left( x \right) \in C^p \left( \left[ 0,\ell \right] \right)\), with \(p \geq 1\), then the following bound holds (cf. Theorem 4.10 in Suetin [35]): \[\label{eq:ritz95gal95suetin} \left\lvert u \left( x \right) - \sum_{i = 0}^{N} \hat{c}_i \widehat{P}_i \left( x \right) \right\rvert \leq \frac{\lgsuetin}{N^{p - \frac{1}{2}}}\,,\quad x \in \left[ 0,\ell \right]\,.\tag{101}\] Here, \(\hat{c}_i\) denotes the coefficient in the Fourier–Legendre series expansion.

Assume that the solution to the problem 97 98 is such that \(w \left( x \right) \in C^2 \left( \left[ 0,\ell \right] \right)\). From estimate 101 , we have: \[\label{eq:ritz95gal95w95diff} \left\lvert w^{\prime} \left( x \right) - S_N \left( x \right) \right\rvert \leq \frac{\lgdiffw}{\sqrt{N}}\,,\quad x \in \left[ 0,\ell \right]\,,\tag{102}\] where \[S_N \left( x \right) = \sum_{i = 0}^{N} \hat{a}_i \widehat{P}_i \left( x \right)\,,\quad \hat{a}_i = \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_i \left( x \right) \mathrm{d}x\,.\] Observe that \[\hat{a}_0 = \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_0 \left( x \right) \mathrm{d}x = \frac{1}{\sqrt{\ell}} \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \mathrm{d}x = \frac{1}{\sqrt{\ell}} \left( w \left( \ell \right) - w \left( 0 \right) \right) = 0\,.\]

Consider the following combination: \[\Phi_N \left( x \right) = \sum_{m = 1}^{N} \hat{a}_m \phi_m \left( x \right) = \sum_{m = 1}^{N} \hat{a}_m \int\limits_{0}^{x} \widehat{P}_m \left( s \right) \mathrm{d}s\,.\] Taking into account that \(\Phi_{N}^{\prime} \left( x \right) = S_N \left( x \right)\), it follows from inequality 102 that: \[\begin{align} \left\lvert w \left( x \right) - \Phi_N \left( x \right) \right\rvert &= \left\lvert \int\limits_{0}^{x} \left( w^{\prime} \left( s \right) - \Phi_{N}^{\prime} \left( s \right) \right) \mathrm{d}s \right\rvert = \left\lvert \int\limits_{0}^{x} \left( w^{\prime} \left( s \right) - S_N \left( s \right) \right) \mathrm{d}s \right\rvert \\ &\leq \int\limits_{0}^{\ell} \left\lvert w^{\prime} \left( x \right) - S_N \left( x \right) \right\rvert \mathrm{d}x \leq \frac{\lgdiffw \ell}{\sqrt{N}}\,. \end{align}\] Thus, we have established the completeness of the system of functions \(\left\{ \phi_m \left( x \right) \right\}_{m = 1}^{\infty}\).

Consider the following approximation: \[w_N \left( x \right) = \sum_{m = 1}^{N} a_m \phi_m \left( x \right)\,,\] where \(a_m\) (\(m = 1,2,\ldots,N\)) represents the solution to the Ritz–Galerkin system.

It is important to note that, in our case, the systems of linear equations emerging from the Ritz–Galerkin and Legendre–Galerkin spectral methods coincide.

On the linear hull spanned by the system of functions \(\phi_1 \left( x \right), \phi_2 \left( x \right),\ldots,\phi_N \left( x \right)\), the functional \(I \left( w \right)\) attains its minimum at \(w \left( x \right) = w_N \left( x \right)\). From this, the following inequality holds: \[I \left( w_N \right) - I \left( w \right) \leq I \left( \Phi_N \right) - I \left( w \right)\,.\] By virtue of 100 , this inequality takes the following form: \[\label{eq:ritz95gal95unif95integ95ineq} \int\limits_{0}^{\ell} \left[ b \left( {w}_{N}^{\prime} - w^{\prime} \right)^{2} + a \left( {w}_{N} - w \right)^2 \right] \mathrm{d}x \leq \int\limits_{0}^{\ell} \left[ b \left( {\Phi}_{N}^{\prime} - w^{\prime} \right)^{2} + a \left( {\Phi}_{N} - w \right)^2 \right] \mathrm{d}x\,.\tag{103}\] It follows immediately from 103 that the following result holds: \[a \int\limits_{0}^{\ell} \left( {w}_{N} - w \right)^2 \mathrm{d}x \leq \int\limits_{0}^{\ell} \left[ b \left( {\Phi}_{N}^{\prime} - w^{\prime} \right)^{2} + a \left( {\Phi}_{N} - w \right)^2 \right] \mathrm{d}x\,.\] or, equivalently: \[\label{eq:ritz95gal95norm95error} a \left\lVert {w}_{N} - w \right\rVert^2 \leq b \left\lVert {\Phi}_{N}^{\prime} - w^{\prime} \right\rVert^2 + a \left\lVert {\Phi}_{N} - w \right\rVert^2\,.\tag{104}\]

We shall now apply the classical Steklov inequality, which states: \[\label{eq:ritz95gal95func95bound95diff} \left\lVert w \right\rVert \leq \frac{\ell}{\pi} \left\lVert w^{\prime} \right\rVert\,.\tag{105}\] Hence, in view of inequality 105 , it follows that \[\label{eq:ritz95gal95diff95ineq95der} \left\lVert \Phi_N - w \right\rVert \leq \frac{\ell}{\pi} \left\lVert \Phi_{N}^{\prime} - w^{\prime} \right\rVert\,.\tag{106}\] Substituting inequality 106 into 104 , we find that: \[\label{eq:ritz95gal95norm95error95refined} \left\lVert {w}_{N} - w \right\rVert \leq b_1 \left\lVert {\Phi}_{N}^{\prime} - w^{\prime} \right\rVert\,,\quad b_1 = \sqrt{\frac{b}{a} + \frac{\ell^2}{\pi^2}}\,.\tag{107}\] Taking into account that \(\Phi_{N}^{\prime} \left( x \right) = S_N \left( x \right)\), the following inequality is obtained from 107 : \[\label{eq:ritz95gal95norm95error95s} \left\lVert {w}_{N} - w \right\rVert \leq b_1 \left\lVert S_N - w^{\prime} \right\rVert\,.\tag{108}\]

Let us assume that the solution to problem 97 98 , \(w \left( x \right) \in C^{p + 1} \left( \left[ 0,\ell \right] \right)\), with \(p \geq 1\). Then, by virtue of bound 101 , inequality 108 leads to the following estimate: \[\label{eq:ritz95gal95estimate95final} \left\lVert {w}_{N} - w \right\rVert \leq \frac{\lgerrorgal}{N^{p - \frac{1}{2}}}\,.\tag{109}\]

Thus, the following theorem is valid.

One should note that the error \(w \left( x \right) - w_N \left( x \right)\) can also be estimated using the uniform norm. To achieve this, additional transformations are required.

It is evident from inequality 103 that the following holds: \[b \left\lVert {w}_{N}^{\prime} - w^{\prime} \right\rVert^2 \leq b \left\lVert {\Phi}_{N}^{\prime} - w^{\prime} \right\rVert^2 + a \left\lVert {\Phi}_{N} - w \right\rVert^2\,.\] From this, taking into account 106 , we deduce: \[\left\lVert {w}_{N}^{\prime} - w^{\prime} \right\rVert^2 \leq \left( 1 + \frac{a {\ell}^2}{b \pi^2} \right) \left\lVert {\Phi}_{N}^{\prime} - w^{\prime} \right\rVert^2\,,\] or, equivalently: \[\label{eq:ritz95gal95unif95diff95error} \left\lVert {w}_{N}^{\prime} - w^{\prime} \right\rVert \leq b_2 \left\lVert S_N - w^{\prime} \right\rVert\,,\quad b_2 = \sqrt{1 + \frac{a {\ell}^2}{b \pi^2}}\,.\tag{110}\]

From the representation \[w \left( x \right) = \int\limits_{0}^{x} w^{\prime} \left( s \right) \mathrm{d}s\,,\] and by standard reasoning, one obtains (cf. Chap. IV, Sec. 4 in Kantorovich and Krylov [34]): \[\max_{0 \leq x \leq \ell} \left\lvert w \left( x \right) \right\rvert \leq \sqrt{\frac{\ell}{2}} \left\lVert w^{\prime} \right\rVert\,.\] By this inequality, we have: \[\max_{0 \leq x \leq \ell} \left\lvert w_N \left( x \right) - w \left( x \right) \right\rvert \leq \sqrt{\frac{\ell}{2}} \left\lVert {w}_{N}^{\prime} - w^{\prime} \right\rVert\,.\] From this, taking 110 into consideration, it follows that: \[\label{eq:ritz95gal95error95unif95norm} \max_{0 \leq x \leq \ell} \left\lvert w_N \left( x \right) - w \left( x \right) \right\rvert \leq b_2 \sqrt{\frac{\ell}{2}} \left\lVert S_N - w^{\prime} \right\rVert\,.\tag{111}\]

If \(w \left( x \right) \in C^{p + 1} \left( \left[ 0,\ell \right] \right)\) with \(p \geq 1\), then from 111 , and in accordance with 101 , the following inequality is obtained: \[\label{eq:ritz95gal95unif95error95final} \max_{0 \leq x \leq \ell} \left\lvert w_N \left( x \right) - w \left( x \right) \right\rvert \leq \sqrt{\frac{\ell}{2}}\frac{b_2 \lgerrunifnorm}{N^{p - \frac{1}{2}}}\,.\tag{112}\]

It is observed that by substituting the values of \(a\) and \(b\) corresponding to equations 82 and 83 into the expression for \(b_2\), the inequality \(b_2 \leq \lgboundbtwo / \tau\), with \(\tau = T / n\), follows immediately. If the number of divisions \(n\) is chosen as the integer part of \(N^s\) (\(n = \left[ N^s \right]\)), where \(p - 1/2 - s >0\), then bound 112 takes the following form: \[\max_{0 \leq x \leq \ell} \left\lvert w_N \left( x \right) - w \left( x \right) \right\rvert \leq \frac{\lgfinalest}{N^{p - \frac{1}{2} - s}}\,,\] where \(p > 1 / 2 + s\), \(p \geq 1\) is a positive integer, and \(s > 0\).

In our opinion, the reader will be interested in an explicit error estimate with a uniform norm in the natural class of solutions for equation ?? .

Let \(\hat{c}_m\) denote the coefficients of the Legendre–Fourier series expansion of the function \(w^{\prime \prime} \left( x \right)\), given by \[\hat{c}_m = \int\limits_{0}^{\ell} w^{\prime \prime} \left( x \right) \widehat{P}_m \left( x \right) \mathrm{d}x \quad \left( m = 0,1,2,\ldots \right)\,.\]

We now establish the relationship between the coefficients \(\hat{c}_m\) and \(\hat{a}_m\), where \(\hat{a}_m\) denotes the coefficients of the Legendre–Fourier series expansion of the function \(w^{\prime} \left( x \right)\). For the orthonormal shifted Legendre polynomials, the following relation holds: \[\label{eq:gener95ritz95gal95orthnorm95leg95diff} \widehat{P}_m \left( x \right) = \frac{\ell}{2} A_{m} \left( A_{m + 1} \widehat{P}_{m + 1}^{\prime} \left( x \right) - A_{m - 1} \widehat{P}_{m - 1}^{\prime} \left( x \right) \right)\,.\tag{113}\] Upon considering 113 , the following result is derived: \[\label{eq:gener95ritz95gal95hat95a95integr} \hat{a}_m = \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_m \left( x \right) \mathrm{d}x = \frac{\ell}{2} A_{m} \left( A_{m + 1} \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_{m + 1}^{\prime} \left( x \right) \mathrm{d}x - A_{m - 1} \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_{m - 1}^{\prime} \left( x \right) \mathrm{d}x \right)\,.\tag{114}\] By applying integration by parts to the first summand of equality 114 , we obtain: \[\begin{align} A_{m + 1} \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_{m + 1}^{\prime} \left( x \right) \mathrm{d}x &= A_{m + 1} \left( \left[ w^{\prime} \left( x \right) \widehat{P}_{m + 1} \left( x \right) \right]_{0}^{\ell} - \int\limits_{0}^{\ell} w^{\prime \prime} \left( x \right) \widehat{P}_{m + 1} \left( x \right) \mathrm{d}x \right) \\ &= \frac{1}{\sqrt{\ell}} \left( w^{\prime} \left( \ell \right) + \left( -1 \right)^m w^{\prime} \left( 0 \right) \right) - A_{m + 1} \hat{c}_{m + 1}\,. \end{align}\] Using the same reasoning, for the second summand of equation 114 , we arrive at the following result: \[A_{m - 1} \int\limits_{0}^{\ell} w^{\prime} \left( x \right) \widehat{P}_{m - 1}^{\prime} \left( x \right) \mathrm{d}x = \frac{1}{\sqrt{\ell}} \left( w^{\prime} \left( \ell \right) + \left( -1 \right)^m w^{\prime} \left( 0 \right) \right) - A_{m - 1} \hat{c}_{m - 1}\,.\] By substituting the last two resulting equalities into 114 , one yields: \[\hat{a}_m = \frac{\ell}{2} A_{m} \left( A_{m - 1} \hat{c}_{m - 1} - A_{m + 1} \hat{c}_{m + 1} \right)\,,\] From this, it follows that: \[\hat{a}_{m}^{2} \leq \frac{\ell^2}{2} A_{m}^{2} \left( A_{m - 1}^{2} \hat{c}_{m - 1}^{2} + A_{m + 1}^{2} \hat{c}_{m + 1}^{2} \right)\,.\] Summing both sides of the preceding inequality from \(N + 1\) to \(\infty\) and invoking Bessel’s inequality leads to the following result: \[\begin{align} \label{eq:gener95ritz95gal95bess95ineq} \left\lVert S_N - w^{\prime} \right\rVert^2 = \sum_{m = N + 1}^{\infty} \hat{a}_{m}^{2} &\leq \frac{\ell^2}{2} \sum_{m = N + 1}^{\infty} A_{m}^{2} \left( A_{m - 1}^{2} \hat{c}_{m - 1}^{2} + A_{m + 1}^{2} \hat{c}_{m + 1}^{2} \right)\nonumber \\ & \leq \frac{\ell^2}{2} A_{N + 1}^{2} \left( A_{N}^{2} \sum_{m = N + 1}^{\infty} \hat{c}_{m - 1}^{2} + A_{N + 2}^{2} \sum_{m = N + 1}^{\infty} \hat{c}_{m + 1}^{2} \right)\nonumber \\ & \leq \frac{\ell^2}{2} A_{N + 1}^{2} \left( A_{N}^{2} + A_{N + 2}^{2} \right) \sum_{m = N}^{\infty} \hat{c}_{m}^{2} \leq \frac{\ell^2}{\left( 2N + 1 \right) \left( 2N + 5 \right)} \left\lVert w^{\prime \prime} \right\rVert^2\,. \end{align}\tag{115}\]

It is observed that, analogous to equation 97 , inequality 111 also holds for equation ?? , where the constant \(b_2\) is determined by equality ?? . Substituting inequality 115 into 111 yields: \[\label{eq:gener95ritz95gal95unif95norm95sec95derr95w} \max_{0 \leq x \leq \ell} \left\lvert w_N \left( x \right) - w \left( x \right) \right\rvert \leq \frac{b_2 \ell \sqrt{\ell}}{\sqrt{2\left( 2N + 1 \right) \left( 2N + 5 \right)}} \left\lVert w^{\prime \prime} \right\rVert\,.\tag{116}\]

For equation ?? , from the condition that the variation of the energy functional vanishes, it follows that: \[\label{eq:gener95ritz95gal95variat95energ95func} \int\limits_{0}^{\ell} \left( p w^{\prime 2} + q w^2 + f w \right) \mathrm{d}x = 0\,.\tag{117}\] By straightforward reasoning, from equation 117 , it can be deduced that: \[\left\lVert w^{\prime} \right\rVert^2 \leq \frac{1}{p_0} \left\lVert f \right\rVert \left\lVert w \right\rVert\,.\] From this, and considering 105 , the following is obtained: \[\label{eq:gener95ritz95gal95inqts95w95diff95w} \left\lVert w^{\prime} \right\rVert \leq \frac{\ell}{\pi p_0} \left\lVert f \right\rVert\,,\quad \left\lVert w \right\rVert \leq \frac{\ell^2}{\pi^2 p_0} \left\lVert f \right\rVert\,.\tag{118}\]

We shall rewrite equation ?? in the following form: \[p w^{\prime \prime} = -p^{\prime} w^{\prime} + q w + f\,.\] Whence, based on the estimates provided in 118 , the following conclusion follows: \[\label{eq:gener95ritz95gal95norm95second95w} \left\lVert w^{\prime \prime} \right\rVert \leq \frac{1}{p_0} \left( \max_{0 \leq x \leq \ell} \left\lvert p^{\prime} \left( x \right) \right\rvert \left\lVert w^{\prime} \right\rVert + \max_{0 \leq x \leq \ell} q \left( x \right) \left\lVert w \right\rVert + \left\lVert f \right\rVert \right) \leq \hat{b}_1 \left\lVert f \right\rVert\,.\tag{119}\] Finally, comparing inequalities 116 and 119 gives the estimate ?? .

8 Illustrative Numerical Results↩︎

In this section, we discuss four benchmark problems corresponding to the initial–boundary value problem 65 68 in order to verify consistency and to validate the performance of the proposed combined numerical schemes against the theoretical findings obtained. For these cases, exact analytical solutions are available. All coefficients appearing in equations 65 66 are set equal to one; i.e., \(\alpha = \beta = \gamma = \delta = a_1 = a_2 = 1\).

In Test [[bm:test4]](#bm:test4){reference-type=“ref” reference=“bm:test4”}, the initial data are prescribed in the form of a wave packet. In this case, an exact analytical solution is not available. In contrast to the previous cases, all coefficients retain the same values, except for the coupling coefficients, which are specified as \(a_1 = a_2 = 0.5\). Furthermore, in all considered settings, the spatial variable satisfies \(x \in \left[ 0, 2 \right]\). The length of the interval is chosen to be \(2\) because it coincides with the extent of the orthogonality interval of the standard Legendre polynomials.

The approximation errors for each temporal layer \(k = 2,3,\ldots n\) (recall that the cases \(k = 0\) and \(k = 1\) correspond to the prescribed initial data and are therefore excluded) are defined as the \(L^2 \left( 0,\ell \right)\)-norms of the differences between the exact solutions and their corresponding numerical approximations. Thus, we set: \[E_{1,k} = \left\lVert u\left( \cdot, t_k \right) - \tilde{u}_{k,N}\left( \cdot \right) \right\rVert\,,\quad \text{and}\quad E_{2,k} = \left\lVert v\left( \cdot, t_k \right) - \tilde{v}_{k,N}\left( \cdot \right) \right\rVert\,.\]

Since the accuracy and computational efficiency of the proposed algorithm are closely related to the effective computation of the integrals involved, these integrals are approximated using a Gauss–Legendre quadrature rule with an error-controlled subdivision of the integration interval.

To visualize the results for all benchmark problems under consideration, we opted to use the Okabe–Ito palette. The provided figures show the solid (green) line depicting the exact analytical solution, whereas the dashed (orange) line represents the corresponding numerical approximation. The approximation errors \(E_{1,k}\) and \(E_{2,k}\) are illustrated using blue and orange lines, respectively, with the temporal grid nodes indicated by circular and square markers.

We now turn our attention to the benchmark problems.

Test ↩︎

[]{#subsec:test\thetest label=“subsec:test\thetest”} Consider the case in which the exact analytical solutions are given by the following functions: \[u\left( x, t \right) = \sin\left( \frac{\pi}{2} t \right) \sin\left( \frac{\lambda \pi}{\ell} x \right)\,, \quad \text{and} \quad v\left( x, t \right) = \sin\left( \frac{\pi}{2} t \right) \sin\left( \frac{\lambda \pi}{\ell} x \right)\,.\]

In this setting, the temporal interval is defined by \(0 \leq t \leq 1\), and the oscillation parameter is taken as \(\lambda = 14\). Moreover, in this benchmark problem, the temporal interval is uniformly partitioned into \(n = 256\) subintervals.

Observe that, in this scenario, the accuracy attained with respect to the time variable is approximately \(10^{-5}\). Owing to the large number of spatial oscillations, it may be presumed that using \(N = 20\) trial functions is insufficient to yield a highly accurate approximation (see Figure 1). Nevertheless, it is important that the combined scheme maintains and replicates the qualitative structure of the exact solution. Thus, it is reasonable to expect that, upon increasing \(N\), the approximate solution should coincide with the exact solution with high accuracy.

Consider the case in which the previously used number of basis functions is increased by \(15\), resulting in \(N = 35\).

As we expected, by increasing the number of basis functions \(N\), the combined numerical scheme captures the oscillatory solution quite well, as illustrated in Figure 2. Moreover, Figure 2 () and Figure 2 () indicate that \(\max_{0 \leq k \leq 256} E_{1,k}\) and \(\max_{0 \leq k \leq 256} E_{2,k}\) are order of \(10^{-6}\).

Test ↩︎

[]{#subsec:test\thetest label=“subsec:test\thetest”} We now consider the setting in which the exact solutions assume the following form: \[u\left( x, t \right) = v\left( x, t \right) = A \left( 1 + \cos\left( \frac{\lambda_1 \pi}{T} t \right) \right) \exp\left( -\frac{\left( 2x - \ell \right)^2}{c^2} \right) \sin\left( \frac{\lambda \pi}{\ell} x \right)\,.\]

As in the preceding benchmark problem, the temporal variable is specified on the interval \(t \in \left[ 0, 1 \right]\), which is uniformly divided into \(n = 256\) subintervals. In contrast to Test [[ subsec:test1]](# subsec:test1){reference-type=“ref” reference=” subsec:test1”}, the amplitude of the sine function now varies with both the spatial and temporal variables. For the numerical computations, the parameters appearing in the given solutions are prescribed as \(A = 0.5\), \(\lambda_1 = 2\), \(c = 1\), and \(\lambda = 19\).

Because the oscillation parameter \(\lambda = 19\) is relatively large, employing \(N = 29\) trial functions does not yield an approximation of high accuracy. On the other hand, it is essential that the numerical solution faithfully replicates the qualitative features of the exact solution (see Figure 3). The adopted temporal discretization (namely, the grid length \(\tau = 2^{-8}\)) is fully adequate for approximating the solution with respect to the temporal variable.

In Figure 4, we observe that increasing the number of approximation basis functions to \(N = 45\) is sufficient to obtain numerical solutions that coincide with the exact analytical solutions with good accuracy. Furthermore, Figure 4 () and Figure 4 () demonstrate that the maximum values of the corresponding errors are approximately of the order of \(10^{-4}\).

Test ↩︎

[]{#subsec:test\thetest label=“subsec:test\thetest”} Let the exact solutions be given by the following functions: \[u\left( x, t \right) = \frac{1}{4} \exp\left( \frac{\pi}{T} t \right) \sin\left( \frac{\lambda \pi}{\ell} x \right)\,, \quad \text{and} \quad v\left( x, t \right) = \frac{1}{4} \exp\left( \frac{\pi}{T} t \right) \sin\left( \frac{\lambda \pi}{\ell} x \right)\,.\]

For this benchmark problem, unlike Test [[ subsec:test1]](# subsec:test1){reference-type=“ref” reference=” subsec:test1”} and Test [[ subsec:test2]](# subsec:test2){reference-type=“ref” reference=” subsec:test2”}, we consider a wider time frame, specifically the temporal interval satisfying \(0 \leq t \leq 4\). The oscillation parameter appearing in the sine functions is prescribed as \(\lambda = 5\). Initially, the uniform temporal grid spacing is chosen as \(\tau = 2^{-7}\), while the number of basis functions is set to \(N = 7\). The numerical results obtained under this configuration are shown in Figure 5.

As illustrated in Figure 5, the chosen configuration is inappropriate for accurately approximating the exact solutions using the proposed combined numerical scheme. Although the oscillation parameter is not particularly large, the time-dependent amplitude of the sine function grows rapidly. Given these facts, it is reasonable to refine both the temporal grid spacing and the number of basis functions, thereby setting \(\tau = 2^{-8}\) and \(N = 15\), respectively.

As depicted in Figure 6, improving both the temporal grid length to \(\tau = 2^{-8}\) and the number of basis functions to \(N = 15\) proved successful, resulting in an evident similarity between the exact and approximate solutions displayed in Figure 6 () and Figure 6 (). It may be observed that the approximation errors associated with each solution are small; refer to Figure 6 () and Figure 6 ().

Test ↩︎

[]{#subsec:test\thetest label=“subsec:test\thetest”} We consider the homogeneous reformulation of the system 65 66 , endowed with the initial-boundary conditions 67 68 , for which no exact analytical solutions are currently known. In this formulation, the spatial and temporal domains are taken to be \(x \in \left[ 0, 2 \right]\) and \(t \in \left[ 0, 1 \right]\), respectively.

The initial data in 67 68 are specified by the following functions: \[\varphi_0\left( x \right) = \psi_0\left( x \right) = A \exp\left( -\frac{\left( 2x - \ell \right)^2}{c^2} \right) \sin\left( \frac{\lambda \pi}{\ell} x \right) \quad \text{and} \quad \varphi_1\left( x \right) = \psi_1\left( x \right) = 0\,.\] In this setup, the prefactor multiplying the sine term is a Gaussian function. For the numerical experiments presented below, the parameters of the Gaussian are chosen as follows: the amplitude is set to \(A = 1\), the shift parameter to \(\ell = 2\), and the width (shape) parameter to \(c = 0.5\). The oscillation parameter appearing in the sine function is fixed at \(\lambda = 10\). For this configuration, the graph of the function \(\varphi_0\left( x \right)\) is depicted in Figure 7.

To describe the approximate solution procedure for the problem under consideration, we introduce the notation required for the subsequent discussion. Let \(n_0\) denote the initial number of subdivisions of the temporal interval \(\left[ 0,T \right]\), and let \(N_0\) denote the initial number of basis functions. At each refinement step, the number of temporal subdivisions is doubled, while at each substep the number of basis functions is increased by one. Thus, if \(n_i\) denotes the number of temporal subdivisions at the \(i\)-th refinement step and \(N_j\) denotes the number of basis functions at the \(j\)-th substep, then: \[\begin{gather} n_i = 2^i n_0\,,\quad t_k^{\left( i \right)} = k \tau_i\,,\quad \tau_i = \frac{T}{n_i}\,,\quad \text{for}\quad i = 0,1,\ldots,\hat{i}\,. \\ N_j = N_0 + j\,,\quad \text{for}\quad j = 0,1,\ldots,\hat{j}\,. \end{gather}\]

Let \(\widetilde{u}_{k,N_j}^{\left( i \right)} \left( x \right)\) and \(\widetilde{v}_{k,N_j}^{\left( i \right)} \left( x \right)\), for \(k = 2,3,\ldots,n_i\), denote the solution corresponding to problem 65 68 . Then, according to relation 88 , we have: \[u \left( x, t_k^{\left( i \right)} \right) \approx \widetilde{u}_{k,N_j}^{\left( i \right)} \left( x \right) = \sum_{m = 1}^{N_j} u_{k,m}^{\left( i,j \right)} \phi_m \left( x \right)\,,\quad v \left( x, t_k^{\left( i \right)} \right) \approx \widetilde{v}_{k,N_j}^{\left( i \right)} \left( x \right) = \sum_{m = 1}^{N_j} v_{k,m}^{\left( i,j \right)} \phi_m \left( x \right)\,,\] where the coefficients \(u_{k,m}^{\left( i,j \right)}\) and \(v_{k,m}^{\left( i,j \right)}\) are obtained by solving the associated Galerkin linear system (see Subsection 7.2) at refinement step \(i\) and substep \(j\).

The underlying idea of the proposed procedure is as follows. For each fixed value of \(n_i\), the computation in the spatial variable is continued (that is, the substep index \(j\) is incremented by one) until convergence is achieved at each temporal level \(k\), in the sense that the \(L^2\)-norms of the differences between two consecutive approximate solutions corresponding to \(u\) and \(v\) both become less than or equal to a prescribed tolerance \(\mathrm{tol} > 0\). More precisely, the process is terminated once the conditions \[\label{eq:num95res95norm95diffs} E_{k,j}^{\left( i \right)} = \left\lVert \widetilde{u}_{k,N_j}^{\left( i \right)} - \widetilde{u}_{k,N_{j - 1}}^{\left( i \right)} \right\rVert \leq \mathrm{tol} \quad \text{and} \quad \hat{E}_{k,j}^{\left( i \right)} = \left\lVert \widetilde{v}_{k,N_j}^{\left( i \right)} - \widetilde{v}_{k,N_{j - 1}}^{\left( i \right)} \right\rVert \leq \mathrm{tol}\tag{120}\] are simultaneously satisfied for each temporal layer \(k = 2,3,\ldots,n_i\).

If the convergence criterion fails to be satisfied during the spatial iteration, namely upon incrementing the iteration index \(j\) by one, then the temporal discretization is refined by doubling the number of time-step subdivisions, and the procedure is restarted from the initial stage. The indices \(\hat{i}\) and \(\hat{j}\) are fixed positive integers, typically chosen in accordance with the constraints imposed by available computational resources, feasible runtime bounds, and the effective stability margins of the proposed combined scheme.

The preceding formula follows immediately from the standard properties of the Legendre polynomials. Indeed, we have: \[\begin{align} \left\lVert \widetilde{u}_{k,N_j}^{\left( i \right)} - \widetilde{u}_{k,N_{j - 1}}^{\left( i \right)} \right\rVert^2 &= \left\lVert \sum_{m = 1}^{N_j} e_{k,m} \phi_m \right\rVert^2 = \sum_{m = 1}^{N_j} e_{k,m} \sum_{s = 1}^{N_j} e_{k,s} \left( \phi_s, \phi_m \right) \\ &= \frac{{\ell}^2}{4} \sum_{m = 1}^{N_j} e_{k,m} \left( - B_{m - 1} e_{k,m - 2} + C_m e_{k,m} - B_{m + 1} e_{k,m + 2} \right) = \frac{{\ell}^2}{4} \boldsymbol{e}_k^{\top} \boldsymbol{\mathcal{H}}_{N_j} \boldsymbol{e}_k\,. \end{align}\]

In the numerical experiments, the computational algorithm is initialized with \(n_0 = 2048\) time steps and \(N_0 = 41\) spatial (Galerkin) modes. Time refinement is controlled by the admissible exponent \(\hat{i} = 1\) (yielding a maximal refinement level of \(n_1 = 4096\)), while spatial refinement is restricted to a maximum index \(\hat{j} = 4\) (corresponding to an upper limit of \(N_4 = 45\)). The stopping criterion of the algorithm is based on the \(L^2\)-norm difference, with the tolerance set to \(\mathrm{tol} = 10^{-4}\).

For a given tolerance, the computational procedure is terminated when the parameters attain the values \(n_0 = 2048\) and \(N_1 = 42\), corresponding to the temporal refinement level \(i = 0\) and the spatial iteration index \(j = 1\), respectively. Numerical experiments conducted for different tolerance values, namely \(\mathrm{tol} = 10^{-s}\) with \(s = 1,2,3,4\), yield results whose qualitative structures are consistent with those illustrated, for instance, in Figure 8 for the case \(\mathrm{tol} = 10^{-4}\) at the time level \(t = 1\). The observed stabilization thus provides strong evidence for the robustness and reliability of the computed numerical solutions.

Note that the functions defining the initial conditions of the problem at hand are symmetric, and the system is considered in the homogeneous setting, i.e., in the absence of source terms. Therefore, one should expect the solution of the system to inherit a certain symmetry. This property is indeed observed in the numerical results displayed in Figure 8.

Initially, the wave packets are localized near the center of the string (see Figure 7). Then, according to linear theory, the initial profile splits into two equal parts that propagate in opposite directions: one to the left and the other to the right. Thus, at the time level \(t = 1\), the two halves of the initial profile have moved away from the center. Therefore, the central part of the string is almost at the equilibrium position (see Figure 8 ()). Furthermore, as depicted in Figure 8 (), when these separated parts reach the fixed endpoints, reflection with a change of sign occurs. Since we are solving a coupled system and have shown that the solution \(v\) corresponding to the linear equation is physically realistic, we may conclude that the solution \(u\), corresponding to the nonlinear equation, is also physically realistic.

Acknowledgments↩︎

The authors sincerely thank Dr. Andreas A. Buchheit and Dr. Daniel Seibel for their careful reading of the initial version of the manuscript and for their insightful comments, which helped improve the paper. We also wish to express our gratitude to Dr. Giorgi Rukhaia for his fruitful remarks during the development of the programming code for the proposed algorithm. Furthermore, the authors would like to thank the anonymous reviewers for their helpful comments, which contributed to improving the final version of the article.

Funding statement↩︎

The second author, Z.V., was supported by the Shota Rustaveli National Science Foundation of Georgia (SRNSFG) under grant number FR-25-215.

Data availability statement↩︎

The research data associated with this work are included in the article itself. The source code used for the implementation of the proposed algorithm is publicly available on GitHub https://github.com/zv1991/abstract_timoshenko_semidiscrete_scheme and in an open-access Zenodo repository [26].

References↩︎