# Numerical Evaluation of Arbitrary Singular Domain Integrals Using Third-Degree B-Spline Basis Functions

Authors: Jin-Xiu Hu; Hai-Feng Peng; Xiao-Wei Gao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A new approach is presented for the numerical evaluation of arbitrary singular domain integrals. In this method, singular domain integrals are transformed into a boundary integral and a radial integral which contains singularities by using the radial integration method. The analytical elimination of singularities condensed in the radial integral formulas can be accomplished by expressing the nonsingular part of the integration kernels as a series of cubic B-spline basis functions of the distancer and using the intrinsic features of the radial integral. In the proposed method, singularities involved in the domain integrals are explicitly transformed to the boundary integrals, so no singularities exist at internal points. A few numerical examples are provided to verify the correctness and robustness of the presented method.


---

## Body

## 1. Introduction

When the boundary element method (BEM) [1, 2] is used to solve complicated engineering problems [3, 4], such as compressible potential flows, heat conduction with heat generation, seepage problems with sources, electromagnetic field problems with electric charge, and elastostatics with body forces, a large number of regular and/or singular domain integrals may appear in the basic integral equations [2]. The accurate evaluation of general regular and/or singular domain integrals is very important and is still a crucial area of research in the BEM [5–11].To compute the domain integrals, the frequently used approach is to apply some special schemes to transform domain integrals into boundary integrals, thus avoiding the discretization of the internal domain. The most extensively used technique is the dual reciprocity method (DRM) [12], which transforms the domain integrals to the boundary by approximating the given body force effect quantities as a series of prescribed basis functions and then employing particular solutions corresponding to these basis functions. As an extension of the idea of DRM, the multiple reciprocity method (MRM) [13, 14] was presented by Nowak and Brebbia [13] and the triple-reciprocity method was developed by Ochiai and Sekiya [15]. The radial integration method (RIM) presented by Gao [16–18] is another powerful method for transforming domain integrals into boundary integrals based on pure mathematical treatments, without using the Laplace operator and particular solutions of the problem.Recently, an efficient approach was presented by Gao and Peng [7] for evaluating arbitrarily high-order singular domain integrals in a unified way based on the radial integration method. In the method, the radial integral can be integrated analytically by expressing the nonsingular part of integrand as polynomials of the global distance r. Therefore, this approach can deal with any type of regular, weakly, and high-order singular domain integrals. However, the convergence performance of the computational result is not very good and is volatility with respect to the order of polynomials for some complicated function over complicated geometry due to the characteristics of high-order polynomial.In this paper, the nonsingular part of the integrand involved in the domain integral is expressed as a series of third-degree B-spline basis functions [19] of the distance r. This approach has two advantages. One is that the coefficients calculation becomes much simpler than that in [7] because of the local support and the endpoint interpolatory properties of B-spline basis functions with open knot vectors [20]. The other is that the integration results are very stable, less dependent on the number of order of polynomials in [7]. A number of singular domain integrals are given to verify the correctness and validity of the presented method.


## 2. Singular Domain Integrals

In this paper, the following singular domain integral is considered:(1)I(xp)=∫Ωf¯(xp,x)rβ(xp,x)dΩ(x),
where f¯(xp,x) is bounded everywhere, being the nonsingular part of the integrand, Ω is the domain of the problem to be considered, and r(xp,x) is the distance between the source point xp and the field point x; that is,
(2)r(xp,x)=x-xp=(x-xp)·(x-xp).The level of the singularity of domain integrals in (1) depends on the value of β. For a finite integration domain Ω, the regular and weakly singular domain integrals [7] always exist (i.e., results are finite). However, for strongly, hyper- and supersingular domain integrals [7], integrals only exist under some conditions, depending heavily on the characteristics of the function f¯(xp,x) [18]. Usually, integrals stemming from real physical problems should exist. In this paper, the singular integrals are evaluated in the Cauchy principal value sense; that is,
(3)I(xp)=∫Ωf¯(xp,x)rβ(xp,x)dΩ(x)=lim⁡ε→0∫Ω-Ωεf¯(xp,x)rβ(xp,x)dΩ(x),
where Ωε is a small circular domain for 2D problems or a small spherical domain for 3D problems with the radius ε around the source point xp.


## 3. The Radial Integration Method [19, 20]

According to the radial integration method, the singular domain integral represented by (1) can be transformed into the equivalent boundary integral as follows:
(4)I(xp)=∫Ωf¯(xp,x)rβ(xp,x)dΩ(x)=∫ΓF(xp,xq)rD-1(xp,xq)∂r∂ndΓ(xq),
where
(5)F(xp,xq)=∫0r(xp,xq)f¯(xp,x)rβ(xp,x)rD-1dr.
In (4) and (5), D=2 for two-dimensional (2D) problems and D=3 for three-dimensional (3D) problems; xq denotes the coordinates at the boundary point q; r(xp,xq) is the distance from the source point xp to the boundary point xq; Γ is the boundary of the integration domain Ω; and n is the unit outward normal vector to the boundary Γ.In order to evaluate the radial integral (5), the coordinate x in the integrand needs to be expressed in terms of the integration variable r; that is,
(6)xi=xip+r,ir,
where xi represents the ith component of Cartesian coordinates at the field point x and r,i is defined as follows:
(7)r,i=∂r∂xiq=xiq-xipr(xp,xq).From (4) and (5), we can see that the key task for the transformation of singular domain integral to the boundary integral is to accurately evaluate the radial integral (5). For simple functions, the radial integral (5) can be integrated analytically. However, for some complicated functions, analytically evaluating the radial integral is not available. In particular, for high-order singular integrals, the numerical evaluation of radial integrals may have infinite value problems. In this paper, a new approach will be presented for evaluating the radial integral (5) analytically in a unified way. This will be described in detail in the next section.


## 4. Evaluation of Singular Domain Integrals

In the section, the analytical evaluation of radial integral (5) is carried out by expressing the nonsingular part f¯ as a series of third-degree B-spline basis functions [19] of the distance r.


### 4.1. Introduction of B-Spline Basis Functions

LetU=u0,…,um be a nondecreasing sequence of real numbers; that is, ui≤ui+1, i=0,1,…,m-1. ui are called knots, and U is the knot vector. The ith B-spline basis function of p-degree, denoted by Ni,pu, is defined as [19, 21] (8a)Ni,0=1ifui≤u<ui+1,0otherwise,(8b)Ni,pu=u-uiui+p-uiNi,p-1u+ui+p+1-uui+p+1-ui+1Ni+1,p-1u,where i=0,1,…,n;  m=n+p+1. If knots are equally spaced in the parametric space, they are said to be uniform [19, 20]. A knot vector is said to be open if its first and last knots appear p+1 times. An initial example of the results of applying (8a) and (8b) to a uniform knot vector U={0,1,2,3,4,5,6} [22] is presented in Figure 1.B-spine basis functions of order 0, 1, 2, 3 for uniform knot vectorU={0,1,2,3,4,5,6}.
(a)
(b)
(c)
(d)Important properties of B-spline basis functions are as follows [19, 20].(1)
They constitute a partition of unity; that is,∀u, ∑i=0nNi,p(u)=1.(2)
The support of eachNi,p is compact and contained in the interval [ui,ui+p+1].(3)
Each basis function is nonnegative; that is,Ni,p(u)≥0,  ∀u.(4)
For an open, nonuniform knot vector, the basis functions are interpolatory at the ends of the interval.In Sections4.2 and 4.3, it will be seen that properties (2) and (3) of B-spline basis functions can bring simplicity to our calculation.


### 4.2. Express Nonsingular Part of the Integration Kernel as a Series of Third-Degree B-Spline Basis Functions

It is assumed that the radial integration equations (4) and (5) are in the Cauchy principle value sense, such that (4) and (5) are written as
(9)Ixp=lim⁡ε→0∫Ω-Ωεf¯xp,xrβxp,xdΩx=∫ΓF(xp,xq)rD-1∂r∂ndΓ(xq),
where
(10)Fxp,xq=lim⁡ε→0∫εrxp,xqf¯xp,xrβrD-1dr=lim⁡ε→0∫εr(xp,xq)f¯(xp,x)rβ-D+1dr.In order to analytically evaluate the singular radial integral (10), the nonsingular part of the integrand is expressed as a series of third-degree B-spline basis functions of the distance r; that is,
(11)f¯(xp,x)≈∑n=0LBnNn,3(r),
where Nn,3(r) are the third-degree B-spline basis functions and Bn are the coefficients to be determined. Once knot vector U={r0,r1,…,rL+3+1} is defined, the third-degree B-spline basis functions Nn,3(r) can be determined. In this paper, the following thirteen types of knot vectors are used, corresponding to L=3,4,…,15, respectively;
(12)Ui=0,0,0,0,r^i,2*r^i,…,i-1*r^i,r^,r^,r^,r^,i=1,2,…,13,
where r^=r(xp,xq).The formulation of third-degree B-spline basis functionsNn,3(r) can be deduced from (8a) and (8b) as follows:
(13)Nn,3r=r-rn3(rn+3-rn)(rn+2-rn)(rn+1-rn)≜Mn1,r∈rn,rn+1,r-rn2(rn+2-r)(rn+3-rn)(rn+2-rn)(rn+2-rn+1)+(r-rn)(rn+3-r)(r-rn+1)(rn+3-rn)(rn+2-rn+1)(rn+3-rn+1)+(rn+4-r)r-rn+12(rn+4-rn+1)(rn+3-rn+1)(rn+2-rn+1)≜Mn2,r∈rn+1,rn+2,(r-rn)(rn+3-r)2(rn+3-rn)(rn+3-rn+2)(rn+3-rn+1)+(rn+3-r)(rn+4-r)(r-rn+1)(rn+3-rn+2)(rn+4-rn+1)(rn+3-rn+1)+rn+4-r2(r-rn+2)(rn+4-rn+1)(rn+3-rn+2)(rn+4-rn+2)≜Mn3,r∈rn+2,rn+3,rn+4-r3(rn+4-rn+1)(rn+4-rn+2)(rn+4-rn+3)≜Mn4,r∈rn+3,rn+4.It should be noted that the open knot vectors listed in (12) are C2-continuous [19, 20] except at the endpoints of the interval, which makes the calculation more accurate. Moreover, the basis functions formed from open knot vectors are interpolatory at the ends of the parametric space interval [0,r(xp,xq)]; that is,(14a)B0=f¯(xp,xp),(14b)BL=f¯(xp,xq).The other coefficientsBn are determined by collocating L+1 distances (0,r1′,r2′,…,rL′) over the integration region [0,r(xp,xq)], in which
(15)rn′=rxp,xn=nrxp,xqL(n=0,1,…,L),
where x0=xp, xL=xq.So the other coefficients can be solved using(16)R(L-1)B=Y,
where R(L-1) is a square matrix of order L-1; that is,
(17)RL-1=N1,3r1′,N2,3r1′,…,NL-1,3r1′N1,3r2′,N2,3r2′,…,NL-1,3r2′⋮⋮⋮N1,3rL-1′,N2,3rL-1′,…,NL-1,3rL-1′
and B and Y are the following vectors:
(18)B=B1B2⋮BL-1,Y=f¯(xp,x1)-B0N0,3(r1′)-BLNL,3(r1′)f¯(xp,x2)-B0N0,3(r2′)-BLNL,3(r2′)⋮f¯(xp,xL-1)-B0N0,3(rL-1′)-BLNL,3(rL-1′).Due to the local support property [19, 20] of B-spline basis functions, the coefficient matrix R(L-1) is banded and (16) is easy to be solved.


### 4.3. Evaluation of Singular Domain Integrals

Once the coefficientsBn are obtained from the above equations, substituting (11) into (10) yields
(19)Fxp,xq=∑n=0LBnlim⁡ε→0∫εrxp,xqNn,3rrβ-D+1dr=∑n=0LBnQn,
where
(20)Qn=lim⁡ε→0∫εr(xp,xq)Nn,3(r)rβ-D+1dr.Substituting (13) into (20) yields
(21)Qn=lim⁡ε→0∫rnrn+1Mn1rβ-D+1dr+∫rn+1rn+2Mn2rβ-D+1dr+∫rn+2rn+3Mn3rβ-D+1dr+∫rn+3rn+4Mn4rβ-D+1dr=∑j=14lim⁡ε→0∫rn+j-1rn+jMnjrβ-D+1dr≜∑j=14Tnj.
From (13), it can be seen that Mnj is cubic polynomial functions of independent variable r. Therefore, Mnj can be expressed as
(22)Mnj=∑k=03Anjkrk,j=1,2,3,4,n=0,1,…,L,
where Anjk can be determined by the polynomial arithmetic of (13).Accordingly,Tnj appearing in (21) can be expressed as
(23)Tnj=∑k=03Anjklim⁡ε→0∫rn+j-1rn+jrk-β+D-1dr=∑k=03AnjkEnjk,
where
(24)Enjk=0,forrn+j-1=rn+j;rn+jk-β+D-rn+j-1k-β+Dk-β+D,forrn+j-1≠rn+j,rn+j-1≠ε;1k-β+D1rn+jβ-k-D-lim⁡ε→01εβ-k-D,forrn+j-1≠rn+j,rn+j-1=ε,k≠β-D;ln⁡rn+j-lim⁡ε→0ln⁡ε,forrn+j-1≠rn+j,rn+j-1=ε,k=β-D.For a physical problem, the infinite terms in (24) should be cancelled out by free terms. That is, if the integral exists, after considering the contributions of all elements around points with the same size of ε, the last term in (24) must be zero [2, 18]. Based on this assumption, the final integral results can be expressed as
(25)Enjk=0,forrn+j-1=rn+j;rn+jk-β+D-rn+j-1k-β+Dk-β+D,forrn+j-1≠rn+j,rn+j-1≠ε;rn+jk-β+Dk-β+D,forrn+j-1≠rn+j,rn+j-1=ε,k≠β-D;ln⁡rn+j,forrn+j-1≠rn+j,rn+j-1=ε,k=β-D.Substituting (25) into (23), the result into (21), and finally applying (19) yield
(26)Fxp,xq=∑n=0LBn∑j=14∑k=03AnjkEnjk.Substituting (26) into (9), it follows that
(27)I(xp)=∫Γ∑n=0LBn∑j=14∑k=03AnjkEnjkrD-1∂r∂ndΓ(xq).Now the singular domain integral has been transformed into boundary integral. As a result, no singularity exists for internal pointsxp and the boundary integral in (27) can be evaluated using Gaussian quadrature [2, 18]. For each Gauss point used in the Gaussian quadrature, the coefficient Bn and quantity Enjk are computed using (16) and (25), respectively. However, for boundary points, singularity still exists and more treatment should be performed (see [22, 23]).


## 4.1. Introduction of B-Spline Basis Functions

LetU=u0,…,um be a nondecreasing sequence of real numbers; that is, ui≤ui+1, i=0,1,…,m-1. ui are called knots, and U is the knot vector. The ith B-spline basis function of p-degree, denoted by Ni,pu, is defined as [19, 21] (8a)Ni,0=1ifui≤u<ui+1,0otherwise,(8b)Ni,pu=u-uiui+p-uiNi,p-1u+ui+p+1-uui+p+1-ui+1Ni+1,p-1u,where i=0,1,…,n;  m=n+p+1. If knots are equally spaced in the parametric space, they are said to be uniform [19, 20]. A knot vector is said to be open if its first and last knots appear p+1 times. An initial example of the results of applying (8a) and (8b) to a uniform knot vector U={0,1,2,3,4,5,6} [22] is presented in Figure 1.B-spine basis functions of order 0, 1, 2, 3 for uniform knot vectorU={0,1,2,3,4,5,6}.
(a)
(b)
(c)
(d)Important properties of B-spline basis functions are as follows [19, 20].(1)
They constitute a partition of unity; that is,∀u, ∑i=0nNi,p(u)=1.(2)
The support of eachNi,p is compact and contained in the interval [ui,ui+p+1].(3)
Each basis function is nonnegative; that is,Ni,p(u)≥0,  ∀u.(4)
For an open, nonuniform knot vector, the basis functions are interpolatory at the ends of the interval.In Sections4.2 and 4.3, it will be seen that properties (2) and (3) of B-spline basis functions can bring simplicity to our calculation.


## 4.2. Express Nonsingular Part of the Integration Kernel as a Series of Third-Degree B-Spline Basis Functions

It is assumed that the radial integration equations (4) and (5) are in the Cauchy principle value sense, such that (4) and (5) are written as
(9)Ixp=lim⁡ε→0∫Ω-Ωεf¯xp,xrβxp,xdΩx=∫ΓF(xp,xq)rD-1∂r∂ndΓ(xq),
where
(10)Fxp,xq=lim⁡ε→0∫εrxp,xqf¯xp,xrβrD-1dr=lim⁡ε→0∫εr(xp,xq)f¯(xp,x)rβ-D+1dr.In order to analytically evaluate the singular radial integral (10), the nonsingular part of the integrand is expressed as a series of third-degree B-spline basis functions of the distance r; that is,
(11)f¯(xp,x)≈∑n=0LBnNn,3(r),
where Nn,3(r) are the third-degree B-spline basis functions and Bn are the coefficients to be determined. Once knot vector U={r0,r1,…,rL+3+1} is defined, the third-degree B-spline basis functions Nn,3(r) can be determined. In this paper, the following thirteen types of knot vectors are used, corresponding to L=3,4,…,15, respectively;
(12)Ui=0,0,0,0,r^i,2*r^i,…,i-1*r^i,r^,r^,r^,r^,i=1,2,…,13,
where r^=r(xp,xq).The formulation of third-degree B-spline basis functionsNn,3(r) can be deduced from (8a) and (8b) as follows:
(13)Nn,3r=r-rn3(rn+3-rn)(rn+2-rn)(rn+1-rn)≜Mn1,r∈rn,rn+1,r-rn2(rn+2-r)(rn+3-rn)(rn+2-rn)(rn+2-rn+1)+(r-rn)(rn+3-r)(r-rn+1)(rn+3-rn)(rn+2-rn+1)(rn+3-rn+1)+(rn+4-r)r-rn+12(rn+4-rn+1)(rn+3-rn+1)(rn+2-rn+1)≜Mn2,r∈rn+1,rn+2,(r-rn)(rn+3-r)2(rn+3-rn)(rn+3-rn+2)(rn+3-rn+1)+(rn+3-r)(rn+4-r)(r-rn+1)(rn+3-rn+2)(rn+4-rn+1)(rn+3-rn+1)+rn+4-r2(r-rn+2)(rn+4-rn+1)(rn+3-rn+2)(rn+4-rn+2)≜Mn3,r∈rn+2,rn+3,rn+4-r3(rn+4-rn+1)(rn+4-rn+2)(rn+4-rn+3)≜Mn4,r∈rn+3,rn+4.It should be noted that the open knot vectors listed in (12) are C2-continuous [19, 20] except at the endpoints of the interval, which makes the calculation more accurate. Moreover, the basis functions formed from open knot vectors are interpolatory at the ends of the parametric space interval [0,r(xp,xq)]; that is,(14a)B0=f¯(xp,xp),(14b)BL=f¯(xp,xq).The other coefficientsBn are determined by collocating L+1 distances (0,r1′,r2′,…,rL′) over the integration region [0,r(xp,xq)], in which
(15)rn′=rxp,xn=nrxp,xqL(n=0,1,…,L),
where x0=xp, xL=xq.So the other coefficients can be solved using(16)R(L-1)B=Y,
where R(L-1) is a square matrix of order L-1; that is,
(17)RL-1=N1,3r1′,N2,3r1′,…,NL-1,3r1′N1,3r2′,N2,3r2′,…,NL-1,3r2′⋮⋮⋮N1,3rL-1′,N2,3rL-1′,…,NL-1,3rL-1′
and B and Y are the following vectors:
(18)B=B1B2⋮BL-1,Y=f¯(xp,x1)-B0N0,3(r1′)-BLNL,3(r1′)f¯(xp,x2)-B0N0,3(r2′)-BLNL,3(r2′)⋮f¯(xp,xL-1)-B0N0,3(rL-1′)-BLNL,3(rL-1′).Due to the local support property [19, 20] of B-spline basis functions, the coefficient matrix R(L-1) is banded and (16) is easy to be solved.


## 4.3. Evaluation of Singular Domain Integrals

Once the coefficientsBn are obtained from the above equations, substituting (11) into (10) yields
(19)Fxp,xq=∑n=0LBnlim⁡ε→0∫εrxp,xqNn,3rrβ-D+1dr=∑n=0LBnQn,
where
(20)Qn=lim⁡ε→0∫εr(xp,xq)Nn,3(r)rβ-D+1dr.Substituting (13) into (20) yields
(21)Qn=lim⁡ε→0∫rnrn+1Mn1rβ-D+1dr+∫rn+1rn+2Mn2rβ-D+1dr+∫rn+2rn+3Mn3rβ-D+1dr+∫rn+3rn+4Mn4rβ-D+1dr=∑j=14lim⁡ε→0∫rn+j-1rn+jMnjrβ-D+1dr≜∑j=14Tnj.
From (13), it can be seen that Mnj is cubic polynomial functions of independent variable r. Therefore, Mnj can be expressed as
(22)Mnj=∑k=03Anjkrk,j=1,2,3,4,n=0,1,…,L,
where Anjk can be determined by the polynomial arithmetic of (13).Accordingly,Tnj appearing in (21) can be expressed as
(23)Tnj=∑k=03Anjklim⁡ε→0∫rn+j-1rn+jrk-β+D-1dr=∑k=03AnjkEnjk,
where
(24)Enjk=0,forrn+j-1=rn+j;rn+jk-β+D-rn+j-1k-β+Dk-β+D,forrn+j-1≠rn+j,rn+j-1≠ε;1k-β+D1rn+jβ-k-D-lim⁡ε→01εβ-k-D,forrn+j-1≠rn+j,rn+j-1=ε,k≠β-D;ln⁡rn+j-lim⁡ε→0ln⁡ε,forrn+j-1≠rn+j,rn+j-1=ε,k=β-D.For a physical problem, the infinite terms in (24) should be cancelled out by free terms. That is, if the integral exists, after considering the contributions of all elements around points with the same size of ε, the last term in (24) must be zero [2, 18]. Based on this assumption, the final integral results can be expressed as
(25)Enjk=0,forrn+j-1=rn+j;rn+jk-β+D-rn+j-1k-β+Dk-β+D,forrn+j-1≠rn+j,rn+j-1≠ε;rn+jk-β+Dk-β+D,forrn+j-1≠rn+j,rn+j-1=ε,k≠β-D;ln⁡rn+j,forrn+j-1≠rn+j,rn+j-1=ε,k=β-D.Substituting (25) into (23), the result into (21), and finally applying (19) yield
(26)Fxp,xq=∑n=0LBn∑j=14∑k=03AnjkEnjk.Substituting (26) into (9), it follows that
(27)I(xp)=∫Γ∑n=0LBn∑j=14∑k=03AnjkEnjkrD-1∂r∂ndΓ(xq).Now the singular domain integral has been transformed into boundary integral. As a result, no singularity exists for internal pointsxp and the boundary integral in (27) can be evaluated using Gaussian quadrature [2, 18]. For each Gauss point used in the Gaussian quadrature, the coefficient Bn and quantity Enjk are computed using (16) and (25), respectively. However, for boundary points, singularity still exists and more treatment should be performed (see [22, 23]).


## 5. Numerical Examples

In order to verify the presented approach in this paper, some 2D and 3D singular domain integrals are analyzed in this section.


### 5.1. 2D Highly Singular Domain Integrals over a Rectangular Region

The first example is aimed at verifying the correctness of the presented formulations. The following domain integral is analyzed:(28)I(xp)=∫04∫-11∂r/∂xrβdxdy.The integration domain considered is a rectangular region withx∈[-1,1] and y∈[0,4] (Figure 2). To evaluate the domain integral, the boundary of the rectangle is discretized into 12 equally spaced three-noded quadratic elements with 24 boundary nodes. Two source points xp (located at nodes 25 and 26 in Figure 2) are considered, the coordinates of which are shown in Table 1. Table 1 lists the computed results for singularity order β from 1 to 6. The computational results for the selected points using L=3,4,…,15 are denoted by RIM_B. For the integral, the analytical solution is easy to be obtained and also listed in Table 1.Table 1
Computational results for points 25 and 26.
Points
25 (0.5, 1)
26 (0.75, 1)
β
Analytical
Rim_B
Analytical
Rim_B
1
−2.0577
−2.057701
−3.12942
−3.129422
2
−1.86663
−1.866635
−3.42229
−3.422286
3
−1.94775
−1.947746
−5.1807
−5.180698
4
−2.29308
−2.293076
−10.341
−10.34099
5
−2.98027
−2.980266
−24.9956
−24.99558
6
−4.18609
−4.186087
−68.1973
−68.19732Computational model of the rectangle.
(a)
(b)From Table1 it can be seen that the results (RIM_B) obtained by the derived formulations in the paper are in excellent agreement with the analytical solutions. It demonstrates that the method described in this paper can indeed evaluate high-order singular domain integrals accurately.It is noted that the current computed results forL=3,4,…,15, which are corresponding to knot vectors defined as (12) for i=1,2,…,13, are the same. That is because the nonsingular part of the integrand is simple. It should also be noted that the more the Gauss points (here the number is 20) used are, the more accurate the results are. And the finer the mesh is, the more accurate the results are.


### 5.2. 2D Domain Integrals with Logarithm and Exponent over a Multiconnected Region

The second example is to demonstrate the capability of the method presented in this paper to evaluate complicated domain integrals over multiconnected regions. The domain considered is a quarter of fan-shaped sector with the inner and outer radiuses being 3 and 12, respectively, which is cut by two cavities (with inner radius of 6 and outer radius of 9, and with a 90/8-degree angle), as shown in Figure3.Figure 3
A multiconnected 2D domain with the discretized boundary.Two domain integrals are examined as follows:(29a)I1(xp)=∬Sln⁡(x+6)rβdxdy,(29b)I2(xp)=∬Sr,2e-x2/10rβdxdy.To compute the singular domain integrals, the boundary of the domain is discretized into 114 three-noded quadratic elements with 228 boundary nodes (Figure3). Three internal points A, B, and C are selected as the source points (see Figure 3), the coordinates of which are shown in Table 2. Computational results for these points using L=15 are listed in Table 3 and denoted by RIM_B. To verify the correctness of the computational results, the two domain integrals are also computed using the method presented in the literature [7]. The corresponding computational results are listed in the columns denoted by [7] in Table 3.Table 2
Coordinates of points A, B, and C.
Point
x
y
A
5.3033
5.3033
B
7.0711
7.0711
C
1.4632
7.3559Table 3
Computational results for the selected points.
Point
A
B
C
β
Ref. [7]
RIM_B
Ref. [7]
RIM_B
Ref. [7]
RIM_B
0
I
1
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
I
2
6.130514
E
+
00
6.130513
E
+
00
-
7.559164
E
-
01
-
7.559170
E
-
01
2.320240
E
+
01
2.320239
E
+
01
1
I
1
7.703283
E
+
01
7.703283
E
+
01
7.612753
E
+
01
7.612753
E
+
01
4.252796
E
+
01
4.252796
E
+
01
I
2
8.426019
E
-
01
8.426019
E
-
01
-
1.194605
E
-
01
-
1.194604
E
-
01
5.074336
E
+
00
5.074332
E
+
00
2
I
1
2.193844
E
+
01
2.193844
E
+
01
2.239734
E
+
01
2.239734
E
+
01
1.109208
E
+
01
1.109208
E
+
01
I
2
7.756043
E
-
02
7.756118
E
-
02
-
1.802284
E
-
02
-
1.802279
E
-
02
1.599750
E
+
00
1.599750
E
+
00
3
I
1
-
4.783002
E
+
00
-
4.783001
E
+
00
-
4.731465
E
+
00
-
4.731465
E
+
00
4.294813
E
+
00
4.294823
E
+
00
I
2
-
1.406200
E
-
02
-
1.406429
E
-
02
-
2.564889
E
-
03
-
2.561466
E
-
03
7.101717
E
-
01
7.101939
E
-
01
4
I
1
-
1.233865
E
+
00
-
1.233880
E
+
00
-
8.742970
E
-
01
-
8.743047
E
-
01
2.256689
E
+
00
2.256630
E
+
00
I
2
-
1.500865
E
-
02
-
1.502554
E
-
02
-
3.388990
E
-
04
-
3.647890
E
-
04
4.033122
E
-
01
4.032076
E
-
01
5
I
1
-
5.133504
E
-
01
-
5.133268
E
-
01
-
2.344087
E
-
01
-
2.343959
E
-
01
1.437597
E
+
00
1.437600
E
+
00
I
2
-
8.750686
E
-
03
-
8.747977
E
-
03
-
4.725454
E
-
05
-
4
.
691195
E
-
05
2.682402
E
-
01
2.682987
E
-
01Table3 shows that when β=0, 1, or 2, the two sets of computational results are in good agreement, while, for high-order singular domain integrals, they are moderately close.To examine the convergence performance of the computational result with respect to the number of B-spline basis functions, the domain integralI1 was evaluated for different values of L. Table 4 gives the computational results of the source point A for L=3,5,7,9,12,15. It can be seen that good results can be obtained using L≥5.Table 4
Computational results of domain integralI1 for point A using different values of L.
L
3
5
7
9
12
15
β
=
0
2.321804
E
+
02
2.321857
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
β
=
1
7.703125
E
+
01
7.703279
E
+
01
7.703283
E
+
01
7.703283
E
+
01
7.703283
E
+
01
7.703283
E
+
01
β
=
2
2.193677
E
+
01
2.193843
E
+
01
2.193844
E
+
01
2.193844
E
+
01
2.193844
E
+
01
2.193844
E
+
01
β
=
3
-
4.781597
E
+
00
-
4.782937
E
+
00
-
4.782993
E
+
00
-
4.782995
E
+
00
-
4.783000
E
+
00
-
4.783001
E
+
00
β
=
4
-
1.232382
E
+
00
-
1.233940
E
+
00
-
1.233920
E
+
00
-
1.233901
E
+
00
-
1.233887
E
+
00
-
1.233880
E
+
00
β
=
5
-
5.142504
E
-
01
-
5.134197
E
-
01
-
5.133521
E
-
01
-
5.133394
E
-
01
-
5.133293
E
-
01
-
5.133268
E
-
01


### 5.3. 3D Singular Domain Integrals over a Hexagonal Prism

The third example is used to examine the performance of the proposed method in solving 3D high-order singular domain integrals. The following two singular domain integrals are considered:(30a)I1=∭Vxy2zrβdxdydz,(30b)I2=∭Vxyezrβdxdydz.The integration domainV is a hexagonal prism with the height and side length being 1 and 2, respectively (Figure 4). To evaluate the integrals I1 and I2, the surfaces of the hexagonal prism are discretized into 90 eight-noded quadratic elements with a total of 272 boundary nodes (Figure 4). Three internal points (273, 274, and 275) with the coordinates of (1,1,0.5), (0.5,0.5,0.25), and (1.5,1.5,0.75) are selected as the source points. The computational results using L=15 are listed in Table 5 and denoted by RIM_B. The two domain integrals are also analyzed by the proposed method in [7] and the computational results are shown in Table 5 in the columns denoted by [7].Table 5
Integration results for integralsI1 and I2 with three source points 273, 274, and 275.
Point
273
274
275
β
Ref. [7]
RIM_B
Ref.  [7]
RIM_B
Ref.  [7]
RIM_B
0
I
1
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
I
2
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
1
I
1
4.410999
E
+
00
4.410999
E
+
00
2.516612
E
+
00
2.516612
E
+
00
6.896879
E
+
00
6.896879
E
+
00
I
2
1.147215
E
+
01
1.147215
E
+
01
6.806745
E
+
00
6.806745
E
+
00
1.521900
E
+
01
1.521900
E
+
01
2
I
1
7.306388
E
+
00
7.306388
E
+
00
2.025978
E
+
00
2.025978
E
+
00
1.948790
E
+
01
1.948790
E
+
01
I
2
2.095453
E
+
01
2.095453
E
+
01
7.064466
E
+
00
7.064466
E
+
00
3.936297
E
+
01
3.936297
E
+
01
3
I
1
9.174757
E
-
01
9.174757
E
-
01
1.297864
E
+
00
1.297864
E
+
00
-
1.805546
E
+
01
-
1.805546
E
+
01
I
2
-
5.447484
E
-
01
-
5.447484
E
-
01
1.803357
E
+
00
1.803357
E
+
00
-
3.043079
E
+
01
-
3.043079
E
+
01
4
I
1
-
4.853797
E
+
00
-
4.853797
E
+
00
9.380380
E
-
01
9.380117
E
-
01
-
6.026606
E
+
01
-
6.026609
E
+
01
I
2
-
2.123922
E
+
01
-
2.123923
E
+
01
-
3.225239
E
+
00
-
3.225267
E
+
00
-
1.080169
E
+
02
-
1.080169
E
+
02
5
I
1
-
4.489810
E
+
00
-
4.489810
E
+
00
5.271243
E
-
01
5.288588
E
-
01
-
7.685096
E
+
01
-
7.684923
E
+
01
I
2
-
1.661141
E
+
01
-
1.660959
E
+
01
-
3.909619
E
+
00
-
3.908771
E
+
00
-
1.376579
E
+
02
-
1.376506
E
+
02Figure 4
Boundary nodes over a hexagonal prism.Comparing the numerical results in Table5, it reveals that the two sets of computational results are in excellent agreement.


### 5.4. 3D Strongly Singular Domain Integral over a Quarter Circular Tube

The fourth example is concerned with the following strongly singular domain integral over a 3D curved circular tube with the radiusr=2 (Figure 5). The curvature of the pipe is determined by a radius of R=8. The following singular domain integral is considered:
(31)Iijklxp=∫ΩEijklxp,xrDxp,xdΩx,
where
(32)Eijklxp,x=14D-1π1-ν×1-2νδikδlj+δjkδli-δijδkl+Dδijr,kr,l+Dν(δlir,jr,k+δjkr,lr,i+δikr,lr,j+δjlr,ir,k)+Dδklr,ir,j-DD+2r,ir,jr,kr,l,
in which the subscripts i,j,k,l=1,2,3 are for the dimensions in the 3D problem, ν is the Poisson ratio (here assumed to be 0.3), r,i is determined by (7), and D=3. Eijkl(xp,x) can be found in the evaluation of the internal stresses in elastoplastic mechanics [2, 24]. Here only I1111, I1122, and I2121 are computed.Boundary nodes over a circular tube.
(a)
(b)To evaluate these three integrals, the surfaces of the tube are discretized into 282 eight-noded quadratic elements with a total of 848 boundary nodes (Figure5).The problem was analyzed by Gao and Peng [7] and the results are used here as comparison. Table 7 gives the computed results using L=9 for three source points (849, 850, and 851) with coordinates listed in Table 6. It can be found that they are in excellent agreement from Table 7.Table 6
Coordinates of selected source points.
Point
x
y
z
849
−1.0711
11.0711
0
850
−3.2
9
0
851
1
12
0Table 7
Computed results for three source points.
Point
849
850
851
Integral
Ref. [7]
RIM_B
Ref. [7]
RIM_B
Ref. [7]
RIM_B
I
1111
-
8.086983
E
-
02
-
8.086983
E
-
02
3.839788
E
-
02
3.839788
E
-
02
-
1.903177
E
-
01
-
1.903177
E
-
01
I
1122
-
2.882396
E
-
02
-
2.882396
E
-
02
-
6.270116
E
-
02
-
6.270116
E
-
02
2.713745
E
-
03
2.713745
E
-
03
I
2121
-
7.232163
E
-
02
-
7.232163
E
-
02
-
6.342540
E
-
02
-
6.342540
E
-
02
-
8.081374
E
-
02
-
8.081374
E
-
02


## 5.1. 2D Highly Singular Domain Integrals over a Rectangular Region

The first example is aimed at verifying the correctness of the presented formulations. The following domain integral is analyzed:(28)I(xp)=∫04∫-11∂r/∂xrβdxdy.The integration domain considered is a rectangular region withx∈[-1,1] and y∈[0,4] (Figure 2). To evaluate the domain integral, the boundary of the rectangle is discretized into 12 equally spaced three-noded quadratic elements with 24 boundary nodes. Two source points xp (located at nodes 25 and 26 in Figure 2) are considered, the coordinates of which are shown in Table 1. Table 1 lists the computed results for singularity order β from 1 to 6. The computational results for the selected points using L=3,4,…,15 are denoted by RIM_B. For the integral, the analytical solution is easy to be obtained and also listed in Table 1.Table 1
Computational results for points 25 and 26.
Points
25 (0.5, 1)
26 (0.75, 1)
β
Analytical
Rim_B
Analytical
Rim_B
1
−2.0577
−2.057701
−3.12942
−3.129422
2
−1.86663
−1.866635
−3.42229
−3.422286
3
−1.94775
−1.947746
−5.1807
−5.180698
4
−2.29308
−2.293076
−10.341
−10.34099
5
−2.98027
−2.980266
−24.9956
−24.99558
6
−4.18609
−4.186087
−68.1973
−68.19732Computational model of the rectangle.
(a)
(b)From Table1 it can be seen that the results (RIM_B) obtained by the derived formulations in the paper are in excellent agreement with the analytical solutions. It demonstrates that the method described in this paper can indeed evaluate high-order singular domain integrals accurately.It is noted that the current computed results forL=3,4,…,15, which are corresponding to knot vectors defined as (12) for i=1,2,…,13, are the same. That is because the nonsingular part of the integrand is simple. It should also be noted that the more the Gauss points (here the number is 20) used are, the more accurate the results are. And the finer the mesh is, the more accurate the results are.


## 5.2. 2D Domain Integrals with Logarithm and Exponent over a Multiconnected Region

The second example is to demonstrate the capability of the method presented in this paper to evaluate complicated domain integrals over multiconnected regions. The domain considered is a quarter of fan-shaped sector with the inner and outer radiuses being 3 and 12, respectively, which is cut by two cavities (with inner radius of 6 and outer radius of 9, and with a 90/8-degree angle), as shown in Figure3.Figure 3
A multiconnected 2D domain with the discretized boundary.Two domain integrals are examined as follows:(29a)I1(xp)=∬Sln⁡(x+6)rβdxdy,(29b)I2(xp)=∬Sr,2e-x2/10rβdxdy.To compute the singular domain integrals, the boundary of the domain is discretized into 114 three-noded quadratic elements with 228 boundary nodes (Figure3). Three internal points A, B, and C are selected as the source points (see Figure 3), the coordinates of which are shown in Table 2. Computational results for these points using L=15 are listed in Table 3 and denoted by RIM_B. To verify the correctness of the computational results, the two domain integrals are also computed using the method presented in the literature [7]. The corresponding computational results are listed in the columns denoted by [7] in Table 3.Table 2
Coordinates of points A, B, and C.
Point
x
y
A
5.3033
5.3033
B
7.0711
7.0711
C
1.4632
7.3559Table 3
Computational results for the selected points.
Point
A
B
C
β
Ref. [7]
RIM_B
Ref. [7]
RIM_B
Ref. [7]
RIM_B
0
I
1
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
I
2
6.130514
E
+
00
6.130513
E
+
00
-
7.559164
E
-
01
-
7.559170
E
-
01
2.320240
E
+
01
2.320239
E
+
01
1
I
1
7.703283
E
+
01
7.703283
E
+
01
7.612753
E
+
01
7.612753
E
+
01
4.252796
E
+
01
4.252796
E
+
01
I
2
8.426019
E
-
01
8.426019
E
-
01
-
1.194605
E
-
01
-
1.194604
E
-
01
5.074336
E
+
00
5.074332
E
+
00
2
I
1
2.193844
E
+
01
2.193844
E
+
01
2.239734
E
+
01
2.239734
E
+
01
1.109208
E
+
01
1.109208
E
+
01
I
2
7.756043
E
-
02
7.756118
E
-
02
-
1.802284
E
-
02
-
1.802279
E
-
02
1.599750
E
+
00
1.599750
E
+
00
3
I
1
-
4.783002
E
+
00
-
4.783001
E
+
00
-
4.731465
E
+
00
-
4.731465
E
+
00
4.294813
E
+
00
4.294823
E
+
00
I
2
-
1.406200
E
-
02
-
1.406429
E
-
02
-
2.564889
E
-
03
-
2.561466
E
-
03
7.101717
E
-
01
7.101939
E
-
01
4
I
1
-
1.233865
E
+
00
-
1.233880
E
+
00
-
8.742970
E
-
01
-
8.743047
E
-
01
2.256689
E
+
00
2.256630
E
+
00
I
2
-
1.500865
E
-
02
-
1.502554
E
-
02
-
3.388990
E
-
04
-
3.647890
E
-
04
4.033122
E
-
01
4.032076
E
-
01
5
I
1
-
5.133504
E
-
01
-
5.133268
E
-
01
-
2.344087
E
-
01
-
2.343959
E
-
01
1.437597
E
+
00
1.437600
E
+
00
I
2
-
8.750686
E
-
03
-
8.747977
E
-
03
-
4.725454
E
-
05
-
4
.
691195
E
-
05
2.682402
E
-
01
2.682987
E
-
01Table3 shows that when β=0, 1, or 2, the two sets of computational results are in good agreement, while, for high-order singular domain integrals, they are moderately close.To examine the convergence performance of the computational result with respect to the number of B-spline basis functions, the domain integralI1 was evaluated for different values of L. Table 4 gives the computational results of the source point A for L=3,5,7,9,12,15. It can be seen that good results can be obtained using L≥5.Table 4
Computational results of domain integralI1 for point A using different values of L.
L
3
5
7
9
12
15
β
=
0
2.321804
E
+
02
2.321857
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
2.321859
E
+
02
β
=
1
7.703125
E
+
01
7.703279
E
+
01
7.703283
E
+
01
7.703283
E
+
01
7.703283
E
+
01
7.703283
E
+
01
β
=
2
2.193677
E
+
01
2.193843
E
+
01
2.193844
E
+
01
2.193844
E
+
01
2.193844
E
+
01
2.193844
E
+
01
β
=
3
-
4.781597
E
+
00
-
4.782937
E
+
00
-
4.782993
E
+
00
-
4.782995
E
+
00
-
4.783000
E
+
00
-
4.783001
E
+
00
β
=
4
-
1.232382
E
+
00
-
1.233940
E
+
00
-
1.233920
E
+
00
-
1.233901
E
+
00
-
1.233887
E
+
00
-
1.233880
E
+
00
β
=
5
-
5.142504
E
-
01
-
5.134197
E
-
01
-
5.133521
E
-
01
-
5.133394
E
-
01
-
5.133293
E
-
01
-
5.133268
E
-
01


## 5.3. 3D Singular Domain Integrals over a Hexagonal Prism

The third example is used to examine the performance of the proposed method in solving 3D high-order singular domain integrals. The following two singular domain integrals are considered:(30a)I1=∭Vxy2zrβdxdydz,(30b)I2=∭Vxyezrβdxdydz.The integration domainV is a hexagonal prism with the height and side length being 1 and 2, respectively (Figure 4). To evaluate the integrals I1 and I2, the surfaces of the hexagonal prism are discretized into 90 eight-noded quadratic elements with a total of 272 boundary nodes (Figure 4). Three internal points (273, 274, and 275) with the coordinates of (1,1,0.5), (0.5,0.5,0.25), and (1.5,1.5,0.75) are selected as the source points. The computational results using L=15 are listed in Table 5 and denoted by RIM_B. The two domain integrals are also analyzed by the proposed method in [7] and the computational results are shown in Table 5 in the columns denoted by [7].Table 5
Integration results for integralsI1 and I2 with three source points 273, 274, and 275.
Point
273
274
275
β
Ref. [7]
RIM_B
Ref.  [7]
RIM_B
Ref.  [7]
RIM_B
0
I
1
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
3.680608
E
+
00
I
2
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
8.928454
E
+
00
1
I
1
4.410999
E
+
00
4.410999
E
+
00
2.516612
E
+
00
2.516612
E
+
00
6.896879
E
+
00
6.896879
E
+
00
I
2
1.147215
E
+
01
1.147215
E
+
01
6.806745
E
+
00
6.806745
E
+
00
1.521900
E
+
01
1.521900
E
+
01
2
I
1
7.306388
E
+
00
7.306388
E
+
00
2.025978
E
+
00
2.025978
E
+
00
1.948790
E
+
01
1.948790
E
+
01
I
2
2.095453
E
+
01
2.095453
E
+
01
7.064466
E
+
00
7.064466
E
+
00
3.936297
E
+
01
3.936297
E
+
01
3
I
1
9.174757
E
-
01
9.174757
E
-
01
1.297864
E
+
00
1.297864
E
+
00
-
1.805546
E
+
01
-
1.805546
E
+
01
I
2
-
5.447484
E
-
01
-
5.447484
E
-
01
1.803357
E
+
00
1.803357
E
+
00
-
3.043079
E
+
01
-
3.043079
E
+
01
4
I
1
-
4.853797
E
+
00
-
4.853797
E
+
00
9.380380
E
-
01
9.380117
E
-
01
-
6.026606
E
+
01
-
6.026609
E
+
01
I
2
-
2.123922
E
+
01
-
2.123923
E
+
01
-
3.225239
E
+
00
-
3.225267
E
+
00
-
1.080169
E
+
02
-
1.080169
E
+
02
5
I
1
-
4.489810
E
+
00
-
4.489810
E
+
00
5.271243
E
-
01
5.288588
E
-
01
-
7.685096
E
+
01
-
7.684923
E
+
01
I
2
-
1.661141
E
+
01
-
1.660959
E
+
01
-
3.909619
E
+
00
-
3.908771
E
+
00
-
1.376579
E
+
02
-
1.376506
E
+
02Figure 4
Boundary nodes over a hexagonal prism.Comparing the numerical results in Table5, it reveals that the two sets of computational results are in excellent agreement.


## 5.4. 3D Strongly Singular Domain Integral over a Quarter Circular Tube

The fourth example is concerned with the following strongly singular domain integral over a 3D curved circular tube with the radiusr=2 (Figure 5). The curvature of the pipe is determined by a radius of R=8. The following singular domain integral is considered:
(31)Iijklxp=∫ΩEijklxp,xrDxp,xdΩx,
where
(32)Eijklxp,x=14D-1π1-ν×1-2νδikδlj+δjkδli-δijδkl+Dδijr,kr,l+Dν(δlir,jr,k+δjkr,lr,i+δikr,lr,j+δjlr,ir,k)+Dδklr,ir,j-DD+2r,ir,jr,kr,l,
in which the subscripts i,j,k,l=1,2,3 are for the dimensions in the 3D problem, ν is the Poisson ratio (here assumed to be 0.3), r,i is determined by (7), and D=3. Eijkl(xp,x) can be found in the evaluation of the internal stresses in elastoplastic mechanics [2, 24]. Here only I1111, I1122, and I2121 are computed.Boundary nodes over a circular tube.
(a)
(b)To evaluate these three integrals, the surfaces of the tube are discretized into 282 eight-noded quadratic elements with a total of 848 boundary nodes (Figure5).The problem was analyzed by Gao and Peng [7] and the results are used here as comparison. Table 7 gives the computed results using L=9 for three source points (849, 850, and 851) with coordinates listed in Table 6. It can be found that they are in excellent agreement from Table 7.Table 6
Coordinates of selected source points.
Point
x
y
z
849
−1.0711
11.0711
0
850
−3.2
9
0
851
1
12
0Table 7
Computed results for three source points.
Point
849
850
851
Integral
Ref. [7]
RIM_B
Ref. [7]
RIM_B
Ref. [7]
RIM_B
I
1111
-
8.086983
E
-
02
-
8.086983
E
-
02
3.839788
E
-
02
3.839788
E
-
02
-
1.903177
E
-
01
-
1.903177
E
-
01
I
1122
-
2.882396
E
-
02
-
2.882396
E
-
02
-
6.270116
E
-
02
-
6.270116
E
-
02
2.713745
E
-
03
2.713745
E
-
03
I
2121
-
7.232163
E
-
02
-
7.232163
E
-
02
-
6.342540
E
-
02
-
6.342540
E
-
02
-
8.081374
E
-
02
-
8.081374
E
-
02


## 6. Conclusions

A new and robust method is described in this paper for evaluating arbitrary singular domain integrals based on RIM and by expressing the nonsingular part of the integration kernel as a series of cubic B-spline basis functions of the global distancer. The distinct feature of the method is that regular, weakly, strongly, hyper-, and supersingular domain integrals are treated in a simple and unified manner. The provided numerical examples demonstrate that the presented method is accurate and robust.



---

*Source: 284106-2014-11-17.xml*