THE METHOD OF INTEGRAL RELATIONS FOR EQUATIONS AND SYSTEMS OF HYPERBOLIC TYPE

(A survey of convergence studies and error estimates)

V. V. BOBKOV, V. I. KRYLOV

With the advent of electronic digital computers, approximate methods for solving partial differential equations began to develop rather intensively. Of great practical importance among these methods are those which approximately reduce partial differential equations to systems of ordinary differential equations or to systems of algebraic or transcendental equations. One such method is the method of integral relations, proposed in the early 1950s by A. A. Dorodnitsyn as a development of the well-known direct method. This method is one of the universal methods for solving partial differential equations and is applicable to equations of various types.

The method of integral relations was originally proposed [1] in the following form:

“Let, in a certain rectangular domain \(a \le x \le b,\ c \le y \le d\), which in special cases may turn into a strip \((a=-\infty,\ b=+\infty)\) or into a half-strip \((b=\infty)\), there be given a system of equations of ‘divergence’ form

\[ \frac{\partial}{\partial x} P_i(x,y;u_1,\ldots,u_m) + \frac{\partial}{\partial y} Q_i(x,y;u_1,\ldots,u_m) = \]

\[ = F_i(x,y;u_1,\ldots,u_m)\quad (i=1,2,\ldots,m). \tag{1} \]

Divide the domain into \(N\) strips by the straight lines
\(y=y_0=c,\ y=y_1,\ldots,\ y=y_N=d\), and integrate each of the equations of system (1) across each strip.

In this way we obtain a system of integral relations

\[ \frac{\partial}{\partial x} \int_{y_n}^{y_{n+1}} P_i\,dy + Q_{i,n+1}-Q_{i,n} = \int_{y_n}^{y_{n+1}} F_i\,dy \quad (i=1,2,\ldots,m;\ n=0,1,\ldots,N-1). \tag{2} \]

Here \(Q_{i,n}\) is the value of the function \(Q_i\) on the line \(y=y_n\). If now some interpolation formula is applied to the functions \(P_i\) and \(F_i\), expressing the values of these functions for any \(y\) through their values on the lines \(y=y_n\), then each of the integrals entering into system (2) will be represented in the form

\[ \int_{y_n}^{y_{n+1}} P_i\,dy \simeq (d-c)\sum_{k=0}^{N} A_{n,k}P_{i,k}, \tag{3} \]

where \(A_{n,k}\) are certain constants depending on the form of the interpolation formula. Substitution of their approximate expressions (3) in place of the integrals transforms the system of integral relations (2) into a system of ordinary differential equations with respect to the unknowns \(u_{i,n}\).

Together with the boundary conditions

\[ \text{for } y=c \qquad \varphi_\nu(x; u_1,\ldots,u_m)=0 \quad (\nu=1,2,\ldots,\alpha), \]

\[ \text{for } y=d \qquad \varphi_\nu(x; u_1,\ldots,u_m)=0 \quad (\nu=\alpha+1,\ldots,m) \]

this system will have the number of unknowns exactly equal to the number of equations*.

The division of the domain of integration into strips in the method of integral relations can, generally speaking, be carried out by curved lines; moreover, they need not necessarily be equidistant, and the method itself can be applied to domains of arbitrary form (including the case of unknown boundaries). In the method, “through” interpolation is usually used with respect to the values of the initial functions on the boundary lines of all the strips. One may also use “piecewise” interpolation. In this case, for example, across each individual strip the functions may be interpolated linearly, across two neighboring strips parabolically, and so on. In this way one can achieve a considerable simplification of the form of the approximating system. In constructing such a system it is also necessary to pay attention to the choice of the direction of approximation and its form. In doing so, one must take into account that the solution obtained here gives, in the direction of integration of the approximating system, a more accurate distribution of the values of the desired functions than in the direction in which they are approximated by interpolation.

A distinctive feature of the method of integral relations is the fact that in this method it is, in essence, the integral that is approximated. As is known, an integral is a smoother function than the integrand expression; consequently, with a small number of interpolation nodes a “good” representation is attained for it more rapidly than for the integrand function. In addition, the integral will admit a continuous representation also in the case when the integrand function has a discontinuity of the first kind.

In the method of integral relations, the original equations are usually taken in the divergent form (1). This form is convenient and useful in that, in this case, integration with respect to one of the variables is carried out exactly.

Substantial difficulties in applying the method of integral relations (especially for nonlinear equations) arise when, for an approximating system of high order, one has to solve a boundary-value problem. In this case the method is effective if it gives a sufficiently accurate solution already at a low order of the approximating system of ordinary differential equations.

To increase the effectiveness of the method, A. A. Dorodnitsyn has recently proposed [2] a generalized method of integral relations. Here, owing to the introduction of smoothing functions, which are chosen in accordance with the expected behavior of the desired functions, a system of ordinary differential equations is obtained, with respect-

significantly smoother than in the simple method. Owing to this, even with a small number of interpolation nodes one can obtain higher accuracy.

Let us set forth the idea of the generalized method of integral relations as applied to system (1) (see [2]).

Multiplying equations (1) by a certain function \(f(y)\) and carrying out integration over the entire interval \([c,d]\) of variation of the variable \(y\), we obtain integral relations which will serve as the starting point for the proposed method:

\[ \frac{d}{dx}\int_c^d P_i f(y)\,dy+Q_{i,d}f(d)-Q_{i,c}f(c)-\int_c^d Q_i f'(y)\,dy = \int_c^d F_i f(y)\,dy . \tag{4} \]

“Let us now choose a system of groups of functions \(f(y)\) \(\{[f_{1,1}], [f_{2,1}, f_{2,2}], \ldots, [f_{k,1}, f_{k,2}, \ldots, f_{k,k}], \ldots\}\) such that the \(k\)-th group contains \(k\) linearly independent functions (functions from different groups may coincide).

Next, in the \(k\)-th approximation, let us divide the whole domain into \(k\) strips and represent the functions \(P_i, Q_i, F_i\) by means of certain interpolation expressions in terms of their values on the boundaries of the strips \(y=y_\nu\) \((\nu=0,1,\ldots,k)\):

\[ P_i[x,y;u_1(x,y),\ldots,u_m(x,y)] \approx \sum_{\nu=0}^{k} P_{i,\nu}\psi_\nu(y), \tag{5} \]

\[ P_{i,\nu}=P_i[x,y_\nu;u_{1,\nu},\ldots,u_{m,\nu}],\qquad u_{\alpha,\nu}=u_\alpha(x,y_\nu). \]

Here \(\psi_\nu(y)\) is an interpolation polynomial, the specific form of which depends on the choice of interpolation formulas (of course, it may contain not only powers of \(y\), but also any transcendental functions).

If we now construct the integral relations (4) for each \(f_{k,n}\) from the \(k\)-th group and substitute their expressions (5) in place of \(P_i, Q_i, F_i\), then a system of \(km\) ordinary differential equations is obtained for \(m(k+1)\) unknown functions \(u_{\alpha,\nu}\), with \(m\) boundary conditions on the boundaries \(y=y_0=c\), \(y=y_k=d\) closing this system.”

Let us also add that if the integrand functions have singularities in the domain of integration, then the functions \(f_{k,n}\), naturally, must be chosen in such a way as to ensure convergence of all integrals in the integral relations.

Another form of the method of integral relations is also possible (see [3]). In two-dimensional problems the integration domain may be divided not into strips, but into subdomains, and integration may be carried out in two directions, with the functions approximated in two variables. In this case the original problem is approximately replaced by a system of algebraic or transcendental equations.

The method of integral relations can also be generalized to the case of three-dimensional partial differential equations. Here two approaches are possible (see [3]). In the first approach, the original system of equations is multiplied by a certain function of two of the variables under consideration and integrated with respect to these two variables. Some quadrature formulas are applied to the resulting double integrals, and one arrives at a system of ordinary differential equations in the third

variable. In the second variant, the original system of equations is first integrated with respect to one of the variables, and the integrand functions are represented with respect to this variable by some interpolating expressions whose coefficients will depend on the two remaining variables. The resulting approximating system will be a system of partial differential equations in these two variables, and it can be solved by the method of integral relations developed for two-dimensional problems.

At present the method has proved itself well in practice, both in our country and abroad, in solving a whole series of gas-dynamic problems. A fairly complete survey of all these works, which appeared chiefly at the Computing Center of the Academy of Sciences of the USSR, may be found in the previously mentioned note [3] by O. M. Belotserkovskii and P. I. Chushkin. In addition to the extensive bibliography on this question given there, we also note here the works [4]—[10], [29]. The method of integral relations has also been successfully applied to other questions, in particular to the approximate solution of the general system of equations for short-range forecasting of the basic meteorological elements [11], to the solution of the problem of nonstationary filtration [12], etc.

A more complete idea of the method itself and of its applications can be obtained by consulting, in addition to the literature mentioned above, the works [13]—[21], [30].

The experience accumulated in the practical application of the method of integral relations has shown that, in most cases, the required accuracy was attained already with a small number of partition lines. However, the conditions for convergence of the method, generally speaking, had not been clarified. In all the works mentioned above only the practical convergence of the method was considered, and this can be judged only after computations have been carried out.

Below we shall dwell on some results obtained by us in the study of the convergence of the method of integral relations.

§ I. APPROXIMATE SOLUTION

OF SECOND-ORDER HYPERBOLIC EQUATIONS

BY REDUCING THEM BY THE METHOD OF INTEGRAL RELATIONS

TO A SYSTEM OF FIRST-ORDER ORDINARY DIFFERENTIAL EQUATIONS

As we have already noted, the method of integral relations has been applied to equations reduced to the form (1). In this section we shall consider only second-order hyperbolic equations and systems of equations of this kind. Taking into account the specific character of the equations, we shall not reduce them to systems of first-order equations of divergent form (1), but shall write them in canonical form, which is always possible under fairly general assumptions concerning the coefficients of these equations. This, in particular, makes it possible to lower the order of the approximating system and to avoid possible difficulties with boundary conditions.

1. Goursat problem. For simplicity of exposition, let us consider the case of a single linear equation of second order. We shall assume it to have been reduced to canonical form. Then the Goursat boundary-value problem for it takes, generally speaking, the following form:

\[ \frac{\partial^{2}u}{\partial x\,\partial y} = a(x,y)\frac{\partial u}{\partial x} + b(x,y)\frac{\partial u}{\partial y} + c(x,y)u + f(x,y); \tag{6} \]

\[ u(x,0)=\varphi(x),\qquad u(0,y)=\psi(y),\qquad \varphi(0)=\psi(0); \tag{7} \]

\[ 0\le x\le l',\qquad 0\le y\le l''. \tag{8} \]

Integrating equation (6) across each of the \(n\) \((1\le n\le N,\; N=[l''/h])\) strips into which, for a chosen step \(h>0\), the rectangle \(D_n(0\le x\le l',\; 0\le y\le y_n)\) is divided by the straight lines \(y=y_m=mh\) \((m=1,2,\ldots,n-1)\), and interpolating the sought function \(u\) linearly by its values on the boundaries of each strip, we obtain, for finding an approximate solution \(u_n=u_n(x)\) of problem (6), (7), (8) on the line \(y=y_n\) \((1\le n\le N,\; 0\le x\le l')\), an approximating Cauchy problem for a system of \(n\) linear ordinary differential equations:

\[ \left. \begin{aligned} &\alpha_{m+1}u'_{m+1}-\beta_m u'_m=\xi_{m+1}u_{m+1}+\eta_m u_m+f_{m+1},\quad m=0,1,\ldots,n-1,\\ &u_0(x)=\varphi(x),\qquad u_i(0)=\psi(y_i),\quad i=1,2,\ldots,n \end{aligned} \right\}. \tag{9} \]

Here

\[ \alpha_{m+1}=\alpha_{m+1}(x)=1-h\int_0^1 a(x,y_m+th)t\,dt, \]

\[ \beta_m=\beta_m(x)=1+h\int_0^1 a(x,y_m+th)(1-t)\,dt, \]

\[ \xi_{m+1}=\xi_{m+1}(x)=h\int_0^1\left[c(x,y_m+th)-\frac{\partial b(x,y_m+th)}{\partial y}\right]t\,dt+b(x,y_{m+1}), \]

\[ \eta_m=\eta_m(x)=h\int_0^1\left[c(x,y_m+th)-\frac{\partial b(x,y_m+th)}{\partial y}\right]\times \]

\[ \times(1-t)\,dt-b(x,y_m), \]

\[ f_{m+1}=f_{m+1}(x)=h\int_0^1 f(x,y_m+th)\,dt,\qquad t=\frac{y-y_m}{h}. \]

The corresponding Cauchy problem for the error of the method \(\gamma_m=\gamma_m(x)=u(x,y_m)-u_m\) has the analogous form:

\[ \left. \begin{aligned} &\alpha_{m+1}\gamma'_{m+1}-\beta_m\gamma'_m=\xi_{m+1}\gamma_{m+1}+\eta_m\gamma_m+r_{m+1},\qquad \gamma_0(x)=0,\\ &\gamma_{m+1}(0)=0,\quad m=0,1,\ldots,n-1 \end{aligned} \right\}, \tag{10} \]

where

\[ r_{m+1}=r_{m+1}(x)=-\frac12 h^3\int_0^1\left\{ \frac{\partial^3 u(x,\tilde y_{m+1})}{\partial x\,\partial y^2}a(x,y_m+th)+ \right. \]

\[ \left. +\frac{\partial^2 u(x,\tilde y_{m+1})}{\partial y^2} \left[c(x,y_m+th)-\frac{\partial b(x,y_m+th)}{\partial y}\right]\right\}\times \]

\[ \times(t-1)t\,dt,\qquad y_m<\tilde y_{m+1}<y_{m+1}. \]

Let us note that if \(a(x,y)\le 0\) in the rectangle \(D_n\), then \(\alpha_{m+1}(x)\ge 1\), \(0\le x\le l'\), \(m=0,1,\ldots,n-1\). For other \(a(x,y)\), it suffices to take \(h<2/A_n\) so that \(\alpha_{m+1}(x)\) \((0\le x\le l',\; m=0,1,\ldots,n-1)\) is positive. Here \(A_n=\max\limits_{D_n}|a(x,y)|\).

Reducing the latter system to a form solved with respect to the derivatives, using matrix notation, and applying the known results of S. M. Lozinskii [22], one can show that, under unlimited decrease of the step \(h\), the approximate solution \(u_n(x)\) in the rectangle \(D\) (8) converges uniformly to the exact solution \(u(x,y_n)\) with rate of order \(h^2\); moreover, the following estimate holds for the error of the method:

\[ |\gamma_n(x)| \le \frac{h^2 y_n \left[M_n^{(1,2)} A_n + M_n^{(2)}(C_n+B_n')\right]\lambda_n \alpha_{(n)}} {12\left[\xi_{(n)}\alpha_{(n)}+(y_n-h)\Delta_n\lambda_n\right]} \times \]

\[ \times \left( e^{\frac{\xi_{(n)}\alpha_{(n)}+(y_n-h)\Delta_n\lambda_n}{\alpha_{(n)}^2}x} -1 \right), \qquad 0\le x\le l',\quad n=1,2,\ldots,N. \tag{11} \]

Here

\[ C_n=\max_{D_n}|c(x,y)|,\qquad B_n=\max_{D_n}|b(x,y)|,\qquad B_n'=\max_{D_n}\left|\frac{\partial b(x,y)}{\partial y}\right|, \]

\[ M_n^{(1,2)}=\max_{D_n}\left|\frac{\partial^3 u(x,y)}{\partial x\,\partial y^2}\right|, \qquad M_n^{(2)}=\max_{D_n}\left|\frac{\partial^2 u(x,y)}{\partial y^2}\right|, \]

\[ \alpha_{(n)}=\min_{1\le i\le n}\ \min_{0\le x\le l'}|\alpha_i(x)|; \]

\[ 1+\frac12 hA_n \ge \alpha_{(n)} \ge \begin{cases} 1+\dfrac12 h a_n, & \text{when } a(x,y)\le 0 \text{ in } D_n,\\[6pt] 1-\dfrac12 hA_n, & \text{for other } a(x,y), \end{cases} \]

\[ a_n=\min_{D_n}|a(x,y)|,\qquad \lambda_n= \begin{cases} 1, & \text{when } a(x,y)\le 0 \text{ in } D_n,\\ \lambda_{(n)}^{\,n-1}, & \text{for other } a(x,y), \end{cases} \]

\[ \lambda_{(n)} = \max_{1\le i\le n}\ \max_{0\le x\le l'} \left|\frac{\beta_i}{\alpha_i}\right| \le \frac{1+\frac12 hA_n}{1-\frac12 hA_n}, \qquad \lambda_n\le e^{A_n y_n}, \]

\[ \Delta_n = \max_{1\le i\le n-1}\ \max_{0\le x\le l'} \frac1h|\alpha_i\eta_i+\beta_i\xi_i| \le A_nB_n+C_n+B_n'+\frac12 hA_n(C_n+B_n'), \]

\[ \xi_{(n)} = \max_{1\le i\le n}\ \max_{0\le x\le l'}|\xi_i(x)|, \qquad \xi_{(n)}\le B_n+\frac12 h(C_n+B_n'). \]

Let us note that for each particular value \(x=x_0\), \(0\le x_0\le l'\), estimate (11) can be somewhat sharpened if the extremal values with respect to \(x\) entering it are taken not on the segment \([0,l']\), but on the segment \([0,x_0]\).

The original second-order differential equations need not necessarily be linear. To obtain results similar to (11), it is sufficient that only one of the first-order partial derivatives of the unknown function enter linearly into the right-hand side of equation (6). In this case the approximating system of ordinary differential equations for the approximate solution, generally speaking, is not

will be linear, and in constructing the corresponding system for the error of the method one will have to use Lagrange’s theorem on finite increments.

The results obtained can also be generalized to systems of equations of the form

\[ \frac{\partial^2 u_i}{\partial x \partial y} = a_i(x,y)\frac{\partial u_i}{\partial x} + \sum_{j=0}^{p} \left[ b_{ij}(x,y)\frac{\partial u_j}{\partial y} + c_{ij}(x,y)u_j \right] + \]

\[ + f_i(x,y),\quad i=1,2,\ldots,p, \tag{12} \]

where the right-hand sides need not be linear (only the partial derivatives with respect to \(y\) must enter linearly).

If equation (6) does not contain first-order partial derivatives of the sought function \(u(x,y)\), i.e. has the form

\[ u_{xy}=c(x,y)u+f(x,y), \tag{13} \]

then estimate (11) is written much more simply:

\[ |\gamma_n(x)| \le \frac{h^2 M_n^{(2)}}{12} \left(e^{C_n y_n x}-1\right), \quad 0\le x\le l', \quad n=1,2,\ldots,N . \tag{14} \]

An analogous estimate in the case of the simplest scheme of the method of straight lines for problem (13), (7), (8) was obtained by induction by B. M. Budak [23]. Comparison of the latter, which guarantees only convergence of order \(h\), with estimate (14) speaks in favor of the method of integral relations.

If, for the representation of the sought function, one uses algebraic interpolation polynomials of degree \(k\) \((1\le k<n)\), then in the case of problem (13), (7), (8) one can in a similar way construct explicit computational schemes converging uniformly in the rectangle \(D\) with order \(h^{k+1}\) (under the natural requirement that the \(k-1\) initial values of the approximate solution \(u_1(x),\ldots,u_{k-1}(x)\) must be found beforehand with sufficient accuracy). In this case the error estimate of the method that interests us takes the following form:

\[ |\gamma_n(x)| \le \frac{ h^{k+1}M_n^{(k+1)}C_{(n)}\beta_k \left[y_n-h(k-1)\right] + hC_k^*\Sigma_0^*\delta+\delta' }{ C_{(n)} \left[ y_n\Sigma-h(2k-1)\Sigma+h\Sigma_1 \right] } \times \]

\[ \times \left\{ e^{C_{(n)}\left[y_n\Sigma-h(2k-1)\Sigma+h\Sigma_1\right]x} -1 \right\}, \quad 0\le x\le l', \quad N\ge n\ge 2k-1 \tag{15} \]

(for \(n<2k-1\) in (15) the bracket \([y_n\Sigma-h(2k-1)\Sigma+h\Sigma_1]\) must be replaced by \(h\Sigma_n^*\)),

where

\[ C_{(n)}= \max_{0\le x\le l'} \max_{(k-1)h\le y\le y_n} |c(x,y)|, \quad C_k^*= \max_{0\le x\le l'} \max_{(k-1)h\le y\le(2k-1)h} |c(x,y)|, \]

\[ M_n^{(k+1)} = \max_{D_n} \left| \frac{\partial^{k+1}u(x,y)}{\partial y^{k+1}} \right|, \quad \Sigma=\sum_{j=0}^{k}|\alpha_j^{(k)}|, \quad \Sigma_1=\sum_{i=0}^{k-1}\sum_{j=0}^{i}|\alpha_j^{(k)}|, \]

\[ \Sigma_n^* = \sum_{i=0}^{n-k}\sum_{j=0}^{i}|\alpha_j^{(k)}|, \quad \Sigma_n^*<\Sigma_1\ (n<2k-1), \quad \Sigma_0^* = \sum_{i=0}^{k-1}\sum_{j=1}^{k-i}|\alpha_{i+j}^{(k)}|, \]

\[ \alpha_i^{(k)} = \frac{(-1)^i}{(k-i)!\,i!} \int_{k-1}^{k} \frac{t(t-1)(t-2)\ldots(t-k)}{t-k+i}\,dt . \]

\[ \beta_k=-\frac{1}{(k+1)!}\int_{k-1}^{k} t(t-1)(t-2)\ldots(t-k)\,dt, \]

\[ \delta=\max_{1\le i\le k-1}\ \max_{0\le x\le l'}|\delta_i(x)|,\qquad \delta'=\max_{0\le x\le l'}|\delta'_{k-1}(x)|, \]
where \(\delta_i(x)\) \((i=1,2,\ldots,k-1)\) are the errors on the initial \(k-1\) straight lines.

The results obtained in the case of interpolation of any fixed order \(k>1\) can also easily be generalized to systems of equations of the form

\[ \frac{\partial^2 u_i}{\partial x\,\partial y} =F_i(x,y;u_1,u_2,\ldots,u_p),\qquad i=1,2,\ldots,p. \]

In conclusion of this subsection we note that, by observing a certain similarity between problems (9) and (10), it is not difficult also to obtain estimates for the approximate solution \(u_n(x)\) \((0\le x\le l',\ n=1,2,\ldots,N)\). Passing to the limit in the estimates obtained as \(h\to 0\), one can also find estimates for the exact solutions of all the linear problems considered above. As an example we give here only the simplest of these estimates, obtained in the case of the Goursat problem for equation (13) under zero boundary conditions:

\[ |u(x,y)|\le Fyx\,e^{Cyx},\qquad 0\le x\le l',\quad 0\le y\le l''. \tag{16} \]

Here
\[ F=\max_{D(x,y)}|f(x,y)|,\qquad C=\max_{D(x,y)}|c(x,y)|, \]
and the rectangle \(D(x,y)\) is determined by the segments \([0,x]\) and \([0,y]\).

In the case of equation (6) and system (12), estimates of type (16) naturally have a somewhat more cumbersome form.

Taking account of the latter results, from the approximating systems and the original equations one can, without any essential difficulty, also obtain estimates for all partial derivatives of \(u(x,y)\). In particular, estimates can be obtained for \(M_n^{(2)}\), \(M_n^{(1,2)}\), which enter into inequalities (11).

2. Other problems. Without first reducing the original second-order differential equations to systems of equations of the form (1), by applying the approach indicated above one can approximately solve, in addition to the Goursat problem, a number of other problems of mathematical physics. We shall dwell briefly on those among them for which it has been possible to prove uniform convergence of the method used and to obtain a priori estimates of its error.

First of all, let us mention the Cauchy problem, which consists in finding the solution \(u=u(x,y)\) of the original second-order equation if, on an arc which at no point assumes a characteristic direction, the values of the function \(u\) and of its derivative in a direction not tangent to the curve carrying the initial values are prescribed. In the case of equation (13), for example, and initial conditions chosen, for definiteness, in the form

\[ u[x,g(x)]=\varphi(x),\qquad \frac{\partial u[x,g(x)]}{\partial y}=\psi(x), \]

\[ -\infty<G\le g'(x)\le g^*<0,\qquad g(0)=0,\qquad -l'\le x\le 0, \]

the error estimate of the method of interest to us was obtained in the following form:

\[ |\gamma_n(x)| \leq \frac{h^2 M_n^{(2)}}{12}\,[e^{C n y_n (x-x_n)}-1],\quad x_n \leq x \leq 0,\quad n=1,\,2,\ldots,N, \]

where \(x_n\) denotes the abscissa of the point of intersection of the line \(y=y_n=nh\) with the curve \(g(x)\), and for the other quantities the notation used in (14) is employed, if by \(D_n\) one understands the curvilinear triangle bounded by the curve \(g(x)\) and the characteristics \(y=y_n\) and \(x=0\).

The last result is generalized to the case of equation (6) and system (12) with the same admissible nonlinearities as in Goursat’s problem. In the linear case, one can also write down here the corresponding estimates for exact solutions and their partial derivatives. As for the explicit computational schemes obtained in Section 1 with the use of interpolation of any fixed order \(k>1\), in the case of the Cauchy problem it is not possible to construct them in the same way as was done above.

Literally the same results can also be achieved in the case of the first mixed problem, which consists in constructing the solution \(u=u(x,y)\) of the original second-order equation, if on the characteristic \(ab\) the value of \(u\) is prescribed, and on the arc \(ac\), which at no point has a characteristic direction, the values of the function \(u\) and of its derivative in a direction not tangent to the curve carrying the initial values are prescribed; here it is assumed that at the common point \(a\) the values of the corresponding functions are compatible and that the characteristic of the second family issuing from the point \(a\) lies inside the angle \(bac\).

More can be achieved along this path in the case of the Picard problem, whose boundary conditions for equation (6), for example, can be prescribed in the following form:

\[ u(x,0)=\varphi(x),\quad -l_1 \leq x \leq l_2,\quad l_1,l_2>0;\quad u[g(y),y]=\psi(y),\quad 0 \leq y \leq Y; \]

\[ \varphi(0)=\psi(0),\quad g(0)=0,\quad |g'(y)| \leq G<\infty,\quad -l_1 \leq g(y) \leq l_2. \]

Here one can write down all results analogous to the corresponding results obtained in the case of Goursat’s problem. We shall now give, as an example, only an estimate of type (16):

\[ |u(x,y)| \leq \frac{F}{C} \left[ e^{Cy\left(|x-g(y)|+\displaystyle\int_0^y |g'(y)|\,dy\right)} -1 \right], \quad -l_1 \leq x \leq l_2,\quad 0 \leq y \leq Y, \]

where for \(F\) and \(C\) the notation used in (16) is employed, with the natural change in the meaning of \(D(x,y)\). All the generalizations made in Section 1 also remain in force.

All these results are somewhat simplified when considering a particular case of the Picard problem, the so-called second mixed problem, which consists in finding the solution of the original equation if compatible values of the solution are prescribed on the characteristic \(ab\) and on the curve \(ac\), which nowhere assumes the given characteristic direction; here it is assumed that the second characteristic issuing from the point \(a\) lies outside the angle \(bac\).

Let us note in conclusion that a similar approach can also be applied to other problems, and in this case the type of the original equations may also be different. However, the convergence of the method in these cases, generally speaking, has not yet been investigated.

§ II. APPROXIMATE SOLUTION OF HYPERBOLIC SYSTEMS OF TWO FIRST-ORDER EQUATIONS BY REDUCING THEM, BY THE METHOD OF INTEGRAL RELATIONS, TO A SYSTEM OF LINEAR ALGEBRAIC EQUATIONS

This section will be devoted primarily to the approximate solution of linear hyperbolic systems of two first-order equations of the general form

\[ \left. \begin{aligned} a_1 u_x + a_2 u_y + a_3 v_x + a_4 v_y + a_5 u + a_6 v + a_7 &= 0,\\ b_1 u_x + b_2 u_y + b_3 v_x + b_4 v_y + b_5 u + b_6 v + b_7 &= 0 \end{aligned} \right\}, \tag{17} \]

where \(a_i\) and \(b_i\) \((i=1,2,\ldots,7)\) are known, and \(u\) and \(v\) are the unknown functions of two independent variables \(x\) and \(y\).

Applying the idea of integration with subsequent interpolation of the integrand functions not to the original system written in divergence form (1), but directly to the conditions on the characteristics of this system, we obtain, as the approximating system, a system of linear algebraic equations, and not a system of ordinary differential equations, as was the case in all the preceding instances.

Let us first consider for system (17) the Goursat problem, which consists in finding the solution \(u, v\) of the original system if on two of its characteristics issuing from one point the values of \(u\) and \(v\) are prescribed (it is understood, of course, that the values of the functions prescribed on the characteristics must be compatible).

The characteristics of system (17) do not depend on the solution of the problem posed and can be found in advance. We take them as coordinate lines in a new coordinate system \(\xi, \eta\). Let \(\varphi_1(\xi)\) and \(\varphi_2(\xi)\) be the prescribed values of \(u\) and \(v\) on the boundary characteristic \(\eta=0\), and \(\psi_1(\eta)\) and \(\psi_2(\eta)\), respectively, on the characteristic \(\xi=0\).

For a given step \(h>0\), construct a grid of characteristics

\[ \xi=\xi^i=ih,\qquad \eta=\eta^j=jh,\qquad i,j=1,2,\ldots . \]

For an approximate solution \(u_{mn}, v_{mn}\) of the problem posed and for the errors of the method \(\gamma_{mn}=u(\xi^m,\eta^n)-u_{mn}\), \(\delta_{mn}=v(\xi^m,\eta^n)-v_{mn}\) at some point \((\xi^m,\eta^n)\) of our grid, integrating the known conditions on the characteristics of system (17), respectively with respect to \(\eta\) from \(\eta=\eta^j\) to \(\eta=\eta^{j+1}\) and with respect to \(\xi\) from \(\xi=\xi^i\) to \(\xi=\xi^{i+1}\), and interpolating the functions \(u(\xi^i,\eta)\) and \(v(\xi^i,\eta)\), \(u(\xi,\eta^j)\) and \(v(\xi,\eta^j)\) linearly according to their values, respectively, on the curves \(\eta=\eta^j\) and \(\eta=\eta^{j+1}\), \(\xi=\xi^i\) and \(\xi=\xi^{i+1}\), we obtain systems of linear algebraic equations for \(u_{mn}, v_{mn}, \gamma_{mn}\), and \(\delta_{mn}\), from which their explicit expressions may be found. For example, for \(\gamma_{mn}\) one obtains

\[ \gamma_{mn} = \frac{1}{q_{mn}} \left[ R_{mn} + \frac{\beta_{mn}}{q_{m-1\,n}}\sigma_{m-1\,n} + \frac{\beta^*_{mn}}{q_{mn-1}} \times \right. \]

\[ \left. \times \sum_{j=1}^{n-1} \left( Q_{mj} - \frac{t_{mj}}{q_{m-1\,j}}\sigma_{m-1\,j} \right) \prod_{i=1}^{\,n-j-1} \frac{p_{mn-i}}{q_{mn-i-1}} \right], \tag{18} \]

where

\[ \sigma_{il} = \frac{t^*_{il}}{q_{il-1}} \sum_{j=1}^{l-1} \left( Q_{ij} - \frac{t_{ij}}{q_{i-1\,j}}\sigma_{i-1\,j} \right) \prod_{k=1}^{\,l-j-1} \frac{p_{il-k}}{q_{il-k-1}} - \frac{p^*_{il}}{q_{i-1\,l}}\sigma_{i-1\,l} - Q^*_{il}, \qquad \sigma_{01}=0, \]

\[ i=1,2,\ldots,m-1;\qquad l=1,2,\ldots,n . \]

(the formulas for \(\delta_{mn}\) are obtained if in (18) we replace \(\gamma_{mn}\) by \(-\delta_{mn}\), \(R_{mn}\) by \(F_{mn}\), \(\beta_{mn}\) by \(\alpha_{mn}\), \(\beta^*_{mn}\) by \(\alpha^*_{mn}\), while the formulas for \(u_{mn}, v_{mn}\) are obtained from the formulas for \(\gamma_{mn}, \delta_{mn}\), respectively, by replacing \(R_{mn}, F_{mn}, \sigma_{ij}, Q_{mj}, Q^*_{in}\) by \(\overline R_{mn}, \overline F_{mn}, \overline\sigma_{ij}, \overline Q_{mj}, \overline Q^*_{in}\)). Here the following notation has been adopted:

\[ a_{ij}=a_i b_j-a_j b_i,\qquad i,j=1,2,\ldots,7;\qquad \Phi=\frac{\partial(\xi,\eta)}{\partial(x,y)}\ne 0, \]

\[ \overline A_1=a_{13}\xi_x\eta_x+a_{14}\xi_y\eta_x+a_{23}\xi_x\eta_y+a_{24}\xi_y\eta_y,\qquad \overline B_1=a_{43}\Phi, \]

\[ \overline C_1=a_{53}\xi_x+a_{54}\xi_y,\qquad \overline D_1=a_{63}\xi_x+a_{64}\xi_y,\qquad \overline F_1=a_{73}\xi_x+a_{74}\xi_y; \]

\[ \overline A_2=a_{13}\eta_x\xi_x+a_{14}\eta_y\xi_x+a_{23}\eta_x\xi_y+a_{24}\eta_y\xi_y,\qquad \overline B_2=a_{34}\Phi, \]

\[ \overline C_2=a_{53}\eta_x+a_{54}\eta_y,\qquad \overline D_2=a_{63}\eta_x+a_{64}\eta_y,\qquad \overline F_2=a_{73}\eta_x+a_{74}\eta_y; \]

\[ A_1^i=\overline A_1(\xi^i,\eta),\qquad A_2^j=\overline A_2(\xi,\eta^j),\qquad A_1^{ij+1}=\overline A_1(\xi^i,\eta^{j+1}) \]

(similarly also for \(B_1^i, C_1^i, D_1^i, F_1^i, B_2^j, C_2^j, D_2^j, F_2^j, B_1^{ij+1}, A_2^{i+1j}, B_2^{i+1j}\));

\[ c_1=C_1^i-\frac{dA_1^i}{d\eta},\qquad d_1=D_1^i-\frac{dB_1^i}{d\eta},\qquad c_2=C_2^j-\frac{dA_2^j}{d\xi},\qquad d_2=D_2^j-\frac{dB_2^j}{d\xi}, \]

\[ \alpha_{ij+1}=A_1^{ij+1}+\frac{1}{h}\int_{\eta^j}^{\eta^{j+1}} c_1(\eta-\eta^j)\,d\eta, \]

\[ \beta_{ij+1}=B_1^{ij+1}+\frac{1}{h}\int_{\eta^j}^{\eta^{j+1}} d_1(\eta-\eta^j)\,d\eta,\qquad \xi_{ij}=A_1^{ij}-\frac{1}{h}\int_{\eta^j}^{\eta^{j+1}} c_1(\eta^{j+1}-\eta)\,d\eta, \]

\[ \eta_{ij}=B_1^{ij}-\frac{1}{h}\int_{\eta^j}^{\eta^{j+1}} d_1(\eta^{j+1}-\eta)\,d\eta,\qquad f_{ij+1}=-\int_{\eta^j}^{\eta^{j+1}} F_1^i\,d\eta, \]

\[ r_{ij+1}=\frac{1}{2}\int_{\eta^j}^{\eta^{j+1}} \left[ c_1\frac{\partial^2 u(\xi^i,\widetilde\eta)}{\partial\eta^2} +d_1\frac{\partial^2 v(\xi^i,\widetilde{\widetilde\eta})}{\partial\eta^2} \right](\eta-\eta^j)(\eta^{j+1}-\eta)\,d\eta, \]

\[ \eta^j<\widetilde\eta,\widetilde{\widetilde\eta}<\eta^{j+1};\qquad \alpha^*_{i+1j}=A_2^{i+1j}+\frac{1}{h}\int_{\xi^i}^{\xi^{i+1}} c_2(\xi-\xi^i)\,d\xi, \]

\[ \beta^*_{i+1j}=B_2^{i+1j}+\frac{1}{h}\int_{\xi^i}^{\xi^{i+1}} d_2(\xi-\xi^i)\,d\xi,\qquad \xi^*_{ij}=A_2^{ij}-\frac{1}{h}\int_{\xi^i}^{\xi^{i+1}} c_2(\xi^{i+1}-\xi)\,d\xi, \]

\[ \eta^*_{ij}=B_2^{ij}-\frac{1}{h}\int_{\xi^i}^{\xi^{i+1}} d_2(\xi^{i+1}-\xi)\,d\xi,\qquad f^*_{i+1j}=-\int_{\xi^i}^{\xi^{i+1}} F_2^j\,d\xi, \]

\[ r^*_{i+1j}=\frac{1}{2}\int_{\xi^i}^{\xi^{i+1}} \left[ c_2\frac{\partial^2 u(\widetilde\xi,\eta^j)}{\partial\xi^2} +d_2\frac{\partial^2 v(\widetilde{\widetilde\xi},\eta^j)}{\partial\xi^2} \right](\xi-\xi^i)(\xi^{i+1}-\xi)\,d\xi, \]

\[ \xi^i<\tilde{\xi}, \qquad \tilde{\xi}<\xi^{i+1}; \]

\[ t_{ij}=\alpha_{ij}\eta_{ij}-\beta_{ij}\xi_{ij}, \quad t_{ij}^{*}=\alpha_{ij}^{*}\eta_{ij}^{*}-\beta_{ij}^{*}\xi_{ij}^{*}, \quad q_{ij}=\alpha_{ij}\beta_{ij}^{*}-\alpha_{ij}^{*}\beta_{ij}, \]

\[ R_{ij}=r_{ij}\beta_{ij}^{*}-r_{ij}^{*}\beta_{ij}, \quad F_{ij}=r_{ij}\alpha_{ij}^{*}-r_{ij}^{*}\alpha_{ij}, \quad \overline{R}_{ij}=\rho_{ij}\beta_{ij}^{*}-\rho_{ij}^{*}\beta_{ij}, \]

\[ \overline{F}_{ij}=\rho_{ij}\alpha_{ij}^{*}-\rho_{ij}^{*}\alpha_{ij}, \quad \rho_{ij}= \begin{cases} f_{ij}, & j\ne 1,\\ f_{i1}+\xi_{i0}\varphi_1(\xi^i)+\eta_{i0}\varphi_2(\xi^i), & j=1, \end{cases} \]

\[ \rho_{ij}^{*}= \begin{cases} f_{ij}^{*}, & i\ne 1,\\ f_{1j}^{*}+\xi_{0j}^{*}\psi_1(\eta^j)+\eta_{0j}^{*}\psi_2(\eta^j), & i=1, \end{cases} \]

\[ p_{ij}=\beta_{ij}^{*}\xi_{ij}-\alpha_{ij}^{*}\eta_{ij}, \quad p_{ij}^{*}=\beta_{ij}\xi_{ij}^{*}-\alpha_{ij}\eta_{ij}^{*}, \quad Q_{ij}=p_{ij}r_{ij}+r_{ij}^{*}t_{ij}, \]

\[ Q_{ij}^{*}=-p_{ij}^{*}r_{ij}^{*}-r_{ij}t_{ij}^{*}, \quad \overline{Q}_{ij}=p_{ij}\rho_{ij}+\rho_{ij}^{*}t_{ij}, \quad \overline{Q}_{ij}^{*}=-p_{ij}^{*}\rho_{ij}^{*}-\rho_{ij}t_{ij}^{*}. \]

Let us note that from the definition of hyperbolicity of the system (17) it follows directly that for sufficiently small \(h>0\) (the restrictions that must be imposed on the choice of the step will be given below) all \(q_{ij}\ne 0\).

From (18) we obtain the following estimate of the error of the method, from which the convergence of the method with rate of order \(h^2\) follows directly:

\[ |\gamma_{mn}|\le h^2 Q\{(l_n-h)PQR_1\lambda_n\beta^{**}+(l_m-h)Q\sigma_n\Lambda_{mn}[\beta+(l_n-h)QT_1\lambda_n\beta^{**}]+ \]
\[ +h[R_1\beta^{**}+R_2\beta+(l_n-h)QR_2T_1\lambda_n\beta^{**}]\} \tag{19} \]

(the estimate for \(|\delta_{mn}|\) is obtained from (19) by replacing \(\beta\) by \(\alpha\), and \(\beta^{**}\) by \(\alpha^{**}\)). Here
\(Q=q^{-1}\), \(0<q=\min\limits_i \min\limits_j |q_{ij}|\), \(q>b_{12}-H_0\) for a step \(h\) satisfying the inequality \(b_{12}-H_{12}>0\),

\[ b_{12}=\min_i \min_j |B_{12}^{ij}|>0, \qquad B_{12}^{ij}=A_1^{ij}B_2^{ij}-A_2^{ij}B_1^{ij}, \]

\[ H_0=\frac{1}{2}h\{A_1(D_2+B_2')+B_1(C_2+A_2')+A_2(D_1+B_1')+B_2(C_1+A_1')+ \]
\[ +\frac{1}{2}h[(C_1+A_1')(D_2+B_2')+(C_2+A_2')(D_1+B_1')]\}, \]

\[ H_{12}=\max_i \max_j |H_{12}^{ij}|; \qquad H_{12}^{ij}=\frac{A_1^{ij}}{h}\int_{\xi^{i-1}}^{\xi^i} d_2(\xi-\xi^{i-1})\,d\xi- \]

\[ -\frac{B_1^{ij}}{h}\int_{\xi^{i-1}}^{\xi^i} c_2(\xi-\xi^{i-1})\,d\xi -\frac{A_2^{ij}}{h}\int_{\eta^{j-1}}^{\eta^j} d_1(\eta-\eta^{j-1})\,d\eta+ \]

\[ +\frac{B_2^{ij}}{h}\int_{\eta^{j-1}}^{\eta^j} c_1(\eta-\eta^{j-1})\,d\eta +\frac{1}{h^2}\int_{\eta^{j-1}}^{\eta^j} c_1(\eta-\eta^{j-1})\,d\eta \times \]

\[ \times \int_{\xi^{i-1}}^{\xi^i} d_2(\xi-\xi^{i-1})\,d\xi -\frac{1}{h^2}\int_{\eta^{j-1}}^{\eta^j} d_1(\eta-\eta^{j-1})\,d\eta \int_{\xi^{i-1}}^{\xi^i} c_2(\xi-\xi^{i-1})\,d\xi . \]

\[ A_1=\max_i\max_\eta |A_1^i|,\quad A_2=\max_j\max_\xi |A_2^j|,\quad A_1'=\max_i\max_\eta \left|\frac{dA_1^i}{d\eta}\right| \]

(similarly for \(B_1, C_1, D_1, F_1, B_2, C_2, D_2, F_2, B_1', A_2', B_2'\)), \(P=B_{12}+H_0\),

\[ B_{12}=\max_i\max_j |B_{12}^{ij}|,\quad B_{12}'=\max_i\max_j \left\{ \frac{\partial B_{12}(\xi^i,\eta^j+\theta_1^j h)}{\partial\eta}, \frac{\partial B_{12}(\xi^i+\theta_2^i h,\eta^j)}{\partial\xi} \right\}, \]

\[ 0<\theta_1^j,\theta_2^i<1,\quad \lambda_n\le \lambda=\exp\frac{l_n(\widetilde B_{12}'+2\widetilde H)}{\widetilde b_{12}} \quad \left( \text{here } l_n=nh,\quad H=\frac{1}{h}H_0, \right. \]

and, by a tilde above, as everywhere below, the limiting values of the corresponding quantities as \(h\to0\) are denoted),

\[ \left. l_m=mh,\quad \Lambda_{mn}\le \Lambda, \right. \]

\[ \Lambda=\exp\frac{l_m}{\widetilde b_{12}} \left[ (\widetilde B_{12}'+2\widetilde H) +\frac{l_n\widetilde T_1\widetilde T_2}{\widetilde b_{12}} \exp\frac{l_n}{\widetilde b_{12}}(\widetilde B_{12}'+2\widetilde H) \right], \]

\[ \widetilde T_i=\widetilde A_i(\widetilde D_i+\widetilde B_i') +\widetilde B_i(\widetilde A_i'+\widetilde C_i),\quad i=1,2, \]

\[ T_i\le \widetilde T_i+\frac12 h(C_i+A_i')(D_i+B_i'),\quad i=1,2, \]

\[ \beta\le B_1+\frac12 h(D_1+B_1'),\quad \beta^*\le B_2+\frac12 h(D_2+B_2'), \]

\[ \alpha\le A_1+\frac12 h(C_1+A_1'), \]

\[ \alpha^*\le A_2+\frac12 h(C_2+A_2'),\quad R_1\le \frac1{12}\,[M(C_1+A_1')+N(D_1+B_1')], \]

\[ R_2\le \frac1{12}\,[\overline M(C_2+A_2')+\overline N(D_2+B_2')], \]

\[ M=\max_i\max_\eta \left|\frac{\partial^2u(\xi^i,\eta)}{\partial\eta^2}\right|,\quad N=\max_i\max_\eta \left|\frac{\partial^2v(\xi^i,\eta)}{\partial\eta^2}\right|, \]

\[ \overline M=\max_j\max_\xi \left|\frac{\partial^2u(\xi,\eta^j)}{\partial\xi^2}\right|,\quad \overline N=\max_j\max_\xi \left|\frac{\partial^2v(\xi,\eta^j)}{\partial\xi^2}\right|, \]

\[ \sigma_n=(l_n-h)T_2Q_1Q\lambda_n+Q_2,\quad Q_1\le PR_1+hT_1R_2,\quad Q_2\le PR_2+hT_2R_1. \]

Let us note that the posed problem is solved analogously also in the case when system (17) is quasilinear.

Moreover, the method applied here also makes it possible to obtain, comparatively simply, an upper estimate for the modulus of the exact solution of the original Goursat problem at the given point \((\xi,\eta)\). We give this estimate, assuming, for simplicity, that the boundary values \(u\) and \(v\) are zero:

\[ |u(\xi,\eta)|\le \frac12\,\widetilde B_{12}\widetilde b_{12}^{-2} \left[ \widetilde B_1\widetilde F_2(\Lambda+\lambda_{(1)})\xi +\widetilde B_2\widetilde F_1(\Lambda_{(1)}+\lambda)\eta +\widetilde b_{12}^{-1}\times \]

\[ \times(\widetilde B_1\widetilde F_1\widetilde T_2+\widetilde B_2\widetilde F_2\widetilde T_1) (\Lambda\lambda+\Lambda_{(1)}\lambda_{(1)})\xi\eta +\widetilde b_{12}^{-2}\widetilde T_1\widetilde T_2\times \]

\[ \times(\widetilde B_1\widetilde F_2\Lambda_{(1)}\lambda_{(1)}^2\xi +\widetilde B_2\widetilde F_1\Lambda\lambda^2\eta)\xi\eta \right], \tag{20} \]

where the values for \(\lambda_{(1)}\) and \(\Lambda_{(1)}\) are obtained, respectively, from the values for \(\lambda\) and \(\Lambda\) by replacing \(\xi\) by \(\eta\), and conversely (the estimate for \(|v(\xi,\eta)|\) is obtained from (20) by replacing \(\widetilde{B}_i\) by \(\widetilde{A}_i\), \(i=1,2\)).

The estimate obtained is attainable. Thus, for example, in the case of the system \(A_1^0 u_\eta-B_1^0 v_\eta-F_1^0=0,\quad A_2^0 u_\xi-B_2^0 v_\xi+F_2^0=0\), where \(A_i^0,\ B_i^0,\ F_i^0\) \((i=1,2)\) are nonnegative constants, and \(A_1^0B_2^0-A_2^0B_1^0>0\), the result given by estimate (20) coincides completely with the exact solution of the corresponding Goursat problem.

In conclusion we note that in a similar manner one can proceed in the case of the Cauchy problem and of both mixed problems for the original system of equations.

A more detailed account of the principal results of both sections may be found in works [24]—[28].

References

1. Dorodnitsyn A. A. Proceedings of the Third All-Union Mathematical Congress, 1956, vol. III, Publishing House of the Academy of Sciences of the USSR, Moscow, 1958, pp. 447—453.

2. Dorodnitsyn A. A. Applied Mechanics and Technical Physics, 1, No. 3, 1960, pp. 111—118.

3. Belotserkovskii O. M., Chushkin P. I. Journal of Computational Mathematics and Mathematical Physics, 2, No. 5, 1962, pp. 731—759.

4. Belotserkovskii O. M. Journal of Computational Mathematics and Mathematical Physics, 2, No. 6, 1962, pp. 1062—1085.

5. Belotserkovskii O. M., Dushin V. K. Journal of Computational Mathematics and Mathematical Physics, 4, No. 1, 1964, pp. 61—77.

6. Belotserkovskii O. M., Golomazov M. M., Shulishnina N. P. Journal of Computational Mathematics and Mathematical Physics, 4, No. 2, 1964, pp. 306—316.

7. Pavlovskii Yu. N. Journal of Computational Mathematics and Mathematical Physics, 4, No. 1, 1964, pp. 178—183.

8. Alikhashkin Ya. I., Favorskii A. P., Chushkin P. I. Journal of Computational Mathematics and Mathematical Physics, 3, No. 6, 1963, pp. 1130—1134.

9. Minailos A. N. Journal of Computational Mathematics and Mathematical Physics, 4, No. 1, 1964, pp. 171—177.

10. Mikheev A. S. Bulletin of Leningrad State University, series Mathematics, Mechanics, Astronomy, issue 1, No. 1, 1964, pp. 108—113.

11. Kadyshnikov V. M. Izvestiya of the Academy of Sciences of the USSR, Geophysics Series, No. 8, 1962, pp. 1083—1092.

12. Tairov M. A. Journal of Computational Mathematics and Mathematical Physics, 2, No. 5, 1962, pp. 938—942.

13. Dorodnitsyn A. A. Arch. Mech. Stosowanej, 12, No. 1, 1960, pp. 13—28.

14. Dorodnitsyn A. A. All-Union Conference “Paths for the Development of Soviet Mathematical Machine-Building and Instrument-Making,” 1956, part 1, VINITI, Moscow, 1956, pp. 44—52.

15. Dorodnicyn A. A. Advances Aeronaut. Sci., 2. London—New-York—Paris—Los Angeles, Pergamon Press, 1959, pp. 832—844.

16. Dorodnitsyn A. A. Proceedings of the Conference on Computational Technology. Publishing House of the Academy of Sciences of the USSR, 1958.

17. Dorodnicyn A. A. Proceedings of 9th Intern. Congress of Applied Mechanics, Brussels, 1, 1957.

18. Dorodnitsyn A. A. Acta mech. sinica, 3, No. 4, 1959, pp. 352—355.

19. Dorodnitsyn A. A. In the collection “Some Problems of Mathematics and Mechanics” (on the 60th anniversary of Academician M. A. Lavrent’ev). Publishing House of the Siberian Branch of the Academy of Sciences of the USSR, Novosibirsk, 1961, pp. 77—83.

20. Babenko K. I. Proceedings of the All-Union Congress on Theoretical and Applied Mechanics, 1960. Publishing House of the Academy of Sciences of the USSR, 1962, pp. 72—76.

21. Hayes U. D., Probstein R. F. Theory of Hypersonic Flows. Foreign Literature Publishing House, Moscow, 1962.

22. Lozinskii S. M. Izvestiya of Higher Educational Institutions, Mathematics Series, No. 5, 1958, pp. 52—90.

23. Budak B. M. Reports of the Academy of Sciences of the USSR, 109, No. 1, 1956, pp. 9—12.

24. Krylov V. I., Bobkov V. V. Reports of the Academy of Sciences of the Byelorussian SSR, 7, No. 7, 1963, pp. 433—438.

25. Bobkov V. V. Izvestiya of the Academy of Sciences of the Byelorussian SSR, Physics-Technical Series, No. 4, 1963, pp. 14—24.

26. Bobkov V. V. Reports of the Academy of Sciences of the Byelorussian SSR, 8, No. 1, 1964, pp. 5—9.

27. Bobkov V. V. Izvestiya of the Academy of Sciences of the Byelorussian SSR, Physics-Technical Series, No. 1, 1964, pp. 18—22.

28. Bobkov V. V. Dissertation. Minsk, 1964.

29. Lun’kin Yu. P., Popov F. D. Journal of Computational Mathematics and Mathematical Physics, 4, No. 5, 1964, pp. 896—904.

30. Dorodnicyn A. A. Inform. Process., 1962. Amsterdam, N. Holland Publ. Co., 1963, pp. 122—125.

Received by the editors
November 20, 1964

Byelorussian State University
named after V. I. Lenin
Institute of Mathematics, Academy of Sciences of the BSSR