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

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 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.


Introduction
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 40-s of the 20th century (see, e.g., [1][2][3][4][5]).
It is known (see e.g., [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 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, e.g., [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 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 to find the best solution (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 e.g., [1]).Besides, it should be noted that in the literature (see,e.g., [4], pp.346-348), while solving the Dirichlet external generalized harmonic simplest problem for a sphere, the existence of 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 singulari-ties 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, e.g., [7,8]); II) A method of conformal mapping (see, e.g., [9]); III) A method of probabilistic solution (see, e.g., [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 even 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 wellposedness 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 is obtained for several cases.
In [13], the Laplace-Dirichlet problem is investigated in three-dimensional case.The Laplace-Dirichlet problem is read in the same words as the Laplace-Neumann equation, excepted for the boundary condition.Corresponding variational formulation is considered.It based on the introduction of an auxiliary unknown by the means of 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].The implementation procedure of the considered analytical-numerical method for solving external boundary problems in three-dimensional semi-space reduces the problem to such a problem that it's possible using traditional techniques and methods of numerical analysis.
In [15] is considered the solution of direct and inverse exterior boundary value problems via the strongly conditioned stochastic method, mainly for exterior harmonic problems.Note that in numerical examples from the Section 4.2, was performed experiments with cones of different thickness and observation points in the near and far field.But, actually the observation points given in Table 2 are "not too far" from the field, when we are considering exterior boundary value problems.
It is known that elliptic problems in external domains arise in many branches of physics.For example, the Laplace equation (∆u(x) = 0, x ∈ D) arises in the studies of thermostatic and electrostatic fields external to given surfaces(see, e. g., [3,16]); the flow of an incompressible irrotational fluid around a body is described by the same equation (see, e. g., [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 three dimensional (3D) case the condition is lim u(x) = 0, for |x| → ∞, where |x| is the distance from a fixed (but chosen arbitrary) origin (see, e.g., [18,19]).
Let D = R 3 \K(O, R) be an infinite domain in the Euclidean space R 3 , where K(O, R) ≡ K R is a kernel with a spherical surface S(O, R) ≡ 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 A. The function g(y) 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(y).It is required to find a function u(x) ) is the Laplace operator, c is a real constant.
On the basis of (1.3), in general, the values of u(y) are not uniquely defined on the curves l k (k = 1, n).In particular, if Problem A concerns the determination of the thermal (or the electric) field, then we must take u(y) = 0 when y ∈ 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 i R (i = 1, m) by the curves l k (k = 1, n), where one of the following conditions holds: n = m, n < m, n > m.Thus, the boundary function g(y) has the following form where: Remark 1.If the domain D is finite with a surface S, then the problem of type A with the conditions (1.1), (1.2), (1.3) has a unique solution depending continuously on the data (see [20,21]).
It should be noted that the condition (1.4) is essential for the uniqueness of solution of the Problem A. Indeed, if out of (1.4) Problem A has any solution u 1 (x), then u 2 (x) = u 1 (x) + k(r − R)/r is its solution as well, where r ≥ R and k is a real constant, i.e., it has an infinite set of solutions.
On the other hand, if Problem A has a solution, then it is unique.Assume that Problem A with the boundary function (1.5) has two solutions u 1 (x) and u 2 (x), then u(x) = u 1 (x) − u 2 (x) is a solution of Problem A with null boundary value and for function u(x) condition (1.4) is fulfilled.It is known (see e.g., [19], p.303) that, in the noted case u(x) ≡ 0 when x ∈ D (or The existence of solution of the Problem A will be shown in Section 2. Based on the above-mentioned, we can investigate by Problem A 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 A 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.
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.

Transition from the infinite domain D to the kernel K R by an inversion and consideration of a new problem on the basis of Kelvin's theorem
In the Problem A the domain D is infinite, therefore it is impossible the direct application of the MPS to its solving.In order to solve the Problem A using MPS, we perform transition from the domain D to the kernel K R by means of the inversion (see e.g., [4,28]).Let a point x(x 1 , x 2 , x 3 ) ∈ D and consider the following inversions: with respect to the sphere S R .The points x and ξ are called symmetric points with respect to the sphere S R .From (2.1) (or (2.2)) we have It is known (see e.g., [4], p.260) that on the basis of (2.1) (or (2.2)) the symmetric points ξ and x with respect to the sphere S R are situated on the ray, whose beginning is at the point |ξ| = 0 (or ξ ≡ O).It is easy to see that by (2.3) 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 = ∞ is transformed into the point ξ = O and vice versa.
On the basis of (2.2) the functions u(x) and g(y) are transformed into the functions respectively, where ξ(ξ It is easy to see that the function u(ξ) is not harmonic in the domain K R \S R .But, 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 is harmonic in the domain K R \{O}.
The function v(ξ) is bounded in the neighborhood of the point ξ 0 = (0, 0, 0), therefore this point for v(ξ) is a removable singular point [4,28].But, actually, on the basis of the Theorem 1, for extension of harmonicity of function v(ξ) at the point ξ 0 we must solve the following Dirichlet generalized harmonic problem.
Problem A*.Find a generalized harmonic function v(ξ) satisfying the conditions: ) ) where in (2.7) It is known (see [20,21]) that the generalized problem (2.6), (2.7), (2.8) is correct, i.e., the solution exists, is unique and depends continuously on the data.Respectively, Problem A is correct.
It is evident that for solving the Problem A * we can apply the MPS.In particular, if we want to find the value of the solution u(x) of Problem A at a point x(x ∈ D), first of all we have to find the image ξ of x by means of (2.1) and then find the solution v(ξ) of the Problem A * at the point ξ.Finally, on the basis of (2.5) we have where Thus, for definition of value of the solution to Problem A at the point x we have the formula (2.9).It is easy to see that for the function u(x), defined by (2.9), the conditions of Problem A are fulfilled.

The method of probabilistic solution and simulation of the Wiener process
In this section we describe the solution of Problem A * by the MPS.It is known(see e.g., [21,24]), that the probabilistic solution of Problem A * at the fixed point ξ ∈ K R \S R has the following form In (3.1) E ξ g(ξ(τ )) is the mathematical expectation of values of the boundary function g(η) at the random intersection points of the trajectory of the Wiener process and the boundary S R ; τ is a random moment of the first exit of the Wiener process ξ(t) = (ξ 1 (t), ξ 2 (t), ξ 3 (t)) from the domain K R .It is assumed that the starting point of the Wiener process is always ξ(t 0 ) = (ξ 1 (t 0 ), ξ 2 (t 0 ), ξ 3 (t 0 )) ∈ K R \S R , where the value of desired function is being determined.If the number N of the random intersection points η j = (η j 1 , η j 2 , η j 3 ) ∈ S R (j = 1, N) is sufficiently large, then according to the law of large numbers, from (3.1) we have or v(ξ) = lim v N (ξ) for N → ∞, in probability.Thus, if we have the Wiener process, the approximate value of the probabilistic solution to the Problem In order to simulate the Wiener process we use the following recursion relations (see e.g., [20,24]): are being determined.In (3.3) γ 1 (t k ), γ 2 (t k ), γ 3 (t k ) 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 (nq) such that 1/nq = √ t k − t k−1 and when nq → ∞, 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 in MATLAB.
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.4. Numerical examples It should be noted that in the 3D case there are no test solutions for generalized problems of type A * , therefore, for the verification of the scheme needed for the numerical solution of Problem A * , the reliability of obtained results can be demonstrated in the following way.
If in boundary conditions (2.7) of Problem A * we take 3 ) ∈ K R , and |η − ξ 0 | denotes the distance between points η and x 0 , then it is evident that the curves l k (k = 1, n) represent removable discontinuity curves for the boundary function g(η).Actually, in the mentioned case instead of generalized problem A * we obtain the following Dirichlet classical harmonic problem.
We solve this problem (by the MPS) using the program used for the Problem A * .It is well-known that the Problem B is well posed, i.e., its solution exists, is unique and depends on data continuously.Evidently, an exact solution of the Problem B is Note that application of the MPS for numerical solution of the Dirichlet classical harmonic problems is interesting and important (see e.g., [10,11,22,23]).In this paper, the Problem B has an auxiliary role.In particular, for the Problem B, verification of the scheme needed for the numerical solution of Problem A * and corresponding program (comparison of the obtained results with exact solution) are carried out first of all, and then actually Problem A * is solved under the boundary conditions (1.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 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 B the absolute errors ∆ i at the points ξ i ∈ K R \S R of v N (ξ) are presented in the MPS approximation, for nq = 100 and various values of N. In particular, we have the approximate solution of Problem B at the point ξ i , which is defined by formula (3.2), and v(ξ 0 , ξ i ) is an exact solution of the test problem is given by (4.1).In tables, for the problems of type A * , the probabilistic solution v N (ξ), defined by (3.2), is presented at the points ξ i ∈ K R \S R .
Example 1.In the first example exterior of the unit sphere S 1 : y 2 1 +y 2 2 +y 2 3 = 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 A the boundary function g(y) has the following form It is evident that in (4.2), 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 A, in the first place, we solve Problem A * with boundary function g(η) (g(η) ≡ g(y)) by the MPS.In order to determine the intersection points η j = (η j 1 , η j 2 , η j 3 )(j = 1, 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 ξ(t k ), defined from (3.3), its location with respect to S 1 is checked, i.e., for the point ξ(t k ) the value is calculated and the following conditions: and the process continues until its trajectory crosses the sphere S 1 , and if d > 1 then ξ(t k )∈K 1 .Let ξ(t k−1 ) ∈ K 1 \S 1 for the moment t = t k−1 and ξ(t k )∈K 1 for the moment t = t k .For determination of the point η j , a parametric equation of a line L passing through the points ξ(t k−1 ) and ξ(t k ) is first obtained: where (ξ 1 , ξ 2 , ξ 3 ) is a point of L, θ is a parameter (−∞ < θ < ∞), and 3).After this , for definition of the intersection points η * and η * * of line L and the sphere S 1 is solved an equation with respect to θ.
It is easy to see that for the parameter θ we obtain an equation Since the discriminant of (4.4) is positive, the equation (4.4) has two solutions θ 1 and θ 2 .Respectively, the points η * and η * * are defined from the (4.3).In the role of the point η j we choose the one from η * and η * * for which |ξ(t k ) − η| is minimal.
In both examples, considered by us for the determination the intersection points η j = (η j 1 , η j 2 , η j 3 )(j = 1, N) of the trajectory of the Wiener process and the surface S 1 is used that scheme, which is above described.As noted above, for verification first we solve the auxiliary Problem B with the program of Problem A * .In the numerical experiments for the example 1, in test Problem B, we took ξ 0 = (0, 5, 0).

Table 4 .
1B. Results for Problem B (in Example 1 )

Table 4 .
1A * the values of the approximate solution v N (ξ) to the Problem A * at the same points ξ i (i = 1, 5) are given.The boundary function (4.2) is symmetric with respect to the plane Ox 1 x 2 .Respectively, in the role of ξ 2 , ξ 3 and ξ 4 , ξ 5 , the points, which are symmetric with respect to the plane Ox 1 x 2 , are taken.The results have sufficient accuracy for many practical problems and agrees with the real physical picture.

Table 4 .
1A. Results for Problem A (in Example 1)

Table 4 .
The values of the numerical solution of Problem A * at the points ξ i ∈ K 1 \S 1 (i = 1, 5) are given in Table4.2A* .Since the boundary function (4.5) is symmetric with respect to the axis Ox 3 , therefore, in the role of ξ i (i = 4, 5), the points which are symmetric with respect to the axis Ox 3 are taken.The obtained results have sufficient accuracy for many practical problems and agrees with the real physical picture.2A.Results for Problem A (in Example 2)