Numerical approximation of Poisson problems in long domains

In this paper, we consider the Poisson equation on a"long"domain which is the Cartesian product of a one-dimensional long interval with a (d-1)-dimensional domain. The right-hand side is assumed to have a rank-1 tensor structure. We will present methods to construct approximations of the solution which have tensor structure and the computational effort is governed by only solving elliptic problems on lower-dimensional domains. A zero-th order tensor approximation is derived by using tools from asymptotic analysis (method 1). The resulting approximation is an elementary tensor and, hence has a fixed error which turns out to be very close to the best possible approximation of zero-th order. This approximation can be used as a starting guess for the derivation of higher-order tensor approximations by an alternating-least-squares (ALS) type method (method 2). Numerical experiments show that the ALS is converging towards the exact solution (although a rigorous and general theoretical framework is missing for our application). Method 3 is based on the derivation of a tensor approximation via exponential sums applied to discretised differential operators and their inverses. It can be proved that this method converges exponentially with respect to the tensor rank. We present numerical experiments which compare the performance and sensitivity of these three methods.


Introduction
In this paper, we consider elliptic partial differential equations on domains which are the Cartesian product of a "long" interval I = (− , ) with a (d − 1)-dimensional domain ω, the cross section -a typical application is the modelling of a flow in long cylinders.As a model problem we consider the Poisson equation with homogeneous Dirichlet boundary conditions and a right-hand side which is an elementary tensor ; i.e., the product of a univariate function (on the long interval) and a (d − 1)variate function on the cross section.Such problems have been studied by using asymptotic analysis, see., e.g., [2].Our first approximation (method 1) is based on this technique and approximates the solution by an elementary tensor where the function on the cross section is the solution of a Poissonproblem on the cross section and the corresponding univariate function is determined afterwards as the best approximation in the Sobolev space H 1 0 on the long interval.In Lemma 2 below, it is shown that this approximation converges exponentially with respect to the length of the cylinder for any subdomain I 0 × ω for fixed 0 < .However, for fixed this is a one-term approximation with a fixed error.
Method 2 uses the result of method 1 as the initial guess for an iterative procedure which is an alternating least squares (ALS) method.Recursively, one assumes that a rank-k tensor approximation of the solution has already been derived and then starts an iteration to compute the k + 1 term: a) one chooses an univariate function on I as an initial guess for this iteration and determines the function on the cross section as the best approximation in H 1 0 of the cross-section.In step b) the iteration is flipped and one fixes the new function on the cross section and determines the corresponding best approximation in H 1 0 of the interval.Steps a) and b) are iterated until a stopping criterion is reached and this gives the k + 1 term in the tensor approximation.We have performed numerical experiments which are reported in Section 4 which show that this method leads to a convergent approximation also for fixed as the tensor rank of the approximation increases.However, it turns out that this method is quite sensitive and requires that the inner iteration a), b) leads to an accurate approximation of the (k + 1) term in order to ensure that the outer iteration is converging.Furthermore, the numerical experiments that we have performed indicate, that the convergence speed can slow down as the number of outer iterations increases.Thus, this method is best suited when a medium approximation accuracy of the Poisson problem is required.An analysis of this method is available in the literature for the case that the higher-dimensional Hilbert space norms are generated from one-dimensional Hilbert space norms In this case, the ALS converges for any initial guess; we refer for this and more details to [10].In our application however, the norms are not generated from one-dimensional norms and convergence is still an open theoretical question.Further results on the alternating least squares approach can be found in [3,8,9].Method 3 is based on a different approach which employs numerical tensor calculus (see [4]).First one defines an exponential sum approximation of the function 1/x.Since the differential operator −∆ is of tensor form, the exponential sum, applied to the inverse of a discretisation of the Laplacian by a matrix which must preserve the tensor format, directly leads to a tensor approximation of the solution u.We emphasize that the explicit computation of the inverse of the discretisation matrix can be avoided by using the hierarchical format for their representation (see [5]).An advantage of this method is that a full theory is available which applies to our application and allows us to choose the tensor rank via an a priori error estimate.It also can be shown that the tensor approximation converges exponentially with respect to the tensor rank (cf.[4]).
The paper is structured as follows.In Section 2 we formulate the problem on the long product domain and introduce the assumptions on the tensor format of the right-hand side.The three different methods for constructing a tensor approximation of the solution are presented in Section 3. The results of numerical experiments are presented in Section 4 where the convergence and sensitivity of the different methods is investigated and compared.For the experiments we consider first the case that the cross section is the one-dimensional unit interval and then the more complicated case that the cross section is an L-shaped polygonal domain.Finally, in the concluding section we summarize the results and give an outlook.

Setting
Let ω be an open, bounded and connected Lipschitz domain in R n−1 , n ≥ 1.In the following we consider Poisson problems on domains of the form where is large.More specifically we are interested in Dirichlet boundary value problems of the form with weak formulation Specifically we are interested in right-hand sides f which have a tensor structure of the form , where g k is a univariate function and f , f k are functions which depend only on the (d − 1)-dimensional variable x ∈ ω.Here, we use the standard tensor notation (g In this paper, we will present and compare methods to approximate u in tensor form.
We consider a right-hand side of the form and derive a first approximation of u as the solution of the (n − 1)-dimensional problem on ω:

Numerical Approximation
In this section we derive three different methods to approximate problem (1).In all three methods we exploit the special structure of the domain Ω and the right-hand side F .Our goal is to reduce the original n-dimensional problem on Ω to one or more (n − 1)-dimensional problems on ω.Compared to standard methods like finite elements methods or finite difference methods, which solve the equations on Ω , this strategy can significantly reduce the computational cost since is considered large and the discretisation in the x 1 direction can be avoided.

Method 1: A one-term approximation based on an asymptotic analysis of problem (1)
Although the right-hand side F in (1) is independent of x 1 , it is easy to see that this is not the case for the solution u , i.e., due to the homogeneous Dirichlet boundary conditions it is clear that u depends on x 1 .However, if is large one can expect that u is approximately constant with respect to x 1 in a subdomain Ω 0 , where 0 < 0 and thus converges locally to a function independent of x 1 for → ∞.The asymptotic behaviour of the solution u when → ∞ has been investigated in [1].It can be shown that where u ∞ is the solution of (3), with an exponential rate of convergence.More precisely, the following theorem holds: Theorem 1 There exist constants c, α > 0 independent of s.t.
For a proof we refer to [1,Theorem 6.6].
Theorem 1 shows that 1 ⊗ u ∞ is a good approximation of u in Ω /2 when is large.This motivates to seek approximations of u in Ω which are of the form where ψ ∈ H 1 0 (− , ).Here, we choose ψ to be the solution of the following best approximation problem: Given u ∈ H 1 0 (Ω ) and u ∞ ∈ H 1 0 (ω), find ψ ∈ H 1 0 (− , ) s.t.
In order to solve problem (4) we define the functional and consider the variational problem of minimizing it with respect to θ ∈ H 1 0 (− , ).A simple computation shows that this is equivalent to finding θ ∈ H 1 0 (I ) such that The strong form of the resulting equation is The solution of this one-dimensional boundary value problem is given by This shows that an approximation of our original problem ( 1) is given by and Note that ψ (a, •) satisfies In Section 4 we report on various numerical experiments that show the approximation properties of this rather simple one-term approximation.Lemma 2 There exist constants c, c > 0 independent of such that, for δ < , (1) ω,δ := 4 and The right-hand side in (8) goes to 0 with an exponential rate of convergence if δ is bounded from below when → ∞.
Proof.For i = 1, 2, . .., let w i be the i-th eigenfunction of −∆ , i.e., and we normalize the eigenfunctions such that (w i , w j ) 2,ω = δ i,j and order them such that (λ i ) i is increasing monotonously.Furthermore let u ,i ∈ H 1 0 (Ω ) be the solution of (∇u ,i , ∇v Then one concludes from ( 7) and ( 10) that This shows that the solutions of ( 3) and ( 1) can be expressed as With ψ as in (5) we get Let δ < .Then, since ω w i w j dx = δ i,j , we get One has for any α > 0 and similarly Since λ 1 ≤ λ i for all i ∈ N and and We employ the estimates (13) and ( 14) in ( 12) and obtain ,ω , which shows the assertion.
Lemma 2 suggests that one cannot expect convergence of the approximation ψ (λ ∞ , •) ⊗ u ∞ on the whole domain Ω .Indeed it can be shown that, in general, Lemma 2 shows that the error on Ω can be estimated as follows: where λ 1 is as in (9).

Method 2: An alternating least squares type iteration
Method 1 can be interpreted as a 2-step algorithm to obtain an approximation u M1 of u .
• Step 1: Solve (3) in order to obtain an approximation of the form 1 ⊗ u ∞ which is nonconforming, i.e., does not belong to H 1 0 (Ω ).
• Step 2: Using u ∞ , find a function ψ that satisfies (4) in order to obtain the conforming In this section we extend this idea and seek approximations of the form by iteratively solving least squares problems similar to (4).We denote by the residual of the approximation and suggest the following iteration to obtain u M2 ,m : q (m) = arg min Then, given q (m) , find p (m) ∈ H 1 0 (I ) s.t.
The algorithm exhibits properties of a greedy algorithm.It is easy to see that in each step of the iteration the error decreases or stays constant.We focus here on its accuracy in comparison with the two other methods via numerical experiments.We emphasize that for tensors of order at least 3, local convergence (under suitable conditions) can be shown for the ALS iteration (see [3,8,9,10]).
One assumption in these papers is that the scalar product of the tensor space is generated by the scalar products of the single spaces -this, however, is not the case in our setting.We note that an analysis of the approximation (15) by the theory of singular value decomposition (SVD) is also not feasible since the function u is unknown.
In each step of the (outer) iteration above we need to solve at least two minimization problems (16) and (17).In the following we derive the strong formulations of these problems.

Resolution of (16)
As before an investigation of the functional shows that q (m) needs to satisfy for all q ∈ H 1 0 (ω), where For the right-hand side we obtain .
In order to obtain the solution of (17) we therefore have to solve in Remark 4 The constants p 1,m−1 , p2,j,m−1 , q 1,m and q2,j,m involve derivatives and Laplace-operators.Note that after solving (18) and ( 19) for q (m) and p (m) , discrete versions of ∆ q (m) and p (m) can be easily obtained via the same equations.Furthermore, since a numerical computation of the gradients can be avoided.

Method 3: Exploiting the tensor product structure of the operator
In this section we exploit the tensor product structure of the Laplace operator and the domain Ω .
Note, that we do not assume that ω has a tensor product structure.Furthermore the Laplace operator in our original problem (1) can be written as We discretise (1) with F as in (2) on a mesh G, e.g., by finite elements or finite differences on a tensor mesh, i.e., each mesh cell has the form (x i−1 , x i ) × τ j , where τ j is an element of the mesh for ω.The essential assumption is that the system matrix for the discrete version of −∆ in (20) is of the tensor form If we discretise with a finite difference scheme on an equidistant grid for I with step size h, then A 1 is the tridiagonal matrix h −2 tridiag [−1, 2, −1] and M 1 is the identity matrix.A finite element discretisation with piecewise linear elements leads as well to It can be shown that the inverse of the matrix A can be efficiently approximated with a sum of matrix exponentials.More precisely the following Theorem holds which is proved in [4], Proposition 9.34.
Theorem 5 Let M (j) , A (j) be positive definite matrices with λ (j) min and λ (j) max being the extreme eigenvalues of the generalized eigenvalue problem A (j) x = λM (j) x and set Then A −1 can be approximated by where the coefficients a ν , α ν > 0 are such that max .The error can be estimated by where M = ⊗ n j=1 M (j) .
Theorem 5 shows how the inverse of matrices of the form (21) can be approximated by sums of matrix exponentials.It is based on the approximability of the function 1/x by sums of exponentials in the interval [a, b].We refer to [4,6] for details how to choose r and the coefficients a ν,[a,b] , α ν, [a,b] in order to reach a given error tolerance ε( 1x , [a, b], r).Note that the interval [a, b] where 1/x needs to be approximated depends on the matrices A (j) and M (j) .Thus, if A changes a and b need to be recomputed which in turn has an influence on the optimal choice of the parameters a ν,[a,b] and α ν, [a,b] .
Numerical methods based on Theorem 5 can only be efficient if the occurring matrix exponential can be evaluated at low cost.In our setting we will need to compute the matrices The evaluation of the first matrix will typically be simpler.In the case where a finite difference scheme is employed and A 1 is a tridiagonal Toeplitz matrix while M 1 is the identity, the matrix exponential can be computed by diagonalizing A 1 , e.g., A 1 = SD 1 S −1 , and using exp The computation of exponentials for general matrices is more involved.We refer to [7] for an overview of different numerical methods.Here, we will make use of the Dunford-Cauchy integral (see [5]).For a matrix M we can write for contour C = ∂D which encircles all eigenvalues of M .We assume here that M is positive definite.Then the spectrum of M satisfies σ( M ) ⊂ (0, M ] and the following (infinite) parabola ds. ( The integrand decays exponentially for s → ±∞.Therefore ( 23) can be efficiently approximated by sinc quadrature, i.e., exp where h > 0 and should be chosen s.t.h = O (N + 1) −2/3 .We refer to [5] for an introduction to sinc quadrature and for error estimates for the approximation in (24).The parameters h and N in our implementation have been chosen such that quadrature errors become negligible compared to the overall discretisation error.For practical computations, the halving rule (cf .[5, §14.2.2.2]) could be faster while the Dunford-Schwartz representation with sinc quadrature is more suited for an error analysis.

The case of a planar cylinder
In this subsection we apply the methods derived in Section 3 to a simple model problem in two dimensions.We consider the planar cylinder and solve (1) for different right-hand sides ) and different lengths .The reduced problem (3) on ω = (−1, 1) is solved using a standard finite difference scheme.We compare the approximations of (1) to a reference solution u ref 2D, that is computed using a finite difference method on sufficiently refined two-dimensional grid.
In Table 1 we state the L 2 Ω 2D -errors of the approximations u M1 2D, for various values of and right-hand sides f .Having in mind that u M1 2D, is a rather simple one-term approximation that only requires the solution of one (n − 1)-dimensional problem (plus some postprocessing), the accuracy of the approximation is satisfactory especially for larger values of .
Figure 2 shows the pointwise, absolute error |u M1 2D, − u ref 2D, | in Ω for = 10 and f (x ) = tanh(4x + 1).As expected the accuracy of the approximation is very high in the interior of the planar cylinder (away from ± ).Lemma 2 (and Figure 2) suggests that the approximation in the interior of the cylinder is significantly better than on the whole domain Ω .Indeed, if the region of interest is only a subdomain Ω 0 ⊂ Ω , where 0 < , the error decreases exponentially as 0 → 0. Figure 3 shows the relative error with respect to 0 for = 20, 50 and the right-hand side f (x ) = tanh(4x + 1).We can see that the exponential convergence sets in almost immediately as l 0 moves away from .
To conclude, Method 1 can be used in applications where • only a limited approximation accuracy is required, • a good starting point for more accurate methods is needed, • the region of interest is a subdomain Ω 0 of Ω with 0 < .

In Method 2 we use u M1
2D, as starting value of the iteration which is then successively refined by approximating the residual in each step with a series of L 2 best approximations.In Table 2 we   2D, ,m for different values of and iterations m.We used f (x ) = tanh(4x + 1) throughout.
state the relative errors of this approach in the case f (x ) = tanh(4x + 1) for different values of and iteration steps.We can see that five iterations are sufficient to reduce the error of the initial approximation u M1 2D, by a factor 100 for all considered values of .However, in this case more iterations do not lead to significantly better results and the convergence seems to flatten.One explanation for this is that the residuals are increasingly difficult to approximate with each step of the iteration.After a few iterations a one-term approximation of these residuals of the form p (m) ⊗ q (m) therefore is not sufficiently accurate which leads to reduced decay of the error in the overall scheme.Note that in the case = 1, Ω cannot be considered as a "long" domain.Therefore, the initial approximation u M1 2D, only exhibits a low accuracy.Nevertheless the error of u M2 2D, ,m decays quickly as m increases and reaches a similar level of accuracy as for larger .This suggests that Method 2 can also be used for more general domains Ω .
In Table 3 we show the relative errors of the approximations u M3 2D, ,r for f (x ) = tanh(4x + 1) and different values of and r.As the theory predicts the error decays exponentially in r and is governed by the approximability of the function 1/x by exponential sums.Note that in this two-dimensional example the arising matrix exponentials could be computed via diagonalization of  2D, ,r for different values of and r.We used f (x ) = tanh(4x + 1) throughout.
the involved finite difference matrices.An approximation of the Dunford-Cauchy integral was not necessary in this case.

A three-dimensional domain with a non-rectangular cross section
In this section we consider the three-dimensional domain where ω is not a rectangle (see Figure 4).As before we solve problem (1) for different right-hand sides f and different values of .The reduced problem (3) on ω is solved using a standard 2D finite difference scheme.As 3D reference solution we use an accurate approximation using method 3, i.e. u M3 3D, ,r for r = 30, which is known to converge exponentially in r.Table 4 shows the relative errors of the approximations u M1 3D, for different values of and righthand sides f .As the theory predicts we cannot observe an exponentially decreasing error as gets large, since we measure the error on the whole domain Ω and not only a subdomain Ω −δ .As before we only have to solve one two-dimensional problem on ω in order to obtain the approximation u M1 3D, .3D, ,m for different values of and iterations m.We used f (x ) = tanh(x 1 x 2 ) throughout.
In Table 5 we show the relative errors of the approximations u M2 3D, ,m for f (x ) = tanh(x 1 x 2 ) and different values of and m (number of iterations).As in the 2D case this method significantly improves the initial approximation u M2 3D, ,1 = u M1 3D, using the alternating least squares type iteration.However, also here we observe that the convergence slows down when a certain accuracy is reached.We remark that a good starting point for the iteration is crucial for this method.In all our experiments u M1 3D, was a good choice which leads to a convergence behaviour similar to the ones in Table 5.Other choices often did not lead to satisfactory results.
In Table 6 we show the relative errors of the approximations u M3 3D, ,r again for f (x ) = tanh(x 1 x 2 ) and different values of and r.As before the error decays exponentially with respect to r.The arising matrix exponentials exp −α ν,[a,b] A x in these experiments were computed using the sinc quadrature approximation (24).The number of quadrature points N was chosen such that the corresponding quadrature error had an negligible effect on the overall approximation.3D, ,r for different values of and r.We used f (x ) = tanh(x 1 x 2 ) throughout.

Conclusion
We have presented three different methods for constructing tensor approximations to the solution of a Poisson equation on a long product domain for a right-hand side which is an elementary tensor.
The construction of a one-term tensor approximation is based on asymptotic analysis.The approximation converges exponentially (on a fixed subdomain) as the length of the cylinder goes to infinity.However, the error is fixed for fixed length since the approximation consists of only one term.The cost for computing this approximation is very low -it consists of solving a Poisson-type problem on the cross section and a cheap post-processing step to find the univariate function in the one-term tensor approximation.
The ALS method uses this elementary tensor and generates step-by-step a rank-k approximation.The computation of the m-th term in the tensor approximation itself requires an inner iteration.If one is interested in only a moderate accuracy (but improved accuracy compared to the initial approximation) this method is still relatively cheap and significantly improves the accuracy.However, the theory for ALS for this application is not fully developed and the definition of a good stopping criterion is based on heuristics and experiments.
Finally the approximation which is based on exponential sums is the method of choice if a higher accuracy is required.A well developed a priori error analysis allows us to choose the tensor rank in the approximation in a very economic way.Since the method is converging exponentially with respect to the tensor rank, the method is also very efficient (but more expensive than the first two methods for the very first terms in the tensor representation).However, its implementation requires the realization of inverses of discretisation matrices in a sparse H-matrix format and a contour quadrature approximation of the Cauchy-Dunford integral by sinc quadrature by using a non-trivial parametrisation of the contour.
We expect that these methods can be further developed and an error analysis which takes into account all error sources (contour quadrature, discretisation, iteration error, asymptotics with respect to the length of the cylinder, H-matrix approximation) seems to be feasible.Also the methods are interesting in the context of a-posteriori error analysis to estimate the error due to the truncation of the tensor representation at a cost which is proportional to the solution of problems on the cross sections.We further expect that more general product domains of the form × d m=1 ω m for some ω m ∈ R dm with dimensions 1 ≤ d m ≤ d such that d m=1 d m = d and domains with outlets can be handled by our methods since also in this case zero-th order tensor approximation can be derived by asymptotic analysis (see [2]).

Figure 1
Figure1shows a plot of ψ (λ ∞ , •) for = 20 and λ ∞ = 2. Since ψ approaches 1 with an exponential rate as x 1 moves away from ± towards the origin, an analogous result to Theorem 1 can be shown for u M1 .

Table 1 :
Relative L 2 Ω 2D -errors of the approximations u M 1 2D, for different values of and f .

Table 2 :
Relative L 2 -errors of the approximations u M 2

Table 3 :
Relative L 2 -errors of the approximations u M 3

Table 4 :
Relative L 2 -errors of the approximations u M 1 3D, for different values of and f .

Table 5 :
Relative L 2 -errors of the approximations u M 2

Table 6 :
Relative L 2 -errors of the approximations u M 3