close

Вход

Забыли?

вход по аккаунту

?

52

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
THE METHOD OF FUNDAMENTAL SOLUTIONS
AND QUASI-MONTE-CARLO METHOD FOR
DIFFUSION EQUATIONS
C. S. CHEN1;∗ , M. A. GOLBERG2 AND Y. C. HON3
1
Department of Mathematical Sciences, University of Nevada, Las Vegas, NV 89154, U.S.A.
2 2025 University Circle, Las Vegas, NV 89119, U.S.A.
3 Department of Mathematics, City University of Hong Kong, Hong Kong
ABSTRACT
The Laplace transform is applied to remove the time-dependent variable in the diusion equation. For nonharmonic initial conditions this gives rise to a non-homogeneous modiÿed Helmholtz equation which we
solve by the method of fundamental solutions. To do this a particular solution must be obtained which we
ÿnd through a method suggested by Atkinson.17 To avoid costly Gaussian quadratures, we approximate the
particular solution using quasi-Monte-Carlo integration which has the advantage of ignoring the singularity in
the integrand. The approximate transformed solution is then inverted numerically using Stehfest’s algorithm.13
Two numerical examples are given to illustrate the simplicity and eectiveness of our approach to solving
diusion equations in 2-D and 3-D. ? 1998 John Wiley & Sons, Ltd.
KEY WORDS: method of fundamental solutions; particular solution; diusion equations; quasi-Monte-Carlo method;
Laplace transform
1. INTRODUCTION
In the past decade Boundary Element Methods (BEM) have emerged as one of the major numerical
techniques for solving Partial Dierential Equations (PDEs). The BEM has successfully avoided
domain discretization and integration for homogeneous PDEs. When the PDE is inhomogeneous,
a domain integration is introduced and diminishes the value of using the BEM. In recent years,
the dual reciprocity method1 and multireciprocity method2 have been developed to convert domain
integrals into equivalent boundary integrals. However, all these methods require more sophisticated
mathematical procedures and the BEM itself involves further numerical integrations with a singular
integrand.
To avoid the boundary integration, the Method of Fundamental Solutions (MFS) has gradually
evolved as an eective boundary free method for solving homogeneous PDEs. The MFS can be
regarded as an indirect BEM or element free BEM with singularities placed outside the boundary.
The MFS retains the same advantage of discretizing only the boundary but requires no numerical
integration at all. If the collocation method is employed, then only a system of linear equations
∗
Correspondence to: C. S. Chen, Department of Mathematical Sciences, University of Nevada, Las Vegas, NV 89154,
U.S.A. E-mail: chen@nevada.edu
Contract / grant sponsor: NATO; Contract / grant number: CRG 970593
CCC 0029–5981/98/081421–15$17.50
? 1998 John Wiley & Sons, Ltd.
Received 23 October 1997
Revised 18 March 1998
1422
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
needs to be solved. When the governing dierential equation is inhomogeneous, two approaches
have been successfully implemented to calculate particular solutions with excellent results:
(i) direct numerical approximation of domain integrals by the quasi-Monte-Carlo method.3; 4
(ii) dual reciprocity method to approximate the forcing term of the governing dierential equation
by using radial basis functions.5 – 7
Coupled with these methods, the MFS has been proven as an eective method in solving Poissontype dierential equations.
It is the purpose of this paper to develop a simple and eective method by extending the MFS
to solve diusion-type equations. In the past, various approaches have been adopted to treat timedependent diusion problems. A detailed review of these approaches is described in References
8 –10. Among these, Laplace transforms have been used to temporarily remove the time dependence; this results in an elliptic dierential equation in the Laplace transform (LT) space. A
standard BEM can then be employed to solve this elliptic dierential equation. Once the solution
in the LT space is found, a numerical inverse transform retrieves the solution in real-time space.
This approach was rst proposed by Rizzo and Shippy11 for certain problems of transient heat
conduction. Since the inverse LT is ill-posed, the results of Rizzo and Shippy11 were hampered
by a poor inversion method. As a result, the applicability of the LT has been limited in solving
diusion equations. In recent years its stability and accuracy has signicantly improved.12 Moridis
and Reddell9 and Zhu et al.10 have successfully implemented Stehfest’s algorithm13 for the numerical inversion of Laplace transform in the context of the BEM and the algorithm has gained
acceptance due to its improved stability and accuracy.
Recently, Zhu et al.10 have employed the dual reciprocity method to deal with the domain integral arising from solving elliptic dierential equations in the LT space. They used the Laplace
operator instead of the modied Helmholtz operator as the main operator.We propose, here, the
use of a modied Helmholtz operator in LT space and then the quasi-Monte-Carlo method introduced by Chen3 to nd a particular solution for the elliptic dierential equation. Gipson14 rst
proposed the Monte-Carlo method to deal with domain integration and singularities for Poissontype equations. It typically provides a simple and direct approach to solving simulation problems. However, the simplicity is often obtained at the expense of accuracy and eciency. So
the Monte-Carlo method has not been seriously considered as an alternative method for solving
PDEs.
In recent years, due to the rapid development of the quasi-Monte-Carlo method in theoretical
number theory, the accuracy of quasi-Monte-Carlo integration has been signicantly increased
over its counterpart Monte-Carlo method. Meanwhile, the increasing power of modern computer
hardware and parallel processing have greatly enhanced its performance and the cost of computing
can be considerably reduced. Thus, the cost of shifting external modelling time and mathematical
analysis to the computer is justied.
2. LAPLACE TRANSFORM FOR DIFFUSION EQUATIONS
Let be a bounded domain in R d ; where d is the dimension of the domain, with boundary
S = S1 ∪ S2 ; S1 ∩ S2 = ∅. We consider the following diusion equation:
1 @u(x; t)
= u(x; t);
k @t
? 1998 John Wiley & Sons, Ltd.
x ∈ ; t¿0
(1)
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
THE METHOD OF FUNDAMENTAL SOLUTIONS
1423
with boundary conditions
u(x; t) = f1 (x; t);
x ∈ S1 ; t¿0
(2)
@u(x; t)
= f2 (x; t);
@n
x ∈ S2 ; t¿0
(3)
x∈
(4)
and initial condition
u(x; 0) = u0 (x);
where n is the unit outward vector normal to S and u0 (x); f1 (x; t) and f2 (x; t) are given functions.
The diusion coecient k is assumed to be constant with respect to space and time.
To approximate the solution to equations (1) – (4), we rst employ the Laplace transform (LT)
to temporarily remove the time variable:
Z ∞
u(x; t)e−st dt
(5)
L[u(x; t)] = U (x; s) =
0
where the transform parameter s is real and positive. By integration by parts, we have
@u(x; t)
= sU (x; s) − u0 (x)
L
@t
By direct substitution of (5) and (6) into (1) – (3), we obtain
u0 (x)
s
U (x; s) = −
; x∈
−
k
k
(6)
(7)
U (x; s) = F1 (x; s);
x ∈ S1
(8)
@U (x; s)
= F2 (x; s);
@n
x ∈ S2
(9)
where Fi (x; s) = L[fi (x; t)]; i = 1; 2: Notice that equation (7) is a modied Helmholtz equation in
the LT space and can be reformulated as
Z
Z
@G(x; ^; s)
@U (x; s)
dS(^) − k G(x; ^; s)
dS(^)
U (x; s) = k U (x; s)
@n
@n
S
S
Z
(10)
+ u0 (x)G(x; ^; s) d
(^)
Here the fundamental solution G(x; ^; s) is given by

r s
1


for 2-D
K

0 r

2
k
G(x; ^; s) =
r 
s
1



for 3-D
exp −r
4r
k
(11)
where K0 is the modied Bessel function of the second kind and of order zero8 and r = k^ − xk
is the Euclidean distance between ^ and x.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1424
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
If the initial condition u0 (x) is harmonic, then the domain integral in (10) can be eliminated
by a change of variable. Let
v(x; t) = u(x; t) − u0 (x)
(12)
Then (1) – (4) become
1 @v(x; t)
= v(x; t);
k @t
x ∈ ; t¿0
v(x; t) = f1 (x; t) − u0 (x);
x ∈ S1 ; t¿0
@u0 (x)
@v(x; t)
= f2 (x; t) −
;
@n
@n
v(x; 0) = 0;
(13)
x ∈ S2 ; t¿0
x∈
By Laplace transformation, we have
s
V (x; s) = 0;
−
k
(15)
(16)
x∈
V (x; s) = F1 (x; s) −
(14)
(17)
u0 (x)
;
s
x ∈ S1
(18)
1 @u0 (x)
@V (x; s)
= F2 (x; s) −
; x ∈ S2
(19)
@n
s @n
where V (x; s) = L[v(x; t)]: Equation (17) is homogeneous and the standard BEM can be applied
without domain integration. After v in (13) – (16) has been computed, the nal solution can be
obtained by u = v + u0 in (12).
Recently, Zhu et al.10 applied the DRM with the Laplace operator as the main operator to
transfer the domain integration to the boundary. As they indicated, the formulation of the DRM
with the modied Helmholtz operator to transfer the domain integral to boundary is quite involved.
However, we believe that by using only the Laplace operator instead of the modied Helmholtz
operator as a whole in the process of the DRM, some information of the governing equation may
be lost and the forcing terms on the right-hand side of the Laplace operator become more dicult
to approximate. The results are expected to be less accurate. In the next section, a direct and
simple process will be implemented to avoid the diculties of domain integration and domain
discretization.
3. THE USE OF PARTICULAR SOLUTIONS IN LAPLACE SPACE
To approximate (7)–(9) numerically, it is convenient to nd a particular solution so that the
original equation can be reduced to an equivalent homogeneous equation. Let Up be a particular
solution of (7) in LT space. This means that Up satises equation (7) without satisfying the
boundary conditions; i.e.,
−u0 (x)
s
Up (x; s) =
(20)
−
k
k
Then
U (x; s) = Uh (x; s) + Up (x; s)
? 1998 John Wiley & Sons, Ltd.
(21)
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1425
THE METHOD OF FUNDAMENTAL SOLUTIONS
where Uh (x; s) satises the following homogeneous equation:
s
Uh (x; s) = 0; x ∈ −
k
Uh (x; s) = F1 (x; s) − Up (x; s);
(22)
x ∈ S1
(23)
@Up (x; s)
@Uh (x; s)
(24)
= F2 (x; s) −
; x ∈ S2
@n
@n
Equations (22)–(24) can be solved by boundary element methods such as the BEM or the MFS.
Notice that the homogeneous solution Uh (x; s) depends on the successful approximation of the particular solution. Hence, our proposed method in this paper relies heavily on the accurate evaluation
of the particular solution.
The classical method for obtaining a particular solution of (20) is to construct the associated
Newtonian potential; i.e.
Z
−u0 (^)
dV (^); x ∈ ∪ @
(25)
Up (x; s) = G(x; ^; s)
k
[Notice that the domain integration term in equation (10) is equivalent to the particular solution
in (25).] For complicated domains and the singularities of the fundamental solution, it is dicult
to evaluate (25) with direct numerical integration, especially for higher-dimensional cases.
Diculties of a complicated domain can be dealt with by domain imbedding and domain transformation based on the idea that the particular solution only satises the governing equation but
poses no restriction on boundary conditions. Atkinson17 developed a general numerical scheme to
evaluate particular solutions by extending the physical domain to a larger and regular ctitious
b under the assumption that the forcing term in (20) can be extended to b with the same
domain order of dierentiability as it had on : Atkinson asserted that
Z
−u0 (^)
dV (^); x ∈ ∪ S
(26)
Up (x; s) = G(x; ^; s)
k
b
is still a particular solution and its normal derivative is given by
Z
@Up (x; s)
@G(x; ^; s) −u0 (^)
=
dV (^); x ∈ ∪ S
@n
@n
k
b
(27)
b properly, the integral in (26) can be evaluated much easier than in (25).
By choosing b
By a change of variables,17 a further domain transformation can be performed by reducing b
to a rectangular box. In the 2-D case let us choose to be an ellipse whose boundary is given
b contains . Using the
by B() = (a cos ; b sin ); 0662; where a and b are chosen so that change of variables,
^ = x + (B() − x);
0661; 0662
and by elementary calculus, we obtain
Z 2 Z 1
ab
K(x; )b
u0 (x; ; ) d d
Up (x; s) =
2k 0
0
? 1998 John Wiley & Sons, Ltd.
(28)
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1426
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
where
K(x; ) = K0
r h
x
i
y
s
1−
cos + sin ;
kB() − xk
k
a
b
x = (x; y) ∈ ∪ S
b
u0 (x; ; ) = u0 (x + (B() − x); ; ):
Note that the integration domain in (28) is mapped into a rectangular box [0; 1] × [0; 2] and
the diculty of a complex domain is further eliminated. For the 3-D case, a similar domain
b be the ellipsoid with boundary
transformation to a rectangular parallelpiped can be achieved. Let B(; ) = (a cos sin ; b sin sin ; c cos );
066; 0662
By a change of variables,
^ = x + (B(; ) − x);
0661; 066; 0662
the integral in (26) becomes
1
Up (x; s) =
4k
where
Z
0
2
Z
0
Z
0
1
K(x; ; ; )b
u0 (x; ; ; ) d d d
(29)
√
e(−kB(; )−xk s=k) sin [abc − abz cos − acy sin sin − bcx sin cos ]
K(x; ; ; ) =
kB(; ) − xk
and b
u0 (x; ; ; ) = u0 (x + (B(; ) − x); ; ; ): As in the 2-D case, the integration domain
becomes [0; 1] × [0; ] × [0; 2] and the quasi-random nodes can be easily generated in this parb as an ellipsoid are sucient for computational purposes in
allepipe. However, (26)–(27) with the higher-dimensional case. In this case, some of the generated quasi-random nodes outside the
ellipsoid have to be discarded.
In contrast to the 2-D case, there is no singularity in the integrand in the 3-D case. To overcome
the diculty of the singularity in the integrand of (28), we employ the idea of quasi-Monte-Carlo
integration which is the subject of the next section.
4. CONSTRUCTION OF QUASI-RANDOM NUMBERS
The Monte-Carlo method has been used in multi-dimensional integration due to its simplicity. It
is a method that is independent of dimensionality, complexity of the integration domain and class
of functions for which it applied. However, the Monte-Carlo method only provides a probabilistic
error bound O(M −1=2 ); where M is the total number of random points.
From the integration point of view, the random distribution of nodes in the integration domain is
not as important as uniform distribution. Based on this fact, the quasi-Monte-Carlo method seeks to
construct a sequence of integration nodes in such a way that they are ‘spread evenly and precisely’
through the domain of integration. In other words, the only dierence between Monte-Carlo and
quasi-Monte-Carlo integration is the procedure for selecting integration nodes by means of sophisticated mathematical analysis so that these sample points maximally avoid each other. Furthermore,
since the integration nodes are predetermined, classical mathematical analysis is possible. As a
result, the quasi-Monte-Carlo integration yields a deterministic error bound O(M −1 (log M ) d−1 )
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1427
THE METHOD OF FUNDAMENTAL SOLUTIONS
where d is the dimension of integration domain.15 Due to its eectiveness, the quasi-Monte-Carlo
method has been developed and has evolved as a sophisticated technique in muti-dimensional integration. Recently, quasi-Monte-Carlo integration has been successfully implemented for solving
Poisson-type equations with high accuracy.3; 4
A common and satisfactory way of generating quasi-random nodes can be obtained by the socalled Hammersley point set.15 For a detailed mathematical analysis of the Hammersley point set,
we refer the reader to the book of Niederreiter.15 Here we simply give the basic construction of
quasi-random nodes.
Let Zb = {0; 1; : : : ; b − 1} be the least residue system mod b. Every integer n¿0 can be expressed as
n=
∞
P
j=0
aj (n)b j
(30)
in base b; where aj (n) ∈ Zb for all suciently large j: Based on this expansion, we dene the
radial-inverse function b as
b (n) =
∞
P
j=0
aj (n)b−j−1
for all integer n¿0
(31)
This leads to the denition of the Hammersley point set.
Denition 1. Let b1 ; b2 ; : : : ; bd−1 be integers ¿2 and relatively prime. Then the Hammersley
point set in the bases b1 ; b2 ; : : : ; bd−1 is the sequence x0 ; x1 ; : : : with
xn =
n
; b1 (n); : : : ; bd−1 (n) ∈ [0; 1) d
M
for all n¿0
(32)
The integration scheme of the quasi-Monte-Carlo method is the same as the Monte-Carlo method
except that the integration nodes are chosen dierently. The basic idea of the Monte-Carlo integrations is to nd the average of the integrand over the domain and then multiply by the area of
the integration domain. Hence, to evaluate the particular solution Up in (26), we have
Up (x;
s) ' UpM (x; s) =
b P
M
(
)
−u0 (xi )
;
G(x; xi ; s)
M i=1
k
b
xi ∈ (33)
b is the area in 2-D and volume in 3-D of b and M is the total number of integration
where (
)
nodes: The particular solution in (28) and (29) can be evaluated as follows:

M
ab P


K(x; i )b
u0 (x; i ; i )i ;
in 2-D

 kM
i=1
M
Up (x; s) ' Up (x; s) =
M

P


K(x; i ; i ; i )b
u0 (x; i ; i ; i )i ; in 3-D

2k i=1
(34)
where (i ; i ) ∈ [0; 1) × [0; 2) in 2-D and (i ; i ; i ) ∈ [0; 1) × [0; 2) × [0; ) in 3-D.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1428
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
5. THE MFS FOR THE MODIFIED HELMHOLTZ EQUATION
After we nd a particular solution of (20), we are ready to solve equations (22) – (24) by either
the BEM or MFS. We prefer the latter.
The MFS can be considered as a modied Tretz method for approximating the solution by
fundamental solutions and boundary data. In the MFS we place the source points outside the
physical domain in order to avoid the singularities of the fundamental solution of (11). Another
attractive feature of the MFS is its spectral convergence.18 For a detailed discussion of the MFS,
we refer the reader to Golberg7 and Bogomolny.19
Let be a ctitious boundary containing and assume that the solution Uh (x; s) of (22) – (24)
has an analytic continuation dened at every point within the region enclosed by and let N be
the number of collocation points on the boundary. The approximate solution UhN (x; s) to Uh can
be represented by
UhN (x; s) =
N
P
i=1
ai G(x; ^; s);
x ∈ S1 ∪ S2 ; ^ ∈
(35)
where {ai }Ni=1 have to be chosen to satisfy the boundary conditions. By collocation, we choose
N
{xj }m
j=1 on S1 and {xj }j=m+1 on S2 such that
N
P
i=1
N
P
i=1
ai
ai G(xj ; ^i ; s) = F1 (xj ; s) − UpM (xj ; s);
{xj }m
j=1 ∈ S1
@UpM (xj ; s)
@G(xj ; ^i ; s)
= F2 (xj ; s) −
;
@n
@n
{xj }Nj=m+1 ∈ S2
(36)
(37)
Equations (36) and (37) are a system of N equations and can be solved by Gaussian elimination.
Once {ai }Ni=1 is found, UhN (x; s) can be evaluated by equation (35) at any point x ∈ ∪ S:
b (x; s) to U (x; s) can be given by
After UhN and UpM are computed, an approximation U
b (x; s) =
U
N
P
i=1
ai G(x; ^; s) + UpM (x; s)
(38)
For a detailed discussion of how to choose source points on the ctitious boundary, we refer
the reader to Bogomonly.19
6. INVERSE LAPLACE TRANSFORM
Notice that the computation of UhN in (35) and UpM in (33) or (34) involves a parameter s in the
LT space. Among all the numerical Laplace transform inverse schemes, Stehfest’s algorithm13 has
been successfully implemented in the context of BEM.9; 10 Based on this algorithm, we numerically
inverted the solution from the LT space to the time domain. In the Stehfest’s algorithm, ns distinct
parameters s have been used to numerically invert the solution from the LT space to the time domain. For detailed implementation, we refer the reader to References 9 and 10 or Stehfest’s paper.13
The accuracy of Stehfest’s algorithm depends on the correct choice of ns , the (even!) number
of terms in approximating the solution in time domain. As ns increases, the accuracy improves
until round-o error becomes a factor and the accuracy declines. This is a common phenomenon
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
THE METHOD OF FUNDAMENTAL SOLUTIONS
1429
in practical numerical computation. The optimal ns has a signicant impact on the nal solution of
our proposed method. According to Stehfest’s test of his algorithm on 50 test functions with known
inverse Laplace transform, he concluded that the optimal value of ns is 10 for single precision
variables and 18 for double precision variables. Moridis and Reddell9 and Zhu et al.10 showed
that there is no signicant dierence in their results for ns between 6 and 16.
7. NUMERICAL EXAMPLES
Comparison of our method with those exploited by Moridis and Reddell9 and Zhu et al.,10 denoted
LTBE and LTDR, respectively, is revealed. Indeed, in the examples that follows, the problems with
chosen constants and test values are identical to those presented by the authors above. Numerical
results were produced on a personal computer (Pentium P5-120) with Microsoftware Fortran Compiler 5·1. Evaluation of K0 in (11) was done with Microsoft IMSL Mathematical and Statistical
Libraries.
Example 1. Consider the ow of heat through the rectangle [−a; a] × [−b; b]; a = b = 0·2; with
a unit initial temperature and boundaries kept at zero temperature. Comparison of the methods is based upon obervation of the absolute errors from the calculation of u(x; 0·025; 9000);
x ∈ [0; a]; k = 5·8 × 10−7 using the analytical solution given by Carslaw and Jaeger.20
For the MFS in Laplace space we choose 16 equally distributed collocation points on the
boundary and the same number of source points on a circle with centre at origin and radius
r = 3·0: Since the convergence rate of quasi-Monte-Carlo method is roughly O(M −1 (log M )); it
has been found that choosing 10006M 63000 to be adequate.4 So 1000 quasi-random nodes have
been generated to evaluate each particular solution.
For our proposed method the agreement of the temperature distribution between the analytical
and numerical solutions is excellent: The eect of increasing ns in Stehfest’s algorithm to at least
10 substantially decreases and stabilizes the absolute error essentially to the order of 10−4 of our
method throughout the interval of construction (see Figure 1). Whereas, it has been reported that
absolute error increases with increasing ns using either LTBE9 and LTDR.10 So far as the accuracy
is concerned, our approach is far superior to that of the LTBE and LTDR.
Since u0 is harmonic, a solution can be obtained without domain integration as shown in
Section 2. Absolute error has, again, been revealed (see Figure 2). Notice the order of magnitude of the absolute error in Figures 1 and 2 is the same.
A comparison of our method and Monte-Carlo methods is obtained by replacing the
Hammerseley points with truly random points. As shown in Figure 1 and Figure 3, our method
is substantially better than that of Monte-Carlo method. This is consistent with the comparison of
the deterministic error, O(M −1 log M ); of the quasi-Monte-Carlo integration and the probabilistic
error, O(M −1=2 ); of the Monte-Carlo integration. We also observed that even the Monte-Carlo
method outperformed the LTBE and LTDR in accuracy.
Figures 4 and 5 show the distribution of 1000 Hammersley points and random points. It is
remarkable that the eect of re-distributing the integration points has such an impact on the
accuracy of numerical integration.
Example 2. We extend Example 1 to a 3-D problem as given in Moridis and Reddell.9 In this
example, we consider the ow of heat in the box [−a; a] × [−b; b] × [−c; c], a = b = c = 0·2; with a
unit initial temperature and boundaries kept at zero temperature. Again, comparison is based upon
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1430
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
Figure 1. Eect of ns on the accuracy of the current method in 2-D
Figure 2. Absolute error of dierent ns without domain integration in 2-D
observation of the absolute error from the calculation of u(x; 0·025; 0·08; 9000), x ∈ [0; a]; using
the analytical solution given by Carslaw and Jaeger.20 For the MFS in Laplace space we choose
62 equally distributed collocation points on the surface of the box (see Figure 6) and the same
number of source points on the box with a = b = c = 3:
The error bound of quasi-Monte-Carlo integration is O(N −1 (log N ) d−1 ) where d is the dimension of integration domain (see Section 5). So to achieve similar accuracy to that in the previous
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
THE METHOD OF FUNDAMENTAL SOLUTIONS
1431
Figure 3. Eect of ns on the accuracy of the Monte-Carlo method
Figure 4. 1000 Hammersley points
example 2000 quasi-random nodes are chosen. Absolute errors is, indeed, on the same order of
magnitude (see Figure 7). Since u0 is, again, harmonic, a solution is obtained without domain
integration and reveals absolute error on the same order of magnitude (see Figure 8). Actually,
this is somewhat surprising, for, no domain integration is involved and the error of the particular
solution is not present. In Figures 7 and 8, we also observe that the absolute errors increase when
inceasing ns for x¡1·2 and behave in a opposite way for x¿1·2: This is in contrast to the results
shown in Figures 1 and 2. However, all of them converge as ns increases. In general, it is dicult
to nd the exact solution. We usually take the converged solutions as the approximate solutions.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1432
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
Figure 5. 1000 random points
Figure 6. Collocation points on each surface of the cube
Table I. Comparison of CPU time for Examples 1 and 2
CPU (s)
ns
2-D
3-D
8
10
12
14
16
11·97
14·87
17·80
20·66
23·50
8·07
10·16
12·19
14·17
16·25
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
THE METHOD OF FUNDAMENTAL SOLUTIONS
1433
Figure 7. Eect of ns on the accuracy of the current method in 3-D
Figure 8. Eect of ns without domain integration in 3-D
In comparing to the previous 2-D case, twice as many quasi-random points in evaluating the
domain integration and three times as many collocation points in the MFS are chosen. However,
CPU time in obtaining solutions at an interior point is less (see Table I). Slower process
p time in
the 2-D case is primarily due to repeated evaluation of the fundamental solution K0 (r s=k ): For
instance, for the case of ns = 10; the total number of evaluations of K0 in computing particular
solutions is 1·7 × 105 and 10×16 2 in the MFS.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
1434
C. S. CHEN, M. A. GOLBERG AND Y. C. HON
8. CONCLUSIONS
There are three major components in our proposed method: the quasi-Monte-Carlo integration for
nding particular solutions, the MFS and the numerical inverse Laplace transform. Each component
has been successfully implemented in the context of the BEM. Doubts about the low accuracy of
the Monte-Carlo method should be re-evaluated by converting to the quasi-Monte-Carlo method.
The only issue left for the Monte or quasi-Monte-Carlo method is their eciency. Since each
particular solution in our proposed method has been computed independently of the others, it can
be easily implemented in parallel. Furthermore, Davies et al.21 have also successfully implemented
a parallel algorithm in the time domain. Hence, our proposed method can be highly parallelized
and eciency is expected to be signicantly improved. This is one of the subjects of our current
research.
ACKNOWLEDGEMENTS
The rst author acknowledges sabbatical support from the University of Nevada, Las Vegas and
partial support from Professor C. A. Brebbia at Wessex Institute of Technology, U.K., during his
sabbatical leave (1996–1997). The rst author also acknowledges partial support of a collaborative
research grant from NATO under reference CRG 970593. The authors would also like to thank
Professor S. P. Zhu for providing computer code of Stehfest’s algorithm and referees’ valuable
comments.
REFERENCES
1. P. W. Partridge, C. A. Brebbia and L. C. Wrobel, The Dual Reciprocity Boundary Element Method, Computational
Mechanics Publications, 1992.
2. A. J. Nowak and C. A. Brebbia, ‘The multiple reciprocity method—a new approach for transforming BEM domain
integrals to the boundary’, Engng. Anal. Boundary Elem., 6, 164 –167 (1989).
3. C. S. Chen, ‘The method of fundamental solutions and the quasi-Monte-Carlo method for Poisson’s equation’, in
H. Niederreiter and P. Shiue (eds.), Lecture Notes in Statistics, Vol. 106, Springer, Berlin, 1995, pp. 158 –167.
4. C. S. Chen and M. A. Golberg, ‘A domain embedding method and the quasi-Monte-Carlo method for Poisson’s
equations’, in C. A. Brebbia, S. Kim, T. A. Osswald and H. Power (eds.), Boundary Elements XVII, Computational
Mechanics Publications, 1995, pp. 115 –122.
5. C. S. Chen, ‘The method of fundamental solutions for nonlinear thermal explosions’, Commun. Numer. Methods
Engr., 11, 675 – 681 (1995).
6. M. A. Golberg and C. S. Chen, ‘The theory of radial basis functions applied to the BEM for inhomogeneous partial
dierential equations’, Boundary Elem. Commun., 5, 57– 61 (1994).
7. M. A. Golberg, ‘The method of fundamental solutions for Poisson’s equation’, Engng. Anal. Boundary Elem., 16,
205 –213 (1995).
8. C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques, Springer, Berlin, 1984.
9. G. L. Moridis and D. L. Reddell, ‘The Laplace Transform boundary element (LTBE) method for the solution of
diusion-type equations’, Boundary Elements XIII, Vol. 1, Springer, Berlin, 1991, pp. 83 –97.
10. S. Zhu, P. Satravaha and X. Lu, ‘Solving linear diusion equations with the dual reciprocity method in Laplace space’,
Engng. Anal. Boundary Elem., 13, 1–10 (1994).
11. F. Rizzo and D. J. A. Shippy, ‘A method of solution for certain problems of transient heat conduction’, AIAA J., 8,
2004 –2009 (1970).
12. B. Davies and B. Martin, ‘Numerical inversion of the Laplace transform: a survey and comparison of methods’,
J. Comput. Phys., 33, 1–32 (1979).
13. H. Stehfest, ‘Algorithm 368: numerical inversion of Laplace transform’, Commun. ACM, 13, 47– 49 (1970).
14. G. S. Gipson, ‘The coupling of Monte-Carlo integration with boundary integral techniques to solve Poisson-type
problems’, Engng. Anal. Boundary Elem., 2, 138 –145 (1985).
15. Harald Niederreiter, Random Number Generation and Quasi-Monte-Carlo Methods, SIAM, CBMS 63, Philadelphia,
PA, 1992.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
THE METHOD OF FUNDAMENTAL SOLUTIONS
1435
16. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes, Cambridge University Press,
Cambridge, 1992.
17. K. E. Atkinson, ‘The numerical evaluation of particular solutions for Poisson’s equation’, IMA J. Numer. Anal., 5,
319 –338 (1985).
18. R. S. C. Cheng, Delta-trigonometric and spline methods using the single-layer potential representation, Ph.D. thesis,
University of Maryland, 1987.
19. A. Bogomolny, ‘Fundamental solutions method for elliptic boundary value problems’, SIAM J. Numer. Anal., 22,
644 – 669 (1985).
20. H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Oxford University Press, London, 1959.
21. A. J. Davies, D. Crann and J. Mushtaq, ‘A Parallel Implementation of the Laplace Transform in BEM’, in
C. A. Brebbia, J. B. Martins, M. H. Aliabadi and N. Haie (eds.), Boundary Elements XVIII, Computational Mechanics
Publications, 1996, pp. 213 –222.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 1421–1435 (1998)
Документ
Категория
Без категории
Просмотров
2
Размер файла
131 Кб
Теги
1/--страниц
Пожаловаться на содержимое документа