Least Squares Solution of a Boundary Value Problem for the Two-Dimensional Poisson Equation. Extension of the Method to Time-Independent Equations

B. F. Krutoi

(Presented by the scientific seminar of the Department of Engineering and Computational Mathematics)

The Poisson equation is given for a plane closed, generally multiply connected domain \(\overline{D}^{(2)} = D^{(2)} + r\):

\[ \Delta u(x,y)=\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=f(x,y), \tag{1} \]

with the boundary condition

\[ \overline{u}=u(x,y)=h(x,y) \tag{2} \]

on the boundary \(R\) of the domain \(\overline{D}^{(2)}\) (Fig. 1).

We replace the unknown function \(u(x,y)\) by its approximate representation \(\varphi(x,y)\) in the form of the expansion

\[ u(x,y)\approx \varphi(x,y)=\sum_{x=0}^{m} c_x \varphi_x(x,y) \tag{3} \]

in terms of certain supporting linearly independent functions \(\varphi_x(x,y)\) chosen by us, where \(c_x\) are unknown coefficients.

To determine the coefficients \(c_x\), we specify \(\mu\) arbitrary points \(\overline{M}_s\) on the boundary \(R\) and \(\nu\) arbitrary points \(M_t\) inside the domain \(\overline{D}^{2}\) (Fig. 1), requiring that

\[ \mu+\nu=n>m+1=\overline{m}, \tag{4} \]

Then for the boundary \(r\) we can write \(\mu\) error equations of the form following from (2), (3):

\[ \sum_{x=0}^{m} c_x \varphi_x(x_s,y_s)-h(x_s,y_s)=\varepsilon_s, \tag{5} \]

\[ s=1,2,\ldots,\mu, \]

where \(\varepsilon_s\) are the corresponding errors of these equations. For the interior domain \(D^{(2)}\), we can write \(\nu\) error equations of the form following from (1), (3):

\[ \sum_{\chi=0}^{m}\left(\frac{\partial^{2}\varphi_{\chi}(x_t,y_t)}{\partial x^{2}}+ \frac{\partial^{2}\varphi_{\chi}(x_t,y_t)}{\partial y^{2}}\right)-f(x_t,y_t)=\varepsilon_t, \tag{6} \]

\[ (t=1,2,\ldots,\nu) \]

where \(\varepsilon_t\) are the errors of equations (6).

Introduce the notation

\[ \begin{aligned} &1)\quad \varphi_{\chi}(x_s,y_s)=a_{s\chi}, \qquad &&2)\quad h(x_s,y_s)=h_s,\\ &3)\quad \frac{\partial^{2}\varphi_{\chi}(x_t,y_t)}{\partial x^{2}}+ \frac{\partial^{2}\varphi_{\chi}(x_t,y_t)}{\partial y^{2}}=b_{t\chi}, \qquad &&4)\quad f(x_t,y_t)=f_t. \end{aligned} \tag{7} \]

Fig. 1. The domain \(\overline{D}^{(2)}, D^{(2)}+\Gamma\), and the interior \(M^t\) and boundary \(M_s\) support points for the function \(u(x,y)=u(M)\).

In the notation (7), the combined system of \(\mu+\nu=n\) error equations (5), (6) is written as

\[ \left\{ \begin{aligned} &\sum_{\chi=0}^{m} a_{s\chi}c_{\chi}-h_s=\varepsilon_s, \quad (s=1,2,\ldots,\mu),\\ &\sum_{\chi=0}^{m} b_{t\chi}c_{\chi}-f_t=\varepsilon_t, \quad (t=1,2,\ldots,\nu) \end{aligned} \right. \qquad \left(\mu+\nu=n<m+1=\overline{m}\right). \tag{8} \]

or, in matrix form,

\[ \begin{cases} a_{(\mu \overline{m})}\,c_{(\overline{m}1)}-h_{(\mu 1)}=\varepsilon_{(\mu 1)},\\[4pt] b_{(\nu \overline{m})}\,c_{(\overline{m}1)}-f_{(\nu 1)}=\varepsilon_{(\nu 1)} \end{cases} \qquad (\mu+\nu=n>m+1=\overline{m}) \tag{8*} \]

For the purpose of further simplifying the notation of the combined system (8), (8) we introduce in (8) the following additional notation:

\[ \begin{aligned} &1)\quad \begin{pmatrix} a_{\mu \overline{m}}\\ b_{\nu \overline{m}} \end{pmatrix} =A_{(\mu+\nu,\overline{m})}=A_{(n\overline{m})}, \qquad &2)\quad \begin{pmatrix} h_{\mu 1}\\ f_{\nu 1} \end{pmatrix} =F_{(\mu+\nu,1)}=F_{(n1)} \tag{9}\\[6pt] &3)\quad \begin{pmatrix} \varepsilon_{\mu 1}\\ \varepsilon_{\nu 1} \end{pmatrix} =E_{(\mu+\nu,1)}=E_{(n1)}, \qquad &4)\quad c_{(\overline{m}1)}=C_{(\overline{m}1)}. \end{aligned} \]

Then the combined system (8), (8*) of \(n\) error equations can be written in one of the following two forms:

\[ \sum_{\chi=0}^{m} A_{s\chi}C_{\chi}-F_s=E_s, \quad (s=1,2,\ldots,n>m+1=\overline{m}). \tag{10} \]

\[ A_{(n\overline{m})}\,C_{(\overline{m}1)}-F_{(n1)}=E_{(n1)}, \quad (n>m+1=\overline{m}). \tag{10*} \]

Since the representation (3) for the desired function \(u(x,y)\) is approximate, the \(\overline{m}\) support functions \(\varphi_\chi(x,y)\) are linearly dependent, and \(n>m\), there exists an infinite set of collections \(c_1^{(\nu)}, c_2^{(\nu)},\ldots,c_m^{(\nu)}\) of the \(\overline{m}\) coefficients \(c_\chi=C_\chi\) satisfying the system of equations (10) with one or another degree of accuracy. Following Gauss, we shall regard as the most reliable that collection \(\overline{c}_0,c_1,\ldots,c_m\) of the coefficients \(c_\chi=C_\chi\) which satisfies the condition of the method of least squares

\[ \sum_{s=1}^{n} E_s^2 = \sum_{s=1}^{n} E_{s1}^2 = \sum_{s=1}^{n} E_{1s}^{*}E_{s1} = E_{(1n)}^{*}E_{(n1)} = Q^2 =\text{min}. \tag{11} \]

But, according to (10),

\[ Q^2=\sum_{s=1}^{n}E_{1s}^{*}E_{s1} = Q^2(C_0,C_1,\ldots,C_m). \]

Therefore condition (11) is equivalent to the \(\overline{m}\) conditions

\[ \frac{\partial}{\partial C_\lambda}\sum_{s=1}^{n}E_s^2 = \frac{\partial}{\partial C_\lambda} Q^2(C_0,C_1,\ldots,C_\chi,\ldots,C_m) =0. \tag{12} \]

\[ (\lambda=0,1,\ldots,m) \]

Carrying out the calculation of conditions (12), we successively find

\[ \begin{aligned} \frac{\partial}{\partial C_\lambda}\sum_{s=1}^{n}E_s^2 &= 2\sum_{s=1}^{n}E_s\frac{\partial E_s}{\partial C_\lambda} = 2\sum_{s=1}^{n}E_s \frac{\partial}{\partial C_\lambda} \left(\sum_{\chi=0}^{m}A_{s\chi}C_\chi-F_s\right) \\ &= 2\sum_{s=1}^{n}E_s\sum_{\chi=0}^{m}A_{s\chi} \frac{\partial C_\chi}{\partial C_\lambda} = 2\sum_{s=1}^{n}E_s\sum_{\chi=0}^{m}A_{s\chi}\delta_\lambda^\chi =0, \end{aligned} \]

where

\[ \delta_x^\lambda = \begin{cases} +1, & x=\lambda,\\ 0, & x\ne\lambda . \end{cases} \]

From this it is seen that

\[ \sum_{x=0}^{m} A_{sx}\delta_x^\lambda = A_{s\lambda}=A_{\lambda s}^{*}. \]

Therefore

\[ \begin{aligned} \frac{1}{2}\frac{\partial}{\partial C_\lambda}\sum_{s=1}^{n} \mathrm{E}_s^2 &=\sum_{s=1}^{n}\left(\sum_{x=0}^{m} A_{sx}C_x-F_s\right)A_{\lambda s}^{*} \\ &=\sum_{x=0}^{m}\left(\sum_{s=1}^{n} A_{\lambda s}^{*}A_{sx}\right)C_x-\sum_{s=1}^{n}A_{\lambda s}^{*}F_s=0, \end{aligned} \tag{*} \]

\[ (\lambda=0,1,\ldots,m) \]

Introduce the notation

\[ \text{1) }\sum_{s=1}^{n} A_{\lambda s}^{*}A_{sx}=B_{\lambda x}, \qquad \text{2) }\sum_{s=1}^{n} A_{\lambda s}^{*}F_s=K_\lambda . \tag{13} \]

or, in matrix form,

\[ \text{1) } A_{(\bar m n)}^{*}A_{(n\bar m)}=B_{(\bar m\bar m)}, \qquad \text{2) } A_{(\bar m n)}^{*}F_{(n1)}=K_{(\bar m1)} . \tag{13*} \]

Then the \(\bar m+1=\bar m\) conditions (12) assume the following final form:

\[ \sum_{x=0}^{m} B_{\lambda x}C_x-K_\lambda=0,\qquad (\lambda=0,1,\ldots,m). \tag{14} \]

or, in matrix form,

\[ B_{(\bar m\bar m)}C_{(\bar m1)}-K_{(\bar m1)}=0_{(\bar m1)} . \tag{14*} \]

Solving the set of \(\bar m\) equations (14), we determine the most advantageous values \(\bar C_x=c_x\) for the \(m\) coefficients \(c_x\) entering into the representation (3) for the sought function \(u(x,y)\). Substitution of the found values of these coefficients into equations (10), or—what is the same—into equations (8), gives the corresponding errors \(\varepsilon_s=e_s,\ \varepsilon_t\) of the indicated equations.

We shall verify the correctness of the calculation of the coefficients \(C_x=c_x\), the errors \(\varepsilon_s=e_s,\ \varepsilon\), and the quantity \(Q^2=\sum_{s=1}^{n} e_s^2\) by multiplying equation (10) on the left by \(A_{(\bar m n)}^{*}\), \(e_{(1n)}^{*}\), and \(F_{(1n)}^{*}\). Then we obtain the following verification equalities, taking into account (13) and (14*):

\[ \text{1) } A_{(\bar m n)}^{*}\mathrm{E}_{(n1)}=0_{(\bar m1)},\qquad \text{2) } Q^2=\mathrm{E}_{(1n)}^{*}\mathrm{E}_{(n1)}=-F_{(1n)}^{*}F_{(n1)}, \tag{15*} \]

\[ \text{3) } F_{(1m)}^{*}F_{(n1)}C_{(1m)}^{*}B_{(\bar m\bar m)}C_{(\bar m1)} =\mathrm{E}_{(1n)}^{*}\mathrm{E}_{(n1)}=Q^2 . \]

Calculation of the approximation \(\varphi(x,y)\) for the sought function \(u(x,y)\) by the least-squares method described above also makes it possible to compute the accuracy of this approximation, specified in the form of the root-mean-square error \(\pm m_\varphi\), or of its square \(m_\varphi^2\). The calculation of \(m_\varphi^2\) is performed as follows [1]:

\[ \begin{gathered} 1)\ \mu^2=\frac{Q^2}{\,n-m\,};\qquad 2)\ m_{c_x}^2=\mu^2\beta_{xx},\ \left(\beta_{\overline{mm}}=B_{\overline{mm}}^{-1}\right); \tag{16}\\[6pt] 3)\ m_\varphi^2=m^2[\varphi(x,y)]= \mu^2\sum_{x,\lambda=1}^{m}\beta_{x\lambda}\varphi_x(x,y)\varphi_\lambda(x,y)=\\[4pt] =\mu^2\varphi_{(1m)}^*(x,y)\beta_{\overline{mm}}\varphi_{(m1)}(x,y). \end{gathered} \]

In conclusion, we note that, in order to construct a set of reference linearly independent functions \(\varphi^x(x,y)\), one may take two series

\[ x^0=1,\quad x^1,x^2,\ldots,x^\mu,\ldots,x^r. \]

\[ y^0=1,\quad y^1,y^1,\ldots,y^\nu,\ldots,y^r. \]

and then form the matrix of products of these types

\[ \begin{gathered} x^0y^0\quad x^1y^0\ldots x^\mu y^0\ldots x^r y^0,\\ x^0y^1\quad x^1y^1\ldots x^\mu y^1\ldots x^r y^1,\\ \ldots\ldots\ldots\ldots\\ x^0y^\nu\quad x^1y^\nu\ldots x^\mu y^\nu\ldots x^r y^\nu,\\ \ldots\ldots\ldots\ldots\\ x^0y^r\quad x^1y^r\ldots x^\mu y^r\ldots x^r y^r. \end{gathered} \]

Then

\[ \varphi_{\mu\nu}(x,y)=x^\mu y^\nu=\varphi_\lambda(x,y). \]

Consequently, in this case \(m=r^2\).

To construct a set of reference functions \(\varphi_x(x,y)\), one may also use the series

\[ T_0^*(x),\ T_1^*(x),\ldots,T_\mu^*(x),\ldots,T_r^*(x). \]

\[ T_0^*(y),\ T_1^*(y),\ldots,T_\nu^*(y),\ldots,T_r^*(y), \]

where \(T_\mu^*(x)\), \(T_\nu^*(y)\) are reduced Chebyshev polynomials of the first kind, defined as follows [2]:

\[ T_0^*(z)=1,\quad T_1^*(z)=z,\quad T_x^*(z)=\frac{1}{2^{x-1}}\cos(x\arccos z); \]

\[ \arccos z=\theta,\quad z=\cos\theta,\quad T_x^*(\theta)=\frac{1}{2^{x-1}}\cos x\theta; \tag{17} \]

\[ T_{x+1}^*(z)=zT_x^*(z)-\frac14 T_{x-1}^*(z). \]

Then

\[ \varphi_{\mu\nu}(x,y)=T_\mu^*(x)T_\nu^*(y)=\varphi_\lambda(x,y). \]

In order that the numbers \(\varphi_{\mu\nu}(x,y)\) not take large values in the domain \(\overline{D}^{\,2}\), it is useful to enclose this domain in the square

\[ -1\le x\le +1 \]

\[ -1\le y\le +1, \]

placing the center of the square \(O\) at the midpoint of the domain \(\overline{D}^{(2)}\) (Fig. 1).

Let us also note that if the number \(\mu\) of points \(M_s\) of the boundary \(R\) is greater than the number \(m+1\) of coefficients \(c_x\) in the representation (3) for \(u(x,y)\) in terms of \(\varphi(x,y)\),

\[ \mu > m+1, \]

then, in order to determine the coefficients \(c_x\) and compute the quadratic residuals \(m_{c_x}^2\), \(m_\varphi^2\) for the found \(c_x\) and \(\varphi(x,y) \approx u(x,y)\), we could make do only with the points \(M_s\) taken on the boundary \(R\).

Generalization. The proposed method for solving the boundary-value problem for the Poisson equation can be applied just as simply to the Helmholtz equation

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

with the boundary condition

\[ \left[\alpha(x,y)+\beta(x,y)\frac{\partial}{\partial n}\right]u(x,y)\big|_r=h(x,y). \tag{18} \]

Moreover, the method set forth here makes it possible to solve, quite simply, even the more difficult boundary-value problem for a linear second-order differential equation with variable coefficients

\[ \left[A^{ij}(x,y)\frac{\partial^2}{\partial x^i\partial y^j}+B^i(x,y)\frac{\partial}{\partial x^i}+C\right]u(x,y)=f(x,y). \tag{19} \]

and with a boundary condition of the form (18).

Finally, it is entirely clear that the proposed method solves, without any complications, the plane boundary-value problem for a linear differential equation of order \(n\) in partial derivatives and with variable coefficients, provided that the sought function \(u\) does not depend on time \(t\), that is, \(u=u(x,y)\).

REFERENCES

1. V. Iordan. A Guide to Geodesy. Moscow, 1939.