1 Introduction

For convenience and without loss of generality, we take into account the following initial boundary value problem of the two-dimensional (2D) viscoelastic wave equation.

Problem 1 Find u=u(x, y, t) to satisfy

(1)

where is a bounded convex polygonal domain with a smooth boundary , , and . ε and γ are two positive constants, T is the final time, the source term f(x, y, t), the boundary value function Q(x, y, t), and the initial value functions G(x, y) and H(x, y) are four smooth enough known functions.

The viscoelastic wave equation (1) has many physical backgrounds. For example, it can be used to describe the vibration wave propagation phenomena through a viscoelastic medium, and the existence and uniqueness of its analytic solution had been proved^[1-4]. However, because the viscoelastic wave equation in the actual engineering applications usually has complex known data or computational domains, though theoretically there exists an analytical solution, it can not be generally solved out so that one has to find its numerical solutions. In more than thirty years, it has been closely watched, and many numerical methods^[5-8] for the viscoelastic wave equation have been developed. Among all numerical methods, the finite difference (FD) scheme^[5] is considered to be one of the simplest and most convenient as well as the most easily programming calculating numerical methods for solving the 2D viscoelastic wave equation. However, the classical FD scheme for the 2D viscoelastic wave equation is a macroscale system of equations including lots of unknowns (i.e., degrees of freedom) so that very large computational loads are required in the real-world engineering applications. Therefore, an important issue is how to lessen the unknowns of the classical FD scheme so as to ease the truncated error amassing and save the CPU time in the numerical computation but hold high enough accuracy of FD solutions.

It has been proved from lots of numerical examples^[9-24] that the proper orthogonal decomposition (POD) method is a very useful tool to lessen the unknowns for numerical models and ease the truncated error amassing in the numerical calculation. But most of existing reduced-order models as mentioned above were established with the POD basis formed from the classical numerical solutions at all time nodes, before computing the reduced-order numerical solutions at the same time nodes, which were some valueless repeated calculations. Since 2014, some reduced-order extrapolating FD schemes^[25-28] based on the POD method for partial differential equations (PDEs) have been established successively by Luo's team in order to avert the valueless repeated computations.

However, as far as we know, there has not been any report that the POD method is used to reduce the unknowns in the classical FD scheme for the 2D viscoelastic wave equation. Therefore, in this article, we intend to study the reduced-order of the classical FD scheme for the 2D viscoelastic wave equation, building an optimized finite difference iterative (OFDI) scheme including very few unknowns but owning high enough accuracy by means of the POD method, analyzing the stability and convergence of the OFDI solutions, and verifying the efficiency and feasibility of the OFDI scheme via numerical simulations.

The main distinction between the OFDI scheme and the existing reduced-order extrapolating FD schemes^[25-28] is twofold, i.e., the viscoelastic wave equation not only involves the second order derivative term about time and those about spacial variables but also includes a mixed derivative term about time of the first order and spacial variables of the second order. Therefore, either the modeling of the OFDI scheme or the demonstration of the stability and convergence of the OFDI solutions faces more difficulties and requires more techniques than the existing reduced-order extrapolating FD schemes as mentioned. However, the OFDI scheme has some specific applications. Fortunately, we adopt the vector and matrix analysis technique to discuss the existence, stability, and convergence of the classical FD and OFDI solutions such that not only the theoretical analysis becomes very simple and convenient but the numerical simulations in computer can also easily implement. Especially, we here only use the classical FD solutions at the initial very few time nodes to formulate the POD basis and build the OFDI scheme so that it has no repeated computation. Consequently, it is a development and an improvement over the existing researches as mentioned above.

The rest of the paper is organized as follows. In Section 2, we first present the classical FD scheme for the 2D viscoelastic wave equation and analyze the existence, stability, and convergence of the classical FD solutions. In Section 3, we establish the OFDI scheme based on the POD method for the 2D viscoelastic wave equation. Next, in Section 4, we analyze the existence, stability, and convergence of the OFDI solutions and provide the error estimates of the OFDI solutions. In Section 5, we use some numerical simulations to verify the efficiency and feasibility of the OFDI scheme. Finally, in Section 6, we summarize some main conclusions.

2 Classical FD scheme for 2D viscoelastic wave equation

Let Δt be the time step increment and Δx and Δ y be, respectively, the spacial step increments in the x-and y-directions. u_{i, j}ⁿ represents the classical FD solutions of u at points (x_i, y_j, t_n), where x_i=a+iΔx, y_j=c+jΔy, t_n=nΔt, 0≤ i≤ I≡[(b-a)/Δx], 0≤ j≤ J≡[(d-c)/Δy], and 0≤ n≤ N≡[T/Δt], where [P] denotes the integer part of the real number P.

If the derivatives of (1) are approximated by the following difference quotients:

we obtain the following classical FD scheme:

(2)

where i=1, 2, …, I-1, j=1, 2, …, J-1, n=1, 2, …, N-1, and f_{i, j}ⁿ=f(x_i, y_j, t_n). Let G_{i, j}=G(x_i, y_j) and H_{i, j}=H(x_i, y_j) as well as Q_{i, j}ⁿ=Q(x_i, y_j, t_n). The boundary conditions of the classical FD scheme are as follows:

(3)

By u_{i, j}¹-u_{i, j}^-1=2ΔtH_{i, j} and the case at n=0 in (2), we can deduce that u_{i, j}¹ and u_{i, j}⁰ satisfy the following equations:

(4)

(5)

Let

with M=(I+1)(J+1). Then, the matrix forms of (4), (5), and (2) are as follows:

(6)

(7)

(8)

where and I is the unit matrix.

Define the norm of matrix denotes the norm of the vector x. For the classical FD scheme, we have the following results.

Theorem 1 The FD scheme (6)-(8) has a unique set of solutions , which is stable and convergent when Δt, Δx, and Δy satisfy

Furthermore, when the analytic solution u∈ C((0, T]; the classical FD solutions Uⁿ (n=1, 2, …, N) have the following error estimates:

(9)

where solutions are formulated by means of the analytic solution of the viscoelastic wave equation (1), and is the positive generic constant independent of Δx, Δy, and Δt.

Proof Because B and C are two positive definite matrices, it is easily known that the matrix is positive definite. Therefore, (6)-(8) have a unique set of solutions and

(10)

Then, the classical scheme (8) is rewritten as follows:

(11)

Because the matrices B and C have the same eigenvalues,

we obtain the eigenvalues of A¹, A, and A-2 I as follows:

Thus, the eigenvalues of AD satisfy

(12)

where j=1, 2, …, M. By the relationships between roots and coefficients in the real coefficient quadric equation (see Theorem 1.5.2 of Ref. [29]), we conclude that the eigenvalues of AD are no more than 1 if and only if

(13)

and

(14)

Further, with the properties^[29] of the matrix norm, we have

and

Hence, from (13) and (14), we obtain that the stable condition of the classical scheme (11), i.e., (8), is as follows:

Thus, by Lax's stability theorem^[29-30], we easily deduce that the FD solutions are conditionally convergent. By Taylor's formulas or the above discrete process, we easily derive the error estimate (9). This finishes the proof of Theorem 1.

Remark 1 In this study, we adopt a simpler explicit FD scheme to discretize the 2D viscoelastic wave equation, but the ideas and approaches here can be easily extended to other FD schemes, for example, the staggered direction FD scheme or the implicit FD scheme.

If f, Q, G, H, Δt, Δx, Δy, and parameters ε and γ are given, then we can gain the set of FD solutions by solving (6)-(8). A subset , called as the snapshots, is extracted from the initial L solution vectors of

3 Establishment of OFDI scheme for 2D viscoelastic wave equation 3.1 Constitution of POD basis

For be the positive eigenvalues of ranked non-increasingly and be the orthonormal eigenvectors of related to the positive eigenvalues. Then, the POD basis is formed from the initial d vectors in U_u and has the following property^[16]:

(15)

Further, there holds

(16)

where are the unit vectors with the nth component being 1. Therefore, is an optimal basis.

Remark 2 The order L of the matrix is far smaller than the order M of the matrix i.e., the number of extracted snapshots L is far smaller than that of the spacial net points M, but their positive eigenvalues are the same. Hence, we may first find the eigenvalues for the matrix and the homologous eigenvectors , and then, by means of the formula , we can gain the eigenvectors corresponding to the positive eigenvalues for the matrix . Thus, we can conveniently compute out the POD basis.

3.2 Formulation of OFDI scheme for 2D viscoelastic wave equation

In Subsection 3.1, we have gained the initial OFDI solutions , where , and Thus, if the solution vectors for the classical FD scheme (6)-(8) are approximated by , i.e., are replaced by , we can obtain the OFDI scheme based on the POD basis as follows:

(17)

where are the known classical FD solution vectors to (8), and the matrices A and A₁ are given in (8). The OFDI scheme (17) is simplified into

(18)

Remark 3 Because the FD scheme (8) includes M=(I+1)(J+1) unknowns at each time node, while the OFDI scheme (18) at the same time node only includes d unknowns (d≪ M), we can clearly realize the superiority of the OFDI scheme (18).

4 Existence, stability, and convergence of OFDI solutions

In the following, we devote ourselves to analyzing the existence, stability, and convergence of the OFDI solutions. We summarize the following main conclusions.

Theorem 2 Under the conditions of Theorem 1, the OFDI scheme (18) has a unique set of solutions . As Δt, Δx, and Δy satisfy the stability condition , the set of solutions is stable and convergent and has the following error estimates:

(19)

where N). Moreover, the errors between the analytic solution vectors for the viscoelastic wave equation (1) and the OFDI solutions have the following estimates:

(20)

Proof By using , the OFDI scheme (17) is reverted into the following equations:

(21)

(22)

Because the classical FD solutions Uⁿ (n=1, 2, …, L) are known and conditionally stable, when n=1, 2, …, L, from (21), we obtain the unique solutions , which are stable since . Furthermore, because the matrix A is positive definite, when n=L+1, L+2, …, N, the OFDI scheme (22) exists a unique set of solutions , and when Δt, Δx, and Δ y satisfy that the stability condition

by using the same approaches as Theorem 1, we easily demonstrate that the OFDI solutions are stable. Thus, the OFDI scheme (22) has a unique set of stable solutions . Furthermore, by using Lax's stability theorem^[29-30], we deduce that the OFDI solutions are convergent.

When n=0, 1, …, L, from (16), we immediately obtain the following error estimates:

(23)

Let . By subtracting (22) from (8) and summing from L to n (n=L, L+1, …, N-1), we obtain

(24)

Because , and , when n=L+1, L+2, …, N, from (8) and (22), we have

(25)

By summing (25) from L to n-1 and using (23) and Gronwall inequality^[29-32], we have

(26)

Combining (23) with (26) yields (19), and combining Theorem 1 with (19) yields (20). This accomplishes the demonstration of Theorem 2.

Remark 4 The error factors and (n=L+1, L+2, …, N) in Theorem 2 are caused for the reduced-order of the classical FD scheme and for the iteration, respectively. They can be used as the suggestion that one chooses the number d of POD bases in the actual numerical computations. In fact, a lot of numerical experiments have shown that the eigenvalues of the symmetrical matrix are usually decreasing quickly to be close to zero. Therefore, if only we choose d such that , we can ensure that the OFDI solutions are convergent and acquire the optimal accuracy.

5 Numerical simulation

In this section, we give some numerical simulations in order to illustrate the superiority of the OFDI scheme for the 2D viscoelastic wave equation.

In the 2D viscoelastic wave equation (1), we choose the computational domain as Ω=. The spatial steps are chosen as Δx=Δy=0.01, and the time step is chosen as Δt=0.000 1.

First, the snapshots are obtained from the classical FD solutions u_{i, j}ⁿ (n=1, 2, …, 20) for the classical FD scheme (8) at the initial twenty steps. And then, the snapshot matrix is compiled, and the eigenvalues and the corresponding eigenvectors of the are computed. It is deduced by computation that when n=20 000, L=20, and γ=1. This implies that as long as choosing the initial four eigenvectors of matrix as the POD basis, the accuracy requirement of the OFDI solution can be met. In the end, the OFDI solutions at t=0.1, 1.0, and 2.0 are computed out via the OFDI scheme (18) with four POD bases, and the classical FD solutions are also obtained by means of the classical FD scheme (8) at the same time nodes, which are exhibited in Figs. 1-6, respectively.

Comparing the numerical conclusions of the classical FD scheme with the OFDI scheme summarizes some results as follows. The charts in Figs. 1 and 2 are basically identical at t=0.1. However, the classical FD solution exhibits a little dispersion (see Fig. 3) at t=1, whereas the OFDI solution is still stable and smooth (see Fig. 4). Especially, the classical FD solution emerges very great dispersion (see Fig. 5) at t=2 due to the truncated error amassing. However, the OFDI solution maintains stable and smooth (see Fig. 6). It follows that the OFDI solutions are much better than the classical FD solutions. It is shown that the OFDI scheme is far more efficient and advanced than the classical FD scheme for solving the 2D viscoelastic wave equation.

Thanks to the OFDI scheme greatly lessening the unknowns, the CPU time of the OFDI scheme is far less than that of the classical FD scheme in the above numerical simulations. For example, the CPU time for the OFDI scheme with four POD bases is only 8 s at t=2, while that of the classical FD scheme is about 120 s on the same PC, but the errors of both solutions do not exceed O(10^-4), which also shows that the numerical computational results are in line with the theoretical ones. The more detailed errors and CPU time are presented in Table 1.

Table 1 shows that the CPU time of the OFDI scheme is far less than that of the classical FD scheme. Even if the errors of the OFDI solutions are larger than those of the classical FD scheme over the initial time interval, as the real time goes on, the errors of the OFDI solutions are smaller than those of the classical FD ones because the OFDI scheme includes fewer degrees of freedom than the classical FD scheme so that the accumulation of the truncated errors of the OFDI solutions is far slower than that of the classical FD solutions.

Furthermore, in the above numerical simulations, we only choose the initial few given classical FD solutions over the very short time span [0, T₀] (T₀≪ T) as the snapshots to formulate the POD basis and establish the OFDI scheme in order to compute the OFDI solutions over the total time span [0, T]. However, for the real-world engineering problems, one may choose the recorded data (over the very short time [0, T₀]) to predict future physical phenomena and changes (over the time span [T₀, T]). Therefore, the OFDI scheme holds a very extensive applied prospect.

6 Conclusions

In this article, we have established the OFDI scheme for the 2D viscoelastic wave equation by means of the POD technique. First, the snapshot ensemble is formed by means of the initial few FD solutions for the 2D viscoelastic wave equation. The POD basis is formulated from the snapshot ensemble by means of the POD technique, and the OFDI scheme with sufficiently high accuracy but very few unknowns is established by replacing the unknown vectors of the FD scheme with the linear-combination of the POD basis. Finally, the stability and convergence of the OFDI solutions are analyzed. The numerical simulation shows that the numerical errors of the classical FD solutions with the OFDI ones correspond to the theoretical results. This implies that the OFDI scheme has high-efficiency and is reliable for solving the 2D viscoelastic wave equation.

Even though we only discuss the OFDI scheme for the 2D viscoelastic wave equation over the domain , the approach here can be extended to more general domains, even to the more complicated engineering problems. Therefore, the present technique has an important applied prospect.