ON DIFFERENCE METHODS FOR THE SOLUTION OF BOUNDARY VALUE PROBLEMS FOR DIFFERENTIAL EQUATIONS OF ELLIPTIC TYPE

A. V. BULEDZA

Reducing a boundary value problem for a differential equation to some difference problem is, at present, apparently the simplest and most effective method for solving boundary value problems for equations of elliptic type.

Although the possibilities for solving the resulting linear algebraic systems of high order are quite considerable with the availability of modern computing devices, they remain limited for a number of reasons.

Let us make several general remarks concerning this question.

1. The matrices of difference operators are often poorly conditioned; moreover, their conditioning worsens as the number of grid nodes increases.

2. In view of 1, the solution of systems of difference equations by direct methods, even when their order is comparatively low, is only slightly effective or altogether impossible because of the accumulation of large computational errors.

3. Iterative methods requiring the computation of a large number of scalar products (the conjugate-gradient method, the method of steepest descent, etc.) are too cumbersome, and their implementation on a computer is not always possible.

4. Iterative methods based on the principle of successive substitutions are the simplest to program and require the least load on the machine memory. This is their great advantage. However, slowly convergent iterative methods often prove practically unsuitable because of loss of stability (in the sense of the definition given in [1]).

5. The use of successive-approximation algorithms with a high rate of convergence, while preserving the simplicity of their implementation on a computer, gives good accuracy of the approximate solution.

Let us illustrate what has been said by an example.

Consider a system of linear algebraic equations of the 20th order

\[ A X = F, \tag{1} \]

where \(A\) is a somewhat transformed matrix of a difference operator approximating the boundary value problem

\[ -\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)=f(x,y), \tag{2} \]

\[ \left.\frac{\partial u}{\partial x}\right|_{x=0} = \left.\frac{\partial u}{\partial x}\right|_{x=a} = \left.\frac{\partial u}{\partial y}\right|_{y=0} =0,\quad \left.u\right|_{y=b}=0 \tag{3} \]

in the rectangle \(\{0\le x\le a,\ 0\le y\le b\}\) (we do not give the values of the elements of the matrix \(A\) and of the components of the vector \(F\)).

The matrix of system (1) is ill-conditioned; for it

\[ \alpha \simeq 281359\cdot 10^{-9}.* \]

System (1) was solved by the Gauss method (compact scheme) on the “Strela” computer. The error of the obtained solution turned out to be approximately equal to 63 (in the first norm), i.e. the result produced by the machine was very far from the truth.

Then system (1) was solved on “Strela” by the formula

\[ X_n=X_{n-1}-\frac{2}{M+m}\,(AX_{n-1}-F), \tag{4} \]

proposed by I. P. Natanson [6].

After 1110 iterations by formula (4), a result was obtained with an error of the same order as in the computation by the Gauss method. The rate of convergence of algorithm (4) has the order of a geometric progression with ratio

\[ q=\frac{M-m}{M+m}\simeq 0.9994 \]

(for our problem), which explains the large error of the solution.

Finally, the solution of system (1) was found by the algorithm proposed by us [2]

\[ x_n=x_{n-1}+\left(\frac{1-\sqrt{\alpha}}{1+\sqrt{\alpha}}\right)^{10}(x_{n-1}-x_{n-2}) +\frac{4}{(1+\sqrt{\alpha})^{10}M}\times \]

\[ \times\sum_{k=0}^{4} c_k A^{4-k}(AX_{n-1}-F), \tag{5} \]

\[ c_0=-\frac{256}{M^4},\quad c_1=\frac{640(1+\alpha)}{M^3},\quad c_2=-\frac{560\alpha^2+1440\alpha+560}{M^2}, \]

\[ c_3=\frac{200\alpha^3+1080\alpha^2+1080\alpha+200}{M}, \]

\[ c_4=-(25\alpha^4+300\alpha^3+630\alpha^2+300\alpha+25), \]

whose rate of convergence has the order

\[ \gamma=\left(\frac{1-\sqrt{\alpha}}{1+\sqrt{\alpha}}\right)^5\simeq 0.8455. \]

After 108 iterations, the error of the obtained solution was (in the first norm) of order \(10^{-6}\). The advantages of algorithm (5) are obvious.

We note that universal algorithms of type (4), (5) are applicable only to systems with symmetric and positive definite matrices. This circumstance is not an essential limitation, since the most commonly used difference analogues of diffe-

* We shall characterize the conditioning of systems of linear algebraic equations (matrices) by the \(\alpha\)-number defined by the equality
\[ \alpha=\frac{\min|\lambda_i|}{\max|\lambda_i|}, \]
where \(\lambda_i\) ranges over the eigenvalues of the matrix.

ential operators of elliptic type are positive definite.

But the algorithms mentioned have another, more serious drawback: for their construction, information about the location of the spectrum of the difference operator is necessary. Nevertheless, in many cases these difficulties too can be successfully overcome.

Below we give formulas for determining the spectra of some difference approximations of the harmonic operator under various boundary conditions, using square grids for the rectangle \(0 \leq x \leq a = Mh\), \(0 \leq y \leq b = Nh\). From these formulas all eigenvalues of the difference operators can be computed exactly (of course, within the capabilities of the computing device). Following V. E. Milne [5], one can determine bounds for the location of the spectrum also for domains of a more general form:

I.

\[ -\Delta u=\lambda u, \]

\[ -\Delta_h u_{m,n}=\frac{1}{h^2}\left(4u_{m,n}-u_{m-1,n}-u_{m+1,n}-u_{m,n-1}-u_{m,n+1}\right). \]

\(1^\circ.\)

\[ u\big|_S=0, \]

\[ \lambda_h=\frac{4}{h^2}\left(\sin^2\frac{\pi m}{2M}+\sin^2\frac{\pi n}{2N}\right), \tag{6} \]

\[ m=1,2,\ldots,M-1,\quad n=1,2,\ldots,N-1\quad \text{(see [7])}. \]

\(2^\circ.\)

\[ \left.\frac{\partial u}{\partial n}\right|_S=0,\qquad \int_S \frac{\partial u}{\partial n}\,dS=0. \]

\[ \lambda_h=\frac{4}{h^2}\left(\sin^2\frac{\pi m}{2(M-1)}+\sin^2\frac{\pi n}{2(N-1)}\right), \tag{7} \]

\[ m=0,1,\ldots,M-2,\quad n=0,1,\ldots,N-2\quad \text{(see [1])}. \]

\(3^\circ.\)

\[ \left[\frac{\partial u}{\partial n}+\sigma u\right]_S=0, \]

\[ \lambda_h=\frac{4}{h^2}\left(\sin^2\frac{\theta_m}{2}+\sin^2\frac{\theta_n}{2}\right), \tag{8} \]

\[ m=1,2,\ldots,M-1,\quad n=1,2,\ldots,N-1; \]

\(\theta_m,\theta_n\) are the nonzero solutions of the equations:

\[ r^2\sin M\theta_m-2r\sin(M-1)\theta_m+\sin(M-2)\theta_m=0, \tag{9} \]

\[ r^2\sin N\theta_n-2r\sin(N-1)\theta_n+\sin(N-2)\theta_n=0, \tag{10} \]

\[ r=1+\sigma h,\quad \sigma=\mathrm{const}\quad \text{(see [3])}. \]

\(4^\circ.\)

\[ \left[\frac{\partial u}{\partial n}+\sigma u\right]_{x=a} = \left[\frac{\partial u}{\partial n}+\sigma u\right]_{y=b} =0, \]

\[ u\big|_{x=0}=u\big|_{y=0}=0, \]

\[ \lambda_h=\frac{4}{h^2}\left(\sin^2\frac{\theta_m}{2}+\sin^2\frac{\theta_n}{2}\right); \tag{11} \]

\[ m=1,2,\ldots,M-1,\quad n=1,2,\ldots,N-1; \]

\(\theta_m, \theta_n\) are nontrivial solutions of the equations:

\[ r \sin M \theta_m - \sin (M - 1)\theta_m = 0, \tag{12} \]

\[ r \sin N \theta_n - \sin (N - 1)\theta_n = 0. \tag{13} \]

\(5^\circ.\)

\[ \left.\frac{\partial u}{\partial n}\right|_{x=0} = \left.\frac{\partial u}{\partial n}\right|_{y=0} =0,\qquad \left.u\right|_{x=a} = \left.u\right|_{y=b} =0, \]

\[ \lambda_h = \frac{4}{h^2} \left( \sin^2 \frac{\pi(2m+1)}{2(2M-1)} + \sin^2 \frac{\pi(2n+1)}{2(2N-1)} \right), \tag{14} \]

\[ m=0,\ 1,\ \ldots,\ M-2,\qquad n=0,\ 1,\ \ldots,\ N-2. \]

\(6^\circ.\)

\[ \left.\frac{\partial u}{\partial n}\right|_{x=0} = \left.\frac{\partial u}{\partial n}\right|_{x=a} = \left.\frac{\partial u}{\partial n}\right|_{y=0} =0,\qquad \left.u\right|_{y=b} =0, \]

\[ \lambda_h = \frac{4}{h^2} \left( \sin^2 \frac{\pi m}{2(M-1)} + \sin^2 \frac{\pi(2n+1)}{2(2N-1)} \right), \tag{15} \]

\[ m=0,\ 1,\ \ldots,\ M-2,\qquad n=0,\ 1,\ \ldots,\ N-2. \]

\(7^\circ.\)

\[ \left[ \frac{\partial u}{\partial n}+\sigma u \right]_{x=0} = \left[ \frac{\partial u}{\partial n}+\sigma u \right]_{y=0} =0,\qquad \left.\frac{\partial u}{\partial n}\right|_{x=a} = \left.\frac{\partial u}{\partial n}\right|_{y=b} =0, \]

\[ \lambda_h = \frac{4}{h^2} \left( \sin^2 \frac{\theta_m}{2} + \sin^2 \frac{\theta_n}{2} \right); \tag{16} \]

\(\theta_m, \theta_n\) are roots of the equations (except the zero roots):

\[ r \sin M\theta_m - (r+1)\sin(M-1)\theta_m + \sin(M-2)\theta_m = 0, \tag{17} \]

\[ r \sin N\theta_n - (r+1)\sin(N-1)\theta_n + \sin(N-2)\theta_n = 0. \tag{18} \]

II.

\[ -\Delta u=\lambda u, \]

\[ -\Delta_h u_{m,n} = \frac{1}{h^2} \left( 4u_{m,n} -u_{m-1,n+1} -u_{m-1,n-1} -u_{m+1,n+1} -u_{m+1,n-1} \right). \]

\(8^\circ.\)

\[ \left.u\right|_S=0, \]

\[ \lambda_h = \frac{4}{h^2} \left( 1-\cos\frac{\pi m}{M}\cos\frac{\pi n}{N} \right), \tag{19} \]

\[ m=1,\ 2,\ \ldots,\ M-1,\qquad n=1,\ 2,\ \ldots,\ N-1. \]

III.

\[ -\Delta u=\lambda u, \]

\[ -\Delta_h u_{m,n} = \frac{1}{6h^2} \left[ 20u_{m,n} - 4\left( u_{m-1,n} + u_{m,n+1} + u_{m+1,n} + u_{m,n-1} \right) - \left( u_{m-1,n+1} + u_{m+1,n+1} + u_{m+1,n-1} + u_{m-1,n-1} \right) \right]. \]

\(9^\circ.\)

\[ \left.u\right|_S=0, \]

\[ \lambda_h = \frac{2}{3h^2} \left[ 9 - \left( 2+\cos\frac{\pi m}{M} \right) \left( 2+\cos\frac{\pi n}{N} \right) \right], \tag{20} \]

\[ m=1,\ 2,\ \ldots,\ M-1,\qquad n=1,\ 2,\ \ldots,\ N-1 \quad (\text{see }[5]), \]

\(10^\circ.\)
\[ \left.\frac{\partial u}{\partial n}\right|_{S}=0,\qquad \int_S \frac{\partial u}{\partial n}\,ds=0, \]

\[ \lambda_h=\frac{2}{3h^2}\left[9-\left(2+\cos\frac{\pi m}{M-1}\right)\left(2+\cos\frac{\pi n}{N-1}\right)\right], \tag{21} \]

\[ m=0,1,\ldots,M-2,\quad n=0,1,\ldots,N-2. \]

\(11^\circ.\)
\[ \left[\frac{\partial u}{\partial n}+\sigma u\right]_S=0, \]

\[ \lambda_h=\frac{2}{3h^2}\left[9-(2+\cos\theta_m)(2+\cos\theta_n)\right], \tag{22} \]

where \(\theta_m,\theta_n\) are determined from equations (9), (10).

\(12^\circ.\)
\[ \left.\frac{\partial u}{\partial n}\right|_{x=0} = \left.\frac{\partial u}{\partial n}\right|_{y=0} =0,\qquad \left.u\right|_{x=a}=\left.u\right|_{y=b}=0, \]

\[ \lambda_h=\frac{2}{3h^2}\left[9-\left(2+\cos\frac{\pi(2m+1)}{2M-1}\right)\left(2+\cos\frac{\pi(2n+1)}{2N-1}\right)\right], \tag{23} \]

\[ m=0,1,\ldots,M-2,\quad n=0,1,\ldots,N-2. \]

References

1. Bondarenko P. S. Studies of computational algorithms for the approximate integration of differential equations by the finite-difference method. Publishing House of Kiev University, 1962.

2. Buledza A. V. Bulletin of Kiev University, ser. mathematics and mechanics, No. 4, issue 1, 1961.

3. Buledza A. V. Reports and Communications of Uzhgorod University, ser. phys.-math. and hist. sciences, No. 5, 1962.

4. Varga R. and Forsythe G. Difference methods for solving differential equations in partial derivatives. Foreign Literature Publishing House, 1963.

5. Milne W. E. Numerical solution of differential equations, Foreign Literature Publishing House, 1955.

6. Natanson I. P. Scientific Notes of the Leningrad State Pedagogical Institute, 64, 1948.

7. Ryabenkii V. S. and Filippov A. F. On the stability of difference equations. GITTL, 1956.

Received by the editors
31 October 1964

Uzhgorod
State University