Fast computation of elastic and hydrodynamic potentials using approximate approximations

We propose fast cubature formulas for the elastic and hydrodynamic potentials based on the approximate approximation of the densities with Gaussian and related functions. For densities with separated representation, we derive a tensor product representation of the integral operator which admits efficient cubature procedures. We obtain high order approximations up to a small saturation error, which is negligible in computations. Results of numerical experiments which show approximation order O(h2M), M = 1,2,3,4, are provided.


Introduction
In this paper we consider volume potentials which arise in the solution of the threedimensional problems in elasticity and hydrodynamics.
Linear isotropic elastic problems in the space R 3 are governed by the threedimensional Lamé system µ ∆ u + (λ + µ)grad div u + f = 0, (1.1) where u = (u 1 , u 2 , u 3 ) is the displacement vector, f = ( f 1 , f 2 , f 3 ) is the volume force, λ and µ are the Lamé constants. The hydrodynamic potentials correspond to the linearized Navier-Stokes equations, the Stokes equations, Here u = (u 1 , u 2 , u 3 ) is the velocity vector, f = ( f 1 , f 2 , f 3 ) is the volume force, P denotes the pressure and ν is the constant viscosity coefficient. The equations (1.2)-(1.3) are considered in the whole space R 3 Suppose that f ∈ S (R 3 ), = 1, 2, 3. By S (R 3 ) we denote the Schwartz space of smooth and rapidly decaying functions. A solution of (1.1) is given by the volume potential Γ k (x − y)g(y)dy, k, = 1, 2, 3.
The concept of approximate approximations and first related results were published by V. Maz'ya in [14], [15]. Various aspects of a general theory of these approximations were systematically investigated (cf. [17] and the references therein). This approximation procedure employs linear combination of smooth and rapidly decaying functions to approximate a given function with high order within a prescribed accuracy. Although in general the method does not converge in a rigorous sense, even simple quasi-interpolants approximate functions very accurately. Additionally, this approximation process has considerable advantages due to the great flexibility in the choice of approximating functions. This flexibility makes it possible to get simple but accurate formulas for the approximation of various integral and pseudo-differential operators. Example of those cubature formulas for volume potentials on R 3 were studied in [16]. By combining semi-analytic cubature formulas for volume potentials based on approximate approximations with the strategy of separated representations, we derive a method for approximating volume potentials, which is accurate and fast also in the multidimensional case and provides approximation formulas of high order. This approach was introduced for the first time in [9] to obtain fast cubature of highdimensional harmonic potentials. Fast cubature of advection-diffusion potentials in the frame of approximate approximations have been obtained in [10]. In [11] the procedure has been extended to parabolic problems of second order and in [12] to the biharmonic operator.
We describe the idea for the elastic potential (1.4). The hydrodynamic potential will be considered in section 4. Our approximation formulas are based on replacing the density f of the integral operator (1.4) by quasi-interpolants of the form with fixed positive parameters h and D and with some generating function η sufficiently smooth and of rapid decay chosen such that the integrals can be effectively determined. Then the sum provides an approximation formula for Γ (k, ) ( f ), k, = 1, 2, 3.
If the generating function η satisfies the moment conditions [17, p.34]) . A proper choice of the parameter D allows to make the saturation error ε 0 (η, D) as small as necessary, e.g., less than the precision machine. The cubature error Therefore it is very efficient to precompute and store the values Γ (k, ) (η 2M )((j − m)/ √ D), which can be used for any gridsize h. However, due to the operation number proportional to h −3 , the convolution (1.12) is not practical. We propose a method which reduces the computational effort. As basis functions, we choose where H k are the Hermite polynomials The functions (1.13) satisfy the moment conditions (1.10) and consequently give rise to approximation formulas of order 2M, plus the saturation error O(h 2 e −Dπ 2 ). Moreover, the action of three-dimensional elastic and hydrodynamic potentials onto these basis functions allows one-dimensional integral representations. For the special case of Gaussian functions, which may be even orthotropic or anisotropic, this was obtained in [17, pp.131-136]. For example, In this paper we obtain new one-dimensional integral representations for the threedimensional integrals Γ (k, ) and Ψ k, applied to the tensor product functions (1.13). The integrand is separated, i.e., it is a product of functions depending only on one of the variables. We show how a suitable quadrature of the one-dimensional integrals combined with the approximation of a separated representation of the density leads to a separated approximation of the integral operator. Since the separated representation of the density can be approximated with high order and only one-dimensional operations are used, the resulting method is very efficient and can provide high order approximations.
The outline of the paper is as follows. In section 2 we derive one-dimensional integral representations of the three-dimensional elastic potential applied to the functions (1.13). In section 3, for densities with separated representation, we derive a tensor product representation of the integral operator which admits efficient onedimensional operations. We provide numerical tests, showing that these formulas are efficient and provide approximation of order O(h 2M ), M = 1, 2, 3, 4. In section 4 we apply the procedure to the hydrodynamic potentials (1.6) and (1.8), and provide approximation up to order O(h 8 ).

Elastic potential
where L denotes the harmonic potential and A one-dimensional integral representation with separable integrand for the harmonic potential is known ([9, p.893]), It remains to determine a one-dimensional integral representation for I k (e −|·| 2 ).
Thus, from the formulas (2.1), (2.4) and (2.5) we obtain Let M > 1. We are looking for a one-dimensional integral representation with separated integrand of the potential Γ (k, ) with the basis functions (1.13). Keeping in mind (2.1), we write A one-dimensional integral representation for the harmonic potential applied to η 2M in (1.13) is known ([9, p.894]) where Q M (x,t) is a polynomial in x whose coefficients depend on t, defined by and, for k = , Proof Using the relation ([17, p.55]) and integrating by parts we get Therefore, by using (2.5), we derive that Then it is easy to verify that Q Hence, For k = , with k, = 1, 2, 3 we have In view of (2.13), (2.14), (2.17) and (2.18), the integrals I (η 2M ) and I k (η 2M ), k = , can be written in the form (2.9) and (2.10), respectively. Direct computation of the left hand sides in (2.15) and (2.16) leads to (2.11) and (2.12), respectively.
3 Implementation and numerical results for the elastic potential In this section we consider fast computation of the elastic potential (1.4) based on (2.19) and (2.20). In view of (1.12), the quadrature of the integrals (2.19) and (2.20) with certain quadrature weights ω p and nodes τ p leads to the approximation formulas M,h (g) are very efficient if the function g has separated representation, i.e., for a given accuracy ε it can be represented as the sum of products of vectors in dimension 1 Then an approximate value of the three-dimensional convolutional sum Γ (k, ) (g)(hj) can be approximated using only one-dimensional operations as follows We use an efficient quadrature rule based on the classical trapezoidal rule, which is exponentially converging for rapidly decaying smooth functions on the real line. We make the substitutions with positive parameters α and β proposed in [20]. Then, the integrals (2.19), (2.20) are transformed to integrals over R with integrands decaying doubly exponentially in the variable u. Thus the classical trapezoidal rule provides very accurate approximations of the integral for a relatively small number of nodes τ p = τ p. After the substitution we have We provide results of some experiments which show the accuracy and convergence rate of the method. We compute the elastic potential of the density f(x) = (e −|x| 2 , 0, 0) which has the exact value u = (Γ (1,1) e −|·| 2 ,Γ (2,1) e −|·| 2 ,Γ (3,1) e −|·| 2 ) with (cf. [17, (5.56), p.110]) In Table 1 we compare the exact values Γ (1,1) (e −|·| 2 ) and the approximate values    Table 3 Absolute error and rate of convergence for Γ (2,1) (e −|·| 2 )(0.8, 0.8, 0.8) using Γ M,h behaves in numerical computations as a usual 2M−order formula. In Table 4 we report on the absolute error and the convergence rate of the elastic potential Γ (1,1) (e −|·| 2 )(0, 1, 0) using formulas Γ  we compute the elastic potential of the density f = ( f 1 , f 2 , f 3 ) with components which provides the exact solution u = (e −|x| 2 /2, 0, 0). In Table 5 we report on the absolute error and the approximation rate for where, rewriting (1.7) as and combining it with (2.2) and (2.3), we have In view of (2.7), (2.9) and (2.10), we can write Ψ (k, ) (η 2M ) in the form with i ∈ {1, 2, 3} \ {k, }. Analogously to sections 2 and 3, formulas (4.1),(4.2) can be used to construct high order cubature formulas for the u = (u 1 , u 2 , u 3 ) in (1.6). The formulas are very efficient if f , = 1, 2, 3, have separated representations. Concerning the approximation of in (1.8), the gradient of L f is approximated by the gradient of L (M h f ). Then, (4.4) The combination of (4.3) with (4.4) gives Theorem 3 [13, Theorem 3.2] The function ∂ ∂ x L (η 2M )(x) admits the following one-dimensional integral representation M given in (2.8), R 1 (x,t) = 0 and, for M > 1, Hence, Based on the mapping properties of the integral operators Ψ (k, ) , estimates of the cubature error can be obtained similar to those of the elastic potential. Indeed, we have ( [17, p.113 which has the exact value u(x) = (curl(0, 0, −e −|x| 2 /2)) = (x 2 e −|x| 2 , −x 1 e −|x| 2 , 0), P(x) = e −|x| 2 /2.
In Table 7 we report on the absolute errors and approximation rates for the value   We conclude this section by showing accuracy and convergence rate of formula (4.5) for the approximation of P in (4.3) (see Table 8). We assume f 1 = e −|x| 2 (3 − 2|x| 2 ), f 2 = f 3 = 0 which provides, keeping in mind that L ( f 1 ) = e −|x| 2 /2, the exact solution P(x) = x 1 e −|x| 2 .
For the calculations we choose the parameter ν = 2, the quadrature rule with α = 5 and β = 6 in (3.1), τ = 0.003 and 250 summands in the quadrature sum. We choose D = 4 to have the saturation error comparable with the double precision rounding errors. The numerical results confirm the h 2M convergence of the cubature formula when M = 1, 2, 3, 4.  Table 8 Absolute error and rate of convergence for P(x) in x = (0.4, 0.4, 0) using (4.5).