The Numerical Solution of the External Dirichlet Generalized Harmonic Problem for a Sphere by the Method of Probabilistic Solution

Abstract

In the present paper, an algorithm for the numerical solution of the external Dirichlet generalized harmonic problem for a sphere by the method of probabilistic solution (MPS) is given, where “generalized” indicates that a boundary function has a finite number of first kind discontinuity curves. The algorithm consists of the following main stages: (1) the transition from an infinite domain to a finite domain by an inversion; (2) the consideration of a new Dirichlet generalized harmonic problem on the basis of Kelvin’s theorem for the obtained finite domain; (3) the numerical solution of the new problem for the finite domain by the MPS, which in turn is based on a computer simulation of the Weiner process; (4) finding the probabilistic solution of the posed generalized problem at any fixed points of the infinite domain by the solution of the new problem. For illustration, numerical examples are considered and results are presented.

1. Introduction

It is well known that there are many problems from the various fields of physics (for example, thermostatics, electrostatics, hydrostatics, etc.) that could be considered as external Dirichlet generalized harmonic problems.

First of all, we should note that the requirement of continuity for the boundary function in Dirichlet’s classical harmonic problem is a very strong constraint, because in practical stationary problems (connected with electric, thermal and other static fields), there are cases when it is necessary to discuss and study 2D or 3D Dirichlet generalized harmonic problems. The problems of this type appeared in the literature mainly from the 1940s (see [1,2,3,4,5]).

It is known (see [1,6]) that the methods used to obtain an approximate solution of the classical boundary-value problems are (a) less suitable or (b) not suitable at all for solving generalized boundary value problems. In the first case, the convergence of corresponding approximate process is very slow in neighborhood of boundary singularities and, consequently, the accuracy of approximate solution of the generalized problem is very low (see [1,2,3,4,5]). In the second case, the process is unstable. For example, a similar phenomenon takes place when solving the 3D Dirichlet generalized harmonic problem by the method of fundamental solutions.

In the above-mentioned literature, simplified, or so-called “solvable” generalized problems (the problems “whose” solutions can be constructed by series with terms, represented by special functions), are considered. The methods of the separation of variables, particular solutions, and heuristic methods are mainly applied for their solution, therefore the accuracy of the solution is rather low. Since heuristic methods do not guarantee that the best solution will be found (moreover, in some cases they may give an incorrect solution), it is necessary to check these solutions in order to establish how well they satisfy all the conditions of a problem (see [1]). In addition, it should be noted that, in the literature (see [4] (pp. 346–348)), while solving the Dirichlet external generalized harmonic simplest problem for a sphere, the existence of a discontinuity curve is ignored.

The study of Dirichlet’s generalized 2D and 3D problems in terms of the existence and uniqueness of solution, and selection of a reliable and effective method for its numerical solution has been intensively carried out since the 21st century in the Department of Computational Methods of the Niko Muskhelishvili Institute of Computational Mathematics.

The choice and construction of computational schemes (algorithms) mainly depend on the problem class, dimension, geometry and location of singularities on the boundary. In particular, for solving the Dirichlet generalized plane harmonic problems the following approaches may be used: (I) a method of reduction of Dirichlet generalized harmonic problems to a classical problem (see [7,8]); (II) a method of conformal mapping (see [9]); (III) a method of probabilistic solution (see [10,11]). It is evident that in the case of 3D Dirichlet harmonic problems, from the above mentioned approaches we can apply only third one.

In order to imagine the difficulties associated with the numerical solution of not only generalized, but also the classical external 3D harmonic Dirichlet problem, we will give below a brief overview of several works related to these topics.

In [12], the boundary conditions for the numerical solution of elliptic equations in exterior regions are considered. Such equations in exterior regions frequently require a boundary condition at infinity to ensure the well-posedness of the problem. Practical applications include examples of the Helmholtz and Laplace’s equations. The constructed algorithm requires the replacement of the condition at infinity by a boundary condition on a finite artificial surface. Direct imposition of the condition at infinity along the finite boundary results in large errors. A sequence of boundary conditions is developed which provides increased accuracy of approximations to the problem in the infinite domain. Estimates of the error are obtained for several cases.

In [13], the Laplace–Dirichlet problem is investigated in a three-dimensional case. The Laplace–Dirichlet problem is investigated similarly to the Laplace–Neumann problem, except for the boundary condition. Corresponding variational formulation is considered. It is based on the introduction of an auxiliary unknown by the means of the decomposition of the function u.

In the paper [14], the analytic–numerical method for solving 3D exterior problems for elliptic equations under Dirichlet and Neumann boundary conditions in half-space is considered. Based on analytical transformation, the external boundary problem is reduced to the internal one; then, the corresponding difference problem is considered based on the grid methods [6,7,8]. By using the considered analytic–numerical method for solving external boundary problems in a three-dimensional semispace, the initial problem is reduced to such a problem, the solution of which is possible by traditional techniques and methods of numerical analysis.

In [15], the solution of direct and inverse exterior boundary value problems via the strongly conditioned stochastic method is considered, mainly for exterior harmonic problems. Note that, in numerical examples from the Section 4.2 of [15], experiments were performed with cones of different thickness and observation points in the near and far fields. However, actually, the observation points given in Table 1 of [15] are “not too far” from the field, when we consider the exterior boundary value problems.

It is known that elliptic problems in external domains arise in many branches of physics. For example, the Laplace equation $(\Delta u(x)=0,x\in D)$ arises in the studies of thermostatic and electrostatic fields external to given surfaces (see [3,16]); the flow of an incompressible irrotational fluid around a body is described by the same equation (see [17]), and so on.

In these cases, infinity can be regarded as a separate boundary. A condition at infinity is required to make the external problem well posed. For the Laplace equation, it is sufficient to impose a condition of regularity at infinity. In the three dimensional (3D) case, the condition is $limu\left(x\right)=0$, for $\left|x\right|\to \infty$, where $\left|x\right|$ is the distance from a fixed (but chosen arbitrary) origin (see [18,19]).

2. Statement of the Problem

Let $D=R^3\backslash K(O,R)$ be an infinite domain in the Euclidean space $R^3$, where $K(O,R)\equiv K_R$ is a kernel with a spherical surface $S(O,R)\equiv S_R$ (with the center at the point O and the radius R, respectively). Since the harmonicity of a function is invariable under the linear transformation of the Cartesian coordinate system, therefore without loss of generality we assume that the origin of coordinates is at the point $O=(0,0,0)$.

Let us formulate the following problem.

Problem 1

(A). The function $g\left(y\right)$ is given on the boundary $S_R$ of the infinite region D and is continuous everywhere, except a finite number of curves $l_1,l_2,...,l_n$, which represent discontinuity curves of the first kind for the function $g\left(y\right)$. It is required to find a function $u\left(x\right)\equiv u(x_1,x_2,x_3)\in C^2\left(D\right)\bigcap C(\overline{D}\backslash {\displaystyle \bigcup _{k=1}^n}{l}_k)$ satisfying the conditions:

(1)$\begin{array}{cc}\hfill \Delta u\left(x\right)=0,& x\in D,\hfill \end{array}$

(2)$\begin{array}{cc}\hfill u\left(y\right)=g\left(y\right),& y\in S_R,y\overline{\in }{l}_k\subset S_R(k=1,\dots ,n),\hfill \end{array}$

(3)$\begin{array}{cc}\hfill \left|u\right(x\left)\right|<c,& x\in \overline{D},\hfill \end{array}$

(4)$\begin{array}{cc}\hfill limu\left(x\right)=0& \mathrm{for}x\to \infty ,\hfill \end{array}$

where $\Delta ={\displaystyle \sum _{i=1}^3}\frac{{\partial }^2}{\partial x_i^2}$ is the Laplace operator, c is a real constant.

On the basis of Equation (3), in general, the values of $u\left(y\right)$ are not uniquely defined on the curves $l_k$ ($k=1,\dots ,n$). In particular, if Problem 1 concerns the determination of the thermal (or the electric) field, then we must take $u\left(y\right)=0$ when $y\in l_k$, respectively. In this case, in the physical sense the curves $l_k$ are non-conductors (or dielectrics). Otherwise, $l_k$ will not be a discontinuity curve.

It is evident that surface $S_R$ is divided into the parts $S_R^i(i=1,\dots ,m)$ by the curves $l_k$ ($k=1,\dots ,n$), where one of the following conditions holds: $n=m$, $n<m$, $n>m$. Thus, the boundary function $g\left(y\right)$ has the following form

(5)$g\left(y\right)=\left\{\begin{array}{cc}{g}_1\left(y\right),\hfill & y\in S_R^1,\hfill \\ g_2\left(y\right),\hfill & y\in S_R^2,\hfill \\ \dots \dots \hfill & \dots \dots \hfill \\ g_m\left(y\right),\hfill & y\in S_R^m,\hfill \\ 0,\hfill & y\in l_k(k=1,\dots ,n),\hfill \end{array}\right.$

where $S_R^i$ ($i=1,\dots ,m$) are the parts of $S_R$ without discontinuity curves, respectively; the functions $g_i\left(y\right)$, $y\in S_R^i$ ($i=1,\dots ,m)$ are continuous on the parts $S_R^i$. It is evident that $S_R=\left({\bigcup }_{j=1}^m{S}_R^j\right)\bigcup \left({\bigcup }_{k=1}^n{l}_k\right).$

Remark 1.

If the domain D is finite with a surface S, then the problem of type A with the conditions (1)–(3) has a unique solution depending continuously on the data (see [20,21]).

It should be noted that the condition (4) is essential for the uniqueness of solution of the Problem 1. Indeed, if out of (4) Problem 1 has any solution $u_1\left(x\right)$, then $u_2\left(x\right)=u_1\left(x\right)+k(r-R)/r$ is its solution as well, where $r\ge R$ and k is a real constant, i.e., it has an infinite set of solutions.

On the other hand, if Problem 1 has a solution, then it is unique. Assume that Problem 1 with the boundary function (5) has two solutions $u_1\left(x\right)$ and $u_2\left(x\right)$, then $u\left(x\right)=u_1\left(x\right)-u_2\left(x\right)$ is a solution of Problem 1 with null boundary value and for function $u\left(x\right)$ condition (4) is fulfilled. It is known (see [19] p. 303) that, in the noted case $u\left(x\right)\equiv 0$ when $x\in \overline{D}$ (or $u_1\left(x\right)=u_2\left(x\right)$, $x\in \overline{D})$. The existence of solution of the Problem 1 will be shown in Section 3. Based on the above-mentioned, we can investigate by Problem 1 only such physical phenomena, which are damped at infinity.

In the case of external 3D Dirichlet generalized harmonic problems, the difficulties become more significant. In particular, there does not exist a standard algorithm which can be applied to a wide class of domains.

On the basis of noted above, for the numerical solution of the Problem 1 we should apply an algorithm which does not require the approximation of a boundary function and in which the existence of discontinuity curves is not ignored. The method of probabilistic solution (MPS) is one such method (see [20,21]), but its direct application in infinite domains is impossible.

Performed investigations showed (see [10,11,22,23,24,25,26,27]) that the MPS is ideally suited for numerical solution of both classical and generalized (2D and 3D) Dirichlet harmonic problems for a wide class of finite domains only.

Therefore, construction of efficient high accuracy computational schemes for the numerical solution of external 3D Dirichlet generalized harmonic problems (which can be applied to a wide class of domains) has both theoretical and practical importance.

3. Transition from the Infinite Domain $\mathit{D}$ to the Kernel ${\mathit{K}}_{\mathit{R}}$ by an Inversion and Consideration of a New Problem on the Basis of Kelvin’s Theorem

In the Problem 1 the domain D is infinite, therefore, the direct application of the MPS to its solving is impossible. In order to solve the Problem 1 using MPS, we perform transition from the domain D to the kernel $K_R$ by means of the inversion (see [4,28]). Let a point $x(x_1,x_2,x_3)\in D$ and consider the following inversions:

(6)$\begin{array}{c}{\xi }_i=\frac{R^2}{{\left|x\right|}^2}{x}_i,x\in D,(i=1,2,3),{\left|x\right|}^2=x_1^2+x_2^2+x_3^2,\end{array}$

(7)$\begin{array}{c}{x}_i=\frac{R^2}{{\left|\xi \right|}^2}{\xi }_i,\xi \in K_R{\backslash \left\{O\right\},(i=1,2,3),\left|\xi \right|}^2={\xi }_1^2+{\xi }_2^2+{\xi }_3^2,\end{array}$

with respect to the sphere $S_R$. The points x and $\xi$ are called symmetric points with respect to the sphere $S_R$. From Equation (6) (or Equation (7)) we have

(8)$\left|x\right|\left|\xi \right|=R^2,\left|\xi \right|\ne 0.$

It is known (see [4] p. 260) that on the basis of Equation (6) (or Equation (7)) the symmetric points $\xi$ and x with respect to the sphere $S_R$ are situated on the ray, whose beginning is at the point $\left|\xi \right|=0$ (or $\xi \equiv O$).

It is easy to see that by Equation (8) the infinite domain D is transformed one-to-one onto the $K_R$. In particular, the points of $S_R$ are transformed into itself, and the point $x=\infty$ is transformed into the point $\xi =O$ and vice versa.

On the basis of Equation (7) the functions $u\left(x\right)$ and $g\left(y\right)$ are transformed into the functions

(9)$u\left(\xi \right)\equiv u\left(\frac{R^2}{{\left|\xi \right|}^2}{\xi }_1,\frac{R^2}{{\left|\xi \right|}^2}{\xi }_2,\frac{R^2}{{\left|\xi \right|}^2}{\xi }_3\right)\mathrm{and}g\left(\eta \right)\equiv g\left(\frac{R^2}{{\left|\eta \right|}^2}{\eta }_1,\frac{R^2}{{\left|\eta \right|}^2}{\eta }_2,\frac{R^2}{{\left|\eta \right|}^2}{\eta }_3\right),$

respectively, where $\xi ({\xi }_1,{\xi }_2,{\xi }_3)\in K_R\backslash S_R,$, $\left|\xi \right|\ne 0$, $\eta ({\eta }_1,{\eta }_2,{\eta }_3)\in S_R$.

It is easy to see that the function $u\left(\xi \right)$ is not harmonic in the domain $K_R\backslash S_R$. However, we can remove the noted defect, if we apply Kelvin’s theorem [4,28].

Theorem 1.

If a function $u(x_1,x_2,x_3)$ is harmonic in the infinite domain D, then the function

(10) $v({\xi }_1,{\xi }_2,{\xi }_3)=\frac{R}{\left|\xi \right|}u\left(\xi \right)$

is harmonic in the domain $K_R\backslash \left\{O\right\}$.

The function $v\left(\xi \right)$ is bounded in the neighborhood of the point ${\xi }^0=(0,0,0)$, therefore this point for $v\left(\xi \right)$ is a removable singular point [4,28]. However, actually, on the basis of the Theorem 1, for extension of harmonicity of function $v\left(\xi \right)$ at the point ${\xi }^0$ we must solve the following Dirichlet generalized harmonic problem.

Problem 2

($A^{*}$). Find a generalized harmonic function $v\left(\xi \right)$ satisfying the conditions:

(11)$\begin{array}{cc}\hfill \Delta v\left(\xi \right)=0,& \xi \in K_R\backslash S_R,\hfill \end{array}$

(12)$\begin{array}{cc}\hfill v\left(\eta \right)=g\left(\eta \right),& \eta \in S_R,\eta \overline{\in }{l}_k\subset S_R,\hfill \end{array}$

(13)$\begin{array}{cc}\hfill \left|v\right(\xi \left)\right|<c,& \xi \in K_R,\hfill \end{array}$

where in Equation (12): $\left|\eta \right|=R$, $\eta \equiv y$ ($y\in S_R$), $l_k$ ($k=1,\dots ,n)$ and c are the same that in Problem 1.

It is known (see [20,21]) that the generalized problem (11)–(13) is well-posed, i.e., the solution exists, is unique and depends continuously on the data. Respectively, Problem 1 is correct.

It is evident that for solving the Problem 2 we can apply the MPS. In particular, if we want to find the value of the solution $u\left(x\right)$ of Problem 1 at a point x ($x\in D$), first of all we have to find the image $\xi$ of x by means of (6) and then find the solution $v\left(\xi \right)$ of the Problem 2 at the point $\xi$.

Finally, on the basis of Equation (10) we have

(14)$u\left(x\right)=\frac{v\left(\xi \right)}{\left|x\right|}R,$

where $x\in D,\xi \in K_R\backslash S_R$.

Thus, for definition of value of the solution to Problem 1 at the point x we have the formula (14). It is easy to see that for the function $u\left(x\right)$, defined by Equation (14), the conditions of Problem 1 are fulfilled.

4. The Method of Probabilistic Solution and Simulation of the Wiener Process

In this section we describe the solution of Problem 2 by the MPS. It is known (see [21,24]) that the probabilistic solution of Problem 2 at the fixed point $\xi \in K_R\backslash S_R$ has the following form

(15)$v\left(\xi \right)=E_{\xi }g\left(\xi \left(\tau \right)\right).$

In Equation (15) $E_{\xi }g\left(\xi \left(\tau \right)\right)$ is the mathematical expectation of values of the boundary function $g\left(\eta \right)$ at the random intersection points of the trajectory of the Wiener process and the boundary $S_R$; $\tau$ is a random moment of the first exit of the Wiener process $\xi \left(t\right)=({\xi }_1\left(t\right),{\xi }_2\left(t\right),{\xi }_3\left(t\right))$ from the domain $K_R$. It is assumed that the starting point of the Wiener process is always $\xi \left(t_0\right)=({\xi }_1\left(t_0\right),{\xi }_2\left(t_0\right),{\xi }_3\left(t_0\right))\in K_R\backslash S_R$, where the value of desired function is being determined. If the number N of the random intersection points ${\eta }^j=({\eta }_1^j,{\eta }_2^j,{\eta }_3^j)\in S_R$ ($j=1,\dots ,N$) is sufficiently large, then according to the law of large numbers, from (15) we have

(16)$v\left(\xi \right)\approx v_N\left(\xi \right)=\frac{1}{N}\sum _{j=1}^{N}g\left({\eta }^j\right)$

or $v\left(\xi \right)=lim{v}_N\left(\xi \right)$ for $N\to \infty$, in probability. Thus, if we have the Wiener process, the approximate value of the probabilistic solution to the Problem 2 at point $\xi \in K_R\backslash S_R$ is calculated by the formula (16).

In order to simulate the Wiener process we use the following recursion relations (see [20,24]):

(17)$\begin{array}{cc}\hfill {\xi }_1\left(t_k\right)=& {\xi }_1\left(t_{k-1}\right)+{\gamma }_1\left(t_k\right)/nq,\hfill \\ \hfill {\xi }_2\left(t_k\right)=& {\xi }_2\left(t_{k-1}\right)+{\gamma }_2\left(t_k\right)/nq,\hfill \\ \hfill {\xi }_3\left(t_k\right)=& {\xi }_3\left(t_{k-1}\right)+{\gamma }_3\left(t_k\right)/nq,,(k=1,2,\cdots ),\hfill \\ \hfill \xi \left(t_0\right)=& \xi ,\hfill \end{array}$

according to which the coordinates of point $\xi \left(t_k\right)=({\xi }_1\left(t_k\right),{\xi }_2\left(t_k\right),{\xi }_3\left(t_k\right))$ are being determined. In Equation (17) ${\gamma }_1\left(t_k\right),{\gamma }_2\left(t_k\right),{\gamma }_3\left(t_k\right)$ are three normally distributed independent random numbers for the k-th step, with zero means and variances equal to one (The above numbers are generated apart); $nq$ is the number of quantification $\left(nq\right)$ such that $1/nq=\sqrt{t_k-t_{k-1}}$ and when $nq\to \infty$, then the discrete process approaches the continuous Wiener process. In the implementation, random process is simulated at each step of the walk and continues until its trajectory crosses the boundary.

In the considered case computations are performed and the random numbers are generated using MATLAB function randn which generates normally distributed random numbers.

Remark 2.

In general, problems of type A can be solved by the MPS for all such locations of discontinuity curves, which give the possibility to establish the part of surface S where the intersection point is located.

5. Numerical Examples

It should be noted that in the 3D case there are no test solutions for generalized problems of type 2, therefore, for the verification of the scheme needed for the numerical solution of Problem 2, the reliability of obtained results can be demonstrated in the following way.

If in boundary conditions (12) of Problem 2 we take $g_i\left(\eta \right)=1/|\eta -{\xi }^0|$, where $\eta \in S_R^i$ ($i=1,\dots ,m$), ${\xi }^0=({\xi }_1^0,{\xi }_2^0,{\xi }_3^0)\overline{\in }{K}_R$, and $|\eta -{\xi }^0|$ denotes the distance between points $\eta$ and $x^0$, then it is evident that the curves $l_k$ ($k=1,\dots ,n$) represent removable discontinuity curves for the boundary function $g\left(\eta \right)$. Actually, in the mentioned case instead of generalized Problem 2 we obtain the following Dirichlet classical harmonic problem.

Problem 3

(B). Find a Function $v\left(\xi \right)\equiv v({\xi }_1,{\xi }_2,{\xi }_3)\in C^2(K_R\backslash S_R)\bigcap C\left(K_R\right)$ satisfying the conditions:

(18)$\begin{array}{cccc}\hfill \Delta v\left(\xi \right)& =0,\hfill & \hfill & \xi \in K_R\backslash S_R,\hfill \end{array}$

(19)$\begin{array}{cccc}\hfill v\left(\eta \right)& =1/|\eta -{\xi }^0|,\hfill & \hfill & \eta \in S_R,{\xi }^0\overline{\in }{K}_R.\hfill \end{array}$

We solve this problem (by the MPS) using the program used for the Problem 2. It is well-known that the Problem 3 is well posed, i.e., its solution exists, is unique and depends on data continuously. Evidently, an exact solution of the Problem 3 is

(20)$v({\xi }^0,\xi )=\frac{1}{|\xi -{\xi }^0|},\xi \in K_R,{\xi }^0\overline{\in }{K}_R.$

Note that application of the MPS for numerical solution of the Dirichlet classical harmonic problems is interesting and important (see [10,11,22,23]). In this paper, the Problem 3 has an auxiliary role. In particular, for the Problem 3, verification of the scheme needed for the numerical solution of Problem 2 and corresponding program (comparison of the obtained results with exact solution) are carried out first of all, and then actually Problem 2 is solved under the boundary conditions (5).

In the present paper, MPS is applied for two examples. In the tables, N is a number of implementation of the Wiener process for the given points ${\xi }^i=({\xi }_1^i,{\xi }_2^i,{\xi }_3^i)\in K_R\backslash S_R$, and $nq$ is a number of quantification. For simplicity, in the considered examples the values of $nq$ and N are the same. In the tables for problems of type 3 the absolute errors ${\Delta }^i$ at the points ${\xi }^i\in K_R\backslash S_R$ of $v_N\left(\xi \right)$ are presented in the MPS approximation, for $nq=100$ and various values of N. In particular, we have ${\Delta }^i=|v_N\left({\xi }^i\right)-v({\xi }^0,{\xi }^i)|$, where $v_N\left({\xi }^i\right)$ is the approximate solution of Problem 3 at the point ${\xi }^i$, which is defined by formula (16), and $v({\xi }^0,{\xi }^i)$ is an exact solution of the test problem is given by Equation (20). In tables, for the problems of type 2, the probabilistic solution $v_N\left(\xi \right)$, defined by Equation (16), is presented at the points ${\xi }^i\in K_R\backslash S_R$.

Remark 3.

The problems of type 2 and 3 for ellipsoidal, spherical, cylindrical, conic, prismatic, pyramidal and axisymmetric finite domains with a cylindrical hole are considered in [24,25,26,27].

Example 1.

In the first example exterior of the unit sphere $S_1:y_1^2+y_2^2+y_3^2=1$, with the center at origin $O(0,0,0)$ and radius $R=1$ is considered in the role of domain D, where $(y_1,y_2,y_3)$ is a point of the surface $S_1$. In the considered case, in Problem 1 the boundary function $g\left(y\right)$ has the following form

(21) $g\left(y\right)=\left\{\begin{array}{cc}1,\hfill & y\in S^1=\{y\in S_1|-1\le y_3<-0.5\},\hfill \\ 2,\hfill & y\in S^2=\{y\in S_1|-0.5<y_3<0.5\},\hfill \\ 1,\hfill & y\in S^3=\{y\in S_1|0.5<y_3\le 1\},\hfill \\ 0,\hfill & y\in l_k(k=1,2).\hfill \end{array}\right.$

It is evident that in Equation (21), the discontinuity curves $l_1$, $l_2$ are the circles, which are obtained by intersection of the planes $x_3=-0.5$, $x_3=0.5$, and the surface $S_1$ (in the physical sense curves $l_1$ and $l_2$ are non-conductors).

According to above noted, for solving the Problem 1, in the first place, we solve Problem 2 with boundary function $g\left(\eta \right)$ ($g\left(\eta \right)\equiv g\left(y\right)$) by the MPS. In order to determine the intersection points ${\eta }^j=({\eta }_1^j,{\eta }_2^j,{\eta }_3^j)$ ($j=1,\dots ,N$) of Wiener process trajectory and the surface $S_1$, we operate the following way.

During the implementation of Wiener process, for each current point $\xi \left(t_k\right)$, defined from Equation (17), its location with respect to $S_1$ is checked, i.e., for the point $\xi \left(t_k\right)$ the value

$d={\xi }_1^2+{\xi }_2^2+{\xi }_3^2$

is calculated and the following conditions: $d=1$, $d<1$, $d>1$ are checked. If $d=1$ then $\xi \left(t_k\right)\in S_1$ and ${\eta }^j=\xi \left(t_k\right)$. If $d<1$ then $\xi \left(t_k\right)\in K_1\backslash S_1$ and the process continues until its trajectory crosses the sphere $S_1$, and if $d>1$ then $\xi \left(t_k\right)\overline{\in }{K}_1$.

Let $\xi \left(t_{k-1}\right)\in K_1\backslash S_1$ for the moment $t=t_{k-1}$ and $\xi \left(t_k\right)\overline{\in }{K}_1$ for the moment $t=t_k$. For determination of the point ${\eta }^j$, a parametric equation of a line L passing through the points $\xi \left(t_{k-1}\right)$ and $\xi \left(t_k\right)$ is first obtained:

(22)$\left\{\begin{array}{cc}{\xi }_1\hfill & ={\xi }_1^{k-1}+({\xi }_1^k-{\xi }_1^{k-1})\theta ,\hfill \\ {\xi }_2\hfill & ={\xi }_2^{k-1}+({\xi }_2^k-{\xi }_2^{k-1})\theta ,\hfill \\ {\xi }_3\hfill & ={\xi }_3^{k-1}+({\xi }_3^k-{\xi }_3^{k-1})\theta ,\hfill \end{array}\right.$

where $({\xi }_1,{\xi }_2,{\xi }_3)$ is a point of L, $\theta$ is a parameter $(-\infty <\theta <\infty )$, and ${\xi }_i^{k-1}\equiv {\xi }_i\left(t_{k-1}\right),{\xi }_i^k\equiv {\xi }_i\left(t_k\right)$ ($i=1,2,3$). After this, for definition of the intersection points ${\eta }^{*}$ and ${\eta }^{**}$ of line L and the sphere $S_1$ is solved an equation with respect to $\theta$.

It is easy to see that for the parameter $\theta$ we obtain an equation

(23)$A{\theta }^2+2B\theta +C=0,$

whose discriminant $d^{*}=B^2-AC>0.$

Since the discriminant of Equation (23) is positive, the Equation (23) has two solutions ${\theta }_1$ and ${\theta }_2$. Respectively, the points ${\eta }^{*}$ and ${\eta }^{**}$ are defined from the Equation (22). In the role of the point ${\eta }^j$ we choose the one from ${\eta }^{*}$ and ${\eta }^{**}$ for which $\left|\xi \right(t_k)-\eta |$ is minimal.

In both examples, considered by us for the determination the intersection points ${\eta }^j=({\eta }_1^j,{\eta }_2^j,{\eta }_3^j)$ ($j=1,\dots ,N$) of the trajectory of the Wiener process and the surface $S_1$ is used that scheme, which is above described. For checking the reliability of calculated results first we solved the auxiliary Problem 3 using the above-described algorithm. In the numerical experiments for the example 1, in test Problem 3, we took ${\xi }^0=(0,5,0)$.

In Table 1 (see also Figure 1), the absolute errors ${\Delta }^i$ of the approximate solution $v_N\left(\xi \right)$ of the test Problem 3 at the points ${\xi }^i\in K_1\backslash S_1$ ($i=1,\dots ,5$) are presented. The notation ($E\pm k$) is used for ${10}^{\pm k}$.

On the basis of Table 1 we see that the obtained results are sufficiently close to the expected results.

We conducted the check experiment. Namely, we calculated the probabilistic solution of Problem 3 at the point $(0,0,0)$ for $N=1\times {10}^5$, $nq=200$ and we obtained ${\Delta }^1=0.18\times {10}^{-4}$ (see Table 1). The result is improved, as expected (see Section 4). In general, increasing $nq$ and N causes the improving of accuracy. For this, it is appropiate to use parallel computing systems and in this case the computational complexity of the algorithm is a matter of interest. Multicore CPUs and/or remote powerful computers can be necessary.

In Table 2 the values of the approximate solution $v_N\left(\xi \right)$ to the Problem 2 at the same points ${\xi }^i$ ($i=1,\dots ,5$) are given. The boundary function (21) is symmetric with respect to the plane $O{x}_1{x}_2$. Respectively, in the role of ${\xi }^2$, ${\xi }^3$ and ${\xi }^4$, ${\xi }^5$, the points, which are symmetric with respect to the plane $O{x}_1{x}_2$, are taken. The results have sufficient accuracy for many practical problems and agrees with the real physical picture.

In Table 3 the values of approximate solution $u_N\left(x\right)$ to Problem 1 at the points $x^i\in D$ ($i=1,\dots ,7$) are given. For definition the values $u_N\left(x^i\right)$ the formula (14) is applied, respectively, in the role of points $x^i$ and ${\xi }^i$ ($i=1,\dots ,7$), the points, which are symmetric with respect to the sphere $S_1$, are taken. It is evident that the points $x^i$ are situated on the ray, whose beginning is at the point $\xi =(0,0,0)$ and located in the plane $O{x}_2{x}_3$. The ray’s inclination angle with plane $O{x}_1{x}_2$ is equal to $\pi /4.$

The points ${\xi }^i$ ($i=1,\dots ,7$) are obtained by the inversion (6). The values $v_N\left({\xi }^i\right)$ are the approximate solution $v_N\left(\xi \right)$ to Problem 2 at the points ${\xi }^i$ for $N=1E+5$ and $nq=100.$ The obtained results agrees with the real physical picture.

Example 2.

Here in a role of infinite domain D we took again the exterior of unit sphere $S_1:y_1^2+y_2^2+y_3^2=1$ with the center at origin $O(0,0,0)$, where $y(y_1,y_2,y_3)$ is a point of the surface $S_1$. In considered case the boundary function $g\left(y\right)=g(y_1,y_2,y_3)$ has the following form

(24) $g\left(y\right)=\left\{\begin{array}{cc}1,\hfill & y\in S^1=\{y\in S_1|y_1>0,y_2>0,y_3>0\},\hfill \\ 0,\hfill & y\in S^2=\{y\in S_1|y_1<0,y_2>0,y_3>0\},\hfill \\ 1,\hfill & y\in S^3=\{y\in S_1|y_1<0,y_2<0,y_3>0\},\hfill \\ 0,\hfill & y\in S^4=\{y\in S_1|y_1>0,y_2<0,y_3>0\},\hfill \\ 2,\hfill & y\in S^5=\{y\in S_1|y_1>0,y_2>0,y_3<0\},\hfill \\ 1,\hfill & y\in S^6=\{y\in S_1|y_1<0,y_2>0,y_3<0\},\hfill \\ 2,\hfill & y\in S^7=\{y\in S_1|y_1<0,y_2<0,y_3<0\},\hfill \\ 1,\hfill & y\in S^8=\{y\in S_1|y_1>0,y_2<0,y_3<0\},\hfill \\ 0,\hfill & y\in l_k(k=1,2,3).\hfill \end{array}\right.$

It is evident that in this case the discontinuity curves $l_k(k=1,2,3)$ are the circles, obtained by intersection of the coordinate planes and the sphere $S_1$. Actually, the sphere $S_1$ is divided into equal parts $S^i$ ($i=1,\dots ,8$) by curves $l_k$. In addition, in the considered case, $l_k$, $S^2$, $S^4$ are non-conductors.

Since in this case the problem domain is the same that in Example 1, therefore for the determination of intersection points ${\eta }^j$ ($j=1,\dots ,N$) of the Wiener process trajectory and the sphere $S_1$, the same algorithm, described in Example 1, is applied.

The values of the numerical solution of Problem 2 at the points ${\xi }^i\in K_1\backslash S_1$ ($i=1,\dots ,5$) are given in Table 4. Since the boundary function (24) is symmetric with respect to the axis $O{x}_3$, therefore, in the role of ${\xi }^i$ ($i=4,5$), the points which are symmetric with respect to the axis $O{x}_3$ are taken. The obtained results have sufficient accuracy for many practical problems and agrees with the real physical picture.

The values of the numerical solution of Problem 1 at the points $x^i\in D$ ($i=1,\dots ,10$) are given in Table 5. The values $u_N\left(x^i\right)$ are obtained by the formula (14), therefore, in the role of points $x^i$ and ${\xi }^i$ ($i=1,\dots ,10$), the points, symmetric to the sphere $S_1$ are taken. The points $x^i$ and ${\xi }^i$ are situated on the axis $O{x}_3$, and the points ${\xi }^i$ are obtained by Equation (6). The values $v_N\left({\xi }^i\right)$ are numerical solution to Problem 2 at the points ${\xi }^i$ ($i=1,\dots ,10$) for $N=1\times {10}^5$ and $nq=100$. The obtained results are agrees with the real physical picture.

In this work we solved problems of type 1 when the boundary functions $g_i\left(y\right)$ ($i=1,\dots ,m$) are constants. This was motivated by our interest to find out how well the obtained results agree with the real physical picture. It is evident that solving Problem 1 under condition (5) is as easy as the Problem 3.

The analysis of results of numerical experiments show that the results obtained by the proposed algorithm are reliable and it is effective for the numerical solution of problems of type 1 and 3. In addition, the algorithm is sufficiently simple for numerical implementation.

6. Concluding Remarks

• 1.. In this work, we have demonstrated that the suggested algorithm is ideally suited for numerically solving of considered problems 1 and 3. It should be noted that it is possible using this algorithm to find the solution of a problem at any point of the domain, unlike other algorithms known in the literature.

• 2.. The algorithm does not require approximation of the boundary function, which is one of its important properties.

• 3.. The algorithm is easy to program, its computational cost is low, it is characterized by an accuracy which is sufficient for many practical problems.

• 4.. In the future we plan to investigate the following:

  • Application of the proposed algorithm to the numerical solution of Dirichlet classical and generalized harmonic problems for the infinite space $R^3$ with a finite number of spherical cavities.

  • Application of the MPS for the same type of problem for iregular pyramidal domains.

  • Application of the MPS for the same type of problem in finite domains which are bounded by several closed surfaces.

  • Application of the MPS for the same problem for infinite 2D domains with a finite number of circular holes.

Author Contributions

Conceptualization, M.Z., Z.T., N.K. and T.D.; methodology, M.Z., Z.T., N.K., T.D., J.M.S. and F.C.-A.; formal analysis, M.Z., Z.T., N.K., T.D. and F.C.-A.; investigation, M.Z., Z.T., N.K., T.D., J.M.S. and F.C.-A.; writing—original draft, M.Z., Z.T., N.K. and T.D.; writing—review and editing, J.M.S. and F.C.-A. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

The research was conducted in the frame of the 5-year scientific plan of the Niko Muskhelishvili Institute of Computational Mathematics.

Conflicts of Interest

The Authors declare no conflict of interests.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figure and Tables

Figure 1. Results from Table 1, for Problem 3 (in Example 1).

Table 1

Results for Problem 3 (in Example 1).

${\mathit{\xi }}^{\mathit{i}}$	$(0,0,0)$	$(0,0,0.5)$	$(0,0,-0.5)$	$(0,0,0.8)$	$(0,0,-0.8)$
$\mathit{N}$	${\Delta }^{\mathbf{1}}$	${\Delta }^{\mathbf{2}}$	${\Delta }^{\mathbf{3}}$	${\Delta }^{\mathbf{4}}$	${\Delta }^{\mathbf{5}}$
5 × 10${}^3$	0.47 × 10${}^{-3}$	0.40 × 10${}^{-3}$	0.85 × 10${}^{-5}$	0.28 × 10${}^{-3}$	0.23 × 10${}^{-3}$
1 × 10${}^4$	0.28 × 10${}^{-3}$	0.25 × 10${}^{-3}$	0.10 × 10${}^{-3}$	0.13 × 10${}^{-3}$	0.41 × 10${}^{-4}$
5 × 10${}^4$	0.62 × 10${}^{-4}$	0.17 × 10${}^{-3}$	0.83 × 10${}^{-4}$	0.41 × 10${}^{-4}$	0.28 × 10${}^{-4}$
1 × 10${}^5$	0.57 × 10${}^{-4}$	0.44 × 10${}^{-4}$	0.45 × 10${}^{-4}$	0.36 × 10${}^{-4}$	0.33 × 10${}^{-4}$
5 × 10${}^5$	0.54 × 10${}^{-4}$	0.56 × 10${}^{-4}$	0.28 × 10${}^{-4}$	0.20 × 10${}^{-5}$	0.35 × 10${}^{-4}$
1 × 10${}^6$	0.29 × 10${}^{-5}$	0.43 × 10${}^{-4}$	0.30 × 10${}^{-4}$	0.51 × 10${}^{-4}$	0.28 × 10${}^{-4}$

Table 2

Results for Problem 2 (in Example 1).

${\mathit{\xi }}^{\mathit{i}}$	$(0,0,0)$	$(0,0,0.5)$	$(0,0,-0.5)$	$(0,0,0.8)$	$(0,0,-0.8)$
$\mathit{N}$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{1}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{2}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{3}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{4}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{5}}\right)$
5 × 10${}^3$	1.50160	1.70600	1.69820	1.89940	1.90000
1 × 10${}^4$	1.49500	1.69480	1.69440	1.89710	1.89430
5 × 10${}^4$	1.50448	1.70062	1.69948	1.89780	1.89870
1 × 10${}^5$	1.49664	1.69822	1.70120	1.89639	1.89609
5 × 10${}^5$	1.50069	1.69935	1.69917	1.89577	1.89609
1 × 10${}^6$	1.50009	1.69903	1.69849	1.89517	1.89587

Table 3

Results for Problem 1 (in Example 1).

${\mathit{x}}^{\mathit{i}}$	${\mathit{\xi }}^{\mathit{i}}$	${\mathit{v}}_{\mathit{N}}\left({\mathit{\xi }}^{\mathit{i}}\right)$	${\mathit{u}}_{\mathit{N}}\left({\mathit{x}}^{\mathit{i}}\right)$
∞	(0,0,0)	1.50069	0
(0, 5000, 5000)	(0, 0.0001, 0.0001)	1.50038	2.12 × 10${}^{-4}$
(0, 1000, 1000)	(0, 0.0005, 0.0005)	1.49971	0.00106
(0, 100, 100)	(0, 0.005, 0.005)	1.49957	0.01060
(0, 50, 50)	(0, 0.01, 0.01)	1.49990	0.02121
(0, 10, 10)	(0, 0.05, 0.05)	1.50153	0.10617
(0, 1, 1)	(0, 0.5, 0.5)	1.67766	1.18628

Table 4

Results for the Problem 2 (in Example 2).

${\mathit{\xi }}^{\mathit{i}}$	$(0,0,0)$	$(0,0,0.8)$	$(0,0,-0.8)$	$(0.5,0.5,0.5)$	$(-0.5,-0.5,0.5)$
$\mathit{N}$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{1}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{2}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{3}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{4}}\right)$	${\mathit{v}}_{\mathit{N}}\left({\xi }^{\mathbf{5}}\right)$
5 × 103	1.00580	0.56320	1.46020	0.97100	0.97180
1 × 104	0.99970	0.54620	1.45060	0.96960	0.06820
5 × 104	0.99992	0.55262	1.44466	0.96930	0.96774
1 × 105	0.99766	0.54962	1.45056	0.96976	0.96990
5 × 105	1.00060	0.55173	1.44925	0.97012	0.96989
1 × 106	1.00002	0.55116	1.44936	0.97000	0.97020

Table 5

Results for Problem 1 (in Example 2).

${\mathit{x}}^{\mathit{i}}$	${\mathit{\xi }}^{\mathit{i}}$	${\mathit{v}}_{\mathit{N}}\left({\mathit{\xi }}^{\mathit{i}}\right)$	${\mathit{u}}_{\mathit{N}}\left({\mathit{x}}^{\mathit{i}}\right)$
∞	(0, 0, 0)	1.00060	0
(0, 0, 10${}^4$)	(0, 0, 10${}^{-4}$)	0.99961	9.9961 × 10${}^{-5}$
(0, 0, 10${}^3$)	(0, 0, 10${}^{-3}$)	1.00176	0.00102
(0, 0, 10${}^2$)	(0, 0, 10${}^{-2}$)	0.99147	0.00915
(0, 0, 10)	(0, 0, 10${}^{-1}$)	0.923912	0.09239
(0, 0, 5)	(0, 0, 0.2)	0.85483	0.17097
(0, 0, 4)	(0, 0, 0.25)	0.81979	0.20495
(0, 0, 2)	(0, 0, 0.5)	0.67068	0.33534
(0, 0, 1.5)	(0, 0, 2/3)	0.59743	0.39829
(0, 0, 1.2)	(0, 0, 5/6)	0.54302	0.45252

References

1. Grinberg, G.A. Selected Questions of the Mathematical Theory of Electrical and Magnetic Phenomena; Izd. AN SSSR: Moscow, Russia, 1948; (In Russian)

2. Smythe, W.B. Static and Dynamic Electricity; Hemisphere Publishing: New York, NY, USA, 1988.

3. Carslaw, H.S.; Jaeger, J.C. Conduction of Heat in Solids; 2nd ed. Oxford Science Publications Oxford University Press: New York, NY, USA, 1988; 522.

4. Košljakov, N.S.; Gliner, E.B.; Smirnov, M.M. The Differential Equations of Mathematical Physics; Košljakov, N.S. Gosudarstv. Izdat. Fiz.-Mat. Lit.: Moscow, Russia, 1962; 767.

5. Budak, B.M.; Samarskiĭ, A.A.; Tihonov, A.N. A Collection of Problems in Mathematical Physics; 3rd ed. Nauka: Moscow, Russia, 1980; 687.(In Russian)

6. Kantorovich, L.V.; Krylov, V.I. Approximate Methods of Higher Analysis; Benster, C.D. Interscience Publishers, Inc.: New York, NY, USA, 1958; xv+681.

7. Koblishvili, N.; Tabagari, Z.; Zakradze, M. On reduction of the Dirighlet generalized boundary value problem to an ordinary problem for harmonic function. Proc. A. Razmadze Math. Inst.; 2003; 132, pp. 93-106.

8. Zakradze, M.; Koblishvili, N.; Karageorghis, A.; Smyrlis, Y. On solving the Dirichlet generalized problem for harmonic function by the method of fundamental solutions. Semin. I. Vekua Inst. Appl. Math. Rep.; 2008; 34, pp. 24-32.

9. Kublashvili, M.; Sanikidze, Z.; Zakradze, M. A method of conformal mapping for solving the generalized Dirichlet problem of Laplace’s equation. Proc. A. Razmadze Math. Inst.; 2012; 160, pp. 71-89.

10. Chaduneli, A.; Zakradze, M.V.; Tabagari, Z. A method of probabilistic solution to the ordinary and generalized plane Dirichlet problem for the Laplace equation. Science and Computing. Proceedings of the Sixth TSTC Scientific Advisory Committee Seminar (Science and Computing); Moscow, Russia, 15–19 September 2003; Volume 2, pp. 361-366.

11. Chaduneli, A.; Zakradze, M.; Tabagari, Z. A computer simulation of probabilistic solution to the Dirichlet plane boundary problem for the Laplace equation in case of an infinite plane with a hole. Bull. Georgian Acad. Sci.; 2005; 171, pp. 437-440.

12. Bayliss, A.; Gunzburger, M.; Turkel, E. Boundary conditions for the numerical solution of elliptic equations in exterior regions. SIAM J. Appl. Math.; 1982; 42, pp. 430-451. [DOI: https://dx.doi.org/10.1137/0142032]

13. Ben Belgacem, F.; Fournié, M.; Gmati, N.; Jelassi, F. On the Schwarz algorithms for the elliptic exterior boundary value problems. M2AN Math. Model. Numer. Anal.; 2005; 39, pp. 693-714. [DOI: https://dx.doi.org/10.1051/m2an:2005030]

14. Kanunnikova, E.A. Analytic–numerical method to solve 3D exterior boundary problems for elliptic equations. St. Petersburg Polytech. Univ. J. Phys. Math.; 2014; 1, pp. 35-38.

15. Charalambopoulos, A.; Gergidis, L.N. A novel stochastic method for the solution of direct and inverse exterior elliptic problems. Quart. Appl. Math.; 2018; 76, pp. 65-111. [DOI: https://dx.doi.org/10.1090/qam/1480]

16. Stratton, J.A. Electromagnetic Theory; McGraw-Hill: New York, NY, USA, 1941.

17. Lamb, H. Hydrodynamics; Cambridge University Press: Cambridge, UK, 1993.

18. Kellogg, O.D. Foundations of Potential Theory; Die Grundlehren der Mathematischen Wissenschaften, Band 31 Springer: Berlin/Heidelberg, Germany, 1967; 384.

19. Tihonov, A.N.; Samarskiĭ, A.A. Equations of Mathematical Physics; Izdat. Nauka: Moscow, Russia, 1966; 724.

20. Zakradze, M.; Kublashvili, M.; Sanikidze, Z.; Koblishvili, N. Investigation and numerical solution of some 3D internal Dirichlet generalized harmonic problems in finite domains. Trans. A. Razmadze Math. Inst.; 2017; 171, pp. 103-110. [DOI: https://dx.doi.org/10.1016/j.trmi.2016.11.001]

21. Dynkin, E.B.; Juškevič, A.A. Theorems and Problems in Markov Processes; Izdat. “Nauka”: Moscow, Russia, 1967; 231.

22. Chaduneli, A.; Tabagari, Z.; Zakradze, M. On solving the internal three-dimensional Dirichlet problem for a harmonic function by the method of probabilistic solution. Bull. Georgian Natl. Acad. Sci.; 2008; 2, pp. 25-28.

23. Zakradze, M.; Sanikidze, Z.; Tabagari, Z. On solving the external three-dimensional Dirichlet problem for a harmonic function by the probabilistic method. Bull. Georgian Natl. Acad. Sci.; 2010; 4, pp. 19-23.

24. Zakradze, M.; Mamporia, B.; Kublashvili, M.; Koblishvili, N. The method of probabilistic solution for 3D Dirichlet ordinary and generalized harmonic problems in finite domains bounded with one surface. Trans. A. Razmadze Math. Inst.; 2018; 172, pp. 453-465. [DOI: https://dx.doi.org/10.1016/j.trmi.2018.08.005]

25. Zakradze, M.; Kublashvili, M.; Koblishvili, N.; Chakhvadze, A. The method of probabilistic solution for determination of electric and thermal stationary fields in conic and prismatic domains. Trans. A. Razmadze Math. Inst.; 2020; 174, pp. 235-246.

26. Zakradze, M.; Kublashvili, M.; Tabagari, Z.; Koblishvili, N. On numerical solving the Dirichlet generalized harmonic problem for regular n-sided pyramidal domains by the probabilistic method. Trans. A. Razmadze Math. Inst.; 2022; 176, pp. 123-132.

27. Zakradze, M.; Tabagari, Z.; Sanikidze, Z.; Abramidze, E. Computer modeling of a probabilistic solution for the Dirichlet generalized harmonic problem in some finite axisymmetric bodies with a cylindrical hole. Bull. TICMI; 2022; 26, pp. 37-52.

28. Vladimirov, V.S. Equations of Mathematical Physics. Pure and Applied Mathematics; Jeffrey, A. Littlewood, A. Marcel Dekker, Inc.: New York, NY, USA, 1971; Volume 3, 418.