Certain issues in the analytical integration of the Boussinesq problem

Mykola Nagirniak, Marek ChaleckiW arsaw University of Life Sciences (SGGW), Polandmykola_nagimiak@sggw.edu.pl, http:// orcid.org/0000-0003-4996-7397marek_chalecki@sggw.edu.pl, http://orcid.org/0000-0003-3451-458X

Introduction

Exact solutions to the equations of the theory of elasticity, both for the 3D and 2D problems, can be obtained only for domains having geometrically simple boundaries and usually only for the space and half-space [1]. The Boussinesq solution [2] constitutes one of the fundamental solutions of the theory of elasticity and refers to an action of concentrated forces in a half-space. If treated as the Green function, it enables to obtain solutions for various shapes of surface loads applied in a plane z = 0 z = 0 z=0z=0z=0 of the elastic half-space. The Boussinesq solution is a basis of a majority of analyses concerning distribution of stresses and displacements in contact issues, especially during determining a settlement of soil under a foundation.
In the related literature, one can find analytical solutions obtained by F. Schleicher [3] and A. Love [4], who, using an integration of the Green functions, determined a displacement and stress field for a half-space loaded in a square domain. Cerruti [5] analyzed normal and tangential loads of a plane limiting a half-space. Terazawa [6] analyzed a number of issues related to various types of loads of an elastic half-space. Becker & Bevis [7] provided expressions for a displacement field (in an implicit form) in an elastic half-space as a result of action of a load on a square domain.
The literature review shows that the problem of displacements and stresses in an elastic half-space due to rectangular loads, known as the Love problem [4], is a classical issue of the engineering mechanics with beginnings reaching over 100 years ago. Despite its rich history and numerous investigations [7-12] which provided solutions for various loads and geometry, many of these results are dispersed and, as Jin et al. [9] have noted, not fully systematized. Moreover, newer research [13-15] focus on more complex models, such as transversally isotropic half-spaces, and that indicates a continuous evolution and significancy of this topic in contemporary investigations.
A relatively frequently applied approach in the contact issues is an application of the Green functions which enable an analytical solution of many problems in mechanics. However, in the practical calculations, many difficulties arise which are connected to an integration of these functions. Computer environments, like Wolfram Mathematica, employed to this analysis, can lead to various, sometimes incorrect results, depending on the software version applied.
The study is aimed on problems related to an analytic integration of the Green functions used to determine displacement and stress fields in an elastic half-space, under loads acting on domains of various shapes. Additionally, it has been performed an evaluation of concordance of results obtained in various versions of Wolfram Mathematica.
The scope of the study encompasses an analysis of problems highlighted above, especially a comparison of results obtained in various versions of the Mathematica software.
The originality of this study lies in identifying and analyzing discrepancies in results of symbolic integration between different versions (releases) of Wolfram Mathematica, rather than in deriving new analytical expressions. These differences lead to mathematically non-equivalent results, which may affect the interpretation and reliability of symbolic computations in engineering applications.

Theoretical Background

The study concerns the problem of an integration of the Green function describing vertical displacements in an elastic half-space. It has been considered a case of a uniform load on a square domain. The analytical solution of this problem is known (e.g. [3], [4], [7]). The integration of the Green functions has been performed in the environment Wolfram Mathematica v. 8.0 [16], 11.3 [17], 12.3 [18], 13.3 [19] and 14.2 [20]. This section focuses on the mathematical aspects of symbolic integration of Green functions, which are essential for the evaluation of displacements and stresses in an elastic half-space.

Green functions

Vertical displacements of an elastic medium loaded with a concentrated force applied in the origin of the coordinate system, for any spatial variables x , y , z x , y , z x,y,zx, y, zx,y,z, are determined from the formula (cf. [21-24]):
(1) w = P 4 π μ R 1 ( z 2 R 1 2 + 2 ( 1 − v ) ) (1) w = P 4 π μ R 1 z 2 R 1 2 + 2 ( 1 − v ) {:(1)w=(P)/(4pi muR_(1))((z^(2))/(R_(1)^(2))+2(1-v)):}\begin{equation*} w=\frac{P}{4 \pi \mu R_{1}}\left(\frac{z^{2}}{R_{1}{ }^{2}}+2(1-v)\right) \tag{1} \end{equation*}(1)w=P4πμR1(z2R12+2(1−v))
where
r = x 2 + y 2 , R 1 = r 2 + z 2 , μ = G = E 2 ( 1 + v ) r = x 2 + y 2 , R 1 = r 2 + z 2 , μ = G = E 2 ( 1 + v ) r=sqrt(x^(2)+y^(2)),quadR_(1)=sqrt(r^(2)+z^(2)),quad mu=G=(E)/(2(1+v))r=\sqrt{x^{2}+y^{2}}, \quad R_{1}=\sqrt{r^{2}+z^{2}}, \quad \mu=G=\frac{E}{2(1+v)}r=x2+y2,R1=r2+z2,μ=G=E2(1+v)
E − E − E-E-E− Young modulus, μ − μ − mu-\mu-μ− Lamé constant, G − G − G-G-G− Kirchhoff modulus, ν − ν − nu-\nu-ν− Poisson ratio. Form. (1) can be treated as a Green function to obtain displacements and stresses in a half-space loaded in any way on the surface z = 0 z = 0 z=0z=0z=0 [22]. This formula enables also a solution for an elastic half-space if the load is applied on domains of various shapes.
In the case of a half-space loaded on a square domain, and treating Form. (1) as the Green function (cf. Fig. 1a), the formula for the displacements directed perpendicularly to the plane closing the half-space w = w ( x , y , z ) w = w ( x , y , z ) w=w(x,y,z)w=w(x, y, z)w=w(x,y,z) can be presented in a form:
(2) w = ( 1 + v ) p 2 π E ∫ − a a ∫ − b b 1 R [ 2 ( 1 − v ) + z R 2 ] d x 0 d y 0 (2) w = ( 1 + v ) p 2 π E ∫ − a a   ∫ − b b   1 R 2 ( 1 − v ) + z R 2 d x 0 d y 0 {:(2)w=((1+v)p)/(2pi E)int_(-a)^(a)int_(-b)^(b)(1)/(R)[2(1-v)+(z)/(R^(2))]dx_(0)dy_(0):}\begin{equation*} w=\frac{(1+v) p}{2 \pi E} \int_{-a}^{a} \int_{-b}^{b} \frac{1}{R}\left[2(1-v)+\frac{z}{R^{2}}\right] d x_{0} d y_{0} \tag{2} \end{equation*}(2)w=(1+v)p2πE∫−aa∫−bb1R[2(1−v)+zR2]dx0dy0
where
R = ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 . R = x − x 0 2 + y − y 0 2 + z 2 . R=sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)).R=\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}} .R=(x−x0)2+(y−y0)2+z2.
Figure 1: Half-space: (a) loaded on its surface with a concentrated force; (b) subjected to a load on a square domain.

Results

Problem of calculation of the definite integral


onsider a problem of loading of an elastic half-space on a square domain (Fig. 1b). In aim to obtain a function of vertical displacements, one has to calculate a following integral (acc. to Form. (2)):
(3) w = ( 1 + v ) p 2 π E ∫ − a a ∫ − b b 1 R [ 2 ( 1 − v ) + z R 2 ] d x 0 d y 0 (3) w = ( 1 + v ) p 2 π E ∫ − a a   ∫ − b b   1 R 2 ( 1 − v ) + z R 2 d x 0 d y 0 {:(3)w=((1+v)p)/(2pi E)int_(-a)^(a)int_(-b)^(b)(1)/(R)[2(1-v)+(z)/(R^(2))]dx_(0)dy_(0):}\begin{equation*} w=\frac{(1+v) p}{2 \pi E} \int_{-a}^{a} \int_{-b}^{b} \frac{1}{R}\left[2(1-v)+\frac{z}{R^{2}}\right] d x_{0} d y_{0} \tag{3} \end{equation*}(3)w=(1+v)p2πE∫−aa∫−bb1R[2(1−v)+zR2]dx0dy0
If the expression (3) depends on three spatial variables, it is impossible to obtain the analytical results of calculation of this expression in a direct way in the Mathematica environment, using a double integral - the software does not return result. The assumption z = 0 z = 0 z=0z=0z=0 (calculation of the displacements on a plane) neither allows to calculate the definite integral (3). This problem can be overcome in a following way:
a) determine the indefinite integral of (3)
w ^ = ( 1 − v 2 ) p E π ( x − x 0 ) ln ⁡ ( y − y 0 + ( x − x 0 ) 2 + ( y − y 0 ) 2 ) + (4) + ( 1 − v 2 ) p E π [ ( y − y 0 ) ln ⁡ ( x − x 0 + ( x − x 0 ) 2 + ( y − y 0 ) 2 ) + y 0 ] w ^ = 1 − v 2 p E π x − x 0 ln ⁡ y − y 0 + x − x 0 2 + y − y 0 2 + (4) + 1 − v 2 p E π y − y 0 ln ⁡ x − x 0 + x − x 0 2 + y − y 0 2 + y 0 {:[ hat(w)=((1-v^(2))p)/(E pi)(x-x_(0))ln(y-y_(0)+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)))+],[(4)+((1-v^(2))p)/(E pi)[(y-y_(0))ln(x-x_(0)+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)))+y_(0)]]:}\begin{align*} & \hat{w}=\frac{\left(1-v^{2}\right) p}{E \pi}\left(x-x_{0}\right) \ln \left(y-y_{0}+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}}\right)+ \\ & +\frac{\left(1-v^{2}\right) p}{E \pi}\left[\left(y-y_{0}\right) \ln \left(x-x_{0}+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}}\right)+y_{0}\right] \tag{4} \end{align*}w^=(1−v2)pEπ(x−x0)ln⁡(y−y0+(x−x0)2+(y−y0)2)+(4)+(1−v2)pEπ[(y−y0)ln⁡(x−x0+(x−x0)2+(y−y0)2)+y0]
b) substitute the integration limits:
(5) w = w ^ | x 0 = a , y 0 = b − w ^ | x 0 = − a , y 0 = − b (5) w = w ^ x 0 = a , y 0 = b − w ^ x 0 = − a , y 0 = − b {:(5)w=( hat(w))|_(x_(0)=a,y_(0)=b)-( hat(w))|_(x_(0)=-a,y_(0)=-b):}\begin{equation*} w=\left.\hat{w}\right|_{x_{0}=a, y_{0}=b}-\left.\hat{w}\right|_{x_{0}=-a, y_{0}=-b} \tag{5} \end{equation*}(5)w=w^|x0=a,y0=b−w^|x0=−a,y0=−b
(for the sake of simplifications of the formulas presented in this study, it has been assumed z = 0 − z = 0 − z=0-z=0-z=0− it does not affect the possibility of calculation of analytical formulas depending on three spatial variables). As a result, the known vertical displacement field w = w ( x , y , 0 ) w = w ( x , y , 0 ) w=w(x,y,0)w=w(x, y, 0)w=w(x,y,0) is obtained which can be presented in a form:
E π ( 1 − v 2 ) p w = ( x + a ) ln ⁡ ( ( a + x ) 2 + ( b + y ) 2 + y + b ( a + x ) 2 + ( b − y ) 2 + y − b ) + ( x − a ) ln ⁡ ( ( a − x ) 2 + ( b − y ) 2 + y − b ( a − x ) 2 + ( b + y ) 2 + y + b ) + (6) + ( y + b ) ln ⁡ ( ( a − x ) 2 + ( b + y ) 2 + a − x ( a + x ) 2 + ( b + y ) 2 − a − x ) + ( y − b ) ln ⁡ ( ( a + x ) 2 + ( b − y ) 2 − a − x ( a − x ) 2 + ( b − y ) 2 + a − x ) . E π 1 − v 2 p w = ( x + a ) ln ⁡ ( a + x ) 2 + ( b + y ) 2 + y + b ( a + x ) 2 + ( b − y ) 2 + y − b + ( x − a ) ln ⁡ ( a − x ) 2 + ( b − y ) 2 + y − b ( a − x ) 2 + ( b + y ) 2 + y + b + (6) + ( y + b ) ln ⁡ ( a − x ) 2 + ( b + y ) 2 + a − x ( a + x ) 2 + ( b + y ) 2 − a − x + ( y − b ) ln ⁡ ( a + x ) 2 + ( b − y ) 2 − a − x ( a − x ) 2 + ( b − y ) 2 + a − x . {:[(E pi)/((1-v^(2))p)w=(x+a)ln((sqrt((a+x)^(2)+(b+y)^(2))+y+b)/(sqrt((a+x)^(2)+(b-y)^(2))+y-b))+(x-a)ln((sqrt((a-x)^(2)+(b-y)^(2))+y-b)/(sqrt((a-x)^(2)+(b+y)^(2))+y+b))+],[(6)+(y+b)ln((sqrt((a-x)^(2)+(b+y)^(2))+a-x)/(sqrt((a+x)^(2)+(b+y)^(2))-a-x))+(y-b)ln((sqrt((a+x)^(2)+(b-y)^(2))-a-x)/(sqrt((a-x)^(2)+(b-y)^(2))+a-x)).]:}\begin{align*} & \frac{E \pi}{\left(1-v^{2}\right) p} w=(x+a) \ln \left(\frac{\sqrt{(a+x)^{2}+(b+y)^{2}}+y+b}{\sqrt{(a+x)^{2}+(b-y)^{2}}+y-b}\right)+(x-a) \ln \left(\frac{\sqrt{(a-x)^{2}+(b-y)^{2}}+y-b}{\sqrt{(a-x)^{2}+(b+y)^{2}}+y+b}\right)+ \\ & +(y+b) \ln \left(\frac{\sqrt{(a-x)^{2}+(b+y)^{2}}+a-x}{\sqrt{(a+x)^{2}+(b+y)^{2}}-a-x}\right)+(y-b) \ln \left(\frac{\sqrt{(a+x)^{2}+(b-y)^{2}}-a-x}{\sqrt{(a-x)^{2}+(b-y)^{2}}+a-x}\right) . \tag{6} \end{align*}Eπ(1−v2)pw=(x+a)ln⁡((a+x)2+(b+y)2+y+b(a+x)2+(b−y)2+y−b)+(x−a)ln⁡((a−x)2+(b−y)2+y−b(a−x)2+(b+y)2+y+b)+(6)+(y+b)ln⁡((a−x)2+(b+y)2+a−x(a+x)2+(b+y)2−a−x)+(y−b)ln⁡((a+x)2+(b−y)2−a−x(a−x)2+(b−y)2+a−x).
Form. (6), firstly obtained by Schleicher [3], allows to determine vertical displacements in any point of the plane z = 0 z = 0 z=0z=0z=0.

Integral of sums and sum of integrals

The integrand in (3) can be presented as a sum of two functions:
(7) 2 ( 1 − v ) + z 2 ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 = 2 ( 1 − v ) ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 + z 2 ( ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) 3 / 2 (7) 2 ( 1 − v ) + z 2 x − x 0 2 + y − y 0 2 + z 2 x − x 0 2 + y − y 0 2 + z 2 = 2 ( 1 − v ) x − x 0 2 + y − y 0 2 + z 2 + z 2 x − x 0 2 + y − y 0 2 + z 2 3 / 2 {:(7)(2(1-v)+(z^(2))/((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))/(sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))=(2(1-v))/(sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))+(z^(2))/(((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))^(3//2)):}\begin{equation*} \frac{2(1-v)+\frac{z^{2}}{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}}{\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}}=\frac{2(1-v)}{\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}}+\frac{z^{2}}{\left(\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}\right)^{3 / 2}} \tag{7} \end{equation*}(7)2(1−v)+z2(x−x0)2+(y−y0)2+z2(x−x0)2+(y−y0)2+z2=2(1−v)(x−x0)2+(y−y0)2+z2+z2((x−x0)2+(y−y0)2+z2)3/2
(the integrand in (3) has been multiplied by 2 π E / p ( 1 + v ) 2 π E / p ( 1 + v ) 2pi E//p(1+v)2 \pi E / p(1+v)2πE/p(1+v) ). Then, the indefinite integrals with respect to the variable x 0 x 0 x_(0)x_{0}x0 in the left and right sides of (7) have been calculated in the Mathematica software. The integral of the left side of (7) is equal to:
a) in the version 8.0, 11.3 and 12.3:
(8) 2 ( 1 − v ) ln ⁡ ( − ( x − x 0 ) + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) − ( x − x 0 ) z 2 ( ( y − y 0 ) 2 + z 2 ) ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 (8) 2 ( 1 − v ) ln ⁡ − x − x 0 + x − x 0 2 + y − y 0 2 + z 2 − x − x 0 z 2 y − y 0 2 + z 2 x − x 0 2 + y − y 0 2 + z 2 {:(8)2(1-v)ln(-(x-x_(0))+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))-((x-x_(0))z^(2))/(((y-y_(0))^(2)+z^(2))sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))):}\begin{equation*} 2(1-v) \ln \left(-\left(x-x_{0}\right)+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)-\frac{\left(x-x_{0}\right) z^{2}}{\left(\left(y-y_{0}\right)^{2}+z^{2}\right) \sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}} \tag{8} \end{equation*}(8)2(1−v)ln⁡(−(x−x0)+(x−x0)2+(y−y0)2+z2)−(x−x0)z2((y−y0)2+z2)(x−x0)2+(y−y0)2+z2
b) in the version 13.3 and 14.2:
(9) − 2 ( 1 − v ) ln ⁡ ( x − x 0 + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) − ( x − x 0 ) z 2 ( ( y − y 0 ) 2 + z 2 ) ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 (9) − 2 ( 1 − v ) ln ⁡ x − x 0 + x − x 0 2 + y − y 0 2 + z 2 − x − x 0 z 2 y − y 0 2 + z 2 x − x 0 2 + y − y 0 2 + z 2 {:(9)-2(1-v)ln(x-x_(0)+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))-((x-x_(0))z^(2))/(((y-y_(0))^(2)+z^(2))sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))):}\begin{equation*} -2(1-v) \ln \left(x-x_{0}+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)-\frac{\left(x-x_{0}\right) z^{2}}{\left(\left(y-y_{0}\right)^{2}+z^{2}\right) \sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}} \tag{9} \end{equation*}(9)−2(1−v)ln⁡(x−x0+(x−x0)2+(y−y0)2+z2)−(x−x0)z2((y−y0)2+z2)(x−x0)2+(y−y0)2+z2
It is evident that the second components of (8) and (9) are equal to each other, however the first ones are different. This difference results from the fact that
ln ⁡ ( − ( x − x 0 ) + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) ≠ ln ⁡ ( x − x 0 + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) ln ⁡ − x − x 0 + x − x 0 2 + y − y 0 2 + z 2 ≠ ln ⁡ x − x 0 + x − x 0 2 + y − y 0 2 + z 2 ln(-(x-x_(0))+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))!=ln(x-x_(0)+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))\ln \left(-\left(x-x_{0}\right)+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right) \neq \ln \left(x-x_{0}+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)ln⁡(−(x−x0)+(x−x0)2+(y−y0)2+z2)≠ln⁡(x−x0+(x−x0)2+(y−y0)2+z2)
The integral of the right side of (7) is equal:
a) in the version 8.0, 13.3 and 14.2:
(10) − 2 ( 1 − v ) ln ⁡ ( x − x 0 + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) − ( x − x 0 ) z 2 ( ( y − y 0 ) 2 + z 2 ) ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 (10) − 2 ( 1 − v ) ln ⁡ x − x 0 + x − x 0 2 + y − y 0 2 + z 2 − x − x 0 z 2 y − y 0 2 + z 2 x − x 0 2 + y − y 0 2 + z 2 {:(10)-2(1-v)ln(x-x_(0)+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))-((x-x_(0))z^(2))/(((y-y_(0))^(2)+z^(2))sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))):}\begin{equation*} -2(1-v) \ln \left(x-x_{0}+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)-\frac{\left(x-x_{0}\right) z^{2}}{\left(\left(y-y_{0}\right)^{2}+z^{2}\right) \sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}} \tag{10} \end{equation*}(10)−2(1−v)ln⁡(x−x0+(x−x0)2+(y−y0)2+z2)−(x−x0)z2((y−y0)2+z2)(x−x0)2+(y−y0)2+z2
b) in the version 11.3 and 12.3:
(11) − 2 ( 1 − v ) arctanh ( x − x 0 ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) − ( x − x 0 ) z 2 ( ( y − y 0 ) 2 + z 2 ) ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 (11) − 2 ( 1 − v ) arctanh x − x 0 x − x 0 2 + y − y 0 2 + z 2 − x − x 0 z 2 y − y 0 2 + z 2 x − x 0 2 + y − y 0 2 + z 2 {:(11)-2(1-v)arctanh((x-x_(0))/(sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))))-((x-x_(0))z^(2))/(((y-y_(0))^(2)+z^(2))sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))):}\begin{equation*} -2(1-v) \operatorname{arctanh}\left(\frac{x-x_{0}}{\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}}\right)-\frac{\left(x-x_{0}\right) z^{2}}{\left(\left(y-y_{0}\right)^{2}+z^{2}\right) \sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}} \tag{11} \end{equation*}(11)−2(1−v)arctanh(x−x0(x−x0)2+(y−y0)2+z2)−(x−x0)z2((y−y0)2+z2)(x−x0)2+(y−y0)2+z2
It is evident that the second components of (10) and (11) are equal to each other, however the first ones are different. Fig. 2a shows the difference between the formulas (8) and (9) and Fig. 2b - between the formulas (10) and (11), both for x 0 = 0 x 0 = 0 x_(0)=0x_{0}=0x0=0, y 0 = 0 , z = 0 , v = 0 y 0 = 0 , z = 0 , v = 0 y_(0)=0,z=0,v=0y_{0}=0, z=0, v=0y0=0,z=0,v=0. The Figures show that the differences mostly are negligible, however for y → 0 y → 0 y rarr0y \rightarrow 0y→0 they are significant and the worst situation is for x → 0 x → 0 x rarr0x \rightarrow 0x→0 and y → 0 y → 0 y rarr0y \rightarrow 0y→0, where the differences go to infinity. It means that, for an assumed area of action of the distribute load, calculations of the displacement may show the highest error for points close to the point x = 0 , y = 0 x = 0 , y = 0 x=0,y=0x=0, y=0x=0,y=0 and it is advisable to select the origin of the coordination system in a certain distance from this area.
Figure 2: Difference between: a) Forms. (8) and (9); b) Forms. (10) and (11).
A related issue is checking if the integral of a sum of integrands is equal to the sum of integrals of individual integrands:
(12) ∫ [ f ( x ) + g ( x ) ] d x = ∫ f ( x ) d x + ∫ g ( x ) d x (12) ∫ [ f ( x ) + g ( x ) ] d x = ∫ f ( x ) d x + ∫ g ( x ) d x {:(12)int[f(x)+g(x)]dx=int f(x)dx+int g(x)dx:}\begin{equation*} \int[f(x)+g(x)] d x=\int f(x) d x+\int g(x) d x \tag{12} \end{equation*}(12)∫[f(x)+g(x)]dx=∫f(x)dx+∫g(x)dx
i.e. if the left and right sides of (7) are equal to each other, despite the differences existing within each side for the individual versions of the Mathematica software. Theorem (12) is fulfilled only for the versions 13.3 and 14.2 of the software (the proof of this statement is also the fact that Forms. (9) and (10), concerning these versions, are equal to each other). However, it is not the case for the versions 8.0, 11.3 and 12.3. A residuum arises which can be written in a form:
a) in the version 8.0 :
(13) 2 ( 1 − v ) ( ln ⁡ ( x − x 0 + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) + ln ⁡ ( − ( x − x 0 ) + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) ) (13) 2 ( 1 − v ) ln ⁡ x − x 0 + x − x 0 2 + y − y 0 2 + z 2 + ln ⁡ − x − x 0 + x − x 0 2 + y − y 0 2 + z 2 {:(13)2(1-v)(ln(x-x_(0)+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))+ln(-(x-x_(0))+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))):}\begin{equation*} 2(1-v)\left(\ln \left(x-x_{0}+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)+\ln \left(-\left(x-x_{0}\right)+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)\right) \tag{13} \end{equation*}(13)2(1−v)(ln⁡(x−x0+(x−x0)2+(y−y0)2+z2)+ln⁡(−(x−x0)+(x−x0)2+(y−y0)2+z2))
b) in the version 11.3 and 12.3:
(14) 2 ( 1 − v ) ( arctanh ( x − x 0 ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) + ln ⁡ ( − ( x − x 0 ) + ( x − x 0 ) 2 + ( y − y 0 ) 2 + z 2 ) ) (14) 2 ( 1 − v ) arctanh x − x 0 x − x 0 2 + y − y 0 2 + z 2 + ln ⁡ − x − x 0 + x − x 0 2 + y − y 0 2 + z 2 {:(14)2(1-v)(arctanh((x-x_(0))/(sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2))))+ln(-(x-x_(0))+sqrt((x-x_(0))^(2)+(y-y_(0))^(2)+z^(2)))):}\begin{equation*} 2(1-v)\left(\operatorname{arctanh}\left(\frac{x-x_{0}}{\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}}\right)+\ln \left(-\left(x-x_{0}\right)+\sqrt{\left(x-x_{0}\right)^{2}+\left(y-y_{0}\right)^{2}+z^{2}}\right)\right) \tag{14} \end{equation*}(14)2(1−v)(arctanh(x−x0(x−x0)2+(y−y0)2+z2)+ln⁡(−(x−x0)+(x−x0)2+(y−y0)2+z2))
Fig. 3 shows the residua (13) and (14) with an assumption that x 0 = 0 , y 0 = 0 , z = 0 , v = 0 x 0 = 0 , y 0 = 0 , z = 0 , v = 0 x_(0)=0,y_(0)=0,z=0,v=0x_{0}=0, y_{0}=0, z=0, v=0x0=0,y0=0,z=0,v=0.
Figure 3: Difference between the integral of sums and sum of integrals resulting from: (a) (13); (b) (14).
Functions arctanh ( x ) arctanh ( x ) arctanh(x)\operatorname{arctanh}(x)arctanh(x) and ln ⁡ ( x ) ln ⁡ ( x ) ln(x)\ln (x)ln⁡(x) in (13) and (14) are related by the following dependence (e.g. [25]): arctanh ( x ) == 1 / 2 ln ⁡ ( ( 1 + x ) / ( 1 − x ) ) arctanh ( x ) == 1 / 2 ln ⁡ ( ( 1 + x ) / ( 1 − x ) ) arctanh(x)==1//2ln((1+x)//(1-x))\operatorname{arctanh}(x)= =1 / 2 \ln ((1+x) /(1-x))arctanh(x)==1/2ln⁡((1+x)/(1−x)). Unfortunately, performing such transformation also does not allow to obtain an equality of the left and right sides of (7).
It can be concluded that one can obtain various functions, not always being equal to each other, depending on the software version used.

Discussion

The results of the performed analysis suggest possible difficulties in the analytical integration of the Green functions in the Wolfram Mathematica environment and can vary depending on the software version used.
The problem of obtaining the explicit version of the results for definite integrals can be overcome by the calculation of an indefinite integral as first and then a substitution of integration limits. This way one can obtain the known displacement field of a half-space subjected to a load on a square domain (e.g. [3], [26]).
During the integration of the Green functions (3), differences between the results obtained in different software versions. In particular, for the versions 8.0, 11.3 and 12.3, the sum of integrals and the integral of sums were not concordant, whereas in the versions 13.3 and 14.2 the basic feature (12) of indefinite integrals was fulfilled.
These discrepancies are most likely caused by modifications in the symbolic integration algorithms introduced in different software versions of the Wolfram Mathematica, particularly in the handling of logarithmic and inverse hyperbolic functions. To confirm this hypothesis, further studies are required to investigate the exact causes of these inconsistencies. Based on the conducted analysis, it is recommended that users of symbolic computation in Wolfram Mathematica:
  • use the most recent versions of the software (13 or newer), which ensure consistency of symbolic integration results and satisfy the identity (12) for the tested functions (Forms. (9) and (10) are equal to each other);
  • verify symbolic results by performing an additional check of the equality between the integral of a sum and the sum of integrals, as this simple test can reveal potential discrepancies in older versions of the software;
  • avoid relying solely on analytical outputs from older versions (such as 8.0 , 11.3 8.0 , 11.3 8.0,11.38.0,11.38.0,11.3, or 12.3 ) without cross-checking, because they may produce mathematically non-equivalent forms of the solution.

Conclusions and resumé

Difficulties related to the analytical integration of the Green functions in the Wolfram Mathematica environment have been analyzed in this study. Discrepancies have been shown between the results obtained in various software versions. The results can be important for practical applications, especially in the cases when integrands are of a complex character. If the integration is performed for other domains than rectangular (e.g. triangular), then the integrands can apply a complex form (after integrating with respect to one of variables and substituting the limits). In such situations, it can arise a need for presenting them as a sum of simpler functions. However, one must draw attention on a version of software being applied because, as it has been shown above, in the versions 8.0 , 11.3 8.0 , 11.3 8.0,11.38.0,11.38.0,11.3 and 12.3 a sum of integrals is not equal to an integral of sums what can lead to incorrect results or their incorrect interpretation.
  • obtaining an explicit form of results for definite integrals,
  • discordance of results obtained in various software versions,
  • differences between a sum of integrals and an integral of sums.
The results of this work can help in seeking benchmarks enabling a verification of solutions obtained with numerical methods, e.g. FEM.
Future research should focus on extending the analysis to other Green functions and more complex load distributions (e.g., over triangular or circular domains). It is also necessary to develop verification procedures for symbolic computations and to investigate the precise mathematical and algorithmic causes of the observed discrepancies in integration results between different versions of computational software

References

[1] Hahn, H.G. (1985). Elastizitätstheorie: Grundlagen der linearen Theorie und Anwendungen auf eindimensionale, ebene und räumliche Probleme. Stuttgart, Vieweg+Teubner Verlag (in German).
[2] Boussinesq, J. (1885). Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques: principalement au calcul des déformations et des pressions que produisent, dans ces solides, des efforts quelconques exercés sur une petite partie de leur surface ou de leur intérieur: mémoire suivi de notes étendues sur divers points de physique, mathematique et d'analyse. Paris, Gauthier-Villars (in French).
[3] Schleicher, F. (1926). Zur Theorie des Baugrundes. Der Bauingenieur, 48, pp. 931-935 (in German).
[4] Love, A.E.H. (1929). The Stress Produced in a Semi-infinite Solid by Pressure on Part of the Boundary. Philos. Trans. R. Soc. Lond. 28, pp. 377-420.
[5] Cerruti, V. (1882). Ricerche intorno all'equilibrio de'corpi elastici isotropi. Memoria del Valentino Cerruti. Roma, Salviucci (in Italian).
[6] Terazawa, K. (1916). On the elastic equilibrium of a semi-infinite solid under given boundary conditions, with some applications. Journal of the College of Science, 37
[7] Becker, J.M. and Bevis, M. (2004). Love's problem. Geophys. J. Int. 156, pp. 171-178, https://doi.org/10.1111/j.1365246X.2003.02150.x.
[8] Dydo, J.R. and Busby, H.R. (1995). Elasticity solutions for constant and linearly varying loads applied to a rectangular surface patch on the elastic half-space. J Elast 38, pp. 153-163, DOI: 10.1007/BF00042496.
[9] Jin, X., Niu, F., Zhang, X., Zhou, Q., Lyu, D., Keer, L. M. and Hu, Y. (2016). Love's rectangular contact problem revisited: A complete solution. Tribology International, 103, pp. 331-342, DOI: 10.1016/j.triboint.2016.07.011.
[10] D'Urso, M.G. and Marmo, F. (2013). On a generalized Love's problem. Comput Geosci 61, pp. 144-151, DOI: 10.1016/j.cageo.2013.09.002.
[11] D'Urso, M.G. and Marmo, F. (2015). Vertical stress distribution in isotropic half-spaces due to surface vertical loadings acting over polygonal domains. Z. Angew. Math. Mech. 95, pp. 91-110, DOI: 10.1002/zamm. 201300034.
[12] Marmo, F. and Rosati, L. (2016). A General Approach to the Solution of Boussinesq's Problem for Polynomial Pressures Acting over Polygonal Domains. J. Elast. 122, pp. 75-112, DOI: 10.1007/s10659-015-9534-5.
[13] Marmo, F., Toraldo, F. and Rosati, L. (2016). Analytical formulas and design charts for transversely isotropic half-spaces subject to linearly distributed pressures. Meccanica 51, pp. 2909-2928, DOI: 10.1007/s11012-016-0443-x.
[14] Marmo, F., Toraldo, F. and Rosati, L. (2017). Transversely isotropic half-spaces subject to surface pressures. Int J Solids Struct 104-105, pp. 35-49, DOI: 10.1016/j.ijsolstr.2016.11.001.
[15] Zhang, Z., Liu, S., Pan, E. and Wang, Q. (2023). Dynamic loading in a transversely isotropic and layered elastic halfspace. Int J Mech Sci 260, 108626, https://doi.org/10.1016/j.ijmecsci.2023.108626.
[16] Wolfram Research, Inc.: Mathematica v8.0 [software], Champaign, IL, 2010.
[17] Wolfram Research, Inc.: Mathematica v11.3 [software], Champaign, IL, 2018.
[18] Wolfram Research, Inc.: Mathematica v12.3 [software], Champaign, IL, 2019.
[19] Wolfram Research, Inc.: Mathematica v13.3 [software], Champaign, IL, 2023.
[20] Wolfram Research, Inc.: Mathematica v14.2 [software], Champaign, IL, 2024.
[21] Nowacki, W. (1970). Teoria sprężystości. Warszawa, PWN (in Polish).
[22] Jemioło, S., Szwed, A. (2017). Zagadnienia statyki sprężystych półprzestrzeni warstwowych. Warszawa, Oficyna Wydawnicza PW (in Polish).
[23] Timoshenko, S. and Goodier, J. (1951). Theory of Elasticity. New York, McGraw-Hill Book Company.
[24] Davis, R. O. and Selvadurai, A. P. S. (1996). Elasticity and Geomechanics. Cambridge, Cambridge University Press.
[25] Ryżyk, I., Gradsztejn, I. (1964). Tablice całek, sum, szeregów i iloczynów. Warszawa, PWN (in Polish).
[26] Florin, V.A. (1959) Fundamentals of Soil Mechanics. 2. Leningrad, State Publishing House of Literature on Construction, Architecture and Building Materials.