Volume 65, 2015, 93–111
Marvin Fleck, Richards Grzhibovskis, and Sergej Rjasanow
A NEW FUNDAMENTAL SOLUTION METHOD BASED ON THE ADAPTIVE CROSS APPROXIMATION
Dedicated to Roland Duduchava on the occasion of his 70th birthday
approximate solution of scalar elliptic boundary value problems is presented.
The construction of the basis functions is based on the Adaptive Cross Approximation (ACA) of the fundamental solutions of the corresponding elliptic operator. An algorithm for an immediate computer implementation of the method is formulated. A series of numerical examples for the Laplace and Helmholtz equations in three dimensions illustrates the efficiency of the method. Extensions of the method to elliptic systems are discussed.
2010 Mathematics Subject Classification. 65N80, 65N12, 65N35.
Key words and phrases. Fundamental solution method, adaptive cross approximation, collocation, condition numbers.
ÒÄÆÉÖÌÄ.
ßÀÒÌÏÃÂÄÍÉËÉÀ ×ÖÍÃÀÌÄÍÔÖÒÉ ÀÌÏÍÀáÓÍÄÁÉÓ ÀáÀËÉ ÀÃÀÐÔÉÖÒÉ ÌÄÈÏÃÉÓ (FSM) ÂÀÌÏÚÄÍÄÁÀ ÓÊÀËÀÒÖËÉ ÄËÉ×ÓÖÒÉ ÓÀÓÀ- ÆÙÅÒÏ ÀÌÏÝÀÍÄÁÉÓ ÌÉÀáËÏÄÁÉÈÉ ÀÌÏáÓÍÉÓÈÅÉÓ. ÓÀÁÀÆÉÓÏ ×ÖÍØÝÉ- ÄÁÉÓ ÀÂÄÁÀ ÄÌÚÀÒÄÁÀ ÛÄÓÀÁÀÌÉÓÉ ÄËÉ×ÓÖÒÉ ÏÐÄÒÀÔÏÒÉÓ ×ÖÍÃÀÌÄÍÔÖ- ÒÉ ÀÌÏÍÀáÓÍÄÁÉÓ ÀÃÀÐÔÉÖÒ ãÅÀÒÄÃÉÍ ÀÐÒÏØÓÉÌÀÝÉÀÓ (ACA). ×ÏÒÌÖ- ËÉÒÄÁÖËÉÀ ÀÌ ÌÄÈÏÃÉÓ ÖÛÖÀËÏ ÊÏÌÐÉÖÔÄÒÖËÉ ÂÀÍáÏÒÝÉÄËÄÁÉÓ ÀËÂÏÒÉÈÌÉ. ÌÄÈÏÃÉÓ Ä×ÄØÔÖÒÏÁÀ ÉËÖÓÔÒÉÒÄÁÖËÉÀ ÓÀÌÂÀÍÆÏÌÉ- ËÄÁÉÀÍÉ ËÀÐËÀÓÉÓÀ ÃÀ äÄËÌäÏËÝÉÓ ÂÀÍÔÏËÄÁÄÁÉÓÈÅÉÓ ÌÏÚÅÀÍÉËÉ ÒÉÝáÅÉÈÉ ÌÀÂÀËÉÈÄÁÉÈ. ÂÀÍáÉËÖËÉÀ ÀÌ ÌÄÈÏÃÉÓ ÂÀÍÆÏÂÀÃÄÁÉÓ ÛÄÓÀÞËÄÁËÏÁÀ ÄËÉ×ÓÖÒÉ ÓÉÓÔÄÌÄÁÉÓÈÅÉÓ.1. Introduction
The Fundamental Solution Method (FSM) is also known as the Method of Fundamental Solutions, Charge Simulation Method or as a special version of the Boundary Collocation Method. It resembles a Trefftz method [7], which means that the solution to a Dirichlet boundary value problem in Ω⊂R3,Γ =∂Ω
Lu(x) = 0 for x∈Ω, u(x) =g(x) for x∈Γ,
is approximated by a linear combination of L-harmonic functions. As the name indicates, the method uses fundamental solutions for basis functions, whose singularities are located outside Ω. It was introduced by Kupradze and Aleksidze [4] in 1963 for treating the Laplace equation. First inves- tigations from a numerical point of view were performed by Mathon and Johnston [5] in 1977. Comprehensive summaries of the attributes of the FSM were written, among others, by Smyrlis [6] and Bogomolny [3].
Two peculiar aspects of the Fundamental Solution Method are an ex- tremely fast convergence, but also a very high condition number of the system matrix, both with respect to a number of collocation points. We address the problem of high condition numbers by adaptively choosing a smaller number of collocation points while keeping the local error below a given threshold, but not necessarily equal to zero, for the remaining collo- cation points. Thus an approximation is obtained, while condition numbers are kept lower due to smaller system matrices. The quality of the approxi- mation is comparable to that of classical FSM. By means of this approach the problems that are too big for classical FSM can be treated. The adap- tive strategy features are new basis functions which vanish at collocation points already treated and thus do not alter the corresponding local approx- imation. The construction of these basis functions uses concepts from the Adaptive Cross Approximation (ACA) [2].
In Section 2 we formulate a model problem and present the classical (collocation-based) Fundamental Solution Method. Section 3 briefly sum- marizes the Adaptive Cross Approximation. The approximation algorithm presented therein leads directly to the construction of basis functions for the Adaptive Fundamental Solution Method in Section 4. In Section 5 we present numerical results for the adaptive method applied to the Laplace and Helmholtz equations, respectively.
2. Formulation of the Problem
We consider the following Dirichlet boundary value problem for an elliptic equation inR3
Lu(x) = 0 for x∈Ω, u(x) =g(x) for x∈Γ,
where L is an elliptic second order differential operator and Ω ⊂ R3 is a Lipschitz domain with the boundaryΓ. In the classical setting, the Dirichlet datumgis assumed to be continuous onΓand the solutionuis assumed to be smooth, i.e.
u∈C2(Ω)∩C(Ω). In this paper, we will consider the Laplace operator
Lu=−∆u and the Helmholtz operator
Lu=−∆u−κ2u .
For these operators, the corresponding fundamental solution u∗, i.e. the solution of the equation
Lu∗=δ (2)
in the distributional sense is known and given by u∗(x) = 1
4π 1
|x| for the Laplace operator, and
u∗(x) = 1 4π
eıκ|x|
|x|
for the Helmholtz operator. In (2)δdenotes the Dirac δ-distribution.
2.1. Fundamental solution method. Let X ⊂ Γ be a discrete set of N pairwise different control (collocation) points on the boundary Γ and Y ⊂R3\Ωa discrete set ofN pairwise different singularity points. Consider a system of basis functions
Φ ={
φ1, . . . , φN
}, φℓ(x) =u∗(x−yℓ), yℓ∈Y, ℓ= 1, . . . , N.
Sinceyℓ̸∈Ω,ℓ= 1, . . . , N, every basis functionφℓ isL-harmonic inΩand the function
uN(x) =
∑N ℓ=1
αℓφℓ(x) = Φ(x)a , a= (α1, . . . , αN)⊤∈RN
can be considered as an approximation of the solution uof the boundary value problem (1). The most simple choice for the coefficients αℓ is the point collocation for the boundary condition
uN(x) =g(x) for x∈X.
This can be equivalently formulated as a linear system for obtaining N coefficientsαl:
∑N ℓ=1
αℓφℓ(xk) =g(xk) for k= 1, . . . , N
or, in a matrix form,
F a=g , (3)
where F =(
φℓ(xk))N
k,ℓ=1∈RN×N, aℓ=αℓ, g
ℓ=g(xℓ), 1≤ℓ≤N.
The main properties of the FSM can be summarised as follows.
1. Since no topology of the discrete point setsX and Y is required, the FSM can be considered as a meshfree numerical method.
2. The entries of the matrixF in (3) are easy to compute as opposed to matrix entries coming from Boundary Element Methods (BEM).
3. The dimension of the matrixF is comparable to those of the BEM (e.g.N ∼104−105for 3D problems).
4. The matrix F is fully populated as in the BEM and M em(F) = O(N2).
5. The condition number of the matrix F grows exponentially, i.e.
cond(F) =O(qN)for some q >1.
For large N the application of a direct solver to the system (3) is expen- sive, while an iterative solver does not converge due to the extremely high condition number of the matrixF.
However, the numerical results for small systems show an exponential convergence of the method not only for the solutionuitself but also for its gradient
gradu= (∂u
∂x1
, ∂u
∂x2
, ∂u
∂x3
)⊤
and even for its Hessian matrix Hu=
( ∂2u
∂xk∂xℓ
)3 k,ℓ=1
, i.e.
O(|u(x)−uN(x)|) =O(grad(u(x)−uN(x)))
=O(H(u(x)−uN(x))
F
)
=O(q−N)
forx∈Ω. Note that the derivatives of the approximate solutionuN can be easily computed analytically.
2.2. Choice of pseudo boundary. In the theoretical analysis of Funda- mental Solution Methods one introduces the concept of pseudo-boundaries, i.e. surfaces where the singularity points are located. Pseudo-boundaries fulfilling the so-called embracing condition provide for the suitability of corresponding fundamental solutions as basis functions [6]. However, one still has great freedom in choosing an actual pseudo-boundary and in the subsequent choice of the location of singularity points.
Here we briefly present the definition and the central theorem for pseudo- boundaries. A thorough overview can be found in [6].
Definition 1 (Segment condition). LetΩ⊂Rd be an open set. Ωfulfills the segment condition, if for everyx∈∂Ωthere exist a neighborhoodU(x) of x and a nonzero vector ξ(x) ∈ Rd such that if y ∈ U(x)∩Ω, then y+tξ(x)∈Ω,∀t∈(0,1).
Definition 2 (Embracing boundary). Let Ω,Ω′ ⊂ Rd be open and con- nected. Ω′ embraces Ω, if:
1. Ω⊂Ω′;
2. For each connected component V of Rd\Ω there is an open con- nected componentV′ ofRd\Ω′ such thatV ′⊂V.
Theorem 1. If Ω⊂Rdfulfills the segment condition andΩ′⊂Rdembraces Ω, then ford≥3andl≥0the spaceXspanned by finite linear combinations of Fundamental solutions
uN(x) =
∑N j=1
αju∗(x−yj) with singularities yj∈∂Ω′ is dense in
Yl={
v∈C2(Ω) : ∆v= 0in Ω}
∩Cl(Ω)
with respect to the norm of Cl(Ω). Ford= 2the density result holds true forX ⊕ {c·1|Ω: c∈R}.
Proof. The proof can be found in [6].
Similar results exist for the operators ∆m, m > 1, and ∆−κ2, κ >0, [6, 3].
One can prove in the two-dimensional case that an increase in the distance between the boundary ∂Ω and the pseudo-boundary ∂Ω′ leads both to a better approximation and to a larger condition number ofF. This can also be observed in three-dimensional settings.
For simple domain shapes the common choice of the singularity points consists in shifting collocation points along the outer normal. This strategy may fail for more complex domains. On the other hand, the construction of pseudo boundaries by means of distance functions may be computationally expensive.
In what follows, we will introduce a new method with the same conver- gence properties but almost without disadvantages of the FSM, i.e. without necessity of numerical solving of big, dense and badly conditioned systems of linear equations. Our main tool is the Adaptive Cross Approximation.
3. Adaptive Cross Approximation
The initial analytical form of the ACA algorithm was been designed to interpolate and, hopefully, to approximate a given functionK:X×Y→R of two variablesxandy by a degenerate functionSn, i.e.
K(x, y)≈Sn(x, y) =
∑n ℓ=1
uℓ(x)vℓ(y),
where ul : X → R, vl : Y → R, l = 1, . . . , n. The construction runs as follows. LetX⊂X⊂R3andY ⊂Y⊂R3be discrete point sets.
Algorithm 1.
1. initialization
1.1 set initial residual and initial approximation R0(x, y) =K(x, y), S0(x, y) = 0;
1.2 choose initial pivot position x0∈X, y0∈Y, R0(x0, y0)̸= 0.
2. recursion for k= 0,1, . . .. 2.1 new residual
Rk+1(x, y) =Rk(x, y)−Rk(x, yk)Rk(xk, y) Rk(xk, yk) ; 2.2 new approximation
Sk+1(x, y) =Sk(x, y) +Rk(x, yk)Rk(xk, y) Rk(xk, yk) ; 2.3 new pivot position
xk+1∈X, yk+1∈Y, Rk+1(xk+1, yk+1)̸= 0.
After n ≥ 1 steps of the ACA-Algorithm 1, we obtain a sequence of residualsR0, . . . , Rn and a sequence of approximationsS0, . . . , Sn with the following properties.
1. Approximation property fork= 0, . . . , n
Rk(x, y) +Sk(x, y) =K(x, y), x∈X, y∈Y; (4) 2. Interpolation property fork= 1, . . . , nandℓ= 0, . . . , k−1
Rk(x, yℓ) =Rk(xℓ, y) = 0, x∈X, y∈Y or
Sk(x, yℓ) =K(x, yℓ), x∈X, Sk(xℓ, y) =K(xℓ, y), y∈Y; 3. Harmonicity property fork= 0, . . . , n.
If
LxK(x, y) = 0, x∈Ω, then
LxRk(x, y) =LxSk(x, y) = 0, x∈Ω;
4. Non-recursive representation fork= 1, . . . , n
Sk(x, y) =u⊤k(x)Vk−1wk(y), Vk ∈Rk×k, uk(x), wk(y)∈Rk (5) with
uk(x) =(
K(x, y0), . . . , K(x, yk−1))⊤ , wk(y) =(
K(x0, y), . . . , K(xk−1, y))⊤ and
Vk =(
K(xi, yj))k−1 i,j=0
The above properties, except the last one, can be easily seen. The proof of the non-recursive representation is more technical and can be found in [2].
4. Adaptive FSM
In this section, we formulate a new adaptive FSM for the boundary value problem (1). Letu∗be the fundamental solution of the differential operator L,X⊂Γa discrete set of the control points,Y ⊂R3\Ωa discrete set of the singularity points andεan upper threshold for the error in the collocation points.
Algorithm 2.
1. initialization
1.1 initial error and initial pivot position Error1=Maxx∈X|g(x)|, x1=ArgMaxx∈X|g(x)|; 1.2 initial residual
R1(x, y) =u∗(x−y);
1.3 first basis function
φ1(x) = R1(x, y1) R1(x1, y1); 1.4 first approximation
u1(x) =α1φ1(x), α1=g(x1).
2. recursion for k= 1,2, . . .
2.1 new error and new pivot position
Errork+1=Maxx∈X|g(x)−uk(x)|, xk+1=ArgMaxx∈X|g(x)−uk(x)|; 2.2 stopping criteria
Stop if Errork+1≤ε or k=# of points in X;
2.3 next residual
Rk+1(x, y) =Rk(x, y)−Rk(x, yk)Rk(xk, y)
Rk(xk, yk) ; (6) 2.4 next basis function
φk+1(x) = Rk+1(x, yk+1)
Rk+1(xk+1, yk+1); (7)
2.5 next approximation
uk+1(x) =uk(x) +αk+1φk+1(x), αk+1=g(xk+1)−uk(xk+1).
After n steps of the above algorithm, we obtain the following approxi- mation:
un(x) =
∑n k=1
αkφk(x). (8)
The basis functionsφk areL-harmonic for allk= 1, . . . , n Lφk(x) = 0 for x∈Ω.
Therefore, the functionun is likewiseL-harmonic Lun(x) = 0 for x∈Ω.
The functionunfulfills the boundary condition pointwise at the pivot points un(xk) =g(xk) for k= 1, . . . , n
and approximates the boundary condition in the other points
|un(x)−g(x)| ≤ε for x∈X\ {x1, . . . , xn}.
Later on, our numerical examples will show that the number n of steps required to obtain a given accuracy is rather small compared to, and seems to be independent of, the number of control points N, i.e. n ≪ N. Due to the ACA interpolation property of the residualsRk, the basis functions φk, k≥2 vanish at all previous pivot points
φk(xℓ) = 0, ℓ= 1, . . . , k−1, k= 2, . . . , n and due to the construction,
φk(xk) = 1, k= 1, . . . , n.
Thus, the coefficientsαk in (8) can be easily computed as in Step 2.5of the Algorithm 2 without the need to solve a system of equations, or more precisely by solving a small system
F a=g , F =(
φℓ(xk))n
k,ℓ=1∈Rn×n, a, g∈Rn with the following triangular matrix
F=
1 0 0 . . . 0
φ1(x2) 1 0 . . . 0
. . . . φ1(xn) φ2(xn) φ3(xn) . . . 1
.
However, the price for the above simple and efficient algorithm is a more complicated evaluation of the basis functionsφk and hence of the approxi- mationun at a given pointx∈Ω. We use the non-recursive representation (5) of the ACA approximationSn, the approximation property (4), and the definition of the basis function inStep 2.4of the Algorithm 2 to obtain
φ1(x) = u∗(x−y1) u∗(x1−y1)
and fork= 2, . . . , n
φk(x) = u∗(x−yk)−u⊤k(x)zk
u∗(xk−yk)−u⊤k(xk)zk
(9) with
uk(x) =(
u∗(x−y1), . . . , u∗(x−yk−1))⊤ ,
zk=Vk−1wk(yk), (10)
wk(yk) =(
u∗(x1−yk), . . . , u∗(xk−1−yk))⊤ and
Vk =(
u∗(xi−yj))k−1 i,j=1.
The vectorszk∈Rk−1, as well as the normalizing constants(u∗(x1−y1))−1 and(u∗(xk−yk)−u⊤k(xk)zk)−1 can be precomputed during the algorithm as follows. Let
V2=L2U2= 1·u∗(x1−y1)
be the LU-decomposition of the1×1-matrixV2. Then, making use of the LU-decomposition of the(k−1)×(k−1)-matrix
Vk =LkUk, k= 2, . . . , n−1, we get for thek×k-matrixVk+1
Vk+1=
(Lk 0 a⊤ 1
) (Uk b 0 c )
with
Lkb=wk(yk), a⊤Uk =u⊤k(xk), c=u∗(xk−yk)−a⊤b.
For the vectorszk, we get
zk=Uk−1L−k1wk(yk) =Uk−1b.
From equations (6) and (7), we can see that
Rk+1(x, y) =Rk(x, y)−Rk(xk, y)R−k1(xk, yk)Rk(x, yk)
=Rk(x, y)−Rk(xk, y)φk(x)
=u∗(x−y)−
∑k j=1
Rj(xj, y)φj(x) ∀k >0 and thus
u∗(x−y) =Rk+1(x, y) +
∑k l=1
Rl(xl, y)φl(x) ∀k≥0.
This leads to
u∗(x1−y1) · · · u∗(x1−yk) ... . .. ... u∗(xk−y1) · · · u∗(xk−yk)
=
1 0 · · · 0
φ1(x2) 1 · · · 0 . . . . φ1(xk) φ2(xk) · · · 1
×
R1(x1, y1) R1(x1, y2) · · · R1(x1, yk) 0 R2(x2, y2) · · · R2(x2, yk) . . . . 0 0 · · · Rk(xk, yk)
.
Due to the uniqueness of the LU decomposition with unit diagonal entries, we see that
Ln =F.
A numerical evaluation of the basis functionφk in (9) requires the scalar productu⊤k(x)zk and, thereforeO(k)arithmetical operations. The approx- imate solutionun will requireO(n2)arithmetical operations for every eval- uation.
The adaptive Fundamental Solution Method allows us to elaborate alter- native strategies for choosing singularity points yi. Instead of introducing fixed pairs (xi, yi) of collocation and singularity points we may equip a simply shaped pseudo-boundary, e.g. an ellipsoid with a large number of uniformly distributed candidate points. The adaptive FSM can be tuned to pick from those candidates a singularity point that maximizes the current basis function’s pivot elementRk+1(xk+1, y)in (7).
5. Numerical Examples
In order to investigate the features of the adaptive Fundamental Solution Method, we perform numerical experiments for the BVP (1) with Laplace or Helmholtz operators. We compare the results of classical FSM, adaptive FSM with the given thresholds as well as a threshold-free adaptive method.
For the latter method we store the maximal local error of an iteration step and terminate, if no improvement is achieved after a given number of further iterations. By dropping the coefficients associated with these additional steps, we restore the currently best result (with respect to local errors).
The indicated condition numbers are calculated by using LAPACK rou- tines [1]. Condition numbers of respective system matrices are labeled condsys, while those of matrices required for evaluation of basis function in the adaptive method are labeled condLU. In the latter case we only in- dicate condition numbers of the respective largest matrix, i.e. the matrix used in the evaluation of the basis function with highest index.
Table 1. Laplace equation in the unit ball,N = 5120.
threshold max.err. rel.err. condsys condLU # nodes classical 8.80·10−10 1.01·10−11 1.42·1021 - 5120
10−4 9.78·10−5 1.08·10−6 1.91·102 4.53·1010 415 10−6 1.01·10−6 1.10·10−8 3.26·102 8.00·1012 711 10−8 1.05·10−8 1.13·10−10 4.94·102 7.26·1014 1069 10−10 1.09·10−10 1.11·10−12 6.09·102 9.81·1016 1531
5.1. Laplace equation. We consider the model problem
−∆v(x) = 0 for x∈Ω, v(x) =g(x) for x∈∂Ω, with the known analytical solution
v(x) =sin(2πx1)x2e−2πx3, g=v
Γ,
in order to display some general observations regarding the adaptive FSM, which are also relevant for other equation types we have already considered.
Performance of the adaptive FSM. As one can see in Table 1, the adaptive method uses only a small number of collocation points which increases upon setting a lower threshold. The accuracy of the full FSM can be achieved even with a relatively small subset of collocation points. This reduction leads to condition numbers of the involved matrices which are significantly lower than those of the full method’s system matrix. Since the approxi- mation space of the adaptive method is always a subset of the full FSM approximation space, the outperformance in the last row in the table can only be explained by a loss of accuracy due to high condition numbers.
Evolution of maximal local error. The strategy of the adaptive method con- sists in eliminating the currently largest residual of all collocation points, while not altering those at collocation points already treated. However, there is no guarantee that after any elimination step the new maximal error is actually smaller than the previous one. In fact, as the maximal error asymptotically decreases during the elimination process, short-term increases are rather typical (cf. Figure 1 for a brick-shaped domain).
For tight thresholds the adaptive method uses a number of collocation points comparable to that of classical FSM. For very large problems this may lead to errors in the evaluation of basis functions (evaluation ofzk in (10)) and ultimately to an asymptotic increase of the maximal local error.
For these cases it is handy to store information about the “best” step so far and restore the corresponding result. Thus, although the threshold is not met, the results in these extreme cases are far better than those of the classical method.
1e-05 0.0001 0.001 0.01 0.1 1 10 100
0 200 400 600 800 1000
max. residual 0.34*0.993^n
Figure 1. Evolution of maximal local error for a Laplace equation on a brick-shaped domain.
Distribution of errors. Both methods, classical and adaptive FSM, control errors at the collocation points only. Therefore there arises the question, how the errors behave inside the domain and on the boundary between the collocation points.
Fundamental solutions and the derived basis functions of the adaptive method are L-harmonic. Due to the maximum principle for the Laplace equations the error assumes its maximum on the boundary ofΩ. This can be illustrated in an error plot along a line segment through the domain (cf. Figure 2). The gradient and the Hessian errors show similar behavior.
Looking at the error on the boundary in case of the adaptive method, one observes a pattern of low error speckles (cf. Figure 3). These correspond to the collocation points where local errors have been eliminated. In this exam- ple, the singularity points were located on an ellipsoidal pseudo-boundary adapted to the domain’s shape. One can observe a higher concentration of speckles in regions located closer to the pseudo-boundary. This is in agree- ment with the theory of classical FSM, where a lower distance between the boundaries leads to higher stability, but to slower convergence [6].
5.2. Helmholtz equation. We perform experiments for the Helmholtz equation
∆v(x) +κ2v(x) = 0, κ= 2n, n= 1, . . . ,5 for x∈Ω, v(x) =g(x) for x∈∂Ω,
1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1
-0.4 -0.2 0 0.2 0.4
error
abs.err.
abs.err. gradient abs.err. Hessian
Figure 2. Absolute values of errors for Laplace problem on a brick-shaped domain measured on a line segment in- tersecting the surface at±12. Downward spikes correspond to the changes of sign.
Figure 3. Absolute error for Laplace problem on the boundary of a crankshaft domain.
Figure 4. Approximations of unit sphere with20,80and 320triangles, respectively.
with the known analytical solution v(x) =exp(
− κ
√2ix1
)
x2sin( κ
√2x3
)
, g=v
Γ.
HereΩis the unit ball; its surfaceΓis approximated by triangulated surface meshes. These meshes are obtained by a quasi-uniform refinement starting from an icosahedron (cf. Figure 4). The collocation points are derived as barycenters of the mesh triangles, singularity points are obtained by shifting the collocation points along the surface normal. Although Fundamental Solution Methods do not require an actual mesh, we will stick to this term, since the collocation points are derived from meshes.
Table 2. Helmholtz equation, varyingκ,N = 20480.
κ rel. err. rel. err. steps required classical FSM ada 10−11 (thres. 10−11) 1 4.57·10−13 1.42·10−11 976 2 2.92·10−13 7.84·10−12 1019 4 6.46·10−13 6.18·10−12 1102 8 1.43·10−12 5.85·10−12 1277 16 3.00·10−12 5.94·10−12 1792 32 1.80·10−11 2.59·10−11 3653(no thres.)
Performance of the adaptive FSM. As one could expect, for both, classi- cal and adaptive FSM, the quality of results gets worse with increasingκ
Table 3. Helmholtz equation,κ= 8,N = 1280.
threshold max.err. rel.err. condsys condLU # nodes classical 3.96·10−10 1.43·10−10 3.11·1013 - 1280
10−8 1.20·10−8 7.40·10−9 8.36·102 2.67·1010 721 10−9 1.43·10−9 8.69·10−10 9.18·102 2.73·1011 872 10−10 5.39·10−10 1.86·10−10 1.01·103 4.94·1012 1042 10−11 3.94·10−10 1.47·10−10 1.08·103 1.44·1014 1200 10−12 3.96·10−10 1.43·10−10 1.11·103 2.64·1014 1276 10−13 3.96·10−10 1.43·10−10 1.12·103 2.66·1014 1278
Table 4. Helmholtz equation,κ= 8,N = 20480.
threshold max.err. rel.err. condsys condLU # nodes classical 2.08·10−12 1.43·10−12 1.26·1021 - 20480
10−10 9.86·10−11 6.06·10−11 3.13·103 2.72·1012 1086 10−11 1.01·10−11 5.85·10−12 3.43·103 1.70·1013 1277 10−12 1.02·10−12 6.15·10−13 3.65·103 1.80·1014 1493 10−13 1.08·10−13 5.97·10−14 3.98·103 2.36·1015 1753 none 7.83·10−14 4.80·10−14 4.00·103 3.99·1015 1781
(cf. Table 2). While classical FSM suffers from a loss of approximation qual- ity, the adaptive method compensates for this by the use of a larger number of basis functions. For strict thresholds, the adaptive method achieves the accuracy of the classical method on coarser meshes (cf. Table 3) and even outperforms it on fine meshes (cf. Table 4). On such meshes, the classical FSM suffers from extremely high condition numbers condsys of the system matrix leading to a loss of accuracy.
Number of required collocation points. Of interest is the observation when the number of available collocation points increases, the number of steps required in the adaptive method to reach a certain threshold does not seem to grow (cf. Figure 5). As is seen from Table 5 on smaller clusters, where the threshold is reached faster, the results are worse. This is due to the fact, that on finer meshes the adaptive FSM has more points to choose from during the error elimination steps. Nevertheless, any threshold can be achieved theoretically by eliminating all (or almost all) errors at the collocation points. In this case, the adaptive method is equivalent to the full FSM.
Effects of large condition numbers. Figure 6 shows the loss of accuracy in the classical Fundamental Solution Method. When the number of colloca- tion points grows beyond a critical value, the error starts to grow slowly.
While this growth does not necessarily lead to very large errors, if a better
600 800 1000 1200 1400 1600 1800
3 3.5 4 4.5 5 5.5 6
# of iteration steps required for given threshold
icosa 2^n e-08
e-09 e-10 e-11 e-12 e-13
Figure 5. Number of iteration steps required to reach a certain threshold for a varying number of collocation points.
Table 5. Helmholtz equation,κ= 8, different geometries, threshold =10−12.
N steps required rel.err. rel.err. classical FSM 1280 1276 1.43·10−10 1.43·10−10 5120 1523 6.47·10−13 5.85·10−13 20480 1493 6.15·10−13 1.43·10−12
81920 1505 5.92·10−13 -
approximation is desired, one has to repeat the calculation with fewer col- location points. In the same figure, the growth of the condition number is indicated. To the same problem the adaptive FSM is applied (cf. Figure 7).
It can be seen that there exists a critical step after which the maximal lo- cal error will grow due to the loss of accuracy in floating point operations.
However, we can still use stored data from the previous steps in order to obtain a better result.
6. Conclusion and Outlook
When applied to large problems, the Fundamental Solution Method fea- tures system matrices with extremely large condition numbers. We have
1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
64 256 1024 4096 16384 65536 1
100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16 1e+18 1e+20 1e+22
relative error condition number
number of collocation points
relative error condition number
Figure 6. Full FSM: Condition number and loss of accu- racy for large numbers of collocation points.
1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000
0 500 1000 1500 2000 1e+06
1e+07 1e+08 1e+09 1e+10 1e+11 1e+12 1e+13 1e+14 1e+15 1e+16 1e+17
max. local error condition number
steps
|max. local error|
condition number
Figure 7. Adaptive FSM: Maximal local error, Helmholtz equation,κ= 8,N = 81920.
presented an adaptive method in which dimensions and the condition num- bers of the matrices involved are reduced by several orders. Numerical results show that the quality of approximations is comparable to that of the classical method. Also, the new method leads to reasonable results even in scenarios, where the classical method fails.
Future work will include the extension of the adaptive method to vector- valued problems with a special focus on elastostatics. We are also planning to investigate the convergence of the method in a theoretical context.
References
1. E. Anderson, Z. Bai, C Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LA- PACK users’ guide. Third edition,Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999.
2. M. Bebendorf, Approximation of boundary element matrices. Numer. Math. 86 (2000), No. 4, 565–589.
3. A. Bogomolny, Fundamental solutions method for elliptic boundary value problems.
SIAM J. Numer. Anal.22(1985), No. 4, 644–669.
4. V. D. Kupradze and M. A. Aleksidze, An approximate method of solving certain boundary-value problems. (Russian)Soobshch. Akad. Nauk Gruzin. SSR30(1963), 529–536.
5. R. Mathon and R. L. Johnston, The approximate solution of elliptic boundary-value problems by fundamental solutions.SIAM J. Numer. Anal.14(1977), No. 4, 638–650.
6. Y.-S. Smyrlis, Applicability and applications of the method of fundamental solutions.
Math. Comp.78(2009), No. 267, 1399–1434.
7. E. Trefftz, Ein Gegenstück zum Ritzschen Verfahren. In2ter Intern. Kongr. für Techn. Mech., pp. 131–137,Zürich, 1926.
(Received 09.06.2015) Authors’ address:
Department of Mathematics, Saarland University, P.O. Box 15 11 50, 66041 Saarbrücken, Germany.
E-mails: [email protected];[email protected];
