UDC 518.1

MATHEMATICAL PHYSICS

Corresponding Member of the USSR Academy of Sciences L. A. VAINSTEIN, M. G. BELKINA

THE METHOD OF DOUBLE REDUCTION AND INFINITE SYSTEMS OF LINEAR EQUATIONS FOR THE EXPANSION COEFFICIENTS OF THE UNKNOWN FUNCTION WITH SINGULARITIES

The numerical solution of many boundary-value problems of mathematical physics can be reduced to the solution of an infinite system of linear equations

\[ \sum_{s=0}^{\infty} A_{rs} X_s = C_r,\qquad r=0,1,2,\ldots, \tag{1} \]

where \(X_s\) are the expansion coefficients of the unknown function

\[ \Psi(x)=\sum_{s=0}^{\infty} X_s \psi_s(x) \tag{2} \]

in a complete system of functions \(\psi_s(x)\). The unknown function \(\Psi(x)\) may be, for example, the current density on an ideally conducting surface constituting part of the field on some auxiliary surface that separates two regions of simple form, etc.

System (1) is usually solved by the reduction method, i.e., it is reduced to a finite system

\[ \sum_{s=0}^{S-1} A_{rs} X_s = C_r,\qquad r=0,1,\ldots,S-1, \tag{1′} \]

and, accordingly, the series (2) is replaced by the finite sum

\[ \Psi(x)=\sum_{s=0}^{S-1} X_s \psi_s(x), \tag{2′} \]

where \(S\) is a sufficiently large number. However, the ordinary reduction method becomes inefficient if the function \(\Psi(x)\) has singularities that are absent from the functions \(\psi_s(x)\); these singularities are conveyed by the slow convergence of the series (2), and, by replacing it with the finite sum (2′), we lose these singularities and do not obtain good results from the system (1′) even for the initial coefficients \(X_0, X_1,\ldots\).

Thus, for example, if the current density on an ideally conducting finite circular cylinder is represented in the form of a trigonometric series and the system for the coefficients of this series (¹) is solved by the reduction method, the convergence proves to be very slow. To accelerate the convergence it was necessary to resort to a new formulation of the problem—expansion of the current density in functions \(\psi_s\) that have at the ends of the cylinder the same singularity as the density itself (²). In the theory of diffraction of a plane wave by an infinite periodic grating, the expansion of the field into plane waves (proper and improper) is natural. However, for a grating of infinitely thin strips such an expansion converges slowly near the grating, and the infinite system (1) obtained for the coefficients of this expansion practically cannot be solved by the reduction method (cf. (³)). This circumstance led to the appearance of works by Agranovich, Marchenko, and Shestopalov

... and Malin(\(^{5}\)), in which a different formulation of the problem is given—the system (1) is transformed into another one, which is then solved numerically.

It turns out, however, that in these and other cases there is no need to reformulate the problem, and one can solve system (1) by applying the method of double reduction, to the exposition of which we now proceed. The point is that the singularities of the function (2) in boundary-value problems arise because the region of space in which we study the field is bounded by a surface with edges and corner points; the singularities of the field near these are known a priori, and therefore all the singularities of the function \(\Psi(x)\) are known, as is the character of the asymptotic series associated with them, determining the behavior of the coefficients \(X_s\) for large \(s\). It may, for example, have the form

\[ X_s=\sum_{j=0}^{\infty}\frac{\rho_j}{s^{\gamma+j}}, \tag{3} \]

where the coefficients \(\rho_0,\rho_1,\ldots\) are unknown, while the exponent \(\gamma>0\) is known in advance.

The method of double reduction consists in the following: the first \(S\) coefficients \(X_s\) (for \(s=0,1,\ldots,S-1\)) of the series (2) are taken into account exactly, as in the usual reduction method, but the remaining coefficients are not set equal to zero; rather, they are replaced by the asymptotic expression

\[ X_s=\sum_{j=0}^{J-1}\frac{\rho_j}{s^{\gamma+j}}, \tag{4} \]

i.e., by the first \(J\) terms of the series (3). Thus, we take into account exactly \(S\) coefficients \(X_s\) of the series (2) and \(J\) coefficients \(\rho_j\) of the series (3), i.e., we perform, as it were, a double reduction. System (1) then takes the form

\[ \sum_{s=0}^{S-1} A_{rs}X_s+\sum_{j=0}^{J-1} B_{rj}\rho_j=C_r,\qquad r=0,1,\ldots,S+J-1, \tag{5} \]

where the elements of the additional matrix

\[ B_{rj}=\sum_{s=S}^{\infty}\frac{A_{rs}}{s^{\gamma+j}},\qquad r=0,1,\ldots,S+J-1,\quad j=0,1,\ldots,J-1, \tag{6} \]

are represented by slowly convergent series; however, all terms of these series are known, and therefore they can be computed, improving convergence by one method or another. We have taken a finite number of coefficients \(X_s\) and \(\rho_j\), and therefore we can satisfy only a finite number of equations (1). The choice of the first \(S+J\) equations is natural, although arbitrary.

In order to obtain an idea of the effectiveness of the method of double reduction, we applied this method to a problem that admits an exact solution(\(^{3}\)). Let the function \(\Phi=\Phi(x,z)\) satisfy the two-dimensional Helmholtz equation for \(z>0\) and the boundary conditions

\[ \begin{gathered} \partial\Phi/\partial z=0 \quad \text{for } (2m-1)d<x<2md,\\ z=0,\quad m=0,\pm1,\pm2,\ldots\\ \Phi=0 \quad \text{for } 2md<x<(2m+1)d, \end{gathered} \tag{7} \]

If a plane wave is normally incident on the plane \(z=0\), then the function \(\Phi\) in the half-space \(z\geqslant 0\) has the form

\[ \Phi(x,z)=e^{-ikz}+\sum_{n=-\infty}^{\infty} R_n e^{ik(z\cos\varphi_n+x\sin\varphi_n)}. \tag{8} \]

Indeed, conditions (7) correspond to a periodic structure (lattice) with period \(2d\) along the \(x\)-axis; the angles \(\varphi_n\) are determined by the relations

\[ \sin\varphi_n = n/2q,\qquad q = kd/2\pi = d/\lambda, \]

\[ \cos\varphi_n = \sqrt{1-\left(\frac{n}{2q}\right)^2} = i\,\frac{|n|}{2q}\sqrt{1-\left(\frac{2q}{n}\right)^2}, \tag{9} \]

and the complex coefficients \(R_n\) are unknown; from symmetry considerations it can be shown that \(R_{-n}=(-1)^n R_n\).

Table 1

For the coefficients \(R_n\) one may write a system of equations (see (3), formula (55.04)), which is readily transformed to the form (1), where

\[ X_s = R_{2s+1},\qquad C_r = -\delta_{0r}, \]

\[ A_{rs}=\frac{i}{\pi}\left(\cos\varphi_{2s+1}+\cos\varphi_{2r}\right) \left(\frac{1}{2s+1-2r}+\frac{1}{2s+1+2r}\right), \tag{10} \]

and the even coefficients \(R_{2m}\) are determined by the formula

\[ R_{2m}=-\delta_{0m}-\frac{2i}{\pi}\left[ \sum_{s=0}^{S-1} \left(\frac{1}{2s+1-2m}+\frac{1}{2s+1+2m}\right)R_{2s+1} +\right. \]

\[ \left. +\sum_{j=0}^{J-1}\rho_j\sum_{s=S}^{\infty} \left(\frac{1}{2s+1-2m}+\frac{1}{2s+1+2m}\right) \frac{1}{s^{\gamma+j}} \right]. \]

The exponent \(\gamma\) can be determined as follows. As is seen from formulas (8) and (9), the coefficients \(R_n\) are the expansion coefficients of the function \(\Phi(x,0)-1\) in the series in the functions \(e^{i\pi nx/d}\), i.e., in a complex Fourier series. The function \(\Phi(x,0)\) itself can be represented in the form

\[ \Phi(x,0)=f(x)p(x), \tag{11} \]

where \(f(x)\) is an analytic function with period \(2d\), and

\[ p(x)=\sqrt{(2md-x)\,[(2m-1)d+x]} \qquad \text{for } (2m-1)d < x < 2md, \]

\[ p(x)=0 \qquad \text{for } 2md < x < (2m+1)d \tag{12} \]

is a “weight function” that takes into account all singularities of \(\Phi(x,0)\). Therefore

\[ R_n=\frac{1}{2}\int_0^1 e^{i\pi n\xi}\,\tilde f(\xi)\sqrt{\xi(1-\xi)}\,d\xi = \frac{i}{2\pi n}\int_0^1 e^{i\pi n\xi}\, \frac{d}{d\xi}\!\left[\tilde f(\xi)\sqrt{\xi(1-\xi)}\right]\,d\xi, \]

where \(\tilde f(\xi)=f(-\xi d)\), and \(n=\pm1,\pm2,\ldots\). For large positive \(n\), the last integral can be transformed into a sum of integrals from \(0\) to \(i\sigma\), from \(i\sigma\) to \(1+i\sigma\), and from \(1+i\sigma\) to \(1\) (\(\sigma>0\)), and we obtain

\[ R_n=\frac{1}{2\pi n} \left\{ \int_0^\sigma e^{-\pi nt} \left[f_0(t)-(-1)^n f_1(t)\right]\frac{dt}{\sqrt{t}} +O\!\left(e^{-\pi\sigma n}\right) \right\}, \]

where \(f_0\) and \(f_1\) are functions expandable, for small \(t\), in power series in \(t\). Hence, for \(n=2s+1\), one obtains the series (3) with exponent \(\gamma=3/2\). An analogous series, but with other coefficients \(\rho_j\), is obtained for \(n=2s\).

In the computations we took \(S=6\), \(J=4\) and obtained (on the Ural-2 machine) the results presented in Table 1. It is seen from the table that these values of \(S\) and \(J\) give 5–6 correct decimal places for \(q \ge 1\) and considerably higher accuracy for small \(q\). Obviously, to obtain the same accuracy for larger \(q\), one must take a higher order of the system.

Table 2

For \(q=0.75\), a study of other variants of the method was carried out (see Table 2). This study shows that decreasing \(J\) (the number of asymptotic coefficients), with the same number of equations \(S+J\), reduces the accuracy of the method.

As this example shows, the method of double reduction is very effective and at the same time does not require cumbersome derivations and transformations; certain difficulties are caused only by the computation of the series (6), which we carried out by direct summation of the initial segment of the series and by applying the Euler–Maclaurin formula to its remainder. At the same time, this method is sufficiently general and, in particular, can be extended to infinite systems that include unknown coefficients of different types, with different asymptotic expansions.

When this work was already nearing completion, we became acquainted with the American paper (6), in which related ideas are developed. Since in paper (6) the system of linear equations for the expansion coefficients is formulated by imposing boundary conditions at a finite number of points, the derivation of the basic relations becomes more cumbersome, and the results, with a much larger number of equations (51), are less accurate (essentially the same problem was solved there). It seems to us more expedient to proceed from system (1), obtained from the boundary conditions on the entire boundary.

We thank N. I. Lesik for his interest in our work.

Institute for Physical Problems named after S. I. Vavilov
Academy of Sciences of the USSR
Moscow

Received
26 VI 1970

CITED LITERATURE

1. P. L. Kapitsa, V. A. Fok, L. A. Vainshtein, ZhTF, 29, No. 10, 1188 (1959).
2. L. A. Vainshtein, ZhTF, 37, No. 7, 1181 (1967).
3. L. A. Vainshtein, Theory of Diffraction and the Factorization Method, Moscow, 1966, §§ 52 and 55.
4. E. S. Agranovich, V. A. Marchenko, V. P. Shestopalov, ZhTF, 32, No. 4, 381 (1962).
5. V. V. Malin, Radio Engineering and Electronics, 8, No. 2, 211 (1963).
6. A. R. Neureuther, K. Zaki, Radio Sci., 3, No. 12, 1158 (1968).