Towards a high-resolution numerical strategy based on separated representations

: Many models in Science and Engineering are deﬁned in spaces (the so-called conformation spaces) of high dimensionality. In kinetic theory, for instance, the micro scale of a ﬂuid evolves in a space whose number of dimensions is much higher than the usual physical space (two or three). Models deﬁned in such a framework suffer from the curse of dimensionality, since the complexity of the problem growths exponentially with the number of dimensions. This curse of dimensionality makes this class of problems nearly intractable if we perform a standard discretization, say, with ﬁnite element methods, for instance. Problems deﬁned in two or three-dimensional spaces, but densely discretized along each spatial dimension are also hardly tractable by ﬁnite element methods. In this paper we present some recent results concerning a method based on the method of separation of variables, originally developed in [1]. We focus on an efﬁcient imposition of essential non-homogeneous boundary conditions and the treatment of problems with a very high number of degrees of freedom.


INTRODUCTION
In problems that involve high dimensional spaces, standard discretization techniques, like the finite element method, fail due to an excessive computation time required to perform accurate numerical simulations.
The main aim of this paper is to present a new numerical strategy able to circumvent some of the difficulties due to the multidimensional character of kinetic theory models, among other fields of science.The results will be presented to illustrate its capability for treating highly dimensional spaces by considering different multidimensional Poisson's problems, and an efficient imposition of essential non-homogeneous boundary conditions.

THE METHOD
For simplicity, we start by considering the Poisson problem defined in a space of dimension N with homogeneous boundary conditions, The problem solution can be written in the form where F kj is the jth basis function, with unit norm, and only depends on the kth coordinate.Note that the solution of numerous problems can be accurately approximated using a finite number Q of approximation functions.
The numerical scheme proposed consist of an iteration procedure that solves at each iteration n the following three steps.
Step 1 A projection step of the solution in a discrete basis, where the coefficients α j are computed.
Step 3 Enrichment of the approximation basis step.
From the coefficients α j just computed the approximation basis can be enriched by adding the new function N k=1 F k(n+1) (x k ) For details of this scheme see [1].

IMPOSITION OF NON-HOMOGENEOUS ES-SENTIAL BOUNDARY CONDITIONS
Consider, for simplicity and without loss of generality, the Poisson problem described before, see Eq. ( 1) subjected to the boundary conditions: Let us assume that we are able to find a function ψ, continuous in Ω, such that −∆ψ ∈ L 2 (Ω) verifying Eq. (3).Then, the solution of the problem given by Eqs. ( 1)-( 3) can be obtained straightforwardly by where we thus face a problem in the z variable solvable by the method of separation of variables before presented.Following this procedure, the weak form of the problem will be: holds.
Several choices for the function ψ exist.In this work we employ R-functions [4], which constitute an efficient means of constructing ψ functions.
The obtained function ψ, which is generally not product or sum separable, is separated approximately by employing the Singular Value Decomposition (SVD) technique.To this end, we employ singular values up to a given precision, which is normally set to 10 −5 .

NUMERICAL EXAMPLES
In this section we describe some simple numerical experiments that show the potential of the method.

Example: 2D problem
Let us consider now the problem subjected to the boundary conditions that in this case do not verify the partial differential equation ( 8).
In this case, using the before mentioned theory of R-functions, we obtain the following function verifying essential boundary conditions: which is represented in Fig. 1.This function is again not sum or product separable, according to [3].Applying SVD, we arrive to the one-dimensional functions depicted in Fig. 2, that approximate the function ψ to the chosen level of precision, 10 −5 in this case.The method before described was applied in meshes composed by 15, 20, 30, 50, 80 and 130 nodes along each spatial direction.The solution was compared to one approximate solution obtained by standard Finite Elements in a mesh composed by 40000 nodes, thus taken as exact.The order of convergence is remarkably high, and has been found to be of the order of 4 (i.e., superconvergent) for all the cases tested.

Extension of the proposed technique to higherdimensional problems
The problem with the straightforward extension of the proposed technique to higher dimensional problems lies in obtaining the appropriate counterpart of the SVD decomposition in high dimensions.
To this end, we have employed the PARAFAC (Parallel Factors) decomposition, also known as CANDECOMP or Canonical Decomposition, provided by the MATLAB tensor toolbox [2].This technique provides a diagonal tensor Σ, suitable for the purpose of this work.The r most relevant entries of this tensor are employed in the decomposition, in order to accomplish with a user-predefined tolerance in the representation of the function ψ.In general, we seek for decompositions of our function ψ, whose nodal values are stored in an order n (n = 3 in the example here considered, for simplicity of graphical representation) tensor A, in the form: where Σ is a diagonal tensor and U , V , W , ... are the components or factors of the decomposition.
Columns of the factors are orthonormal in standard, two-dimensional SVD, but this is not possible, in general, in tensors of order three or higher.The products involved in the previous decomposition are defined as follows: is such that To analyze the behavior of the technique, we have extended the previous problem to three-dimensional settings: subjected to the boundary conditions The higher-order SVD decomposition of the function ψ up to the tolerance 10 −6 provides the functions depicted in Fig. 4.Only four functions along each spatial direction were necessary.With these functions thus calculated, we let the method run to obtain the solution to the problem as depicted along each direction in Fig. 5.The global, three-dimensional, solution to the problem is depicted in Fig. 6.In this case, 18 functions were necessary along each spatial direction to obtain the before mentioned accuracy.The figure depicts the intersection of the solution with the plane x = 0.

A problem with 10 9 degrees of freedom
In order to show the behavior of the proposed method we have analyzed the problem (the same analyzed in [1]), subjected to the boundary conditions T = 0 on ∂Ω, that has been discretized in a regular lattice of 1000 3 elements, thus giving 1001 3 degrees of freedom.This value is on the limit of current capabilities of nowadays computers, and is analyzed here only for the purpose of showing the performance of the method.
The analytical solution to this problem is given by the function The functions by the method in each spatial direction are depicted in Fig. 7.As can be noticed, the method obtains a very accurate solution in only one iteration, as was the case in [1], employing fewer degrees of freedom.The solution obtained is represented in Fig. 8(top), where a downsampling has been applied to the results by employing a coarser mesh.For obvious reasons, it is difficult to depict the results on a mesh of 10 9 nodes.The error obtained is depicted in  In this example we do not pursue accuracy, but testing the possibility of employing a huge number of degrees of freedom.What is noticeable about this solution, easily reachable with traditional Finite Elements, is that, despite the extremely high number of degrees of freedom involved, the Matlab code ran in around 1 hour on a PC with four processors (only one was employed, no parallel calculation was allowed) with 1.96 Gb of RAM memory each, running under Linux.This time comprised data file writing, generation of figures, etc.

Figure 3 :
Figure 3: Solution convergence obtained by the proposed technique.

Figure 4 :
Figure 4: Functions obtained along x (a), y (b) and z (c) spatial directions for the boundary conditions of the problem stated in Eq. (13).

Figure 5 :
Figure5: Functions obtained along x, y and z spatial directions for the problem stated in Eq. (13).

Figure 6 :
Figure 6: Solution obtained by the proposed technique for the problem with non-homogeneous boundary conditions in three dimensions.Mesh of 40 elements along each direction.

Figure 7 :
Figure 7: Functions obtained along x, y and z spatial directions for the problem stated in Eq. (14).

Fig. 8 (
Fig. 8(bottom), where a slice on an arbitrary section of the domain has been performed for clarity.

Figure 8 :
Figure 8: Solution of the problem stated in Eq. (7) and a slice on the error map.