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)

1/--страниц