Mathematics

Academician S. L. SOBOLEV

ON ONE METHOD FOR COMPUTING COEFFICIENTS FOR MECHANICAL CUBATURE FORMULAS

The coefficients of mechanical cubature formulas, exact for all polynomials of a given degree at given nodes, can often be successfully constructed by using the Fourier transform.

Suppose that we need to find \(C_k\) such that the functional

\[ (l, f)=\int_{\Omega} f\,dx-\sum_{k=1}^{N} C_k f(x^{(k)})=\int_{-\infty}^{+\infty} f(x) f(x)\,dx, \tag{1} \]

where \(l(x)\) is the generalized function

\[ l(x)=\mathcal E_{\Omega}(x)-\sum_{k=1}^{N} C_k \delta(x-x^{(k)}), \tag{2} \]

vanishes on all polynomials of degree \(m-1\). Here, obviously, \(\mathcal E_{\Omega}(x)\) is the characteristic function of the domain \(\Omega\), and \(\delta(x-x^{(k)})\) is the generalized Dirac function.

Theorem 1. In order that the functional \((l,\varphi)\) be equal to zero on all polynomials of degree \((m-1)\), it is necessary and sufficient that its Fourier transform \(\tilde l(p)\) have a zero of multiplicity \(m\) at the origin.

From the condition of the theorem it follows that for all \(\alpha\) such that \(|\alpha|\le m-1\), we have

\[ l*x^\alpha=0, \tag{3} \]

where \(\alpha\) is a vector, and \(x^\alpha\) denotes: \(x^\alpha=x_1^{\alpha_1}\cdots x_n^{\alpha_n}\).

Since the Fourier transform of \(x_\alpha\) is \((2\pi)^{n/2}D^\alpha\delta(p)\), and convolution is transformed into multiplication, we obtain

\[ D^\alpha \tilde l(p)=0. \tag{4} \]

Formula (4) means the existence of a zero of multiplicity \(m\) in the function \(\tilde l(p)\). The theorem is proved.

It is convenient to use this theorem in the case when the Fourier transform of the characteristic function \(\mathcal E_{\Omega}(x)\) is computed in finite form. Such a property is possessed, for example, by ellipsoids or polyhedra. In this case the use of the theorem leads to a system of linear equations for determining \(C_k\). This system will, generally speaking, be underdetermined, with the number of equations less than the number of unknowns. By finding its solutions, we construct the required coefficients.

We formulated Theorem 1 for finite domains. If the domain \(\Omega\) is infinite, then the integrals (2) may, generally speaking, not exist.

However, the Fourier transform of the characteristic function \(\mathcal E_{\Omega}(x)\) sometimes remains meaningful. We shall call the generalized function \(l(x)\) a function of order \(m\) if it has a generalized Fourier transform that is \(m\) times continuously differentiable in a neighborhood of the origin, and if the origin is a zero of multiplicity \(m\) for \(\tilde l(p)\).

Any linear combination of functions of order \(m\) will in turn be a function of order \(m\). The use of generalized functions of order \(m\) makes it possible, in some cases, conveniently to compute the coefficients for cubature formulas of order \(m\). We indicate an example.

Let us introduce notation for some functions of one variable, defined on an infinite interval. Let

\[ \psi_0(x)=1,\qquad \Phi_0(x)=\sum_{k=-\infty}^{+\infty}\delta(x-k), \]

\[ \psi_1(x)= \begin{cases} 1, & x>0,\\ 0, & x\leqslant 0; \end{cases} \qquad \Phi_1(x)=\frac12\delta(x)+\sum_{k=1}^{\infty}\delta(x-k), \]

\[ \chi_0(x)=\sum_{k=-\infty}^{+\infty}\delta\left(x-k-\frac12\right),\qquad \chi_1(x)=\sum_{k=0}^{\infty}\delta\left(x-k-\frac12\right). \]

The characteristic functions for coordinate \(s\)-hedral angles \(\Omega^{(i_1,i_2,\ldots,i_s)}\) in \(n\)-dimensional space are represented in the form

\[ \mathcal E^{(i_1,i_2,\ldots,i_s)} = \psi_1(x_{i_1})\psi_1(x_{i_2})\cdots \psi_1(x_{i_s}). \tag{5} \]

A linear transformation makes it possible to consider \(s\)-hedral angles between hyperplanes of arbitrary direction. One can express the characteristic function of some parallelepiped \(0<x_j<a_j\) in the form of a sum of functions of the form (5).

For two variables we obtain, for example:

\[ \mathcal E_{\Omega}(x,y) = \mathcal E^{(1,2)}(x,y) + \mathcal E^{(1,2)}(a-x,y) + \mathcal E^{(1,2)}(x,b-y) + \]

\[ + \mathcal E^{(1,2)}(a-x,b-y) - \mathcal E^{(1)}(x,y) - \mathcal E^{(1)}(a-x,y) - \mathcal E^{(2)}(x,y) - \]

\[ - \mathcal E^{(2)}(x,b-y) + 1. \]

For arbitrary \(n\) the formulas are analogous.

If for all \(\mathcal E^{(i_1,i_2,\ldots,i_s)}\) we have found \(C_k\) such that the differences

\[ \mathcal E^{(i_1,i_2,\ldots,i_s)} - \sum_{k\in K} C_k\delta\bigl(x-x^{(k)}\bigr) \]

are functions of order \(m\), then the corresponding linear combination of them gives the error functional for a certain cubature formula in a finite parallelepiped. This formula will be exact for all polynomials of degree \(m-1\). Instead of a parallelepiped one may take any convex polyhedron with rational faces. We shall seek a cubature formula of order \(m\) for infinite domains \(\Omega^{(i_1,i_2,\ldots,i_s)}\) in the form

\[ \mathcal E^{(i_1,i_2,\ldots,i_s)} - \sum_{k\in K_1} C\delta\bigl(x-x^{(k)}\bigr) - \sum_{k\in K_2} C_k\delta\bigl(x-x^{(k)}\bigr), \]

where \(x^{(k)}\) are elements of the set of nodes of some parallelepipedal system of points; moreover, \(x^{(k)}\) for \(k\in K_1\) are all grid points lying inside the domain \(\Omega^{(i_1,\ldots,i_s)}\), while for \(k\in K_2\) they are points of a certain “boundary layer,” situated at a finite distance from the boundary in this domain. We shall establish in subsequent publications that such a choice can be made so that it is, in a certain sense, close to the most advantageous one.

We shall determine the coefficients in the boundary layer, distinguishing them by orders. Let us call the set \(K_2^{(r)}\) the boundary layer of the \(r\)-th order if it consists of points lying at a finite distance from \(r\) coordinate planes.

We have:

\[ K_2=K_2^{(1)}\supset K_2^{(2)}\supset \cdots \supset K_2^{(n)}. \]

We shall determine the coefficients \(C_k\) in the layer \(K_2^{(r)}\) so that they are common to all functions

\[ \Omega^{(i_1,i_2,\ldots,i_r)},\qquad \Omega^{(i_1,i_2,\ldots,i_{r+1})},\ldots,\qquad \Omega^{(1,2,\ldots,n)}. \]

The Fourier transform makes it possible to compute all coefficients of boundary layers of order \(r\). The boundary layers obtained will be correct in a certain sense, which we shall indicate below.

Theorem 2. For any finite convex polyhedron with rational faces, in each of whose \((s+1)\)-dimensional edges exactly \(s\) \((n-1)\)-dimensional faces intersect, the error functional

\[ l(x)=\mathscr{E}_{\Omega}(x)-C\sum_{k\in K_1}\delta\bigl(x-x^{(k)}\bigr) -\sum_{k\in K_2} C_k\delta\bigl(x-x^{(k)}\bigr) \]

will vanish on all polynomials of degree \(m-1\), if the coefficients in the boundary layer \(K_2\) are chosen the same as in the infinite boundary layers for \(s\)-faced angles formed by the faces of this polyhedron.

We have carried out computations of the simplest boundary layers in the three-dimensional case for the optimal periodic lattice.

In the article \((^1)\) it was established that, for a lattice with fundamental matrix \(H\), with \(|H|=1\), the error in computing the integral of periodic functions from \(L_2^{(m)}\) has the estimate

\[ |(l,\varphi)|^2 \leqslant \sum_{\gamma^*\ne 0}\frac{1}{(\gamma H^{-1})^{2m}}\, \|\varphi\|_{L_2^{(m)}}^2 . \]

The constant

\[ \sum_{\gamma^*\ne 0}\frac{1}{(\gamma H^{-1})^{2m}}=B_{n,m}(H) \tag{6} \]

gives a characteristic of the quality of the lattice.

We shall prove in one of the following publications that the same constant estimates the quality of the lattice also in the nonperiodic case. From the form of \(B_{n,m}(H)\) it follows that, for large \(m\), of the entire sum (6) only its first term has significance,

\[ \frac{1}{r_{\min}^{2m}}, \]

where \(r_{\min}\) is the shortest distance between points of the lattice \(H^{-1}\), dual to the given one, i.e. the maximum diameter of nonintersecting balls with centers at the nodes of this lattice.

It follows from this that the most advantageous lattice will be one for which the diameters, and hence also the volumes, of such balls are the largest.

Since the number of balls in a large domain for various \(H\) with \(|H|=1\) is constant and numerically equal to \(V\), the volume of the domain, we see that, in the optimal case, the lattice \(H^{-1}\) must be the lattice of the densest packing of balls in \(n\)-dimensional space.

It follows from this that for \(n=3\) the optimal \(H^{-1}\) will be the face-centered cubic lattice, and the optimal \(H\) the body-centered cubic lattice. We have constructed formulas for the coordinate planes of two-faced and three-faced nodes in such a lattice.

The computation of the coefficients of such formulas by the method indicated above is based on the Fourier transforms of the functions \(\psi_0,\psi_1,\Phi_0,\Phi_1,\chi_0\) and \(\chi_1\) given above.

The formulas for such a transform are:

\[ \widetilde{\psi}_0(p)=\sqrt{2\pi}\,\delta(p),\qquad \widetilde{\psi}_1(p)=\sqrt{\frac{\pi}{2}}\,\delta(p)+\frac{i}{\sqrt{2\pi}\,p}, \]

\[ \widetilde{\Phi}_0(p)=\sqrt{2\pi}\sum_{k=-\infty}^{+\infty}\delta(p-2k\pi), \]

\[ \widetilde{\Phi}_1(p)=\sqrt{\frac{\pi}{2}}\sum_{k=-\infty}^{+\infty}\delta(p-2k\pi) +\frac{i}{2\sqrt{2\pi}}\operatorname{ctg}\frac{p}{2}, \]

\[ \widetilde{\chi}_0(p)=\sqrt{2\pi}\sum_{k=-\infty}^{+\infty}(-1)^k\delta(p-2k), \]

\[ \widetilde{\chi}_1(p)=\sqrt{\frac{\pi}{2}}\sum_{k=-\infty}^{+\infty}(-1)^k\delta(p-2k\pi) +\frac{i}{2\sqrt{2\pi}}\frac{1}{\sin p/2}. \]

L. V. Voitishek carried out a numerical calculation of the boundary-layer coefficients in such formulas.

For the trihedral angle \(x_1>0,\ x_2>0,\ x_3>0\), the formula for the error functional can be written as:

\[ \begin{aligned} l(x)={}&\varepsilon_0^{(1,2,3)} -\frac12 \Phi_1(x_1)\Phi_1(x_2)\Phi_1(x_3) -\frac12 \chi_1(x_1)\chi_1(x_2)\chi_1(x_3)\\ &-\Phi_1(x_1)\Phi_1(x_2)\Phi_1(x_3) \sum_{j_3=1}^{3}\sum_{k=0}^{2} \alpha_{2k}\frac{\delta(x_{j_3}-k)}{\Phi(x_{j_3})} -\chi_1(x_1)\chi_1(x_2)\chi_1(x_3)\times\\ &\quad \times \sum_{j_3=1}^{3}\sum_{k=0}^{1} \alpha_{2k+1}\frac{\delta(x_{j_3}-k-1/2)}{\chi_1(x_{j_3})} \\ &-\Phi_1(x_1)\Phi_1(x_2)\Phi_1(x_3)\times \sum_{\substack{j_1=1\\ j_1>j_2}}^{3} \sum_{j_2=1}^{3}\sum_{k=0}^{2} \alpha_{2k,\,2l} \frac{\delta(x_{j_1}-k)\delta(x_{j_2}-l)} {\Phi_1(x_{j_1})\Phi_1(x_{j_2})} \\ &-\chi_1(x_1)\chi_1(x_2)\chi_1(x_3) \sum_{\substack{j_1=1\\ j_1>j_2}}^{3} \sum_{j_2=1}^{3}\sum_{k=0}^{2} \alpha_{2k+1,\,2l+1} \frac{\delta(x_{j_1}-k-1/2)\delta(x_{j_2}-l-1/2)} {\chi_1(x_{j_1})\chi(x_{j_2})} \\ &-\sum_{k=0}^{2}\sum_{l=0}^{2}\sum_{j=0}^{2} \alpha_{2k,\,2l,\,2j}\delta(x_1-k)\delta(x_2-l)\delta(x_3-j) \\ &-\sum_{k=0}^{1}\sum_{l=0}^{1}\sum_{j=0}^{1} \alpha_{2k+1,\,2l+1,\,2j+1} \delta(x_1-k-1/2)\delta(x_2-l-1/2)\delta(x_3-j-1/2). \end{aligned} \]

Here, as it turns out, one may take the following values of the coefficient \(\alpha\):

\[ \begin{aligned} \alpha_0&=-0.850694444\cdot 10^{-1},\\ \alpha_2&=-0.116666667,\\ \alpha_4&=-0.93750000\cdot 10^{-2},\\ \alpha_1&=0.160416667,\\ \alpha_3&=0.506944444\cdot 10^{-1},\\ \alpha_{00}&=0.318612557\cdot 10^{-1},\\ \alpha_{02}&=\alpha_{20}=0.633680555\cdot 10^{-1},\\ \alpha_{22}&=0.361111111\cdot 10^{-1},\\ \alpha_{04}&=\alpha_{40}=0.372902200\cdot 10^{-1},\\ \alpha_{24}&=\alpha_{42}=0.607638888\cdot 10^{-1},\\ \alpha_{44}&=0.320818866\cdot 10^{-2},\\ \alpha_1&=-0.127097801,\\ \alpha_{13}&=\alpha_{31}=-0.444299768\cdot 10^{-1},\\ \alpha_{33}&=0.484664352\cdot 10^{-2}, \end{aligned} \qquad \begin{aligned} \alpha_{13}&=\alpha_{31}=-0.444299768\cdot 10^{-1},\\ \alpha_{33}&=0.484664352\cdot 10^{-2},\\ \alpha_{000}&=-0.335894549\cdot 10^{-2},\\ \alpha_{002}&=\alpha_{020}=\alpha_{200}=0.309787330\cdot 10^{-2},\\ \alpha_{022}&=\alpha_{202}=\alpha_{220}=0.527777778\cdot 10^{-1},\\ \alpha_{222}&=\alpha_{224}=\alpha_{242}=\alpha_{442}=\alpha_{333}=0,\\ \alpha_{004}&=\alpha_{040}=\alpha_{400}=0.969939054\cdot 10^{-2},\\ \alpha_{024}&=\alpha_{042}=\alpha_{204}=\alpha_{402}=\alpha_{240}=\alpha_{420}\\ &=0.63964836\cdot 10^{-2},\\ \alpha_{440}&=\alpha_{404}=\alpha_{044}=0.210458262\cdot 10^{-3},\\ \alpha_{442}&=\alpha_{424}=\alpha_{244}=0.969509547\cdot 10^{-2},\\ \alpha_{444}&=0.927847402\cdot 10^{-2},\\ \alpha_{111}&=0.114626736,\\ \alpha_{113}&=\alpha_{131}=\alpha_{311}=0.247395833\cdot 10^{-1},\\ \alpha_{133}&=\alpha_{313}=\alpha_{331}=0.742187498\cdot 10^{-2},\\ \alpha_{333}&=0. \end{aligned} \]

By the same method one computes the boundary layer also near the boundary separating regions with different grid densities, if such regions exist.

Institute of Mathematics with Computing Center
Siberian Branch of the Academy of Sciences of the USSR

Received
21 III 1963

References

1. S. L. Sobolev, DAN, 137, No. 3, 527 (1961).