UDC 517.946.6 : 517.947.42

ON THE QUESTION OF A PRACTICAL APPLICATION OF A NEW APPROXIMATION METHOD

M. A. Aleksidze

In the book by V. D. Kupradze [1] a new method is presented for the approximate solution of boundary-value problems, called the method of generalized Fourier series. The principal point in the justification of this method is the proof of the linear independence and completeness, on any closed surface \(S\), of the system of functions \(\{\omega_i(M)\}\), where \(\omega_i(M)\) are so-called potential functions [3], chosen in a definite way, fundamental solutions of the corresponding boundary-value problem. In the case of the Laplace equation and the Dirichlet problem such a system of functions will be \(\{\omega_i(M)\}=\{r^{-1}(x_i,M)\}\), where \(r(x_i,M)\) is the distance between the points \(x_i\) and \(M\), and \(x_i\) are elements of a countable set of points situated everywhere densely on an auxiliary surface \(S_1\). The domain with surface \(S_1\) either entirely contains the domain with surface \(S\) (for interior boundary-value problems), or is entirely contained in the domain with surface \(S\) (for exterior boundary-value problems). In both cases the surfaces \(S\) and \(S_1\) do not touch one another.

In the present note we shall make two remarks concerning the practical application of the method of generalized Fourier series.

Let us first recall some properties of systems of orthonormal functions.

Let \(\{\omega_i(M)\}\) be any linearly independent system of functions. Then the functions

\[ \varphi_n(M)= \frac{ \left| \begin{array}{ccccc} \int \omega_1^2\,dS, & \int \omega_1\omega_2\,dS, & \cdots, & \int \omega_1\omega_{n-1}\,dS, & \omega_1\\ \int \omega_1\omega_2\,dS, & \int \omega_2^2\,dS, & \cdots, & \int \omega_2\omega_{n-1}\,dS, & \omega_2\\ \cdots & \cdots & \cdots & \cdots & \cdots\\ \cdots & \cdots & \cdots & \cdots & \cdots\\ \int \omega_n\omega_1\,dS, & \int \omega_n\omega_2\,dS, & \cdots, & \int \omega_n\omega_{n-1}\,dS, & \omega_n \end{array} \right| }{ \sqrt{G_{n-1}\cdot G_n} }, \tag{1} \]

\[ n=1,2,\ldots \]

form an orthonormal set. Here \(G_n=G_n(\omega_1,\omega_2,\ldots,\omega_n)\) is the Gram determinant,

\[ G_n= \left| \begin{array}{ccccc} \int \omega_1^2\,dS, & \int \omega_1\omega_2\,dS, & \cdots, & \int \omega_1\omega_n\,dS\\ \int \omega_1\omega_2\,dS, & \int \omega_2^2\,dS, & \cdots, & \int \omega_2\omega_n\,dS\\ \cdots & \cdots & \cdots & \cdots\\ \cdots & \cdots & \cdots & \cdots\\ \cdots & \cdots & \cdots & \cdots\\ \int \omega_n\omega_1\,dS, & \int \omega_n\omega_2\,dS, & \cdots, & \int \omega_n^2\,dS \end{array} \right|, \tag{2} \]

for which, as is known, the inequality [4]

\[ G_n(\omega_1,\omega_2,\ldots,\omega_n)\leq G_m(\omega_1,\omega_2,\ldots,\omega_m)G_{n-m}(\omega_{m+1},\omega_{m+2},\ldots,\omega_n). \tag{3} \]

is valid. Following [2], we shall call a system of linearly independent functions reliable if \(\lim_{n\to\infty}G_n=G>0\); for \(G=0\) we shall call the system unreliable.

Let us prove the proposition. For any \(\varepsilon>0\) there exists an \(N\) such that for \(n>N\) the inequality

\[ G_n(\omega_1,\omega_2,\ldots,\omega_n)<\varepsilon \]

will hold.

Obviously, this means that the system of potential functions \(\{\omega_i(M)\}\) is unreliable in the sense of [2]. We shall assume that the points \(x_i\) \((i=1,2,\ldots,n)\) are renumbered so that \(r(x_i,x_{i+1})<h\) and \(h\to0\) as \(n\to\infty\). It is clear that, under the assumptions made, this is always possible. From (3) we obtain, for even \(n=2k\),

\[ G_n(\omega_1,\omega_2,\ldots,\omega_n)\leq \prod_{i=1}^{k}G_2(\omega_{2i-1},\omega_{2i}) \]

and for odd \(n=2k+1\)

\[ G_n(\omega_1,\omega_2,\ldots,\omega_n)\leq \iint_S \omega_n^2\,dS\prod_{i=1}^{k}G_2(\omega_{2i-1},\omega_{2i}). \]

For \(G_2(\omega_{2i-1},\omega_{2i})\) we obtain

\[ G_2(\omega_{2i-1},\omega_{2i})= \iint_S \omega_{2i-1}^2\,dS \iint_S \omega_{2i}^2\,dS -\left(\iint_S \omega_{2i-1}\omega_{2i}\,dS\right)^2. \]

Taking into account that

\[ \omega_{2i-1}=\frac{1}{r(x_{2i-1},M)} =\frac{1}{r(x_{2i},M)+\xi(M)}, \tag{4} \]

where \(|\xi(M)|=|r(x_{2i-1},M)-r(x_{2i},M)|\leq r(x_{2i-1},x_{2i})<h\), for \(G_2(\omega_{2i-1},\omega_{2i})\) we obtain, up to terms of higher order of smallness in \(h\),

\[ G_2(\omega_{2i-1},\omega_{2i})\approx O(h) \tag{5} \]

and for \(G_n\) we shall have \(G_n\approx O(h^k)\).

Since as \(n\to\infty\), \(h\to0\), from the last approximate equality there follows the unreliability of the system \(\{\omega_i(M)\}\).

Thus, the Gram determinant \(G_n\) tends to zero as \(n\to\infty\). To check the rate of convergence, the following numerical experiments were carried out for the plane case. \(S\) is a circle of radius 1, on which it is required to orthonormalize the system of functions \(\{\ln r(x_i,M)\}\) \((i=1,2,\ldots,28)\), where \(x_i\in S_1^{(1)}\) and \(S_1^{(1)}\) is a concentric circle of radius 2. The points \(x_i\) \((i=1,2,\ldots,28)\) are distributed uniformly on \(S_1^{(1)}\), with step \(\pi/14\). The Gram determinant of the system was computed. The integration was performed with an error not exceeding \(10^{-6}\). The rank of the determinant with accuracy up to \(10^{-9}\) is equal to 9, i.e., after going through all possible determinants of the tenth and higher orders, we were unable to find among them one different from machine zero (all computations were carried out on the BESM-2 machine).

Below are the maximum values \(\bar G_i\)—determinants of order \(i\), for \(i=1,2,\ldots,9\):

\[ \bar G_1=3.8;\qquad \bar G_2=9.7;\qquad \bar G_3=9.5;\qquad \bar G_4=5.9;\qquad \bar G_5=0.4; \]

\[ \overline{G}_6 = 2.7 \cdot 10^{-2}; \qquad \overline{G}_7 = 9.2 \cdot 10^{-4}; \qquad \overline{G}_8 = 2.8 \cdot 10^{-5}; \qquad \overline{G}_9 = 3.7 \cdot 10^{-8}. \]

Then the points \(x_i\) were taken uniformly on the concentric circle \(S_1^{(2)}\) of radius equal to \(1.1\), with step \(\pi/12\); the number of points is 24. Below are the values of the Gram determinants \(G_i\) \((i = 1, 2, \ldots, 24)\) for this case:

\[ \begin{aligned} G_1 &= 3.6 \quad &G_5 &= 80 \quad &G_9 &= 95 \quad &G_{13} &= 1.8 \\ G_2 &= 8.4 \quad &G_6 &= 130 \quad &G_{10} &= 38 \quad &G_{14} &= 0.7 \\ G_3 &= 28 \quad &G_7 &= 184 \quad &G_{11} &= 13 \quad &G_{15} &= 0.2 \\ G_4 &= 48 \quad &G_8 &= 236 \quad &G_{12} &= 4.7 \quad &G_{16} &= 8.9 \cdot 10^{-2} \end{aligned} \]

\[ \begin{aligned} G_{17} &= 1.6 \cdot 10^{-2} \quad &G_{21} &= 1.5 \cdot 10^{-5} \\ G_{18} &= 3.0 \cdot 10^{-3} \quad &G_{22} &= 2.9 \cdot 10^{-6} \\ G_{19} &= 5.1 \cdot 10^{-4} \quad &G_{23} &= 4.9 \cdot 10^{-7} \\ G_{20} &= 8.6 \cdot 10^{-5} \quad &G_{24} &= 4.9 \cdot 10^{-8}. \end{aligned} \]

We see that in this case the Gram determinant is considerably better conditioned.

Finally, the points \(x_i\) were taken uniformly distributed on the concentric circle \(S_1^{(3)}\) of radius \(1.05\), with the same step \(\pi/12\). The number \(n\) of points is 24. We give the corresponding values of the Gram determinants \(G_i\) \((i = 1, 2, \ldots, 24)\):

\[ \begin{aligned} G_1 &= 4.1 \quad &G_7 &= 4.7 \quad &G_{13} &= 2.9 \quad &G_{19} &= 0.5 \\ G_2 &= 4.7 \quad &G_8 &= 4.6 \quad &G_{14} &= 2.5 \quad &G_{20} &= 0.3 \\ G_3 &= 4.8 \quad &G_9 &= 4.3 \quad &G_{15} &= 2.0 \quad &G_{21} &= 0.2 \\ G_4 &= 4.9 \quad &G_{10} &= 4.0 \quad &G_{16} &= 1.6 \quad &G_{22} &= 0.1 \\ G_5 &= 4.9 \quad &G_{11} &= 3.7 \quad &G_{17} &= 1.2 \quad &G_{23} &= 0.08 \\ G_6 &= 4.8 \quad &G_{12} &= 3.3 \quad &G_{18} &= 0.8 \quad &G_{24} &= 0.03. \end{aligned} \]

The results of these numerical experiments show that, when the auxiliary contour \(S_1\) approaches the principal contour \(S\), the corresponding Gram determinant increases. Let us note, however, that the proof given above of the unreliability of systems of potential functions relies essentially on the constancy of the auxiliary contour \(S_1\). Otherwise, the functions \(r(x_{2i}, M)\) and \(\xi(M)\) may have the same order of smallness.

The second remark*) concerns the choice of the algorithm for orthonormalizing the system \(\{\omega_i(M)\}\). We shall show that, abandoning the method of orthonormalization by formula (1) and using the method according to which, if the orthonormalized elements \(\varphi_1, \varphi_2, \ldots, \varphi_{n-1}\) have been constructed, then \(\varphi_n\) is found from the relation

\[ \varphi_n = \frac{\overline{\varphi}_n}{\iint\limits_S \overline{\varphi}_n^{\,2}\, dS}, \tag{6} \]

where

\[ \overline{\varphi}_n = \omega_n + \sum_{k=1}^{n-1} (\omega_n, \varphi_k)\varphi_k, \qquad (\omega_n, \varphi_k) = \iint \omega_n \varphi_k\, dS, \]

we achieve a significant increase in the stability of the computational scheme. It is clear that if computations by both algorithms are carried out absolutely

*) This remark follows directly from [2], where it is proved that preliminary orthonormalization by formula (6) has a number of essential advantages over the direct solution of the Ritz system.

exact (with an infinite number of digits), then the results will also be identical. The results will also be sufficiently close in the case when the system being orthonormalized is reliable [2]. But, as was shown above, the system \(\{\omega_i(M)\}\) is not reliable, and for the practical realization of the approximate method it is extremely important to choose, from the two methods of orthonormalization, the one that gives the more stable computational scheme (ensures a larger number of correct digits).

Although computations by formula (1) are easier to organize on a universal machine, since ready-made standard subroutines exist for computing determinants, nevertheless, as will be shown below, the orthonormalization should be carried out by formula (6), since the corresponding algorithm is considerably more stable.

We shall regard the normalized functions \(\omega_i(M)\) as vectors (with initial point at 0) in the space \(L_2\), and denote by \(\alpha_k\) the angle between the vector \(\omega_k(M)\) and the hyperplane passing through the vectors \(\omega_1, \omega_2, \ldots, \omega_{k-1}\). It is known [2] that determinant (2) is equal to the square of the volume of the parallelepiped constructed on the vectors \(\omega_1, \omega_2, \ldots, \omega_n\)

\[ G_n=\prod_{i=1}^{n}\sin^2 \alpha_i. \]

Thus, for the denominator in formula (1) we obtain

\[ d_n=\sqrt{G_{n-1}\cdot G_n}=\sin \alpha_n\prod_{i=1}^{n-1}\sin^2 \alpha_i. \]

As for the sum

\[ \sum_{k=1}^{n-1}(\omega_n,\varphi_k)\varphi_k, \tag{7} \]

it represents [2] the projection of the element \(\omega_n\) onto the subspace of the vectors \(\varphi_1,\varphi_2,\ldots,\varphi_{n-1}\), or, in view of the equivalence of the subspaces of the vectors \(\varphi_1,\varphi_2,\ldots,\varphi_{n-1}\) and \(\omega_1,\omega_2,\ldots,\omega_{n-1}\), (7) represents the projection of \(\omega_n\) onto the subspace of the vectors \(\omega_1,\omega_2,\ldots,\omega_{n-1}\). Hence it is obvious that the denominator in formula (6) is equal to

\[ \iint_S \varphi_n^{-2}\,dS=\sin \alpha_n, \]

and, when orthonormalizing \(n\) elements, in all it will be necessary to divide by

\[ \bar d_n=\prod_{i=1}^{n}\sin \alpha_i. \]

For the ratio \(d_n:\bar d_n\) we obtain

\[ d_n:\bar d_n=\prod_{i=1}^{n-1}\sin \alpha_i=\bar d_{n-1}. \tag{8} \]

It follows that formula (1) gives a considerably less stable computational scheme in comparison with formula (6).

We attempted*) to orthonormalize the system \(\{\ln r(x_i,M)\}\), \(i=1,2,\ldots,28\), \(x_i\in S_1^{(1)}\), by means of formula (6); however, an emergency stop of the machine occurred. It turned out that this happened during division by

*) All numerical calculations were carried out by the scientific staff members of the Computing Center of the Academy of Sciences of the Georgian SSR, N. Lekishvili and N. Arveladze.

\((\overline{\varphi}_{21}, \overline{\varphi}_{21})\). This fact (taking into account that all Gram determinants of tenth order are equal to zero) agrees well with formula (8), from which it follows that formula (6) makes it possible, for a fixed number of digits with which the numerical calculations are performed, to orthonormalize approximately twice as many functions as formula (1). It is clear that if the orthonormalization of a sufficiently unreliable system can be carried out both by formula (6) and by formula (1), then the latter will give a considerably cruder result. To confirm this, the system \(\{\ln r(x_i,M)\}\), \(x_i \in S_1^{(2)}\) \((i=1,2,\ldots,24)\), was orthonormalized. The table below gives the orthonormalization coefficients \(A_{22,i}^{(1)}\) for the function \(\varphi_{22}\), obtained by formula (1), and the same coefficients \(A_{22,i}^{(2)}\), obtained by formula (6). In the third and fourth columns are given

\[ a_{22,i}^{(1)}=\int_{S_1^{(2)}} \varphi_{22}^{(1)} \varphi_i^{(1)}\,dS \quad\text{and}\quad a_{22,i}^{(2)}=\int_{S_1^{(2)}} \varphi_{22}^{(2)} \varphi_i^{(2)}\,dS, \]

where \(\varphi_i^{(1)}\) and \(\varphi_i^{(2)}\) are the orthonormalized functions obtained respectively with the aid of formulas (1) and (6); we see that formula (6) gives a considerably more accurate orthonormalization than formula (1).

References

1. Kupradze V. D. Potential Methods in the Theory of Elasticity. Moscow, Fizmatgiz, 1963.
2. Mikhlin S. G. Izv. Vuzov, Mathematics, No. 5 (6), 1958.
3. Kupradze V. D. On the completeness of certain classes of functions. Reports of the Academy of Sciences of the Georgian SSR, vol. XXXVII: 2, 1965.
4. Akhiezer N. I. Lectures on Approximation Theory. Moscow, Fizmatgiz, 1965.

Received by the editors
March 22, 1966

Computing Center of the Academy of Sciences of the Georgian SSR