We consider an iteration method for solving an elliptic type boundary value problem , where a positive definite operator is generated by a quasi-periodic structure with rapidly changing coefficients (a typical period is characterized by a small parameter ϵ). The method is based on using a simpler operator (inversion of is much simpler than inversion of ), which can be viewed as a preconditioner for . We prove contraction of the iteration method and establish explicit estimates of the contraction factor q. Certainly the value of q depends on the difference between and . For typical quasi-periodic structures, we establish simple relations that suggest an optimal (in a selected set of “simple” structures) and compute the corresponding contraction factor. Further, this allows us to deduce fully computable two-sided a posteriori estimates able to control numerical solutions on any iteration. The method is especially efficient if the coefficients of admit low-rank representations and if algebraic operations are performed in tensor structured formats. Under moderate assumptions the storage and solution complexity of our approach depends only weakly (merely linear-logarithmically) on the frequency parameter .
Problems with periodic and quasi-periodic structures arise in various natural sciences models and technical applications. Quantitative analysis of such problems requires special methods oriented towards their specific features. For perfectly periodic structures, efficient methods are developed within the framework of the homogenization theory (see, e.g., [1, 3, 8] and other literature cited therein). However, classical homogenization methods cover only one class of problems (all cells are self-similar and the amount of cells is very large). In this paper, we use a different idea and suggest another modus operandi for quantitative analysis of boundary value problems with periodic and quasi-periodic coefficients. It generates approximations converging (in the energy space) to the exact solution and provides guaranteed and computable error estimates. The approach is applicable to (see, e.g., Figures 1, 2)
periodic structures, in which the amount of cell is considerable (e.g., – ) but not large enough to neglect the error generated by the respective homogenized model;
quasi-periodic structures that contain cells with defects and deformations;
multi-periodic structures where the coefficients reflect the combined effect of several functions with different periodicity.
In general terms, the idea of the method is as follows. We consider the problem :
where V is a reflexive Banach space with the norm , is the space conjugate to V (the respective duality pairing is denoted by ), and is a bounded linear operator. It is assumed that the operator is positive definite and invertible, so that problem (1.1) is well posed. However, is viewed as a very difficult problem because is generated by a complicated physical structure, which may contain a huge amount of details. Therefore, attempts to solve (1.1) numerically by standard methods may lead to enormous expenditures. Similar difficulties arise if we wish to verify the quality of a numerical solution.
Assume that the operator is approximated by a simplified positive definite operator and the inversion of is much simpler than the inversion of . By means of , we construct an iteration method based on solving a “simple” problem : . In other words, the method is based on the operation . It also includes the operation , which can be performed very efficiently by tensor-type decomposition methods provided that physical structures generated have low-rank representations. We prove that iterations generate a sequence of functions converging to the exact solution of (1.1) with a geometrical rate. Furthermore, we deduce explicitly computable and guaranteed a posteriori error estimates adapted to this class of problems. They evaluate the accuracy of approximations computed on each step of the iteration algorithm. These estimates also use only inversion of and operations of the type . In the iteration methods and error estimates inversion of the operator is avoided.
In this paper, we consider one class of problems associated with divergent type elliptic equations where and . Here is a bounded operator induced by a complicated quasi-periodic structure while and are conjugate operators, i.e.,
where Y is a Hilbert space with the scalar product and the norm . The operators Q and are induced by differential operators or certain finite-dimensional approximations of them. Henceforth, it is assumed that , where is a Hilbert space with the scalar product . This space is intermediate between V and , i.e., .
The operator contains the operator generated by a simplified structure. We assume that the operators Λ and are Hermitian (i.e., and ) and satisfy the conditions
Then, the structural operators Λ and are spectrally equivalent:
where the constants are the minimal and maximal eigenvalues of the generalized spectral problem . Obviously, they satisfy the estimates and (which may be rather coarse).
Concerning the operator Q, we assume that there exists a positive constant c such that
Generalized solutions of the problems and are defined by the variational identities
In Section 2, we show that a sequence converging to u in V can be constructed by solving problems (1.4) with specially constructed right-hand sides generated by the residual of (1.3). In proving convergence, the key issue is analysis of the spectral radius of the operator
and selection of such relaxation parameter ρ that provides the best convergence rate. Moreover, iteration procedures of such a type become contracting if the iteration parameter is properly selected. This fact is often used in proving analytical results (see, e.g., , where classical results on existence and uniqueness of a variational inequality are established by contraction arguments). Also, these ideas were used in the construction of various numerical methods (see, e.g., ). However, achieving our goals requires more than the fact of contraction. We need explicit and realistic estimates of the contraction factor (which are used in error analysis) and a practical method of finding with minimal q. The latter task leads to a special optimization problem that defines the most efficient “simplified” operator among a certain class of “admissible” operators. This question is studied in Section 3. In general, Λ and can be induced by scalar, vector, and tensors functions. We show that selection of the optimal structural operator is reduced to a special interpolation type problem, which is purely algebraical and does not require solving a differential problem (therefore a suitable can be found a priori). We discuss several examples and suggest the corresponding optimal (or quasi-optimal) , which guarantees convergence of the iteration sequence with explicitly known contraction factor.
Now, it is worth discussing the main differences between our approach and the classical homogenization method developed for regular periodic structures. This method operates with a homogenized boundary value problem , where is defined by means of an auxiliary problem with periodical boundary conditions in the cell of periodicity. The respective solution contains an irremovable (modeling) error depending on the cell diameter ϵ. Moreover, if ϵ tends to zero, then typically converges to u only weakly (e.g., in ). Getting a better convergence (e.g., in ) requires certain corrections, which lead to other (more complicated) boundary value problems in the cell of periodicity. The respective “corrected” solution also contains an error. Typically, the error is proportional to and can be neglected only if the amount of cells is very large. If our method is applied to perfectly periodical structures then setting is one possible option. In this case, the homogenized operator (defined without correction procedures) is used for a different purpose: construction of a suitable preconditioning operator. The latter operator generates numerical solutions converging to the exact solution in the energy norm (i.e., the method is free from irremovable errors) and can be applied for a rather wide range of ϵ. In addition, the theory suggests other simpler ways of selecting suitable . In this context, it is interesting to know whether or not the choice always yields the minimal value of the contraction factor. In Section 3, we briefly discuss this question and present an example of that the best may differ from .
In Section 4, we deduce a posteriori estimates that provide fully computable and guaranteed estimates of the distance to the exact solution u for any numerical approximation computed for an approximation subspace . These estimates are established by combining functional type a posteriori estimates (see [25, 22, 26] and references cited therein) and estimates generated by the contraction property of the iteration method (see [24, 29]).
The second part of the paper is devoted to a fast solution method for the basic iteration problem (2.1). The key idea consists of using tensor-type representations for approximations, what is quite natural if both coefficients of the respective quasi-periodic structure and the right-hand side admit low-rank tensor-type representations. We notice that the amount of structures representable in terms of low-rank formats is much larger than the amount of periodic structures covered by the homogenization method. The idea of tensor-type approximations of partial differential equations traces back to . In computational mechanics this method is known as the Kantorovich–Krylov (or extended Kantorovich) method. However, it is rarely used in modern numerical technologies. In part, this is due to restrictions on the shape of the domain imposed by the Kantorovich method. Henceforth, we assume that the domain Ω satisfies these restrictions, i.e., it is a tensor-type domain (e.g., rectangular) or a union of tensor-type domains. This assumption induces certain geometrical limitations. However, they can be bypassed by such methods as coordinate transformation and domain decomposition, which are widely used in modern computational mathematics (e.g., in iso-geometric analysis).
The recent tensor numerical methods (for steady state and dynamical problems) based on the advanced nonlinear tensor approximation algorithms have been developed in the last ten years. Literature survey on the modern tensor numerical methods for multi-dimensional PDEs can be found in [14, 16, 13]. In the context of problems considered in the paper, we are mainly concerned with another specific feature: very complicated material structure. In this case, direct application of standard finite element methods suffers from the necessity to account huge information encompassed in coefficients (especially in multi-dimensional problems). We show that tensor-type methods allow us to reduce computations to a collection of one-dimensional problems, which can be solved very efficiently using low-rank representations with the small storage requests. Similar ideas are applied for computing a posteriori error estimates.
Section 5 discusses numerical aspects of the method and exposes several examples. Typical behavior of quasi-periodic coefficients is described by oscillation around a constant, modulated oscillation around a given smooth function, or oscillation around a piecewise constant function.
Figure 1 (1D case) represents examples of highly oscillating (left) and modulated periodic coefficients (right) functions. Figure 2 (2D case) illustrates the well-separable equation coefficient obtained by a sum of step-type and uniformly oscillating functions, namely,
where , , , and .
We show that specially constructed FEM type approximations of PDEs with slightly perturbed or regularly modulated periodic coefficients on d-fold tensor grids in may lead to the discretized algebraic equations with the low Kronecker rank stiffness matrix of size , where is proportional to the large frequency parameter . In this case, the rank decomposition with respect to the d spacial variables is applied, such that the discrete solution can be calculated in the low-rank separable form, which requires storage size of only instead of complexity representations which are mandatory for the traditional FEM techniques (the latter quickly leads to the bottleneck in case of small parameter ).
The arising linear system of equations can be solved by preconditioned iteration with the simple preconditioner , such that the storage and numerical costs scale almost linearly in the univariate discrete problem size n. Numerical examples in Section 5 demonstrate the stable geometric convergence of the preconditioned CG (PCG) iteration with the preconditioner and confirm the low-rank approximate separable representation to the solution with respect to d spacial variables even in the case of complicated quasi-periodic coefficients.
This approach is well suited for applying the quantized-TT (QTT) tensor approximation  to functions discretized on large tensor grids of size proportional to the frequency parameter, i.e., , as it was demonstrated in the previous paper  for the case . The use of tensor-structured preconditioned iteration with the adaptive QTT rank truncation may lead to the logarithmic complexity in the grid size, , see [14, 16, 23] for the rank-truncated iterative methods, [11, 12, 4, 10] for various examples of the QTT tensor approximation to lattice structured systems, and  for tensor approximation of complicated functions with multiple cusps in .
In Section 6, we conclude with the discussion on further perspectives of the presented approach for 2D and 3D elliptic PDEs with quasi-periodic coefficients.
2 The Iteration Method
Let and . Consider the problem: find such that
Obviously, the right-hand side of (2.1) is a bounded linear functional on V, so that this problem has a unique solution . Thus, we have a mapping , which becomes a contraction if the parameter ρ is properly selected. Indeed, for any and in V, we obtain
where , , , and . Hence
From (2.2) we find that
where is defined by (1.5). If ρ is selected such that
then (2.3) shows that is a contractive mapping.
Since Λ and are invertible with trivial kernels, μ and are an eigenvalue and the respective eigenfunction of if and only if they are an eigenvalue and the eigenfunction of the problem . This means that
and (2.5) implies
The minimum of the expression in round brackets is attained if . For , we find that
Hence is a contractive mapping with explicitly known contraction factor . Well-known results in the theory of fixed points (see, e.g., ) yield the following theorem.
For any and , the sequence of functions satisfying the relation
converges to u in V and as .
From (2.3) we obtain
This relation yields a simple (but not very sharp) estimate of the contraction factor.
For further analysis, it is convenient to estimate the right-hand side of (2.3) by a different method. Let denote the operator norm
Hence (2.3) yields the estimate
which shows that is a contraction provided that
In applications is a self-adjoint bounded operator acting in a finite-dimensional space, so that verification of this condition amounts to finding ρ which yields the respective spectral radius of (see Section 4).
3 Selection of
In this section, we discuss how to select in order to minimize q, which is crucial for two major aspects of quantitative analysis: convergence of the iteration method and guaranteed a posteriori estimates. We assume that V, , and Y are spaces of functions defined in a Lipschitz bounded domain Ω (namely for a.e. where may coincide with , , or ) and the operators Λ and are generated by bounded scalar functions, matrices or tensors. In this case,
whose computation is reduced to solving algebraic problems at a.e. , i.e.,
Let be a certain set of “simple” operators defined a priori (e.g., it can be a finite-dimensional set formed by piecewise constant or polynomial functions). Then, finding the best “simplified” operator amounts to solving the following problem: find such that is minimal. In other words, the optimal is defined by the problem
Notice that (3.1) is an algebraic problem, which should be solved (analytically or numerically) before computations. The respective solution defines the best operator to be used in the iteration method (2.6) and yields the respective contraction factor. Below we discuss some particular cases, where analysis of this problem generates an optimal (or almost optimal) .
Problem (3.1) is explicitly solvable if and Λ have a special structure, namely,
where is the unit operator and and are positive bounded functions defined in Ω. Then,
Define and . It is not difficult to show that
Minimization with respect to ρ yields the best value and the respective value
In accordance with (3.1), the identification of the optimal simplified problem is reduced to the problem
where is a given set of functions.
We illustrate the above relations by means of several examples.
Example 3.1 (Constant Coefficients).
From Theorem 2.1, it follows that
Example 3.2 (Oscillation Around a Given Function).
Consider a somewhat different example. Let be a function oscillating around a certain mean function so that
If g is a relatively simple function, then it is natural to set . By (3.2), we find that , , and . Hence the method is very efficient for small ϵ (i.e., if a oscillates around g with a relatively small amplitude). Figures 1 and 2 illustrate three examples of quasi-periodic coefficients a and respective corresponding to the case of oscillation around a constant with smooth modulation, oscillation around a given smooth function, or oscillation around a piecewise constant function.
Example 3.3 (Piecewise Constant Coefficients).
Consider a more complicated case, where Ω is divided into N nonoverlapping subdomains and if . Define the numbers
Since the constants are defined up to a common multiplier, we can without a loss of generality assume that
In accordance with (3.3), the maximum of is attained if
It is interesting to compare these results with those generated by homogenized models in the case of perfectly periodic structures. For this purpose, we consider a simple one-dimensional problem
where is a perfectly periodical function attaining only two values and . The Lebesgue measure of the set where is , . Analogously, is a periodical function attaining only two values and . The Lebesgue measure of the set where is , . We assume that , and the amount of periods is very large. Then, the homogenization method can be successfully applied. The corresponding homogenized problem has the coefficients (see, e.g., )
It is easy to see that
and . Therefore, the coefficients and may not generate the best piecewise constant , which produces the smallest contraction factor in the iteration procedure (2.6).
4 Error Estimates
4.1 General Estimate
This estimate cannot be directly applied because is generally unknown (it is the exact solution of a boundary value problem). Instead, we must use a numerical approximation (in our analysis, we impose no restrictions on the method by which the function was constructed). Thus, the difference is a known function and the quantity is directly computable. It is easy to see that
For any and , we have
We estimate the first term on the right-hand side of (4.4) as follows:
where . The second term meets the estimate
where is the dual norm. Hence
It is not difficult to show that the last term of can be estimated via an explicitly computable quantity provided that y has the same regularity as the true flux (see ). However, in our subsequent analysis and these advanced forms of the majorant are not required. In this case, the majorant has a simpler form:
It is important that the computation of the majorant does not require inversion of the operator Λ associated with a complicated quasi-periodic problem.
The error is subject to the estimate
where is defined by (4.5), , and y is a function in Y.
Here and are directly computable and is defined in accordance with relations presented in the previous section. Hence the cost of (4.6) is mainly related to the quantity , which is an a posteriori error majorant of the functional type. The derivation of such estimates is performed by purely functional methods and does not exploit special features of approximations (e.g., Galerkin orthogonality), numerical method, and exact solution (e.g., extra regularity). Properties of the majorants are well studied (see [25, 26] and the literature cited therein). Numerous tests performed for different boundary value problems have confirmed high practical efficiency of error majorants of the functional type. It was shown that is a guaranteed and efficient majorant of the global error and generates good indicators of local errors if y is replaced by a certain numerical reconstruction of the exact dual solution. There are many different ways to obtain suitable reconstructions with minimal expenditures (concerning this point we refer to  where the reader will find a systematic discussion of computational aspects in the context of various boundary value problems). Error majorants of this type can be also used for the evaluation of modeling errors (see [28, 27]).
Usually, the cost of a good estimate (with the efficiency index between 1 and 2) is comparable with the cost of a numerical solution. However, the proportion essentially depends on the numerical method used. For the classical FEM schemes the expenditures are maximal (because this method generates rather coarse approximations of fluxes). For the dual mixed method, finite volume method, isogeometric approximations, and other methods producing locally equilibrated fluxes, the expenditures may be two to three times smaller than for the numerical solution.
Now we shortly discuss applications of Theorem 4.2 to problems where Q and are defined by the operators and , respectively, , , , and .
In order to apply Theorem 4.2, we set , where and μ is a constant. Then and . The best constant μ is defined by minimization of , which has the form
Since , the problem is reduced to minimization of the second term and the best μ satisfies the equation
and (4.6) yields the estimate
Here v and are two consequent numerical approximations (e.g., finite element approximations and computed on a mesh ). Then
are directly computable. Since is a “simple” function, the integrals
are easy to compute. Other integrals
contain a highly oscillating coefficient a multiplied by piecewise polynomial mesh functions. If a and f have low QTT rank tensor representations , then the integrals can be efficiently computed by tensor-type methods discussed in  (see also Section 5 below). We have
Notice that the quantity is the value of the majorant, where the flux has been selected in the best way. It is not difficult to show that
In other words, this term coincides with the error of the Galerkin solution related to the simplified boundary value problem (4.7), where . In accordance with Section 3, we set and find that . Now (4.8) yields easily computable lower and upper bounds of the error encompassed in :
It is worth adding comments on convergence properties of the quantities and entering (4.9). If the mesh is fixed, then tends to the Galerkin approximation of problem (1.1) on this mesh. This fact follows from Theorem 2.1 applied to the case where the iterations are performed on the respective finite-dimensional space (considered as the space V). Then, and for any h the term tends to zero with the geometric rate. The quantity is equal to the error of the Galerkin solution to the simplified problem. It has different asymptotic properties. It mainly depends on , and for a given mesh it does not tend to zero when . However, for any given v (which in our example is defined by ) this term goes to zero if provided that the mesh satisfies the standard regularity conditions. The problem with is assumed to be much more regular than the problem with a. Therefore, in terms of h the approximation error (associated with ) will decrease faster than the analogous error in the original problem (e.g., for , the term is proportional to h).
Since both quantities and are explicitly known, estimate (4.9) (and other analogous estimates) contains a very useful information unavailable in the context of purely asymptotic error analysis. Using this information, we can organize the computational process in the best possible way by comparing iteration and discretization errors. In this process, rapidly converging iterations with respect to k should be continued until . If , then further iterations on the mesh are unable to essentially improve the numerical solution. Instead, we should refine , project on the refined mesh, and use it as a starting point for a new series of iterations generated by problem (4.7). For each step, we compute the right-hand side of (4.9) and stop the process when it becomes smaller than the desired tolerance.
The computation of for 2D problems can be also reduced to the computation of one-dimensional integrals. Certainly, on the multidimensional case the amount of integrals is much larger. However, the basic tensor decomposition methods remain the same. Below we briefly discuss them with the paradigm of a simple case where
Assume that approximations are represented in the form of series formed by one-dimensional functions and (which may be supported locally or globally), so that
In this case,
We define another set of one-dimensional functions and , which form the vector function
Here is a given function, which can be defined in different ways. In particular, we set
The functions must satisfy the usual linear independence conditions in order to guarantee unique solvability of the respective approximation problem. For any smooth function w vanishing on , we have
Thus, and we can use the simplified form of .
In the simplest case , where is a constant. The best y minimizes the quantity
which shows that y must satisfy the relation . We select that defines the Galerkin approximation of this function, and we arrive at the system
Introduce the following matrices:
It is not difficult to see that
where is the fourth-order tensor. Hence the left-hand side of the system (4.12) has the form . In the right-hand side we have the term
where , . Another term is
where , .
Now (4.12) implies
5 Low-Rank Solution of the Discrete Equation
We consider the following elliptic diffusion equation with quasi-periodic coefficient (whose oscillations are characterized by the parameter ϵ):
where the function f corresponds to the modified right-hand side in problem (4.3). In this case , , , and .
In what follows we assume that f and a admit low-rank representation i.e.,
where the parameters and are called the separation rank. Then one may assume that the exact FEM solution can be well approximated by , where K depends on the separation rank of f and a. In some cases this important property can be rigorously proven (e.g., for the Laplacian and other closely related operators). Similar low-rank approximations can be observed for the QTT tensor approximations (see ). Existence of a low-rank solution means that for some moderate K we have up to the rank truncation threshold.
First, we sketch the rank-structured computational scheme. In our set of examples the original problem is to find u such that
It is approximated by the following Galerkin problem for the low-rank representation :
where is a subset of formed by functions of the type
Therefore, in terms of the general scheme exposed in the introduction, the problem is now problem (5.2) and we solve it by iterations with the help of the simplified (preconditioned) problem
where depends on and is a simple function (i.e., it is representable by a sum of terms with very simple multipliers).
Given the right-hand side, problem (5.3) is much simpler than the initial problem and the stiffness matrix associated with (5.3) is computed much easier and has a simpler (low Kronecker rank) form that allows a rank-structured representation of its inverse.
5.1 Kronecker Product Representation of the Stiffness Matrix
Figure 3 illustrates a 2D example of the periodic coefficient with corresponding to the choice . In this example, the scalar coefficient is represented by the separable function
where the generating univariate function has the shape of six uniformly distributed bumps of height 1 as shown in Figure 3, right. Figure 3, left, presents the oscillating part of a 2D coefficients function, which is . Here, the coefficient bumps are displaced on the coarse grid of size in such a way that bumps occupy the central box in each of the cells, which compose the whole lattice-type decomposition of Ω (we have in Figure 3, i.e., the size of the coarse grid is , while in Figures 4 and 5, i.e., the size of the coarse grid is ). Hence the axis scale denotes the coarse grid in both and that describes the construction of coefficient in detail.
We apply the FEM Galerkin discretization of equation (5.1) by means of tensor-product piecewise affine basis functions (instead of “linear finite elements”)
where are 1D finite element basis functions (say, piecewise linear hat functions).
We associate the univariate basis functions with the uniform grid , on with the mesh size . In this construction we have basis functions . Notice that the univariate grid size is of the order of designating the total problem size .
For ease of exposition we first consider the case , and further assume that the scalar diffusion coefficient can be represented in the form
with a small rank parameter R.
The stiffness matrix is constructed by the standard mapping of the multi-index into the N-long univariate index i representing all degrees of freedom. For instance, we use the so-called big-endian convention for and :
respectively. Hence all matrices and vectors are defined on the long index i as usual, however, the special Kronecker structure allows the low-storage and low-complexity matrix vector multiplications when appropriate, i.e., when a vector also admits the low-rank Kronecker form representation. In particular, the basis function is designated via the long index, i.e., .
First, we consider the simplest case and let . We construct the Galerkin stiffness matrix in the form of a sum of Kronecker products of small “univariate” matrices. Recall that given matrix A and matrix B, their Kronecker product is defined as a matrix C via the block representation
We say that the Kronecker rank of the matrix A in the representation above equals 1. Now the elements of the Galerkin stiffness matrix take the form
which leads to the rank-2 Kronecker product representation
where denotes the conventional Kronecker product of matrices. Here and denote the univariate stiffness matrices and and define the corresponding weighted mass matrices, e.g.,
By simple algebraic transformations (e.g., by lamping of the tri-diagonal mass matrices, which does not effect the approximation order of the FEM discretization) the matrix can be simplified to the form
where are the diagonal matrices. The matrix A corresponds to the FEM discretization of the initial elliptic PDE with complicated highly oscillating coefficients.
The simple choice of the spectrally equivalent preconditioner corresponds to the operator Laplacian. In this case the representation in (5.4) is simplified to the discrete Laplacian matrix in the form of rank-2 Kronecker sum
where and denote the identity matrices of the corresponding size. Here the simple tree-diagonal matrices and represent the FEM/FDM Laplacian in 1D. This matrix will be used in what follows as a prototype preconditioner for solving the linear system of equations
The matrix A is constructed in general for the R-term separable coefficient with which leads to the rank- Kronecker sum representation
with matrices of the respective size.
5.2 On Existence of the Low-Rank Solution
In this paper we discuss the approach based on the low-rank separable ϵ-approximation of the solution to equation (5.6) that is considered as the d-dimensional real-valued array . In general, for the case this favorable property is not guaranteed by the low Kronecker rank representation to the Galerkin system matrix A, discussed in Section 5.1.
Let and . The existence of the low-rank approximation to the solution of equation (5.6) with the low-rank right-hand side
taking into account that the matrices and commute with and , respectively. Hence equation (5.7) represents the accurate rank- Kronecker product approximation to the preconditioner which can be applied directly to the right-hand side to obtain
The numerical efficiency of the representation (5.7) can be explained by the fact that the quadrature parameters can be chosen in such a way that the low Kronecker rank approximation converges to exponentially fast in M. For example, under the choice , with there holds 
in the Frobenius norm, which means that the approximation error can be achieved with the number of terms of the order of .
Figures 4 and 5 demonstrate the singular values of the discrete solution on the grid for and 191, indicating very moderate dependence of the ϵ-rank on the grid size n. As in the case of Figure 3, in Figures 4 and 5 we only represent the oscillating part of the coefficients and omit the small constant .
Further enhancement of the tensor approximation can be based on the application of the quantized-TT (QTT) tensor approximation which has been already applied in  to the 1D equations with quasi-periodic coefficients. The power of the QTT approximation method is due to the perfect low-rank decompositions applied to the wide class of function-related tensors . See  for a more detailed discussion and a number of numerical examples.
One can apply QTT approximations to problems with quasi-periodic coefficients, which can be described by oscillation with smooth modulation around a constant value, oscillation around a given smooth function, or oscillation around a piecewise constant function, see Figure 1 and examples in .
Let the vector , , be obtained by sampling a continuous function (or even piecewise smooth functions), on the uniform grid of size N. For the following examples of univariate functions the explicit QTT-rank estimates of the corresponding QTT tensor representations are valid uniformly in the vector size N, see :
for complex exponentials, , .
for trigonometric functions, , , .
for polynomials of degree m.
For a function f with the QTT-rank modulated by another function g with the QTT-rank r (say, step-type function, plain wave, polynomial) the QTT rank of a product fg is bounded by a multiple of r and ,
Furthermore, the following result holds : the QTT rank for the periodic amplification of a reference function on a unit cell to a rectangular lattice is of the same order as that for the reference function.
5.3 Numerical Test on the Rank Decomposition of u
The PCG solver for the system of equations (5.6) with the discrete Laplacian inverse as the preconditioner demonstrates robust convergence with the rate . The next example demonstrates the rank behavior in the singular value decomposition (SVD) of a matrix representing the solution vector to equation (5.6) with the periodic coefficient shown in Figure 4, left. Figure 7 represents the rank behavior in the SVD decomposition of the solution in the case of the periodic coefficient.
Comparing Figures 4 and 7 indicates that the exponential decay of the approximation error in the rank parameter is stable with respect to the size of the lattice structure of the coefficient, i.e., the behavior of the singular values remains almost the same for different parameters .
Our iterative scheme includes only the matrix-vector multiplication with the stiffness matrix A that has the small Kronecker rank , and the action of the preconditioner defined by the approximate inverse to the Laplacian type matrix. The latter has low Kronecker rank of order as shown above.
Given rank-1 vector , the standard property of the Kronecker product matrices
indicates that the matrix-vector multiplication enlarges the initial rank by the factor of 2, and similar with the action of the preconditioner. Hence each iterative step should be supplemented with certain rank truncation procedure which can be implemented adaptively to the chosen approximation threshold or fixed bound on the rank parameter.
Notice that for the transformed matrix takes the form
and it obeys the d-term Kronecker sum representation. Hence in the general case of and the Kronecker rank of the matrix is given by
We present a preconditioned iteration method for solving an elliptic type boundary value problem in with the operator generated by a quasi-periodic structure with rapidly changing coefficients characterized by a small length parameter ϵ. We use tensor product FEM discretization that allows to approximate the stiffness matrix A in the form of a low-rank Kronecker sum. The preconditioner is constructed based on certain averaging (homogenization) procedure of the initial equation coefficients such that the inversion of is much simpler than the inversion of A. We prove contraction of the iteration method and establish explicit estimates of the contraction factor . For typical quasi-periodic structures we deduce fully computable two-sided a posteriori estimates which are able to control numerical solutions on any iteration.
We apply the tensor-structured approximation which is especially efficient if the equation coefficients admit low-rank representations and algebraic operations are performed in tensor structured formats. Under moderate assumptions the storage and solution complexity of our approach depends only weakly (merely linear-logarithmically) on the frequency parameter . Numerical tests demonstrate that the FEM solution allows the accurate low-rank separable approximation which is the basic prerequisite for application of the tensor numerical methods to the problems of geometric homogenization.
The approach allows further enhancement based on the quantized-TT (QTT) tensor approximation which is the topic for future research work. Another direction is related to fully tensor structured implementation of the computable two-sided a posteriori error estimates. The interesting question arises how far the presented approach can be extended to the numerical analysis of elliptic equations with rather unstructured jumping coefficients arising in stochastic homogenization, see, e.g., .
SR appreciates the support provided by the Max-Planck Institute for Mathematics in the Sciences (Leipzig, Germany) during his scientific visit in 2016. The authors are thankful to Dr. V. Khoromskaia (MPI MIS, Leipzig) for the help with numerical experiments.
 N. S. Bakhvalov and G. Panasenko, Homogenisation: Averaging Processes in Periodic Media. Mathematical Problems in the Mechanics of Composite Materials, Springer, Berlin, 1989. Search in Google Scholar
 P. Benner, V. Khoromskaia and B. N. Khoromskij, Range-separated tensor formats for numerical modeling of many-particle interaction potentials, preprint (2016), http://arxiv.org/abs/1606.09218. Search in Google Scholar
 A. Bensoussan, J.-L. Lions and G. Papanicolaou, Asymptotic Analysis for Periodic Structures, North-Holland, Amsterdam, 1978. Search in Google Scholar
 S. Dolgov, V. Kazeev and B. N. Khoromskij, The tensor-structured solution of one-dimensional elliptic differential equations with high-dimensional parameters, Preprint 51/2012, MPI MiS, Leipzig, 2012. Search in Google Scholar
 I. P. Gavrilyuk, W. Hackbusch and B. N. Khoromskij, Hierarchical tensor-product approximation to the inverse and related operators in high-dimensional elliptic problems, Computing 74 (2005), 131–157. Search in Google Scholar
 A. Gloria and F. Otto, Quantitative estimates on the periodic approximation of the corrector in stochastic homogenization, ESAIM Proc. 48 (2015), 80–97. Search in Google Scholar
 R. Glowinski, J.-L. Lions and R. Trémolierés, Analyse Numérique des Inéquations Variationnelles, Dunod, Paris, 1976. Search in Google Scholar
 V. V. Jikov, S. M. Kozlov and O. A. Oleinik, Homogenization of Differential Operators and Integral Functionals, Springer, Berlin, 1994. Search in Google Scholar
 L. V. Kantorovich and V. L. Krylov, Approximate Methods of Higher Analysis, Interscience, New York, 1958. Search in Google Scholar
 V. Kazeev, O. Reichmann and C. Schwab, Low-rank tensor structure of linear diffusion operators in the TT and QTT formats, Linear Algebra Appl. 438 (2013), no. 11, 4204–4221. Search in Google Scholar
 V. Khoromskaia and B. N. Khoromskij, Grid-based lattice summation of electrostatic potentials by assembled rank-structured tensor approximation, Comp. Phys. Commun. 185 (2014), no. 12, 3162–3174. Search in Google Scholar
 V. Khoromskaia and B. N. Khoromskij, Tensor approach to linearized Hartree–Fock equation for lattice-type and periodic systems, preprint (2014), https://arxiv.org/abs/1408.3839. Search in Google Scholar
 V. Khoromskaia and B. N. Khoromskij, Tensor numerical methods in quantum chemistry: From Hartree–Fock to excitation energies, Phys. Chem. Chem. Phys. 17 (2015), 31491–31509. Search in Google Scholar
 B. N. Khoromskij, Tensor-structured preconditioners and approximate inverse of elliptic operators in , J. Constr. Approx. 30 (2009), 599–620. Search in Google Scholar
 B. N. Khoromskij, -quantics approximation of N-d tensors in high-dimensional numerical modeling, Constr. Approx. 34 (2011), 257–280. Search in Google Scholar
 B. N. Khoromskij, Tensors-structured numerical methods in scientific computing: Survey on recent advances, Chemometr. Intell. Lab. Syst. 110 (2012), 1–19. Search in Google Scholar
 B. N. Khoromskij and S. Repin, A fast iteration method for solving elliptic problems with quasiperiodic coefficients, Russian J. Numer. Anal. Math. Modelling 30 (2015), no. 6, 329–344. Search in Google Scholar
 B. N. Khoromskij, S. Sauter and A. Veit, Fast quadrature techniques for retarded potentials based on TT/QTT tensor approximation, Comput. Methods Appl. Math. 11 (2011), no. 3, 342–362. Search in Google Scholar
 B. N. Khoromskij and G. Wittum, Numerical Solution of Elliptic Differential Equations by Reduction to the Interface, Lect. Notes Comput. Sci. Eng. 36, Springer, Berlin, 2004. Search in Google Scholar
 J.-L. Lions and G. Stampacchia, Variational inequalities, Comm. Pure Appl. Math. 20 (1967), 493–519. Search in Google Scholar
 O. Mali, P. Neittaanmaki and S. Repin, Accuracy Verification Methods. Theory and Algorithms, Springer, New York, 2014. Search in Google Scholar
 P. Neittaanmaki and S. Repin, Reliable Methods for Computer Simulation. Error Control and a Posteriori Estimates, Elsevier, Amsterdam, 2004. Search in Google Scholar
 I. V. Oseledets and S. V. Dolgov, Solution of linear systems and matrix inversion in the TT-format, SIAM J. Sci. Comput. 34 (2012), no. 5, A2718–A2739. Search in Google Scholar
 A. Ostrowski, Les estimations des erreurs a posteriori dans les procédés itératifs, C. R. Acad. Sci Paris Sér. A–B 275 (1972), A275–A278. Search in Google Scholar
 S. Repin, A posteriori error estimation for variational problems with uniformly convex functionals, Math. Comp. 69 (2000), no. 230, 481–500. Search in Google Scholar
 S. Repin, A Posteriori Estimates for Partial Differential Equations, Walter de Gruyter, Berlin, 2008. Search in Google Scholar
 S. Repin, T. Samrowski and S. Sauter, Combined a posteriori modeling-discretization error estimate for elliptic problems with complicated interfaces, ESAIM Math. Model. Numer. Anal. 46 (2012), no. 6, 1389–1405. Search in Google Scholar
 S. Repin, S. Sauter and A. Smolianski, A posteriori estimation of dimension reduction errors for elliptic problems on thin domains, SIAM J. Numer. Anal. 42 (2004), no. 4, 1435–1451. Search in Google Scholar
 E. Zeidler, Nonlinear Functional Analysis and Its Applications. I: Fixed-Point Theorems, Springer, New York, 1986. Search in Google Scholar
© 2017 Walter de Gruyter GmbH, Berlin/Boston