RANK STRUCTURED APPROXIMATION METHOD FOR QUASIPERIODIC ELLIPTIC PROBLEMS
B. KHOROMSKIJ AND S. REPIN

arXiv:1701.00039v1 [math.NA] 31 Dec 2016

Abstract. We consider an iteration method for solving an elliptic type boundary value problem Au = f , where a positive definite operator A is generated by a quasiperiodic structure with rapidly changing coefficients (typical period is characterized by a small parameter ) . The method is based on using a simpler operator A0 (inversion of A0 is much simpler than inversion of A), which can be viewed as a preconditioner for A. 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 A and A0. For typical quasiperiodic structures, we establish simple relations that suggest an optimal A0 (in a selected set of "simple" structures) and compute the corresponding contraction factor. Further, this allows us to deduce fully computable twosided a posteriori estimates able to control numerical solutions on any iteration. The method is especially efficient if the coefficients of A 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 1/ , providing the FEM approximation of the order of O( 1+p), p > 0.
AMS Subject Classification: 65F30, 65F50, 65N35, 65F10 Key words: elliptic problems with periodic and quasiperiodic coefficients, precondition methods, tensor type methods, guaranteed error bounds

1. Introduction

Problems with periodic and quasiperiodic 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, 6] 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 quasiperiodic 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., Fig. 1.1, 1.2)
(1) periodic structures, in which the amount of cell is considerable (e.g. 103104) but not large enough to neglect the error generated by the respective homogenized model;
(2) quasiperiodic structures that contain cells with defects and deformations; (3) multiperiodic structures where the coefficients reflect combined effect of several func-
tions with different periodicity.
In general terms, the idea of the method is as follows. We consider the problem P

(1.1)

Au = f, f  V ,
1

2

B. KHOROMSKIJ AND S. REPIN

where V is a reflexive Banach space with the norm  V , V  is the space conjugate to V (the respective duality pairing is denoted by < v, v >), and A : V  V  is a bounded linear operator. It is assumed that the operator A is positive definite and invertible, so that the problem (1.1) is well posed. However, P is viewed as a very difficult problem because A is generated by a complicated physical structure, which may contain a huge amount 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 A is approximated by a simplified positive definite operator A and the inversion of A is much simpler than inversion of A. By means of A, we construct an iteration method based on solving a "simple" problem P0: Au = g. In other words, the method is based on the operation g  A- 1g. It also includes the operation v  Av, which can be performed very efficiently by tensor type decomposition methods provided that physical structures generated A 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 A and operations of the type v  Av. In the iteration methods and error estimates inversion of the operator A is avoided.
In the paper, we consider one class of problems associated with divergent type elliptic equations where A = QQ and A = QQ. Here  : Y  Y is a bounded operator induced by a complicated quasiperiodic structure while Q : V  Y and Q : Y  V  are conjugate operators, i.e.,

(1.2)

(y, Qw) =< Qy, w > y  Y and w  V,

where Y is a Hilbert space with the scalar product (, ) and the norm  . The operators Q and Q are induced by differential operators or certain finite dimensional approximations
of them. Henceforth, it is assumed that f  V, where V is a Hilbert space with the scalar product (, )V. This space is intermediate between V and V , i.e., V  V  V .
The operator A = QQ contains the operator  generated by a simplified structure. We assume that the operators  and  are Hermitiam (i.e., (y, z) = (y, z) and (y, z) = (y, z)) and satisfy the conditions

(1.3) (1.4)

 y 2  (y, y)   y 2 y  Y,  y 2  (y, y)   y 2,  < .

Then, the structural operators  and  are spectrally equivalent

(1.5)

c1(y, y)  (y, y)  c2(y, y),

where the constants are the minimal and maximal eigenvalues of the generalized spectral

problem

y - y

= 0.

Obviously,

they

satisfy

the

estimates

c1



 

and

c2



 

(which

may be rather coarse).

Concerning the operator Q, we assume that there exists a positive constant c such that

(1.6)

Qw  c w V w  V.

Generalized solutions of the problems P and P0 are defined by the variational identities

(1.7)

(Qu, Qw) =< f, w >  w  V,

RANK STRUCTURED METHOD

3

and

(1.8)

(Qu0, Qw) =< f , w >  w  V.

In Sect. 2, we show that a sequence {uk} converging to u in V can be constructed by solving
problems (1.8) with specially constructed right hand sides fk generated by the residual of (1.7). In proving convergence, the key issue is analysis of the spectral radius of the operator

(1.9)

B := I - - 1,

and selection of such relaxation parameter  that provides the best convergence rate. More-

over, iteration procedures of such a type become contracting if the iteration parameter is

properly selected. This fact is often used in proving analytical results (e.g., see [25], where

classical results on existence and uniqueness of a variational inequality has been established

by contraction arguments ). Also, these ideas were used in construction of various numer-

ical methods (see, e.g., [10]). 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" . These questions are studied in Sect. 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

selection of a suitable  can be done 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 saying about 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 QHQ uH = f , where H is defined by means of an auxiliary problem with periodical boundary conditions in the cell of periodicity. The

respective solution uH contains an irremovable (modeling) error depending on the cell diam-
eter . Moreover, if tends to zero, then typically uH converges to u only weakly (e.g., in L2). Getting a better convergence (e.g., in H1) requires certain corrections, which lead to

other (more "corrected"

complicated) boundary value solution ucH also contains an

problems in the cell of periodicity. The respective 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  := H 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 weather or not the choice  := H always yields minimal value of the contraction factor. In Sect. 3, we briefly discuss this question and present an example of that the best  may differ from H .
In Sect. 4, we deduce a posteriori estimates that provide fully computable and guaran-

teed estimates of the distance to the exact solution u for any numerical approximation uk,h

4

B. KHOROMSKIJ AND S. REPIN

computed for an approximation subspace Vh. These estimates are established by combining functional type a posteriori estimates (see [31, 29, 32] and references cited therein) and estimates generated by the contraction property of the iteration method (see [30, 37]).
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 quasiperiodic 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 [11]. In computational mechanics this method is known as the KantorovichKrylov (or extended Kantorovich) method. However, it is rarely used in modern numerical technologies, which are mainly based upon various finite element technologies. In part, this is due 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. Certainly, this fact induces some limitations, which however could be bypassed by known methods (coordinate transformation, domain decomposition, iso-geometric analysis, etc.).
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 [19, 21, 18]. 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 constant, modulated oscillation around given smooth function, or oscillation around piecewise constant function.

1.5 2

1.5 1
1

0.5

0.5

500

1000

1500

2000

500

1000

1500

2000

Figure 1.1. Examples of periodic and modulated periodic coefficients in 1D.

RANK STRUCTURED METHOD

5

Figure 1.1 (1D case) represents examples of highly oscillating (left) and modulated periodic coefficients (right) functions.
Figure 1.2 (2D case) illustrates the well separable equation coefficient obtained by a sum of step-type and uniformly oscillating functions.

Figure 1.2. An example of modulated piecewise periodic coefficients in 2D.

We show that specially constructed FEM type approximations of PDEs with slightly perturbed or regularly modulated periodic coefficients on d-fold n      n tensor grids in

Rd may lead to the discretized algebraic equations with the low Kronecker rank stiffness matrix of size nd  nd, where n = O( 1 ) is proportional to the large frequency parameter 1/ .

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

the

only

O(dn)

storage

size

instead

of

O(nd)

=

O(

1
d

)

complexity

representations

which

are

mandatory for the traditional FEM techniques (the latter quickly leads to the bottleneck in

case of small parameter > 0).

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, i.e., they are estimated by

O(dn logp( 1 ))

1 O( ), p > 0,

d

where d is the spatial dimension. Numerical examples in Section 5 demonstrate the stable geometric convergence of the preconditioned CG (PCG) iteration with the preconditioner  and confirm the 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 [20] to functions discretized on large tensor grids of size proportional to the frequency parameter, i.e. n = O(1/ ), as it was demonstrated in the previous paper [23] for the case d = 1. The use of tensor-structured preconditioned iteration with the adaptive QTT rank truncation may lead to the logarithmic complexity in the grid size, O(logp n), see [19, 21, 26] for the rank-truncated iterative methods, [15, 16, 14, 13] for various examples of the QTT tensor approximation to lattice structured systems, and [2] for tensor approximation of complicated functions with multiple cusps in Rd.
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.

6

B. KHOROMSKIJ AND S. REPIN

2. The iteration method

Let v  V and   R+. Consider the problem: find uv such that

(2.1)

(Quv, Qw) = v(w) -  v(w) w  V,

where

v(w) := (Qv, Qw)- < f, w >

and

 v

(w)

:=

(Qv,

Qw).

Obviously, the right hand side of (2.1) is a bounded linear functional on V , so that this problem has a unique solution uv. Thus, we have a mapping T : V  V , which becomes a contraction if the parameter  is properly selected. Indeed, for any v1 and v2 in V , we obtain

(2.2)

(Q, Qw) = (Q - Q, Qw)) w  V,

where u1 = Tv1, u2 = Tv2,  := v1 - v2, and  := u1 - u2. Hence

(2.3)



2 

:=

(Q,

Q)

=

( Q ,

Q)

-

(Q ,

Q)

= (Q, Q) - (- 1Q, Q) = (Q - - 1Q, Q)

   Q -  Q, Q -  - 1Q 1/2 .

From (2.3) we find that

(2.4)



2 



(Q, Q) - 2(Q, Q) + 2(- 1Q, Q)

= (Q, Q) - 2(- 1Q, Q) + 2(- 1- 1Q, Q)

= ((I - 2- 1 + 2- 1- 1)Q, Q) = (B2Q, Q)

 (B2Q, B2Q)1/2(Q, Q)1/2,

where B is defined by (1.2). If  is selected such that

(2.5)

(B2Q, Q)  q2  2, for some q < 1,

then (2.4) shows that T is a contractive mapping. It is not difficult to show that  satisfying (2.5) can be always found. Indeed, in view of
(1.5)

(2.6) (B2Q, Q) = (Q, Q) - 2(Q, Q) + 2(- 1Q, Q)  (1 - 2c1)(Q, Q) + 2(- 1Q, Q).

Since  and  are invertible with trivial kernels,  and y are an eigenvalue and the respective eigenfunction of y = y if and only if they are an eigenvalue and the eigenfunction of the problem - 1y = y. This means that
c1(y, y)  (- 1y, y)  c2(y, y)  c22(y, y).

Hence

(- 1Q, Q)  c22



2 

and (2.6) implies

(2.7)

(B2Q, Q)  1 - 2c1 + 2c22  2.

RANK STRUCTURED METHOD

7

Minimum

of

the

expression

in

round

brackets

is

attained

if



=



:=

. c1
c22

For



=

,

we

find

that

(2.8)

q2

:=

1

-

c21 c22



q^2

:=

1

-

2 2 22



[0, 1).

Hence, T is a contractive mapping with explicitly known contraction factor q. Well known results in the theory of fixed points (e.g., see [37]) yield the following result.

Theorem 2.1. For any u0  V and  =  the sequence {uk}  V of functions satisfying the relation

(2.9) (Quk, Qw) = (Quk-1, Qw) -  (Quk-1, Qw)- < f, w >

w  V

converges to u in V and uk - u   qk u0 - u  as k  +.

Remark 2.2. From (2.4) we obtain



2 



1 0,min

B2Q







B2 20,min



2.

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.4) by a different

method. Let |B| denote the operator norm

(2.10)

|B| := sup yY

By  . y

Then By   |B| y  and

(B2y, y)  |B|2

y

2 

.

Hence, (2.4) yields the estimate

(2.11)

   |B|  ,

which shows that T is a contraction provided that

(2.12)

|B| < 1.

In applications B is a self adjoint bounded operator acting in a finite dimensional space, so that verification of this condition amounts finding  which yields the respective spectral radius of B (see Section 4).

3. Selection of 
In this section, we discuss how to select  in order to minimize q what is crucial for two major aspects of quantitative analysis: convergence of the iteration method and guaranteed a posteriori estimates. We assume that V , V, and Y are spaces of functions defined in a Lipschitz bounded domain  (namely y(x)  T for a.e. x   where T may coincide with R, Rd, or Mdd) and the operators  and  are generated by bounded scalar functions, matrices or tensors. In this case,

(y, y) := (x)y y dx, and (y, y) := (x)y y dx,





8

B. KHOROMSKIJ AND S. REPIN

where denotes the respective product of scalar, vector, or tensor functions. In view of

(2.10)

and

(2.12),

the

value

of



should

minimize

the

quantity

sup
yY

. (By,By) (y,y)

This

procedure

yields the contraction factor

(3.1)

(x)B(x)y

q2 = Q(, ) := inf sup   yY

(x)y



B(x)y dx ,
y dx

which computation is reduced to B(x) to solving algebraic problems at a.e. x  , i.e.,

(3.2)

Q(, )

:=

inf 

sup
x

sup
 T

(x)B(x) (x)

B(x) 

Let S be a certain set of "simple" operators defined a priori (e.g., it can be a finite dimensional set formed by piece vise constant or polynomial functions). Then, finding the
best "simplified" operator amounts solving the problem: find   S such that Q(, ) is minimal. In other words, optimal  is defined by the problem

(3.3)

inf
0 S,

sup
x

(x)B(x) (x)

R  T

B(x) 

= q2.

Notice that (3.3) 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.9) and yields the respective contraction factor. Below we discuss some particular cases, where analysis of this problem generates optimal (or almost optimal) .
Problem (3.3) is explicitly solvable if  and  have a special structure, namely,

 = a(x)I,  = a(x)I,

where I is the unit operator and a(x) and a(x) are positive bounded functions defined in . Then,

B(x) = (1 - h(x))I,

a(x) h(x) :=
a(x)

and

(1 - h(x))2 

sup

 T

| |2

 = |1 - h(x)|2

 x  .

Define h := min h(x) and h := max h(x). It is not difficult to show that

x

x

sup |1 - h(x)| = max{|1 - h |, |1 - h|}.
x

Minimization with respect to  yields the best value  =

h

2 +h

and the respective value

(3.4)

Q(, ) =

h - h h + h

2
=

1 - J (a, a)

2
< 1,

1 + J (a, a)

h

J

(a, a)

=

. h

RANK STRUCTURED METHOD

9

In accordance with (3.3) identification of the optimal simplified problem is reduced to the problem

(3.5)

sup J (a, a).
a0S

where S is a given set of functions.

We illustrate the above relations by means of several examples.

Example 1. Constant coefficients. In the simplest case, we set S = P 0, i.e., a0 is a constant.

From

(3.5)

it

follows

that

q

=

a-a ,
a+a

where

a

:=

min a(x)
x

and

a

:=

max a(x).
x

Then



=

2a0 a+a

and the iteration procedure (2.9) with  =  has the form

(3.6)

Quk


Qw dx =


2a

1- a+a

Quk-1

2

Qw dx +

f w dx

a+a



From Theorem 2.1, it follows that

|Q(uk - u)|2 dx  C

a - a 2k .
a+a



Example 2. Oscillation around a given function. Consider a somewhat different example. Let a(x) be a function oscillating around a certain mean function g(x) so that

a(x)  [1 - , 1 + ],  (0, 1).
g(x)

If g is a relatively simple function, then it is natural to set a(x) = g(x). By (3.4), we find that h = 1 + , h = 1 - , and q = . Hence the method is very efficient for small

(i.e., if a oscillates around g with a relatively small amplitude). Figures 1.1 and 1.2 illustrate

three examples of quasi-periodic coefficients a and respective a corresponding to the case of oscillation around constant with smooth modulation, oscillation around given smooth

function, or oscillation around piecewise constant function.

Example 3. Piecewise constant coefficients. Consider a more complicated case, where 

is divided into N nonoverlapping subdomains i and (x) = ciI if x  i. Define the

numbers

a(i)

:=

max
xi

a(x),

a(i)

:=

min
xi

a(x),

h = min

a(1) a(2)

a(N )

, , ...,

c1 c2

cN

,

and

h = max

a(1) , a(2) , ..., a(N)

c1 c2

cN

.

Since the constantans ci are defined up to a common multiplier, we can without a loss of generality assume that

(3.7)

(N )
i = 1,
i=1

where

1

i

=

. ci

In accordance with (3.5), maximum of Q(, ) is attained if

(3.8)

min 1a(1), 2a(2), ..., N a(N) max 1a(1), 2a(2), ..., N a(N)

 max,

10

B. KHOROMSKIJ AND S. REPIN

where i > 0 and satisfy (3.7). If N = 2, then the problem (3.8) has a simple solution, which

shows

that

the

ratio

1 2

(i.e.,

) c2
c1

can

be

any

in

the

interval

[1, 2],

where

1

=

min{

a(2) a(1)

,

} a(2)
a(1)

and

2

=

max{

a(2) a(1)

,

}. a(2)
a(1)

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 1-dimensional

problem

(au ) - f = 0 in (0, 1)

with

a(x) = a(1)(x) a(x) = a(2)(x)

in 1 = (0, ), in 2 = (, 1),

  (0, 1),

where a(1)(x) is a perfectly periodical function attaining only two values a(1) (Lebesgue measure of this set is 1|1|, 1  (0, 1)) and a(1) (Lebesgue measure of this set is (1-1)|1|). Similarly, a(2)(x) is a perfectly periodical function attaining only two values a(2) (Lebesgue measure of this set is 2|2|, 2  (0, 1)) and a(2) (Lebesgue measure of this set is (1-2)|2|).
Assume that the amount of periods is very large and, therefore, the homogenization method

can be successfully applied. The corresponding homogenized problem has the following

coefficients



-1

a(1) :=  1 

1 a(1)(x) dx

0

in 1



1

-1

and a(2) :=  1 1-

1 a(2)(x) dx



in 2.

It is easy to see that

a(1)

=

1a(1)

a(1)a(1) + (1 - 1)a(1)



(a(1), a(1)),

a(2)

=

2a(2)

a(2)a(2) + (1 - 2)a(2)



(a(2), a(2))

Hence

a(2)

min{a(1), a(2)}

max{a(1), a(2)}

a(1)  (1, 2), where 1 := max{a(1), a(2)} , 2 := min{a(1), a(2)} .

It is clear that 1  1 and 2  2. Therefore, homogenized coefficients may not generate the best piece wise constant a, which produces the smallest contraction factor q.

4. Error estimates
4.1. General estimate. Since T is a contractive mapping, we can use the Ostrowski estimates (see [30, 37, 32]), which yield the estimate of the distance between v  V and the fixed point:

(4.1)

v-u  

, 1 + q() 1 - q()

,

where := Tv - v .

The is estimate cannot be directly applied because v := Tv is generally unknown (it is the exact solution of a boundary value problem). Instead, we must use a numerical approximation v (in our analysis, we impose no restrictions on the method by which the

RANK STRUCTURED METHOD

11

function v  V was constructed). Thus, the difference  := v - v is a known function and the quantity  =   is directly computable. It is easy to see that

(4.2)

 - v - v   v - v    + v - v .

To deduce a fully computable majorant of the norm v - v  we use the method suggested in [31, 32]. First, we rewrite (2.1) in the form

(4.3)

(Qv, Qw) = (Qv, Qw) -  (Qv, Qw)- < f, w > .

For any y  Y and w  V0, we have

(4.4) (Q(v - v), Qw) = (Q(v - v), Qw)) -  (Qv, Qw)- < f, w > = (Q(v - v) - Qv + y, Qw))- < Qy + f, w > .
We estimate the first term in the right hand side of (4.4) as follows:

(Q(v - v) - Qv + y, Q(v - v))) = (Q(v - v) - - 1Qv + - 1y, Q(v - v)))  Q(v - v) + , Q(v - v) + - 1 1/2 v - v ,

where  := y - Qv. The second term meets the estimate

< Qy + f, v - v >  |Qy + f | v - v



1 ( )1/2

|Qy

+

f

|

v - v

,

where

|w|

=

sup
wV

<w,w> w

is

the

dual

norm.

Hence,

(4.5)

v - v   Q + , Q + - 1 1/2 +

1 

|Qy

+

f |

=:

M(,

 ).

Notice that

inf M(,  ) = v - v .
yY
Indeed, set y = Q(v - v) + Qv. Then,  = Q(v - v). In view of (4.3), Qy + f = 0, and the majorant is equal to v - v 2. Hence, the estimate (4.5) has no gap.
It is worth noting that computation of the majorant M does not require inversion of the operator  associated with a complicated quasiperiodic problem.

Remark 4.1. M(,  ) is an a posteriori error majorant of the functional type (its derivation is performed by purely functional methods based on generalized formulation of the boundary value problem and special properties of approximations or numerical method are not used). Properties of such type error majorants are well studied (see [31, 32] and the literature cited therein). It is not difficult to show that the last term of M(,  ) can be estimated via an explicitly computable quantity provided that y has the same regularity as the true flux. However, in our subsequent analysis these advanced forms of the majorant are not required. Therefore we omit this discussion (interested reader can find the respective analysis in [32]). Numerous tests performed for different boundary value problems have confirmed high practical efficiency of error majorants of the functional type. It was shown that M 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

12

B. KHOROMSKIJ AND S. REPIN

different ways to obtain suitable reconstructions (see [27] for a systematic discussion of computational aspects of this error estimation method). Error majorants of this type can be also used for the evaluation of modeling errors (see [35, 34]).

Now, (4.1), (4.2), and (4.5) yield the following result

Theorem 4.2. The error e = v - u is subject to the estimate

(4.6)

e

max

0,  - M(,  ) 1 + q()

,  + M(,  ) 1 - q()

,

where  := y - Qv and y is a function in Y and M is defied by (4.5). If Qy + f = 0 then M2 (,  ) = (Q, Q) + (- 1,  ) - 2(Q,  ).

4.2. Examples. Now we shortly discuss applications of Theorem 4.2 to problems, where Q
and Q are defined by the operators  and div, respectively,  = a(x)I,  = a(x)I, x  ,

and V =H 1().

4.2.1. d = 1. Let  = (0, 1). The equation (1.1) has the form (a(x)u ) - f = 0. In this case, Qw = w , Qy = -y , and (4.3) is reduced to

1

1

(4.7)

a(v - v) w dx +  (av w + f w) dx = 0.

0

0

In order to apply Theorem 4.2, we set y = (g(x) + ), where g(x) = -

x 0

f dx

and



is

a

constant. Then -y - f = 0 and  = (g(x) + ) - av = ( + g - av ). The best constant

 is defined by minimization of M2 (,  ), which has the form

1

(a()2 + a- 12( + g - av )2 - 2( + g - av )dx

0

Since

1 0

dx

=

0,

the

problem

is

reduced

to

minimization

of

the

second

term

and

the

best

1

 satisfies the equation a- 1( + g(x) - av )dx = 0. Hence

0

 =  :=

1 0

a- 1(av

1 0

a- 1

-g dx

)dx

,

and (4.6) yields the estimate

(4.8)

e

max

0,  - I(v, v) 1 + q()

,  + I(v, v) , 1 - q()

where

1 2
I2 (v, v) = a- 1 a(v - v) - ( + g - av ) dx.

0
Here v and v are two consequent numerical approximations (e.g., finite element approximations vhk and vhk+1 computed on a mesh Ih. Then
 = hk := vhk - vhk+1 and  = k := vhk - vhk+1 

RANK STRUCTURED METHOD

13

are directly computable. Since a is a "simple" function, the integrals

1

1

1

1

F1 = a- 1 dx, F2 = a- 1 g dx, F3 = a hk 2 dx, F4 = a ( + g)2 dx, F5 =

0

0

0

0

are easy to compute. Other integrals

1
f hk dx
0

1

1

1

1

G1 = a- 1a vhk dx, G2 = a (vhk) hk dx, G3 = ( + g)a- 1a vhk dx, G4 = a- 1 a2 vhk ) 2 dx

0

0

0

0

contain highly oscillating coefficient a multiplied by piece wise polynomial mesh functions. If a has a low QTT rank tensor representation [20], then the integrals can be efficiently computed by tensor type methods already discussed in [23]. We have

I2 (v, v) = F3 + 2G2 + 2F5 + 2(F4 - 2G3 + G4) =: k,

 = G1 - F2 . F1

Here  =

h

2 +h

is selected in accordance

with Section

3.

The respective contraction factor

is

q

=

. h-h
h +h

Now (4.8) yields easily computable lower and upper bounds of the error

encompassed in vhk:

k - k 
1+q

vhk - u





k + k 1-q

4.2.2. d = 2. Computation of M 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

f = f (1)(x1)f (2)(x2) and a = a(1)(x1)a(2)(x2).

Assume that approximations are represented in the form of series formed by one dimensional functions (i1) and (j2) (which may be supported locally or globally), so that

n1 n2

v=

ij (i1) (x1 )(j2) (x2 ),

i=1 j=1

n1 n2

v =

ij (i1) (x1 )(j2) (x2 ).

i=1 j=1

In this case,

 =

n1 i=1

n2 j=1

ij

(i1) x1

(j2)

,

n1 i=1

n2 j=1

ij (i1)

(j2) x2

,

where ij = ij - ij.

We define another set of one dimensional functions Wk(1)(x1) and Wl(2)(x2), which form the vector function

(4.9)

m1 m2

y = 0 +

klkl,

k=1 l=1

kl =

Wk(1)

Wl(2) x2

;

-

Wk(1) x1

Wl(2)

.

Here 0 is a given function, which can be defined in different ways. In particular, we set

0 =

W0(1)(x1)W0(2)(x2) ; 0

, W0(1)(x1) =

x1 0

f (1)dx1

and

W0(2)

=

-f (2).

The

functions

kl

must satisfy the usual linear independence conditions in order to guarantee unique solvability

14

B. KHOROMSKIJ AND S. REPIN

of the respective approximation problem. For any smooth function w vanishing on , we have

(0  w - f w)dx1dx2 = 0 and kl  wdx1dx2 = 0.





Thus, |Qy + f| = |divy - f| = 0 and we can use the simplified form of M. In the simplest case  = aI, where a is a constant. The best y minimizes the quantity

(4.10) M2 (,  ) = a  dx + a- 1y  ydx + 2 a- 1a2v  vdx







- 2 (a- 1av + )  y dx + 2 a   dx,





which shows that y must satisfy the relation y = av + a. We select kl that defines Galerkin approximation of this function and arrive at the system

m1 m2

(4.11)

kl kl  stdx1dx2 + 0  stdx1dx2

k=1 l=1





n1 n2

=

(aij + aij)

i=1 j=1 

(i1) x1

(j2)

,

(i1)

(j2) x2

Introduce the following matrixes

 stdx1dx2

D(1) = Dk(1l) , D(2) = Dk(2l) ,

Dk(1l) = Dk(2l) =

a 0

Wk(1) x1

Wl(1) x1

dx1,

b 0

Wk(2) x2

Wl(2) x2

dx2,

W(1) = W(2) =

Wk(l1) Wk(l2)

a

, Wk(l1) =

Wk(1)Wl(1) dx1,

0

b

, Wk(l2) =

Wk(2)Wl(2) dx2,

0

F(1) = Fi(k1) ,

Fi(k1) =

a 0

(i1) x1

Wk(1)dx1

,

F(2) = Fj(l2) ,

Fj(l2) =

b 0

(j2)

Wl(2) x2

dx2,

G(1) = G(ik1) , G(2) = G(j2l ) ,

G(ik1) =

a

(i1)

Wk(1) x1

dx1,

0

G(j2l ) =

b

(j2) x2

Wl(1)

dx2,

0

F(1) =

Fi(k1)

, Fi(k1) =

0

a

a1(x1)

(i1) x1

Wk(1)dx1

,

G(1) =

G(ik1)

, G(ik1) =

a 0

a1(x1)(i1)

Wk(1) x1

dx1,

F(2) =

Fj(l2)

, Fj(l2) =

0

b

a2(x2)(j2)

Wl(2) x2

dx2

,

G(2) =

G(j2l )

, G(j2l ) =

b 0

a2(x2)

(j2) x2

Wl(1)

dx2.

and vectors

g(1) = {gk(1)},

a

gk(1) =

W0(1)Wk(1) dx1,

0

g(2) = {gl(2)},

gl(2) =

b 0

W0(2)

Wl(2) x2

dx2.

Notice that all coefficients are presented by one dimensional integrals, which can be efficiently

computed with the help of special (tensor type) methods (see, e.g., [20]-[24]).

RANK STRUCTURED METHOD

15

It is not difficult to see that

Yklst := kl  st dx = Wk(s1)Dl(t2) + Dk(1s)Wl(t2)



and

0  stdx1dx2 =

W0(1)

Ws(1)W0(2)

Wt(2) x2

dx1dx2

=

gs(1)gt(2),





where Y = {Yklst} is the fourth order tensor. Hence the left hand side of the system (4.11) has the form Y + g(1)  g(2). In the right hand side we have the term

aij


(i1) x1

(j2),

(i1)

(j2) x2

 stdx1dx2 = aH,

where H = {Hijst}, Hstij = Fi(s1)Fj(t2) - G(is1)G(j2t). Another term is

aij


(i1) x1

(j2),

(i1)

(j2) x2

 stdx1dx2 = H,

where H = {Hijst}, Hstij = Fi(s1)Fj(t2) - G(is1)G(j2t). Now (4.11) implies  = Y-1(H + aH - g(1)  g(2)) and the value of M is obtained by
(4.6), (4.9), and (4.10).

5. Low-rank solution of the discrete equation

In what follows we assume that f and a admit low rank representation (e.g.,

f=

Rf i=1

f1i

(x1

)f2i(x2

),

a

=

Ra j=1

aj1(x1

)aj2(x2

)).

Then one may assume that the ex-

act FEM solution can be well approximated by uK(x) =

K j=1

uj1(x1)uj2(x2),

where

K

depends on the separation rank of f and a. In some cases this important property can be

rigorously proven (say, for Laplacian like operators). The similar low rank approximation

can be observed for the QTT tensor approximation (see [23]). Existence of low rank solution

means that for some K we have uK  u up to the rank truncation threshold. Here we sketch the rank-structured computational scheme. In our set of examples the

original problem: find u such that

(5.1)

a(x)u  w dx = f wdx w  V0 := H01





is replaced by the Galerkin problem for low rank representations

(5.2)

a(x)uK  wK dx = f wdx wK  V0K,





where V0K is a subset of V0 formed by functions of the type

K
wK(x) = j1(x1)j2(x2).
j=1

16

B. KHOROMSKIJ AND S. REPIN

Therefore, in terms of the general scheme exposed in the introduction, the Problem P is now the problem (5.2) and we solve it by iterations with the help of simplified (preconditioned) problem

(5.3)

a(x)uk  wK dx = fk-1wdx wK  V0K,





where a is a simple (mean) function and fk-1 depends on uk-1. Given the right-hand side, the problem (5.3) is much simpler than the initial equation since
the matrix , generated by the coefficient a is easily invertible. Moreover, the coefficient a may be rather complicated and admits a representation with rank R, i.e.,

R
a(x) = a(1s)(x1)...a(ds)(xd),
s=1
where R is a small integer. When we construct the low-rank Kronecker representation of
stiffness matrix for this a, which is presented by elements of 4R matrices computed by only 1D integrals containing oscillating functions a(is)(xi).
If we use (5.3), then a is a simple function, it may be a even a constant, or a function representable in the form a1(x1)...ad(xd) with very simple multipliers. Then, the respective Kronecker stiffness matrix  is computed much easier and has a simple (low rank) form that allows the low rank representation of its inverse.

5.1. Kronecker product representation of the stiffness matrix. We consider the elliptic diffusion equation with quasi-periodic coefficient a(x) > 0 (whose oscillations are characterized by the parameter )
(5.4) Au = -div(a(x)u) = f (x), x = (x1, . . . , xd)   = (0, 1)d, u| = 0,
where the function f corresponds to the modified right hand side in the problem (4.3),  = , and the right-hand side f (x1, . . . , xd) can be represented with a low separation rank.
Figure 5.1 illustrates a 2D example of L  L periodic coefficient with L = 6 corresponding to the choice = 1/L. In this example, the scalar coefficient is represented by the separable

Figure 5.1. Example of the 2D periodic oscillating coefficients (left) and the 1D
factor a1(x1).
function a(x) = C + a1(x1)a1(x2), C > 0, where the generating univariate function a1(x1)

RANK STRUCTURED METHOD

17

has the shape of six uniformly distributed bumps of hight 1 as shown in Figure 5.1, right. Figure 5.1, left, presents the oscillating part of 2D coefficients function, a1(x1)a1(x2).
The examples of other possible shapes of the equation coefficient corresponding to the cases (1), (2) and (3) specified in Introduction are presented in Figures 1.1 and 1.2.
We apply the FEM Galerkin discretization of equation (5.4) by means of tensor-product piecewise affine basis functions (instead of "linear finite elements")

{i(x) := i1(x1)    id(xd)}, i = (i1, . . . , id), i  I = {1, . . . , n }, = 1, . . . , d,
where ik are 1D finite element basis functions (say, piecewise linear hat functions). We associate the univariate basis functions with the uniform grid {j}, j = 1, . . . , n , on
[0, 1] with the mesh size h = 1/(n + 1). In this construction we have N = n1n2...nd basis functions i. Notice that the univariate grid size n is of the order of n = O(1/ ) designating the total problem size N = O(1/ d).
For ease of exposition we, first, consider the case d = 2, and further assume that the scalar diffusion coefficient a(x1, x2) can be represented in the form
R
a(x1, x2) = a(k1)(x1)a(k2)(x2) > 0
k=1
with a small rank parameter R. The N  N stiffness matrix is constructed by the standard mapping of the multi-index i
into the N -long univariate index i representing all degrees of freedom. For instance, we use the so-called big-endian convention for d = 3 and d = 2

i  i := i3 + (i2 - 1)n3 + (i1 - 1)n2n3, i  i := i2 + (i1 - 1)n2,

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 i is designated via the long index, i.e. i = i.
First, we consider the simplest case R = 1 and let d = 2. We construct the Galerkin stiffness matrix A = [aij]  RNN in the form of a sum of Kronecker products of small "univariate" matrices. Recall that given p1  q1 matrix A and p2  q2 matrix B, their Kronecker product is defined as a p1p2  q1q2 matrix C via the block representation

C = A  B = [aijB], i = 1, . . . , p1, j = 1, . . . , q1.

We say that the Kronecker rank of the matrix A in the representation above equals to 1. Now the elements of Galerkin stiffness matrix take a form

18

B. KHOROMSKIJ AND S. REPIN

(5.5) aij = Ai, j = a(1)(x1)a(2)(x2)i(x) j(x)dx



1

1

=

a(1)

(x1

)



i1 (x1 ) x1



j1 (x1 ) x1

dx1

a(2)(x2)i2 (x2)j2 (x2)dx2

0

0

1

1

+

a(1)(x1 )i1 (x1)j1 (x1 )dx1

a(2)(x2

)



i2 (x2 x2

)



j2 (x2 x2

)

dx2

,

0

0

which leads to the rank-2 Kronecker product representation

A = [aij] = A1  M2 + M1  A2,

where  denotes the conventional Kronecker product of matrices. Here A1 = [ai1j1]  Rn1n1 and A2 = [ai2j2]  Rn2n2 denote the univariate stiffness matrices and M1 = [mi1j1]  Rn1n1 and M2 = [mi2j2]  Rn2n2 define the corresponding weighted mass matrices, e.g.,

1

ai1j1 =

a(1)

(x1)



i1 (x1 x1

)



j1 (x1 x1

)

dx1,

0

1
mi1j1 = a(1)(x1)i1 (x1)j1 (x1)dx1.
0

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 A can be simplified to the form

(5.6)

A  A = A1  D2 + D1  A2,

where D1, D2 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 A corresponds to the operator Laplacian. In this case the representation in (5.6) is simplified to the discrete Laplacian matrix in the form of rank-2 Kronecker sum

(5.7)

A   = A1  I2 + I1  A2,

where I1 and I2 denote the identity matrices of the corresponding size. This matrix will be used in what follows as a prototype preconditioner for solving the linear system of equations

(5.8)

Au = f .

The matrix A is constructed in general for the R-term separable coefficient a(x1, x2) with R  1 which leads to the rank-2R Kronecker sum representation

R
A = [A1,k  D2,k + D1,k  A2,k],
k=1
with matrices of the respective size.

RANK STRUCTURED METHOD

19

5.2. Existence of the low-rank solution. In this paper we discuss the approach based on the low rank separable -approximation of the solution to the equation (5.8) that is considered as the d-dimensional real valued array, u  Rn1׷nd. In general, for the case R > 1 this favorable property is not guaranteed by the low Kronecker rank representation to the Galerkin system matrix A, discussed in the previous section.
Let R = 1 and d = 2, the existence of the low rank approximation to the solution of the equation (5.8) with the low-rank right-hand side

Rf

f=

fk(1)  fk(2),

k=1

fk( )  Rn ,

and with the system matrix in the form (5.7) can be justified by plugging the representation (5.7) in the sinc-quadrature approximation to the Laplace integral transform [8]

(5.9)

- 1 =

M

M

e-tdt  BM :=

cke-tk =

cke-tkA1  e-tkA2 ,

R+

k=-M

k=-M

taking into account that the matrices A1 and A2 commute with I1 and I2, respectively. Hence, the equation (5.9) represents the accurate rank-(2M + 1) Kronecker product approximation to - 1 which can be applied directly to the right-hand side to obtain

M

Rf

u = - 1f  BM f =

ck

e-tkA1 fm(1)  e-tkA2 fm(2).

k=-M m=1

The numerical efficiency of the representation (5.9) can be explained by the fact that the quadrature parameters tk, ck can be chosen in such a way that the low Kronecker rank approximation BM converges to - 1 exponentially fast in M . For example, under the choice tk = ekh, ck = htk with h = / M there holds [8]

- 1 - BM  Ce- M - 1 ,
which means that the approximation error > 0 can be achieved with the number of terms RB = 2M + 1 of the order of RB = O(| log |2).
Figures 5.2 and 5.3 demonstrate the singular values of the discrete solution on the n  n grid for n = 95, 143, 191 indicating very moderate dependence of the -rank on the grid size n. As in the case of Figure 5.1, in above figures we represent the only oscillating part of the coefficients and omit the small constant C > 0.

L=12

n=95

1

0
10

n=143

n=191

0

-1 1
0.5

00

10-5

1

-10

0.5

10

5

10

15

Figure 5.2. Rank decomposition of the solution for 12  12 periodic coefficient.

20

B. KHOROMSKIJ AND S. REPIN

1

0

-1 1
0.5

00

n=95

n=143

10 0

n=191

10 -5

1

0.5

10 -10

5

10

15

20

Figure 5.3. Rank decomposition of the solution for 12  12 modulated periodic coefficient.
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 [23] to the 1D equations with quasi-periodic coefficients. The power of QTT approximation method is due to the perfect low rank decompositions applied to the wide class of function-related tensors [20], see [23] for the 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 piecewise constant function, see Figure 1.1 and examples in [23].
Let the vector x  CN , N = 2L, be obtained by sampling a continuous function f  C[0, 1] (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 [20]: (A) r = 1 for complex exponentials, f (x) = eix,   R. (B) r = 2 for trigonometric functions, f (x) = sin x, f (x) = cos x,   R. (C) r  m + 1 for polynomials of degree m. (D) For a function f with the QTT-rank r0 modulated by another function g with the QTTrank r (say, step-type function, plain wave, polynomial) the QTT rank of a product f g is bounded by a multiple of r and r0,
rankQT T (f g)  rankQT T (f )rankQT T (g).
(E) Furthermore, the following result holds ([15]): 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.
The rank of the QTT tensor representation to the 1D Galerkin FEM matrix in the case of oscillating coefficients was discussed in [14, 23].
5.3. Numerical test on the rank decomposition of u. Figure 5.4 represents the righthand side f1(x1, x2) and the respective solution for the discretization to equation (5.4) (with the coefficient depicted in Figure 5.1) on 400  400-grid, where
f1(x1, x2) = sin(2x1) sin(2x2).
The PCG solver for the system of equations (5.8) with the discrete Laplacian inverse as the preconditioner demonstrates robust converges with the rate q 1. Next example demonstrates the rank behavior in the singular value decomposition (SVD) of a matrix representing the solution vector u  Rn1n2 to the equation (5.8) with 12  12 periodic

RANK STRUCTURED METHOD

21

1
0.5
0
-0.5
-1 1
0.8 0.6 0.4 0.2 00

1 0.8 0.6 0.4 0.2

4
2
0
-2
-4 1
0.8 0.6 0.4 0.2 00

1 0.8 0.6 0.4 0.2

Figure 5.4. The right-hand side and solution for periodic oscillating coefficients
shown in Figure 5.1.

coefficient shown in Figure 5.2, left. Figure 5.5 represents the rank behavior in the SVD decomposition of the solution in the case of 8  8 periodic coefficient.
It is worth to observe that comparison of Figures 5.2 and 5.5 indicates that the exponential decay of the approximation error in the rank parameter is stable with respect to the size of L  L lattice structure of the coefficient, i.e. the behavior of the singular values remains almost the same for different parameters = 1/L.

1
0.5
0
-0.5
-1 1 0.5

n=63

n=95

0

10

n=127

n=197

-5
10

00

1

-10

0.5

10

5

10

15

20

25

Figure 5.5. Accuracy of the rank decomposition of the solution vs. rank param-
eter for 8  8 periodic coefficient and grid size n  n.

Our iterative scheme includes only the matrix-vector multiplication with the stiffness matrix A that has the small Kronecker rank 2R, and the action of the preconditioner defined by the approximate inverse to the Laplacian type matrix. The latter has low Kronecker rank of order RB = O(| log |2) as shown above.
Given rank-1 vector u = u1 u2, the standard property of the Kronecker product matrices
Au = A1u1  M2u2 + M1u1  A2u2,
indicates that the matrix-vector multiplication enlarges the initial rank by the factor of 2 and similar with action of 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.

22

B. KHOROMSKIJ AND S. REPIN

Remark 5.1. Notice that for d = 3 the transformed matrix A takes a form A = A1  I2  I3 + I1  A2  I3 + I1  I2  A3,
and it obeys the d-term Kronecker sum representation in the general. Hence, in the general case of d  2 and R  1 the Kronecker rank of the matrix A is given by
rankKron(A) = d R.

6. Conclusions
We present a preconditioned iteration method for solving an elliptic type boundary value problem in Rd with the operator generated by a quasiperiodic 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 low-rank Kronecker sum. The preconditioner A0 is constructed based on certain averaging (homogenization) procedure of the initial equation coefficients such that inversion of A0 is much simpler than inversion of A. We prove contraction of the iteration method and establish explicit estimates of the contraction factor q < 1. For typical quasiperiodic structures we deduce fully computable twosided 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 1/ . 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 twosided a posteriori error estimates.
Acknowledgements. 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 numerical experiments.

References
[1] Bakhvalov, N. S., Panasenko, G. Homogenisation: Averaging Processes In Periodic Media: Mathematical Problems In The Mechanics Of Composite Materials. Springer, 1989.
[2] P. Benner, V. Khoromskaia and B. N. Khoromskij. Range-separated tensor formats for numerical modeling of many-particle interaction potentials. E-preprint, http://arxiv.org/abs/1606.09218, 2016.
[3] Bensoussan, A., Lions, J.-L., Papanicolaou, G. (1978): Asymptotic analysis for periodic structures. Amsterdam: North-Holland
[4] S. Brenner and R. Scott. The mathematical theory of finite element methods. Springer, 1994. [5] J. P. Davis. Circulant matrices. New York. John Wiley & Sons, 1979. [6] Jikov, V.V., Kozlov, S.M., Oleinik, O.A. (1994): Homogenization of differential operators and integral
functionals. Berlin: Springer [7] Friedman, A. (1976): Partial Differential Equations. R. E. Krieger Pub. Co., Huntington, NY [8] 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.

RANK STRUCTURED METHOD

23

[9] Antoine Gloria and Felix Otto. Quantitative estimates on the periodic approximation of the corrector in stochastic homogenization In: ESAIM / Proceedings, 48 (2015), p. 80-97. MIS-Preprint 12/2015, DOI: 10.1051/proc/201448003.
[10] R. Glowinski, J.-L. Lions, R. Tremolieres. Analyse numerique des inequations variationnelles. Dunod, Paris, 1976.
[11] Kantorovich L. V. and Krylov V. L., Approximate Methods of Higher Analysis. Interscience, New York, 1958.
[12] V. Kazeev, I. Oseledets, M. Rakhuba, and Ch. Schwab. QTT-finite-element approximation for multiscale problems I: model problems in one dimension. Adv. Comput. Math., 2016. DOI: 10.1007/s10444-0169491-y.
[13] V. Kazeev, O. Reichmann, and Ch. Schwab. Low-rank tensor structure of linear diffusion operators in the TT and QTT formats. Linear Algebra and its Applications, v. 438(11), 2013, 4204-4221.
[14] 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.
[15] V. Khoromskaia and B. N. Khoromskij. Grid-based lattice summation of electrostatic potentials by assembled rank-structured tensor approximation. Comp. Phys. Commun., 185 (12), 2014, pp. 3162-3174.
[16] V. Khoromskaia, and B.N. Khoromskij. Tensor Approach to Linearized Hartree-Fock Equation for Lattice-type and Periodic Systems. E-preprint arXiv:1408.3839, 2014.
[17] V. Khoromskaia and B.N. Khoromskij. Fast tensor method for summation of long-range potentials on 3D lattices with defects. Numerical Linear Algebra with Applications, 2016, v. 23: 249-271.
[18] V. Khoromskaia and B.N. Khoromskij. Tensor numerical methods in quantum chemistry: from HartreeFock to excitation energies. Phys. Chem. Chem. Phys., 17:31491 - 31509, 2015.
[19] B.N. Khoromskij. Tensor-Structured Preconditioners and Approximate Inverse of Elliptic Operators in Rd. J. Constr. Approx. 30 (2009) 599-620.
[20] B.N. Khoromskij. O(d log N )-Quantics Approximation of N -d Tensors in High-Dimensional Numerical Modeling. Constr. Approx. 34 (2011) 257280.
[21] B.N. Khoromskij. Tensors-structured Numerical Methods in Scientific Computing: Survey on Recent Advances. Chemometr. Intell. Lab. Syst. 110 (2012), 1-19.
[22] B.N. Khoromskij and G. Wittum. Numerical Solution of Elliptic Differential Equations by Reduction to the Interface. Research monograph, LNCSE, No. 36, Springer-Verlag, 2004.
[23] B.N. Khoromskij and S. Repin. A fast iteration method for solving elliptic problems with quasiperiodic coefficients. Russ. J. Numer. Anal. Math. Modelling 2015; 30 (6):329-344. E-preprint arXiv:1510.00284, 2015.
[24] B.N. Khoromskij, S. Sauter, and A. Veit. Fast Quadrature Techniques for Retarded Potentials Based on TT/QTT Tensor Approximation. Comp. Meth. in Applied Math., v.11 (2011), No. 3, 342 - 362.
[25] J.-L. Lions and G. Stampacchia. Variational inequalities. Comm. Pure Appl. Math. 20 1967 493519. [26] Ivan V Oseledets, and S.V. Dolgov. Solution of linear systems and matrix inversion in the TT-format.
SIAM Journal on Scientific Computing, v. 34(5), 2012, A2718-A2739. [27] O. Mali, P. Neittaanmaki, S. Repin. Accuracy verification methods. Theory and algorithms. Springer,
2014 [28] G. I. Marchuk and V. V. Shaidurov. Difference methods and their extrapolations. Applications of Math-
ematics, New York: Springer, 1983. [29] P. Neittaanmaki and S. Repin. Reliable methods for computer simulation. Error control and a posteriori
estimates. Elsevier, 2004. [30] A. Ostrowski. Les estimations des erreurs a posteriori dans les procedes iteratifs, C. R. Acad. Sci, Paris,
Ser. AB 275 (1972), pp. A275A278. [31] S. Repin. A posteriori error estimation for variational problems with uniformly convex functionals,
Math. Comput., 69(2000), 230, 481500. [32] S. Repin. A Posteriori Estimates for Partial Differential Equations. Walter de Gruyter, Berlin, 2008. [33] S. Repin, T. Samrowski, and S. Sauter. A posteriori error majorants of the modeling errors for elliptic
homogenization problems. C. R. Math. Acad. Sci. Paris 351 (2013), no. 23-24, 877-882 [34] 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.

24

B. KHOROMSKIJ AND S. REPIN

[35] 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, 14351451.
[36] U. Schollwock. The density-matrix renormalization group in the age of matrix product states, Ann.Phys. 326 (1) (2011) 96-192.
[37] E. Zeidler. Nonlinear functional analysis and its applications. I. Fixed-point theorems, Springer-Verlag, New York, 1986.

Max Planck Institute for Mathematics in the Sciences, Inselstr. 22-26, 04103, Leipzig, Germany; E-mail: bokh@mis.mpg.de
V.A. Steklov Institute of Mathematics, Fontanka 27, 191 011 St. Petersburg, Russia, and University of Jyvaskyla, Finland ; E-mail: repin@pdmi.ras.ru; serepin@jyu.fi