- Email: [email protected]
A Trefftz type method for time-dependent problems
Engineering Analysis with Boundary Elements 28 (2004) 13–21 www.elsevier.com/locate/enganabound
A Trefftz type method for time-dependent problems S.Yu. Reutskiy Laboratory of Magnetohydrodynamics, P.O. Box 136, Moskovski ave. 199, Kharkov 61037, Ukraine Received 29 October 2002; revised 27 May 2003; accepted 30 May 2003
Abstract This paper presents the application of a Trefftz type method with approximate trial functions for time-dependent problems of the parabolic type. At first it was suggested for elliptic equations only. The method was developed for PDEs with variable coefficients when a total differential operator of a problem can be represented as a sum of the one-dimensional ones. A comparison of the numerical solutions with analytic solutions is performed. q 2003 Elsevier Ltd. All rights reserved. Keywords: Trefftz type method; Time differencing; Particular solution; Approximate trial functions; Eigenfunctions
1. Introduction In this paper the numerical method proposed in Ref. [1] is extended onto the initial and boundary value problems (IBVP) for partial differential equations (PDEs) of the type:
›uðx;tÞ ¼ LðxÞ½uðx;tÞþhðx;tÞ; x [ V , Rd ; t $ 0; d ¼ 1;2;3: ›t ð1Þ Here V is a simply connected domain bounded by a simple closed curve ›V: There are two general approaches for solving such equations in the framework of boundary methods: (i) taking the Laplace transform in t and (ii) finite differencing in t. The first case results in the elliptic equation for the Laplace transform u^ ðx; sÞ : ^ sÞ; ðLðxÞ 2 sÞ½^uðx; sÞ ¼ 2uðx; 0Þ 2 hðx;
ð2Þ
which should be solved for a sequence ðsk ; k ¼ 1; …; Ns Þ: The parameters sk (complex in general) depend on the method of inversion used. Finite differencing in time transforms Eq. (1) to a sequence of coupled stationary equations. For example, if the Crank – Nicholson scheme is used, one gets the sequence: ðLðxÞ 2 pÞ½ujþ1 ðxÞ ¼ 2ðLðxÞ þ pÞ½uj ðxÞ 2 2hjþ1=2 ðxÞ:
ð3Þ
E-mail address: [email protected] (S.Yu. Reutskiy). 0955-7997/04/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00115-2
Here uj ðxÞ ¼ uðx;tj Þ; hjþ1=2 ðxÞ ¼ hðx; ðtj þ tjþ1 Þ=2Þ; tj ¼ jDt; Dt is a time step and p ¼ 2=Dt: These techniques were developed in [2 – 7]. In both cases on gets equations with non-zero right hand sides even for a homogeneous initial PDE Eq. (1). So, a problem of finding a particular solution arises when a boundary method is used to solve Eqs. (2) or (3). When L ¼ D the Laplacian, the classical method for obtaining a particular solution is to construct the associated Newtonian potential. Difficulties of a complicated solution domain can be dealt with e.g. by the Atkinson’s embedding technique [9]. The key idea of his method is to embed the solution domain into a larger and simpler Cartesian region V0 instead of the original integration domain V: Combination of this idea with the quasi-Monte Carlo method is performed in Refs. [3,10 –12]. Next approach—the dual reciprocity method (DRM) was developed in order to transform domain integrals to boundary ones [13 – 15]. The key idea is to approximate the right hand side f ðxÞ P of the PDE under consideration LðuÞ ¼ f by a series Nj¼1 aj wj with an appropriate basis functions so that L½Cj ¼ wj can be solved analytically. In the context of the problem considered significant results have been achieved in Refs. [4,5,16]. In the first work analytic formulae were given for Cj when wj was thin plate splines. These formulae proved developing efficient meshfree algorithms for solving diffusion equation in Ref. [4]. In Ref. [5] these results were generalized by using higher order splines as basis functions wj : The corresponding particular
14
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
solutions for D ^ x can be obtained without numerical integration. This is important from the point of view of practical calculations because search for the particular solution should be performed for each sk in Eq. (2) or repeated in each time-layer in Eq. (3). When the initial PDE is reduced to a homogeneous form any boundary technique can be applied. Recently the method of fundamental solutions (MFS) has gained great popularity among researchers. Recall that according to this technique an approximate solution of Eq. (3) is looked for in the form of a linear combination: ujþ1 ðxlq1 ; …; qK Þ ¼ vjþ1 ðxÞ þ
K X
qkjþ1 Fk ðxÞ:
ð4Þ
k¼1
Here vjþ1 ðxÞ is the particular integral in the j þ 1th time layer and the trial functions Fk ðxÞ satisfy exactly the homogeneous PDE ðLðxÞ 2 pÞ½Fk ¼ 0 but do not necessary satisfy boundary conditions. They are used to determine the unknown q kjþ1 : Similar to BEM, the MFS employs fundamental solutions as trial functions but their singular points are situated outside a solution domain. So, the MFS can be treated as a Trefftz-type method [17 – 21]. The techniques of applying the MFS to diffusion-type equations are described in Refs. [3,4,6] and in reviews [7,8]. Note that similar techniques have been developed recently to other branches of applied mathematics under other names: the charge simulation method, the superposition method, the source function method, the modified Trefftz method, etc. As an example of such technique is the method of discrete sources—a very effective technique for scattering problems in homogeneous mediums developed last few years [22]. Other references on this subject can be found in Ref. [8,23]. Through the paper we deal with the differencing in time technique. However, as it will be seen below the same approach also can be used together with the Laplace transform procedure. Similar to Atkinson’s approach, the method presented in this paper includes the embedding procedure. More precisely, the solution domain V is assumed to be embedded into a simple Cartesian region V0 ¼ ½21; þ1d ; d ¼ 1; 2; 3: If this is not the case originally, then appropriate translation and scaling operations may be performed to make it so. The sequence of PDE Eq. (3) is replaced by a new one:
from zero only inside some neighborhood of the source point jk : These functions are analogous in some sense to Dirac’s functions dðx 2 jk Þ which in fact are placed in the right hand sides of given PDEs when the MFS procedure is applied. Free parameters qk ; as usual, are determined from boundary conditions on ›V: The source functions are taken in the form of truncated series IðxljÞ ¼
þ
qjþ1 k Iðxljk Þ;
ð5Þ
k¼1
which differs from the original in an additional sum that is added to the P right hand side of the equation. The sum Kk¼1 qkjþ1 Iðxljk Þ appears due to the boundary nature of the method. It contains free parameters qk and dshaped source functions Iðxljk Þ which essentially differ
def
rn ðM; xÞwn ðjÞwn ðxÞ ¼
n¼1
M X
cn ðM; x; jÞwn ðxÞ;
ð6Þ
n¼1
over some orthogonal basis system wn ðxÞ: Note that n is the multiindex, e.g. in the 2D case n ¼ ðn1 ; n2 Þ; wn ðxÞ ¼ wn1 ðxÞwn2 ðyÞ; cn ðM; x; jÞ ¼ cn1 ðM; x; jÞcn2 ðM; x; hÞ; and it is a double series. It is easy to see that series Eq. (6) coincides with the expansion of Dirac’s functions dðx 2 jÞ over the orthogonal system wn ðxÞ up to the coefficients rn ðM; xÞ ¼ rn1 ðM; xÞrn2 ðM; xÞ: The regularizing coefficients rn ðM; xÞ depend on the particular choice of the system wn ðxÞ [24,25]. Let us dwell in brief on the problem of the sources location. We use the technique in which the coordinates of the sources jk are fixed. Another approach leads to an adaptive version of the MFS when jk are considered as unknown values and should be determined together with the coefficients qjþ1 by the use of a non-linear k minimization process to match the boundary conditions [8]. Bogomolny [27] showed that the fixed sources can be distributed equally on an auxiliary circle in 2D case (sphere in 3D). As indicated by Bogomolny, the larger the radius of the circle (sphere), the better the approximation can be expected. Other results and references on this subject can be found in Ref. [23]. In all the calculations presented in this work the source points jk are placed inside V0 on an auxiliary equidistant contour as far as possible from the boundary ›V: In the 1D case (Examples 1 and 2 in Section 2, Examples 6 in Section 3) the source points j1 ; j2 are placed near the endpoints of the interval V0 ¼ ½21; þ1: The method of this paper is developed for PDEs with differential operators which can be represented as a sum:
ðLðxÞ 2 pÞ½ujþ1 ¼ 2 ðLðxÞ þ pÞ½uj 2 2hjþ1=2 ðxÞ K X
M X
LðxÞ ¼
d X
li ðxi Þ;
i¼1
where each one-dimensional operator li has the form: 1 › › pi ðxi Þ 2 qi ðxi Þ : li ðxi Þ ¼ ›x i ri ðxi Þ ›xi
ð7Þ
The equations like Eq. (1) with such L often arise and are widely used in applied researches and in mathematical
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
modeling. When the solution domain V is a regular region (e.g. the rectangular one) and there exists an exactly known complete system of eigenfunctions for each one-dimensional li ðxi Þ; the problem can be solved using an analytical procedure, e.g. by the method of separation of variables, the methods of integral transforms, etc. See e.g. Ref. [28] for more details. In a general case it can be solved only numerically. We assume that li ðxi Þ is a negative defined in the solution domain, i.e.
ri ðxi Þ; pi ðxi Þ . 0; qi ðxi Þ $ 0;
ð8Þ
when
Examples: (1) For l1 ½w ¼ d2 w=dx2 and the boundary conditions wð21Þ ¼ wðþ1Þ ¼ 0; one gets: xþ1 np 2 ð1Þ wð1Þ ðxÞ ¼ sin np m ¼ ; ; n n 2 2 ð11Þ p ð1Þ g ¼ : 2 (2) The second system w xþ1 x sin np wð2Þ ðx; wÞ ¼ exp ; n 2 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi w2 np 2 p ð2Þ þ ; n ¼ 1; …; 1; gð2Þ ¼ ; mn ¼ 2 2 4
21 , minðxi lx [ VÞ ¼ ai # xi # bi ¼ maxðxi lx [ VÞ , þ1:
corresponds to the 1D convection – diffusion operator
Go back to the choice of the regularizing coefficients rn ðM; xÞ: Let {wn ðxÞ; mn } be a solution of a Sturm –Liouville problem
l2 ðx; wÞ ¼
l½w¼2mw; ða1 ›x þ b1 Þwð21Þ¼0; ða2 ›x þ b2 Þwð1Þ¼0;
ð9Þ
a21;2 þ b21;2 –0; a1 b1 #0; a2 b2 $0; where the operator lðxÞ is given in Eq. (7). Under the assumptions (8) it satisfies the following conditions:
15
›2 › 2w ; 2 ›x ›x
ð12Þ
ð13Þ
ð14Þ
pðxÞ ¼ rðxÞ ¼ expð2wxÞ; qðxÞ ¼ 0; ð2Þ and the boundary conditions wð2Þ n ð21Þ ¼ wn ðþ1Þ ¼ 0: (3) The Bessel’s operator
l3 ðxÞ ¼
›2 1 › ; rðxÞ ¼ pðxÞ ¼ x; qðxÞ ¼ 0; þ 2 x ›x ›x
ð15Þ
together with conditions 0, m1 , m2 ,…, mn !þ1;
wn ð0Þ is bounded; wn ð1Þ ¼ 0;
and there exists the limit:
forms the system of the eigenfunctions: pffiffi 2J0 ðlð3Þ n xÞ ð3Þ ð3Þ 2 wn ðxÞ ¼ ; mð3Þ n ¼ ðln Þ ; g ¼ p; ð3Þ J1 ðln Þ
pffiffiffiffi mn : n!1 n
g ¼ lim
The eigenfunctions wn ðxÞ form an orthogonal system on [2 1, þ 1] with a scalar product: ðþ1 rðxÞwn ðxÞwm ðxÞdx¼gn dn;m : 21
In this case the regularizing coefficients are [24]: 1 sin nðn;MÞ ; ½sn ðMÞx ; sn ðMÞ¼ gn nðn;MÞ pffiffiffiffi p mn : nðn;MÞ¼ gðMþ1Þ rn ðM; xÞ¼
ð10Þ
where mn ; gn ; g correspond to the Sturm –Liouville problem considered and x is the parameter of regularization. For the sake of simplicity, further we omit the arguments M and x in notation of the source functions and their coefficients. A multi-dimensional source function Iðx;jÞ like Eq. (6) can be taken as a product of the one-dimensional ones. We begin with the simplest case when each onedimensional operator (7) possesses a set of eigenfunctions which are known in an explicit analytic form.
ð16Þ
where lð3Þ n is nth zero of the Bessel’s function J0 ðxÞ: These three one-dimensional operators and their combinations will be used in all the numerical examples in the next section to form a total differential operator of the PDE under consideration LðxÞ: The corresponding eigenfunctions wðiÞ n ðxÞ and their products are used as a basis system for an approximate solution. However, in most of practical problems, the corresponding differential operators usually do not have a known spectrum. As it will be shown in Section 3, an approximate solution of the Strum – Liouville problems Eq. (9) can be used in the algorithm. In each section we present results of the numerical experiments justifying the method presented.
2. Differential operators with a known spectrum As it mentioned above the functions wðiÞ n ðxÞ and their products form a natural basis in the problems considered. The trial functions in Eq. (4) satisfy the equation: ðL 2 pÞ½Fk ðxÞ ¼ Iðx; jk Þ; jk [ V0 \ V:
ð17Þ
16
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Due to the special choice of the basis system the solution can be written in the explicit form:
Fk ðxÞ ¼
M X
DðkÞ n wn ðxÞ;
DðkÞ n ¼2
n¼1
cn ðM; x; jk Þ ðjÞ mðiÞ n1 þ mn2 þ p
:
ð18Þ
Suppose a two dimensional case L ¼ li þ lj is considered here. The similar formulae can be written in the 1D and 3D cases. To get a particular solution the function hjþ1=2 ðxÞ in the right hand side of Eq. (5) is approximated by the truncated series over the same basis: hjþ1=2 ðxÞ ¼
M X
Hnjþ1=2 wn ðxÞ:
ð19Þ
the expansion for all the space functions fi ðxÞ: Next,Pto get Hnjþ1=2 at each time-layer as a sum: Hnjþ1=2 ¼ Ni¼1 ai ðtjþ1=2 ÞFi;n : To demonstrate the effectiveness of algorithm we give some numerical examples. Example 1 Consider the following problem:
›u ›2 u ¼ 2 þ hðx; tÞ; ›t ›x uðaÞ ¼ Ua ðtÞ;
uj ðxÞ ¼
M X
Unj wn ðxÞ:
ð20Þ
n¼1
The first expansion of such type can be obtained using the same procedure as the one for hjþ1=2 ðxÞ: Under these conditions: the particular solution vjþ1 ðxÞ (see Eq. (4)) can be written in the same form: vjþ1 ðxÞ ¼
M X
Vnjþ1 wn ðxÞ;
n¼1
Vnjþ1 ¼
ð21Þ j jþ1=2 2 mðjÞ n2 ÞUn þ 2Hn ðjÞ ðiÞ mn1 þ mn2 þ p
mðiÞ n1
ðp 2
;
after determining the free parameters qkjþ1 the solution on the new j þ 1th time layer ujþ1 ðxÞ can be written in the same form as uj ðxÞ : ujþ1 ðxÞ ¼
M X
Unjþ1 wn ðxÞ;
n¼1
Unjþ1
¼
Vnjþ1
þ
K X
ð22Þ qkjþ1 DðkÞ n :
k¼1
So, the code is closed inside a finite set of the basis functions ½wn ðxÞM n¼1 : Some remarks to the practical use of the algorithm. When hðx; tÞ ; 0; i.e. a homogeneous PDE is considered, then the algorithm requires a single Fourier expansion for the initial function uðx; 0Þ to get particular solutions on the all time layers. P When the function hðx; tÞ is a sum of products hðx; tÞ ¼ Ni¼1 ai ðtÞfi ðxÞ; or can be approximated well by such sum, it is more suitable to perform the algorithm in the following way. First, to get the coefficients Fi;n of
uðbÞ ¼ Ub ðtÞ;
uðx; 0Þ ¼ u0 ðxÞ;
t $ 0;
ð23Þ
t $ 0;
ð24Þ
x [ ½a; b , ½21; þ1:
ð25Þ
The source functions take the form:
n¼1
This series can be obtained using the Atkinson’s technique [9] or in another appropriate way. In particular we use the collocation method described in Refs. [1,25]. We also suppose the solution at the jth time layer uj ðxÞ already is represented in the same form of expansion over wn ðxÞ :
x [ ½a; b;
Iðx; jÞ ¼
M X
ð1Þ rn ðM; xÞwð1Þ n ðjÞwn ðxÞ;
n¼1
where the basis functions are given in Eq. (11). In this case the regularizing coefficients are: rn ðM; xÞ ¼ ½sn ðMÞx ; vðn; MÞ ¼
sn ðMÞ ¼
sin½vðn; MÞ ; vðn; MÞ
np : Mþ1
Here sn ðMÞ are the Lanczos sigma-factors which are used to overcome the so-called Gibb’s phenomenon in the Fourier expansion of non-smooth functions [26]. Let the solution domain be the interval ½a; b ¼ ½20:5; þ0:5: The right hand side hðx; tÞ; the boundary functions Ua ðtÞ; Ub ðtÞ and the initial function u0 ðxÞ correspond to the exact solution. uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ:
ð26Þ
Following the rule formulated in Section 1, the two source functions Iðx; jk Þ have parameters j1 ¼ 20:95; j2 ¼ 0:95: The maximal absolute error: ea ¼ max luex ðxi ; tj Þ 2 uðjÞ M ðxi Þl; {xi ;tj }
ð27Þ
as a function of the parameters M and x is presented in Table 1. Here: tj ¼ 0:5; 1:0; …; 10 are the 20 checking time layers and xi ¼ 0:01ði 2 0:5Þ 2 0:5; i ¼ 1; …; 100 are Table 1 The maximal absolute error in solution of the equation: ›u=›u ¼ ›2 u= ›x2 þ hðx; tÞ; uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:05
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
1 2 3 4 5 6 7
Div Div 3.7 £ 1023 2.2 £ 1022 4.7 £ 1022 7.5 £ 1022 0.1
Div Div Div 4.0 £ 1024 1.3 £ 1025 5.6 £ 1025 4.4 £ 1024
Div Div Div 4.7 £ 1023 3.5 £ 1025 1.3 £ 1025 1.3 £ 1025
Div Div Div 1.2 £ 1023 2.0 £ 1025 1.3 £ 1025 1.3 £ 1025
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
the checking points. All the calculations presented in the table were performed with the time step Dt ¼ 0:05: One can see that for each M there exists a minimal xdiv such that the solution process diverges for x , xdiv and it converges for x $ xdiv : The minimal value of ea ðM; xÞ is achieved for xm , xdiv and ea ðM; xÞ grows slightly when x increases. As a function of M the error ea ðM; xm Þ first decreases from ea ð10; 3Þ ¼ 3:7 £ 1023 to ea ð20; 5Þ ¼ 1:3 £ 1025 but then it does not change. This can be explained by the fact that the grows of M can reduce the errors due to approximation of the right hand side hjþ1=2 ðxÞ and the errors due to approximation of the trial functions Fk ðxÞ but is does not affect the error of the time discretization. So, it seems that 1.3 £ 1025 is the error introduced by the Crank – Nicholson scheme which is dominant for these values of M and x: This error should be a square function of the time step Dt: To verify this supposition the calculations were repeated with the time step Dt ¼ 0:01: The results are placed in Table 2. One can see that the behavior of ea ðM; xÞ here is the same as in the previous example but the minimal value of error now is: ea ð30; 8Þ ¼ 5:2 £ 1027 : The ratio 1:3 £ 1025 =5:2 £ 1027 ¼ 25 ¼ ð0:05=0:01Þ2 confirms the nature of the error. So, 5.2 £ 1027 is a new error due to the discretization in time and it cannot be reduced by the increase of the number of harmonics M: Example 2. The similar calculations were performed for the 1D convection –diffusion PDE.
›u ›2 u ›u ¼ 2w þhðx;tÞ; x[½20:5;þ0:5; t $0; ›t ›x 2 ›x
ð28Þ
with the same exact solution (26). The basis functions (12) are used here. The results correspond to w¼5 and the time step Dt ¼0:05 are placed in Table 3. Behavior of the error is the same as above. Here the errors value 1.5 £ 1024 is caused by the Crank – Nicholson time scheme. It decreases to 6.0 £ 1028 when one takes ðM; xÞ¼ð40;11Þ and Dt ¼0:001: Table 2 The maximal absolute error in solution of the equation: ›u=›t ¼ ›2 u= ›x2 þ hðx; tÞ; uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:01
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
3 4 5 6 7 8 9 10 11
3.7 £ 1023 2.2 £ 1022 4.7 £ 1022 7.5 £ 1022 0.1 0.1 0.2 0.2 0.2
Div Div Div 5.5 £ 1025 4.4 £ 1024 1.4 £ 1023 3.1 £ 1023 5.5 £ 1023 8.7 £ 1023
Div Div Div Div Div 5.2 £ 1027 1.5 £ 1026 1.1 £ 1025 4.0 £ 1025
Div Div Div Div Div Div Div 5.2 £ 1027 5.2 £ 1027
17
Table 3 The maximal absolute error in solution of the equation: ›u=›t ¼ ›2 u= ›x2 2 w›u=›x þ hðx; tÞ; uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:05
x
M ¼ 10
M ¼ 20
M ¼ 30
1 2 3 4 5 6 7
Div 3.8 £ 1021 1.5 £ 1023 8.9 £ 1023 2.0 £ 1022 3.5 £ 1022 6.8 £ 1022
Div Div Div 2.4 £ 1024 1.5 £ 1024 1.6 £ 1024 3.0 £ 1024
Div Div Div Div Div 1.5 £ 1024 1.5 £ 1024
Example 3. Consider the heat diffusion in an axially symmetric body with a distributed heat source ›u 1 › ›u ›2 u r 2 þ z2 ¼ : þ 2 þ12 r ›r ›r ›t ›z a2 pffiffiffiffiffiffiffiffiffi Here r ¼ x2 þ y2 ; z are cylindrical coordinates, the solution domain in rz-plane is half of a disk: V ¼ {ðr; zÞ : r 2 þ z2 # a2 ; r $ 0} So, the body is a solid sphere with the radius a: The boundary and initial condition are: uðr; z; tÞ ¼ 0;
r 2 þ z2 ¼ a2 ; t $ 0;
uðr; z; 0Þ ¼ 0;
i.e. the initial temperature and the temperature of the surface is zero. The differential operator of the problem is a sum of the one-dimensional ones l3 ðrÞ and l1 ðzÞ: So the twoð1Þ dimensional basis functions are: wn ðr; zÞ ¼ wð3Þ n1 ðrÞwn2 ðzÞ (see Eqs. (11) and (16)). The problem has an analytic solution [29]. ! ! 1 r2 r2 12 2 uðr; z; tÞ ¼ 723 2 60 a a 1 12a3 X ð21Þnþ1 sin an r 2a2n t e ; ð29Þ rp5 n¼1 n5 pffiffiffiffiffiffiffiffiffi where r ¼ r 2 þ z2 ; an ¼ np=a: In Table 4 this solution is compared with the one obtained by the algorithm described above. All the calculations were performed with Dt ¼ 0:05:
2
Table 4 The maximal absolute error in solution of the 2D equation: ›u 1 › ›u ›2 u r 2 þ z2 ¼ ; þ 2 þ12 r ›r ›r ›t ›z a2 the exact solution is given in Eq. (29), the time step is Dt ¼ 0:05 t
M ¼ 10; x ¼ 2
M ¼ 20; x ¼ 6
M ¼ 30; x ¼ 8
0.1 0.5 1.0 5.0 10.0
3.3 £ 1022 5.0 £ 1022 5.0 £ 1022 5.5 £ 1022 5.5 £ 1022
1.1 £ 1023 8.1 £ 1025 8.3 £ 1025 8.0 £ 1025 8.1 £ 1025
1.1 £ 1023 1.1 £ 1025 6.6 £ 1026 6.1 £ 1026 6.1 £ 1026
18
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Table 5 The maximal absolute error in solution of the 3D PDF:
3. General case
›u ›2 u ›2 u ›2 u x2 þ y2 þ z2 ¼ 2 þ 2 þ 2 þ12 ; Dt ¼ 0:05 ›t ›x ›y ›z a2
For the sake of simplicity consider a one-dimensional equation
t
M ¼ 10; x ¼ 2
M ¼ 20; x ¼ 6
M ¼ 30; x ¼ 8
0.1 0.5 1.0 5.0 10.0
1.1 £ 1023 2.2 £ 1024 2.2 £ 1024 2.3 £ 1024 2.3 £ 1024
9.8 £ 1024 1.7 £ 1025 1.3 £ 1025 1.2 £ 1025 1.2 £ 1025
9.7 £ 1024 9.2 £ 1026 1.8 £ 1026 1.3 £ 1026 1.2 £ 1026
›u ¼ lðxÞ½u þ f ðx; tÞ; ›t
ð30Þ
with lðxÞ given in Eqs. (7) and (8). The two basis systems
cn ðxÞ ¼ cosð0:5ðn 2 1Þpðx þ 1ÞÞ;
ð31Þ
wn ðxÞ ¼ sinð0:5npðx þ 1ÞÞ; Example 4. Consider the same problem as a 3D one in Cartesian coordinates x; y; z :
›u ›2 u ›2 u ›2 u x2 þ y2 þ z2 ¼ 2 þ 2 þ 2 þ12 : ›t ›x ›y ›z a2 The basis system is formed by the products ð1Þ ð1Þ wð1Þ n1 ðxÞwn2 ðyÞwn3 ðzÞ; K ¼ 160 source points are distributed uniformly on the auxiliary sphere of the radius 0.95, Dt ¼ 0:05: The results are placed in Table 5. Example 5. Consider 2D time-dependent convection– diffusion equation:
are used in this section to approximate: (i) the space operator lðxÞ : 2cn ðxÞ and (ii) solutions of the corresponding Sturm – Liouville problems: 2wn ðxÞ: The approximation of lðxÞ pursues the following two goals: (1) to extend the initial PDE from region V to V0 : Note that in one-dimensional case this means a simple extension from the interval ½a; b to the bigger one [2 1, þ 1]. However, for a multi-dimensional PDE the initial irregular solution domain is replaced by a simple Cartesian one. (2) The second goal is to transform the initial PDE into the form convenient for applying the Galerkin projection procedure. The approximations are sought in the form of series: r ð0Þ ðxÞ ¼
›u ›2 u ›2 u ›u ›u ¼ 2 þ 2 2 wx 2 wy þ hðx; y; tÞ: ›t ›x ›y ›x ›y
N X
Ri ci ðxÞ;
i¼1
qð0Þ ðxÞ ¼
ð2Þ The basis system is wð2Þ n1 ðx; wx Þwn2 ðy; wy Þ: The solution domain in xy-plane is a disk: V ¼ {ðx; yÞ : x2 þ y2 # a2 }: The boundary condition of the Dirichlet type, the initial condition and the right hand side hðx; y; tÞ correspond to the exact solution: uex ðx; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ expð2yÞÞ cos t: Table 6 displays the maximal absolute error for the two cases: (1) wx ¼ wy ¼ 5; (2) wx ¼ 15; wy ¼ 215: The pairs ðM; xÞ are placed in parentheses. The time step Dt ¼ 0:01 and K ¼ 40 source points are used in both cases.
N X
pð0Þ ðxÞ ¼
N X
Pi ci ðxÞ;
i¼1
Qi ci ðxÞ:
i¼1
To get the unknown coefficients Ri ; Pi ; Qi we write the collocation conditions P for some inner points xj [ ½a; b; j ¼ 1; …; N $ N; e.g. Ni¼1 Ri ci ðxj Þ ¼ rðxj Þ and then solve the system by an appropriate numerical procedure. Note that the Atkinson’s technique mentioned above also can be used there. As a result one gets a new differential operator 1 › ð0Þ › ð0Þ l ðxÞ ¼ ð0Þ p ðxÞ 2 q ðxÞ ; ›x r ðxÞ ›x ð0Þ
Table 6 The maximal absolute error in solution of the 2D convection-diffusion equation: 2
2
›u › u › u ›u ›u ¼ 2 þ 2 2 wx 2 wy þ hðx; y; tÞ; ›t ›x ›y ›x ›y exact solution Dt ¼ 0:001 t
0.5 1.0 1.5 2.0
uex ðx; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ expð2yÞÞcos t;
wx ¼ 5; wy ¼ 5
wx ¼ 15; wy ¼ 215
(20,5)
(30,8)
(10,3)
(20,5)
(30,8)
3 £ 1022 2 £ 1022 3 £ 1023 1 £ 1022
3 £ 1023 2 £ 1023 3 £ 1024 1 £ 1023
5 £ 1026 3 £ 1026 6 £ 1027 2 £ 1026
1 £ 10þ3 7 £ 10þ2 1 £ 10þ2 5 £ 10þ2
3 £ 1023 2 £ 1023 3 £ 1024 1 £ 1023
3 £ 1024 1 £ 1024 2 £ 1025 1 £ 1024
ð33Þ
defined in [2 1, þ 1] which approximates lðxÞ when x [ ½a; b: We assume that lð0Þ ðxÞ continues to be a negative defined operator, i.e. the conditions (8) are true for r ð0Þ ðxÞ; pð0Þ ðxÞ; qð0Þ ðxÞ: Let us consider the Sturm – Liouville problem on [2 1, þ 1]: lð0Þ ðxÞ½v ¼ 2mv;
(10,3)
ð32Þ
vð21Þ ¼ vðþ1Þ ¼ 0:
ð34Þ
An approximate solution of the problem can be looked for in the form of Mth term expansion:
vðxÞ ¼
M X n¼1
W n wn ðxÞ:
ð35Þ
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Projecting with the weight r ð0Þ ðxÞ on the finite basis w1 ðxÞ; …; wM ðxÞ one gets the sequence of equations: klð0Þ ðxÞ½v þ mv; wn ðxÞl ¼ 0;
ð36Þ
where the scalar product is: ðþ1 r ð0Þ ðxÞuðxÞvðxÞdx: ku; vl ¼
ð37Þ
19
when lm 2 nl . N; i.e., when the number of harmonics M exceeds the number of terms in the expansion of the coefficients, then the matrices have a band-like form. It is well-known that under the assumptions mentioned m }; above problem (38) has exactly M solutions {mm ; w m ¼ 1; …; M such that the eigenvectors w m and corresponding functions
21 M X
Substituting Eq. (35) in Eq. (36) one gets the eigenvalue problem:
vm ðxÞ ¼
^w ^ w; S· ¼ mR·
satisfy the conditions [30]:
ð38Þ
kw k; w m lR ¼
w ¼ ½W ; …; W ; M T
ðþ1 21
r ð0Þ ðxÞwn ðxÞwm ðxÞdx:
ð40Þ
Substitution Eq. (32) in Eqs. (39) and (40) yields: Sn;m ¼
N X
½0:25nmp2 Pi Z1 ði 2 1; n; mÞ þ Qi Z2 ði 2 1; n; mÞ;
i¼1
N X
Ri;j Wki Wmj ¼ dkm ;
kvk ; vm l ¼ dkm :
ð46Þ
The first condition means orthogonality of the vectors w k with the weight matrix Ri;j and the second one means orthogonality of the functions vk ðxÞ in the sense of the scalar product (37). In particular orthogonality of the system vm ðxÞ; m ¼ 1; …;M means that a series f ðxÞ ¼ P M wn ðxÞ can be n¼1 Fn wn ðxÞ over the basis functions P v further re-expanded over vk ðxÞ : f ðxÞ ¼ M m¼1 Fm vm ðxÞ: The coefficients of both expansion Fn and Fmv are connected by the relations: Fmv ¼
M X
Tm;n Fn ;
n¼1
ð41Þ Rn;m ¼
M X i;j¼1
^ are positive defined matrices with the entries: S^ and R ðþ1 dw ðxÞ dwm ðxÞ Sn;m ¼Sm;n ¼ pð0Þ ðxÞ n dx dx 21 þ qð0Þ ðxÞwn ðxÞwm ðxÞ dx; ð39Þ Rn;m ¼ Rm;n ¼
ð45Þ
n¼1
where 1
Wmn wn ðxÞ;
Ri Z2 ði 2 1; n; mÞ;
ð42Þ
Tm;n ¼
M X
Fn ¼
M X
Wmn Fmv ;
m¼1
ð47Þ
Wmk Rn;k :
k¼1
i¼1
The following property of the v-basis plays a central role in the algorithm [24]:
where Z1 ði; n; mÞ ¼ ¼
ðþ1 21
ci ðxÞcn ðxÞcm ðxÞdx
klð0Þ ðxÞ½vm ðxÞ; vn ðxÞl ¼ 2mn dmn :
1 ½Dði þ n 2 mÞ þ Dði þ m 2 nÞ þ Dðn 2
þ m 2 iÞ þ Dði þ n þ mÞ; ðþ1 Z2 ði; n; mÞ ¼ ci ðxÞwn ðxÞwm ðxÞdx
ð43Þ
(
DðkÞ ¼
where Iðx; jÞ is given in Eq. (6). Re-expanding it over the v-basis and projecting to each vn ðxÞ one gets
1 ½Dði þ n 2 mÞ þ Dði þ m 2 nÞ 2 Dðn 2 þ m 2 iÞ 2 Dði þ n þ mÞ;
1
if k ¼ 0
0
otherwise
In particular it is possible to write an approximation of the trial function Fk ðxÞ using the v-basis. Let us consider the PDE for Fk ðxÞ : ðlð0Þ ðxÞ 2 pÞ½Fk ðxÞ ¼ Iðx; jk Þ;
21
¼
ð48Þ
ð44Þ
:
^ are It is important to note that entries of the matrices S^ and R calculated in an analytic form without any numerical integration due to the special approximation (32). In general case M 2 numerical integrations are needed to get the matrices ^ Moreover, it easy to prove that Sn;m ¼ 0; Rn;m ¼ 0 S^ and R:
Fk ðxÞ ¼
M X n¼1
dnv ðjk ;M; xÞvn ðxÞ;
dnv ðjk ;M; xÞ ¼ 2
cvn ðj;M; xÞ ; mn þ p
where cvn ðj;M; xÞ are obtained from cn ðj;M; xÞ using linear transformation (47). In the similar way in the 2D case one gets an analog of Eq. (18) where cn ðM; x; jk Þ is replaced by cvn ðM; x; jÞ ¼ cvn1 ðM; x; jÞcvn2 ðM; x; hÞ: In general, all the formulae derived in the previous section for differential operators with a known spectrum are valid when one takes the expansion over the v-basis.
20
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Table 7 The maximal absolute error in solution of the equation: ›u › ›u ¼ ex cosðxÞ 2 x2 u þ hðx; tÞ; ›x ›t ›x
Table 8 The maximal absolute error in solution of the equation: ›u › ›u ›2 u ¼ ex cosðxÞ 2 x2 u þ 2 þ hðx; y; tÞ; ›x ›t ›x ›y
uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:05
uex ðx; y; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ e2y Þcos t
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
0 1 2 3 4 5 6 7 8 9 10 11
Div Div Div 4.5 £ 1023 2.2 £ 1022 4.6 £ 1022 7.3 £ 1022 0.1 0.1 0.2 0.2 0.2
Div Div Div Div 3.3 £ 1022 2.2 £ 1025 5.5 £ 1025 4.4 £ 1024 1.4 £ 1023 3.0 £ 1023 5.4 £ 1023 8.6 £ 1023
Div Div Div Div Div 1.1 £ 1023 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025 1.9 £ 1025 4.5 £ 1025
Div Div Div Div Div Div Div 8.5 £ 1025 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025
N ¼ 5; Dt ¼ 0:05; K ¼ 10 4 5 6 7 8
0.1 0.2 0.3 0.4 0.5
Div 2.8 £ 1022 2.4 £ 1022 2.1 £ 1022 2.0 £ 1022
Div Div 4.4 £ 1022 4.1 £ 1022 3.9 £ 1022
Div Div 9.9 £ 1022 5.0 £ 1022 4.9 £ 1022
0.3 0.4 0.7 1.0 1.0 1.0
1.9 £ 1024 8.9 £ 1024 3.8 £ 1023 9.4 £ 1023 1.9 £ 1022 3.2 £ 1022
Div 3.6 £ 1023 6.7 £ 1024 2.1 £ 1024 2.0 £ 1024 2.1 £ 1024
Div Div 1.7 £ 1022 1.7 £ 1023 2.1 £ 1024 2.1 £ 1024
0.3 0.4 0.7 1.0 1.0 1.0
2.2 £ 1024 9.1 £ 1024 3.6 £ 1023 9.5 £ 1023 1.9 £ 1022 3.2 £ 1022
Div Div Div 1.6 £ 1024 1.8 £ 1024 2.6 £ 1024
Div Div Div Div 0.9 1.5 £ 1024
0.2 0.4 0.6 1.0 1.0 1.0
7.6 £ 1025 8.0 £ 1024 3.6 £ 1023 9.6 £ 1023 1.9 £ 1022 3.2 £ 1022
Div Div Div 3.3 £ 1026 2.6 £ 1025 1.3 £ 1024
Div Div Div Div 1.6 £ 1022 5.1 £ 1026
Example 6. Consider the following PDE: ›u › ›u ¼ ex cosðxÞ 2 x2 u þ hðx; tÞ; ›x ›t ›x with hðx; tÞ; boundary and initial data corresponding to the exact solution (26). Here the coefficients of the operator lðxÞ are: pðxÞ ¼ cosðxÞ;
qðxÞ ¼ x2 e2x ;
rðxÞ ¼ e2x :
To approximate these coefficients N ¼ 10 terms are used in expansion (32). This provides the maximal absolute errors in approximation of lðxÞ :8.1 £ 1026, 5.1 £ 1025 , 1.6 £ 1025 for pðxÞ; qðxÞ and rðxÞ correspondingly. The maximal absolute error in solution is shown in Table 7 and corresponds to the time step Dt ¼ 0:05: Example 7. As an example of 2D-problems consider the PDE:
›u › ›u ›2 u ¼ ex cosðxÞ 2 x2 e2x u þ 2 þ hðx; y; tÞ: ›x ›t ›x ›y Here the solution domain V is the ellipse bounded by the curve x ¼ 0:5 cos u; y ¼ 0:3 sin u; 0 # u # 2p: The function hðx; y; tÞ; initial and boundary conditions, are taken in accordance with the exact solution uex ðx; y; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ e2y Þcos t: The basis system is formed by the products vn ðxÞwð1Þ m ðxÞ: The source points ðjk ; hk Þ are placed on the equidistant ellipse x ¼ 0:95 cos u; y ¼ 0:75 sin u: To evaluate the maximal absolute error 100 checking points are taken inside V on each time-layer tj : In the upper part of Table 8 the results of the calculations with: N ¼ 5 2 number of terms in the approximation (35) of the coefficients of the differential operator; K ¼ 102number of terms in an approximation of the boundary data by the linear combinations like (4) and Dt ¼ 0:05 are presented. One can see that the increase of M
N ¼ 5; Dt ¼ 0:05; K ¼ 30 5 6 7 8 9 10 N ¼ 5; Dt ¼ 0:01; K ¼ 30 5 6 7 8 9 10 N ¼ 10; Dt ¼ 0:01; K ¼ 30 5 6 7 8 9 10
cannot reduce an error to a level below 4 4 5 £ 1022 : When K grows up to 30 the error decreases to a level 2 £ 1024. These data are placed in the second part of the table. An attempt to improve exactness of the calculations by diminishing Dt in five times from 0.05 to 0.01 is shown in the third part of the table. The error has not been reduced. It is due to low approximation of the coefficients of the differential operator. The increase of N from 5 to 10 allows to reduce the error ea approximately in 100 times. These data are presented in the last part of the table.
4. Conclusions A Trefftz type method has been proposed to solve timedependent PDEs. The method suggested follows the general scheme of the MFS with trial functions which satisfy the initial PDE only approximately. Using these trial functions we introduce an additional error in the approximate solution. But as the numerical experiments have shown,
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
this error can be done of the same level as the one due to approximate satisfaction of the boundary condition.
References [1] Reutskiy S. A boundary method of Trefftz type with approximate trial functions. Engng Anal Bound Elem 2002;26(4):341–53. [2] Zhu S, Satravaha P, Lu X. Solving linear diffusion equations with the dual reciprocity method in Laplace space. Engng Anal Bound Elem 1995;13(1):1 –10. [3] Chen CS, Golberg MA, Hon YC. The method of fundamental solutions and quasi-Monte Carlo method for diffusion equations. Int J Numer Meth Engng 1998;43:1421–35. [4] Chen CS, Rashed YF, Golberg MA. A mesh free method for linear diffusion equations. Numer Heat Transfer, Part B 1998;33:469– 86. [5] Muleshkov AS, Golberg MA, Chen CS. Particular solutions of Helmholtz-type operators using higher order polyharmonic splines. Computat Mech 1999;23:411– 9. [6] Golberg MA, Chen CS, Muleshkov AS. The method of fundamental solutions for time-dependent problems. In: Chen CS, Brebbia CA, Pepper DW, editors. Boundary element technology, vol. XIII. Boston: WIT Press; 1999. p. 377– 86. [7] Golberg MA, Chen CS. The method of fundamental solutions for potential, Helmholtz and diffusion problems. In: Golberg MA, editor. Boundary inetgral methods—numerical and mathematical aspects. Southampton, UK: Computational Mechanics Publications; 1998. p. 103–76. [8] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comp Math 1998;9:69–95. [9] Atkinson KE. The numerical evaluation of particular solutions for Poisson’s equation. IMA J Numer Anal 1985;5:319–38. [10] Chen CS. The method of fundamental solutions and the quasi-Monte Carlo method for Poisson’s equation. In: Niederreiter H, Skiue P, editors. Lecture notes in statistics, vol. 106. Berlin: Springer; 1995. p. 158–67. [11] Chen CS, Golberg MA. A domain embedding method and the quasiMonte Carlo method for Poisson’s equation. In: Brebbia CA, Kim S, Osswald TA, Power H, editors. Boundary elements, vol. XVII. Southampton: Computational Mechanics Publications; 1995. p. 158–67. [12] Chen CS, Golberg MA, Hon YC. Numerical justification of the method of fundamental solutions and quasi-Monte Carlo method for Poisson-type equations. Engng Anal Bound Elem 1998;22:37–50.
21
[13] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary elements methods. Southampton, UK: Computational Mechanics Publications; 1992. [14] Chen CS, Brebbia CA, Power H. Dual reciprocity method using compactly supported radial basis functions. Commun Numer Meth Engng 1999;15:137– 50. [15] Chen CS, Golberg MA, Schaback RS. Recent developments of the dual reciprocity method using compactly supported radial basis functions. In: Rashed YF, editor. Transformation of domain effects to the boundary. Boston: WIT Press; 2002. [16] Chen CS, Rashed YF. Evaluation of thin plate spline based particular solutions of Helmholtz-type operators for the DRM. Mech Res Commun 1998;25:195 –201. [17] Herrera I. Trefftz method. In: Brebbia CA, editor. Progress in boundary element methods, vol. 3. New York: Wiley; 1983. [18] Zielinski AP, Zienkiewicz OC. Generalized finite element analysis with T-complete boundary solution functions. Int J Numer Meth Engng 1985;21:509– 28. [19] Cheung YK, Jin WG, Zienkiewicz OC. Direct solution procedure for solution of harmonic problems using complete, non-singular, Trefftz functions. Commun Appl Numer Meth 1989;5:159–69. [20] Plinter R. Recent development in the Trefftz method for finite element and boundary element applications. Adv Engng Software 1995;24: 107–15. [21] Portela A, Charafi A. Programming Trefftz boundary elements. Adv Engng Software 1997;28:509–23. [22] Doicu A, Eremin Y, Wriedt T. Acoustic and electromagnetic scattering analysis using discrete sources. San Diego: Academic Press; 2000. [23] Zielinski AP. On trial functions applied in the generalized Trefftz method. Adv Engng Software 1995;24:147 –55. [24] Reutskiy S, Tirozzi B. Quasi Trefftz spectral method for separable linear elliptic equations. Comp Math Appl 1999;37:47 –74. [25] Reutskiy S, Tirozzi B. A new boundary method for electromagnetic scattering from inhomogeneous bodies. J Quant Spectrosc Radiative Transfer 2002;72:837–52. [26] Lanczos C. Applied analysis. New York: Prentice Hall; 1956. [27] Bogomolny A. Fundamental solutions method for elliptic boundary value problems. SIAM J Numer Anal 1985;22:644 –69. [28] Morse PM, Feshbach H, Methods of theoretical physics, vol. 1 –2. New York: McGraw-Hill; 1953. [29] Carslaw HS, Jaeger JC. Conduction of heat in solids. London: Oxford University Press; 1995. [30] Strang G. Linear algebra and its applications. New York: Academic Press; 1976.
A Trefftz type method for time-dependent problems S.Yu. Reutskiy Laboratory of Magnetohydrodynamics, P.O. Box 136, Moskovski ave. 199, Kharkov 61037, Ukraine Received 29 October 2002; revised 27 May 2003; accepted 30 May 2003
Abstract This paper presents the application of a Trefftz type method with approximate trial functions for time-dependent problems of the parabolic type. At first it was suggested for elliptic equations only. The method was developed for PDEs with variable coefficients when a total differential operator of a problem can be represented as a sum of the one-dimensional ones. A comparison of the numerical solutions with analytic solutions is performed. q 2003 Elsevier Ltd. All rights reserved. Keywords: Trefftz type method; Time differencing; Particular solution; Approximate trial functions; Eigenfunctions
1. Introduction In this paper the numerical method proposed in Ref. [1] is extended onto the initial and boundary value problems (IBVP) for partial differential equations (PDEs) of the type:
›uðx;tÞ ¼ LðxÞ½uðx;tÞþhðx;tÞ; x [ V , Rd ; t $ 0; d ¼ 1;2;3: ›t ð1Þ Here V is a simply connected domain bounded by a simple closed curve ›V: There are two general approaches for solving such equations in the framework of boundary methods: (i) taking the Laplace transform in t and (ii) finite differencing in t. The first case results in the elliptic equation for the Laplace transform u^ ðx; sÞ : ^ sÞ; ðLðxÞ 2 sÞ½^uðx; sÞ ¼ 2uðx; 0Þ 2 hðx;
ð2Þ
which should be solved for a sequence ðsk ; k ¼ 1; …; Ns Þ: The parameters sk (complex in general) depend on the method of inversion used. Finite differencing in time transforms Eq. (1) to a sequence of coupled stationary equations. For example, if the Crank – Nicholson scheme is used, one gets the sequence: ðLðxÞ 2 pÞ½ujþ1 ðxÞ ¼ 2ðLðxÞ þ pÞ½uj ðxÞ 2 2hjþ1=2 ðxÞ:
ð3Þ
E-mail address: [email protected] (S.Yu. Reutskiy). 0955-7997/04/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00115-2
Here uj ðxÞ ¼ uðx;tj Þ; hjþ1=2 ðxÞ ¼ hðx; ðtj þ tjþ1 Þ=2Þ; tj ¼ jDt; Dt is a time step and p ¼ 2=Dt: These techniques were developed in [2 – 7]. In both cases on gets equations with non-zero right hand sides even for a homogeneous initial PDE Eq. (1). So, a problem of finding a particular solution arises when a boundary method is used to solve Eqs. (2) or (3). When L ¼ D the Laplacian, the classical method for obtaining a particular solution is to construct the associated Newtonian potential. Difficulties of a complicated solution domain can be dealt with e.g. by the Atkinson’s embedding technique [9]. The key idea of his method is to embed the solution domain into a larger and simpler Cartesian region V0 instead of the original integration domain V: Combination of this idea with the quasi-Monte Carlo method is performed in Refs. [3,10 –12]. Next approach—the dual reciprocity method (DRM) was developed in order to transform domain integrals to boundary ones [13 – 15]. The key idea is to approximate the right hand side f ðxÞ P of the PDE under consideration LðuÞ ¼ f by a series Nj¼1 aj wj with an appropriate basis functions so that L½Cj ¼ wj can be solved analytically. In the context of the problem considered significant results have been achieved in Refs. [4,5,16]. In the first work analytic formulae were given for Cj when wj was thin plate splines. These formulae proved developing efficient meshfree algorithms for solving diffusion equation in Ref. [4]. In Ref. [5] these results were generalized by using higher order splines as basis functions wj : The corresponding particular
14
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
solutions for D ^ x can be obtained without numerical integration. This is important from the point of view of practical calculations because search for the particular solution should be performed for each sk in Eq. (2) or repeated in each time-layer in Eq. (3). When the initial PDE is reduced to a homogeneous form any boundary technique can be applied. Recently the method of fundamental solutions (MFS) has gained great popularity among researchers. Recall that according to this technique an approximate solution of Eq. (3) is looked for in the form of a linear combination: ujþ1 ðxlq1 ; …; qK Þ ¼ vjþ1 ðxÞ þ
K X
qkjþ1 Fk ðxÞ:
ð4Þ
k¼1
Here vjþ1 ðxÞ is the particular integral in the j þ 1th time layer and the trial functions Fk ðxÞ satisfy exactly the homogeneous PDE ðLðxÞ 2 pÞ½Fk ¼ 0 but do not necessary satisfy boundary conditions. They are used to determine the unknown q kjþ1 : Similar to BEM, the MFS employs fundamental solutions as trial functions but their singular points are situated outside a solution domain. So, the MFS can be treated as a Trefftz-type method [17 – 21]. The techniques of applying the MFS to diffusion-type equations are described in Refs. [3,4,6] and in reviews [7,8]. Note that similar techniques have been developed recently to other branches of applied mathematics under other names: the charge simulation method, the superposition method, the source function method, the modified Trefftz method, etc. As an example of such technique is the method of discrete sources—a very effective technique for scattering problems in homogeneous mediums developed last few years [22]. Other references on this subject can be found in Ref. [8,23]. Through the paper we deal with the differencing in time technique. However, as it will be seen below the same approach also can be used together with the Laplace transform procedure. Similar to Atkinson’s approach, the method presented in this paper includes the embedding procedure. More precisely, the solution domain V is assumed to be embedded into a simple Cartesian region V0 ¼ ½21; þ1d ; d ¼ 1; 2; 3: If this is not the case originally, then appropriate translation and scaling operations may be performed to make it so. The sequence of PDE Eq. (3) is replaced by a new one:
from zero only inside some neighborhood of the source point jk : These functions are analogous in some sense to Dirac’s functions dðx 2 jk Þ which in fact are placed in the right hand sides of given PDEs when the MFS procedure is applied. Free parameters qk ; as usual, are determined from boundary conditions on ›V: The source functions are taken in the form of truncated series IðxljÞ ¼
þ
qjþ1 k Iðxljk Þ;
ð5Þ
k¼1
which differs from the original in an additional sum that is added to the P right hand side of the equation. The sum Kk¼1 qkjþ1 Iðxljk Þ appears due to the boundary nature of the method. It contains free parameters qk and dshaped source functions Iðxljk Þ which essentially differ
def
rn ðM; xÞwn ðjÞwn ðxÞ ¼
n¼1
M X
cn ðM; x; jÞwn ðxÞ;
ð6Þ
n¼1
over some orthogonal basis system wn ðxÞ: Note that n is the multiindex, e.g. in the 2D case n ¼ ðn1 ; n2 Þ; wn ðxÞ ¼ wn1 ðxÞwn2 ðyÞ; cn ðM; x; jÞ ¼ cn1 ðM; x; jÞcn2 ðM; x; hÞ; and it is a double series. It is easy to see that series Eq. (6) coincides with the expansion of Dirac’s functions dðx 2 jÞ over the orthogonal system wn ðxÞ up to the coefficients rn ðM; xÞ ¼ rn1 ðM; xÞrn2 ðM; xÞ: The regularizing coefficients rn ðM; xÞ depend on the particular choice of the system wn ðxÞ [24,25]. Let us dwell in brief on the problem of the sources location. We use the technique in which the coordinates of the sources jk are fixed. Another approach leads to an adaptive version of the MFS when jk are considered as unknown values and should be determined together with the coefficients qjþ1 by the use of a non-linear k minimization process to match the boundary conditions [8]. Bogomolny [27] showed that the fixed sources can be distributed equally on an auxiliary circle in 2D case (sphere in 3D). As indicated by Bogomolny, the larger the radius of the circle (sphere), the better the approximation can be expected. Other results and references on this subject can be found in Ref. [23]. In all the calculations presented in this work the source points jk are placed inside V0 on an auxiliary equidistant contour as far as possible from the boundary ›V: In the 1D case (Examples 1 and 2 in Section 2, Examples 6 in Section 3) the source points j1 ; j2 are placed near the endpoints of the interval V0 ¼ ½21; þ1: The method of this paper is developed for PDEs with differential operators which can be represented as a sum:
ðLðxÞ 2 pÞ½ujþ1 ¼ 2 ðLðxÞ þ pÞ½uj 2 2hjþ1=2 ðxÞ K X
M X
LðxÞ ¼
d X
li ðxi Þ;
i¼1
where each one-dimensional operator li has the form: 1 › › pi ðxi Þ 2 qi ðxi Þ : li ðxi Þ ¼ ›x i ri ðxi Þ ›xi
ð7Þ
The equations like Eq. (1) with such L often arise and are widely used in applied researches and in mathematical
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
modeling. When the solution domain V is a regular region (e.g. the rectangular one) and there exists an exactly known complete system of eigenfunctions for each one-dimensional li ðxi Þ; the problem can be solved using an analytical procedure, e.g. by the method of separation of variables, the methods of integral transforms, etc. See e.g. Ref. [28] for more details. In a general case it can be solved only numerically. We assume that li ðxi Þ is a negative defined in the solution domain, i.e.
ri ðxi Þ; pi ðxi Þ . 0; qi ðxi Þ $ 0;
ð8Þ
when
Examples: (1) For l1 ½w ¼ d2 w=dx2 and the boundary conditions wð21Þ ¼ wðþ1Þ ¼ 0; one gets: xþ1 np 2 ð1Þ wð1Þ ðxÞ ¼ sin np m ¼ ; ; n n 2 2 ð11Þ p ð1Þ g ¼ : 2 (2) The second system w xþ1 x sin np wð2Þ ðx; wÞ ¼ exp ; n 2 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi w2 np 2 p ð2Þ þ ; n ¼ 1; …; 1; gð2Þ ¼ ; mn ¼ 2 2 4
21 , minðxi lx [ VÞ ¼ ai # xi # bi ¼ maxðxi lx [ VÞ , þ1:
corresponds to the 1D convection – diffusion operator
Go back to the choice of the regularizing coefficients rn ðM; xÞ: Let {wn ðxÞ; mn } be a solution of a Sturm –Liouville problem
l2 ðx; wÞ ¼
l½w¼2mw; ða1 ›x þ b1 Þwð21Þ¼0; ða2 ›x þ b2 Þwð1Þ¼0;
ð9Þ
a21;2 þ b21;2 –0; a1 b1 #0; a2 b2 $0; where the operator lðxÞ is given in Eq. (7). Under the assumptions (8) it satisfies the following conditions:
15
›2 › 2w ; 2 ›x ›x
ð12Þ
ð13Þ
ð14Þ
pðxÞ ¼ rðxÞ ¼ expð2wxÞ; qðxÞ ¼ 0; ð2Þ and the boundary conditions wð2Þ n ð21Þ ¼ wn ðþ1Þ ¼ 0: (3) The Bessel’s operator
l3 ðxÞ ¼
›2 1 › ; rðxÞ ¼ pðxÞ ¼ x; qðxÞ ¼ 0; þ 2 x ›x ›x
ð15Þ
together with conditions 0, m1 , m2 ,…, mn !þ1;
wn ð0Þ is bounded; wn ð1Þ ¼ 0;
and there exists the limit:
forms the system of the eigenfunctions: pffiffi 2J0 ðlð3Þ n xÞ ð3Þ ð3Þ 2 wn ðxÞ ¼ ; mð3Þ n ¼ ðln Þ ; g ¼ p; ð3Þ J1 ðln Þ
pffiffiffiffi mn : n!1 n
g ¼ lim
The eigenfunctions wn ðxÞ form an orthogonal system on [2 1, þ 1] with a scalar product: ðþ1 rðxÞwn ðxÞwm ðxÞdx¼gn dn;m : 21
In this case the regularizing coefficients are [24]: 1 sin nðn;MÞ ; ½sn ðMÞx ; sn ðMÞ¼ gn nðn;MÞ pffiffiffiffi p mn : nðn;MÞ¼ gðMþ1Þ rn ðM; xÞ¼
ð10Þ
where mn ; gn ; g correspond to the Sturm –Liouville problem considered and x is the parameter of regularization. For the sake of simplicity, further we omit the arguments M and x in notation of the source functions and their coefficients. A multi-dimensional source function Iðx;jÞ like Eq. (6) can be taken as a product of the one-dimensional ones. We begin with the simplest case when each onedimensional operator (7) possesses a set of eigenfunctions which are known in an explicit analytic form.
ð16Þ
where lð3Þ n is nth zero of the Bessel’s function J0 ðxÞ: These three one-dimensional operators and their combinations will be used in all the numerical examples in the next section to form a total differential operator of the PDE under consideration LðxÞ: The corresponding eigenfunctions wðiÞ n ðxÞ and their products are used as a basis system for an approximate solution. However, in most of practical problems, the corresponding differential operators usually do not have a known spectrum. As it will be shown in Section 3, an approximate solution of the Strum – Liouville problems Eq. (9) can be used in the algorithm. In each section we present results of the numerical experiments justifying the method presented.
2. Differential operators with a known spectrum As it mentioned above the functions wðiÞ n ðxÞ and their products form a natural basis in the problems considered. The trial functions in Eq. (4) satisfy the equation: ðL 2 pÞ½Fk ðxÞ ¼ Iðx; jk Þ; jk [ V0 \ V:
ð17Þ
16
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Due to the special choice of the basis system the solution can be written in the explicit form:
Fk ðxÞ ¼
M X
DðkÞ n wn ðxÞ;
DðkÞ n ¼2
n¼1
cn ðM; x; jk Þ ðjÞ mðiÞ n1 þ mn2 þ p
:
ð18Þ
Suppose a two dimensional case L ¼ li þ lj is considered here. The similar formulae can be written in the 1D and 3D cases. To get a particular solution the function hjþ1=2 ðxÞ in the right hand side of Eq. (5) is approximated by the truncated series over the same basis: hjþ1=2 ðxÞ ¼
M X
Hnjþ1=2 wn ðxÞ:
ð19Þ
the expansion for all the space functions fi ðxÞ: Next,Pto get Hnjþ1=2 at each time-layer as a sum: Hnjþ1=2 ¼ Ni¼1 ai ðtjþ1=2 ÞFi;n : To demonstrate the effectiveness of algorithm we give some numerical examples. Example 1 Consider the following problem:
›u ›2 u ¼ 2 þ hðx; tÞ; ›t ›x uðaÞ ¼ Ua ðtÞ;
uj ðxÞ ¼
M X
Unj wn ðxÞ:
ð20Þ
n¼1
The first expansion of such type can be obtained using the same procedure as the one for hjþ1=2 ðxÞ: Under these conditions: the particular solution vjþ1 ðxÞ (see Eq. (4)) can be written in the same form: vjþ1 ðxÞ ¼
M X
Vnjþ1 wn ðxÞ;
n¼1
Vnjþ1 ¼
ð21Þ j jþ1=2 2 mðjÞ n2 ÞUn þ 2Hn ðjÞ ðiÞ mn1 þ mn2 þ p
mðiÞ n1
ðp 2
;
after determining the free parameters qkjþ1 the solution on the new j þ 1th time layer ujþ1 ðxÞ can be written in the same form as uj ðxÞ : ujþ1 ðxÞ ¼
M X
Unjþ1 wn ðxÞ;
n¼1
Unjþ1
¼
Vnjþ1
þ
K X
ð22Þ qkjþ1 DðkÞ n :
k¼1
So, the code is closed inside a finite set of the basis functions ½wn ðxÞM n¼1 : Some remarks to the practical use of the algorithm. When hðx; tÞ ; 0; i.e. a homogeneous PDE is considered, then the algorithm requires a single Fourier expansion for the initial function uðx; 0Þ to get particular solutions on the all time layers. P When the function hðx; tÞ is a sum of products hðx; tÞ ¼ Ni¼1 ai ðtÞfi ðxÞ; or can be approximated well by such sum, it is more suitable to perform the algorithm in the following way. First, to get the coefficients Fi;n of
uðbÞ ¼ Ub ðtÞ;
uðx; 0Þ ¼ u0 ðxÞ;
t $ 0;
ð23Þ
t $ 0;
ð24Þ
x [ ½a; b , ½21; þ1:
ð25Þ
The source functions take the form:
n¼1
This series can be obtained using the Atkinson’s technique [9] or in another appropriate way. In particular we use the collocation method described in Refs. [1,25]. We also suppose the solution at the jth time layer uj ðxÞ already is represented in the same form of expansion over wn ðxÞ :
x [ ½a; b;
Iðx; jÞ ¼
M X
ð1Þ rn ðM; xÞwð1Þ n ðjÞwn ðxÞ;
n¼1
where the basis functions are given in Eq. (11). In this case the regularizing coefficients are: rn ðM; xÞ ¼ ½sn ðMÞx ; vðn; MÞ ¼
sn ðMÞ ¼
sin½vðn; MÞ ; vðn; MÞ
np : Mþ1
Here sn ðMÞ are the Lanczos sigma-factors which are used to overcome the so-called Gibb’s phenomenon in the Fourier expansion of non-smooth functions [26]. Let the solution domain be the interval ½a; b ¼ ½20:5; þ0:5: The right hand side hðx; tÞ; the boundary functions Ua ðtÞ; Ub ðtÞ and the initial function u0 ðxÞ correspond to the exact solution. uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ:
ð26Þ
Following the rule formulated in Section 1, the two source functions Iðx; jk Þ have parameters j1 ¼ 20:95; j2 ¼ 0:95: The maximal absolute error: ea ¼ max luex ðxi ; tj Þ 2 uðjÞ M ðxi Þl; {xi ;tj }
ð27Þ
as a function of the parameters M and x is presented in Table 1. Here: tj ¼ 0:5; 1:0; …; 10 are the 20 checking time layers and xi ¼ 0:01ði 2 0:5Þ 2 0:5; i ¼ 1; …; 100 are Table 1 The maximal absolute error in solution of the equation: ›u=›u ¼ ›2 u= ›x2 þ hðx; tÞ; uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:05
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
1 2 3 4 5 6 7
Div Div 3.7 £ 1023 2.2 £ 1022 4.7 £ 1022 7.5 £ 1022 0.1
Div Div Div 4.0 £ 1024 1.3 £ 1025 5.6 £ 1025 4.4 £ 1024
Div Div Div 4.7 £ 1023 3.5 £ 1025 1.3 £ 1025 1.3 £ 1025
Div Div Div 1.2 £ 1023 2.0 £ 1025 1.3 £ 1025 1.3 £ 1025
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
the checking points. All the calculations presented in the table were performed with the time step Dt ¼ 0:05: One can see that for each M there exists a minimal xdiv such that the solution process diverges for x , xdiv and it converges for x $ xdiv : The minimal value of ea ðM; xÞ is achieved for xm , xdiv and ea ðM; xÞ grows slightly when x increases. As a function of M the error ea ðM; xm Þ first decreases from ea ð10; 3Þ ¼ 3:7 £ 1023 to ea ð20; 5Þ ¼ 1:3 £ 1025 but then it does not change. This can be explained by the fact that the grows of M can reduce the errors due to approximation of the right hand side hjþ1=2 ðxÞ and the errors due to approximation of the trial functions Fk ðxÞ but is does not affect the error of the time discretization. So, it seems that 1.3 £ 1025 is the error introduced by the Crank – Nicholson scheme which is dominant for these values of M and x: This error should be a square function of the time step Dt: To verify this supposition the calculations were repeated with the time step Dt ¼ 0:01: The results are placed in Table 2. One can see that the behavior of ea ðM; xÞ here is the same as in the previous example but the minimal value of error now is: ea ð30; 8Þ ¼ 5:2 £ 1027 : The ratio 1:3 £ 1025 =5:2 £ 1027 ¼ 25 ¼ ð0:05=0:01Þ2 confirms the nature of the error. So, 5.2 £ 1027 is a new error due to the discretization in time and it cannot be reduced by the increase of the number of harmonics M: Example 2. The similar calculations were performed for the 1D convection –diffusion PDE.
›u ›2 u ›u ¼ 2w þhðx;tÞ; x[½20:5;þ0:5; t $0; ›t ›x 2 ›x
ð28Þ
with the same exact solution (26). The basis functions (12) are used here. The results correspond to w¼5 and the time step Dt ¼0:05 are placed in Table 3. Behavior of the error is the same as above. Here the errors value 1.5 £ 1024 is caused by the Crank – Nicholson time scheme. It decreases to 6.0 £ 1028 when one takes ðM; xÞ¼ð40;11Þ and Dt ¼0:001: Table 2 The maximal absolute error in solution of the equation: ›u=›t ¼ ›2 u= ›x2 þ hðx; tÞ; uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:01
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
3 4 5 6 7 8 9 10 11
3.7 £ 1023 2.2 £ 1022 4.7 £ 1022 7.5 £ 1022 0.1 0.1 0.2 0.2 0.2
Div Div Div 5.5 £ 1025 4.4 £ 1024 1.4 £ 1023 3.1 £ 1023 5.5 £ 1023 8.7 £ 1023
Div Div Div Div Div 5.2 £ 1027 1.5 £ 1026 1.1 £ 1025 4.0 £ 1025
Div Div Div Div Div Div Div 5.2 £ 1027 5.2 £ 1027
17
Table 3 The maximal absolute error in solution of the equation: ›u=›t ¼ ›2 u= ›x2 2 w›u=›x þ hðx; tÞ; uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:05
x
M ¼ 10
M ¼ 20
M ¼ 30
1 2 3 4 5 6 7
Div 3.8 £ 1021 1.5 £ 1023 8.9 £ 1023 2.0 £ 1022 3.5 £ 1022 6.8 £ 1022
Div Div Div 2.4 £ 1024 1.5 £ 1024 1.6 £ 1024 3.0 £ 1024
Div Div Div Div Div 1.5 £ 1024 1.5 £ 1024
Example 3. Consider the heat diffusion in an axially symmetric body with a distributed heat source ›u 1 › ›u ›2 u r 2 þ z2 ¼ : þ 2 þ12 r ›r ›r ›t ›z a2 pffiffiffiffiffiffiffiffiffi Here r ¼ x2 þ y2 ; z are cylindrical coordinates, the solution domain in rz-plane is half of a disk: V ¼ {ðr; zÞ : r 2 þ z2 # a2 ; r $ 0} So, the body is a solid sphere with the radius a: The boundary and initial condition are: uðr; z; tÞ ¼ 0;
r 2 þ z2 ¼ a2 ; t $ 0;
uðr; z; 0Þ ¼ 0;
i.e. the initial temperature and the temperature of the surface is zero. The differential operator of the problem is a sum of the one-dimensional ones l3 ðrÞ and l1 ðzÞ: So the twoð1Þ dimensional basis functions are: wn ðr; zÞ ¼ wð3Þ n1 ðrÞwn2 ðzÞ (see Eqs. (11) and (16)). The problem has an analytic solution [29]. ! ! 1 r2 r2 12 2 uðr; z; tÞ ¼ 723 2 60 a a 1 12a3 X ð21Þnþ1 sin an r 2a2n t e ; ð29Þ rp5 n¼1 n5 pffiffiffiffiffiffiffiffiffi where r ¼ r 2 þ z2 ; an ¼ np=a: In Table 4 this solution is compared with the one obtained by the algorithm described above. All the calculations were performed with Dt ¼ 0:05:
2
Table 4 The maximal absolute error in solution of the 2D equation: ›u 1 › ›u ›2 u r 2 þ z2 ¼ ; þ 2 þ12 r ›r ›r ›t ›z a2 the exact solution is given in Eq. (29), the time step is Dt ¼ 0:05 t
M ¼ 10; x ¼ 2
M ¼ 20; x ¼ 6
M ¼ 30; x ¼ 8
0.1 0.5 1.0 5.0 10.0
3.3 £ 1022 5.0 £ 1022 5.0 £ 1022 5.5 £ 1022 5.5 £ 1022
1.1 £ 1023 8.1 £ 1025 8.3 £ 1025 8.0 £ 1025 8.1 £ 1025
1.1 £ 1023 1.1 £ 1025 6.6 £ 1026 6.1 £ 1026 6.1 £ 1026
18
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Table 5 The maximal absolute error in solution of the 3D PDF:
3. General case
›u ›2 u ›2 u ›2 u x2 þ y2 þ z2 ¼ 2 þ 2 þ 2 þ12 ; Dt ¼ 0:05 ›t ›x ›y ›z a2
For the sake of simplicity consider a one-dimensional equation
t
M ¼ 10; x ¼ 2
M ¼ 20; x ¼ 6
M ¼ 30; x ¼ 8
0.1 0.5 1.0 5.0 10.0
1.1 £ 1023 2.2 £ 1024 2.2 £ 1024 2.3 £ 1024 2.3 £ 1024
9.8 £ 1024 1.7 £ 1025 1.3 £ 1025 1.2 £ 1025 1.2 £ 1025
9.7 £ 1024 9.2 £ 1026 1.8 £ 1026 1.3 £ 1026 1.2 £ 1026
›u ¼ lðxÞ½u þ f ðx; tÞ; ›t
ð30Þ
with lðxÞ given in Eqs. (7) and (8). The two basis systems
cn ðxÞ ¼ cosð0:5ðn 2 1Þpðx þ 1ÞÞ;
ð31Þ
wn ðxÞ ¼ sinð0:5npðx þ 1ÞÞ; Example 4. Consider the same problem as a 3D one in Cartesian coordinates x; y; z :
›u ›2 u ›2 u ›2 u x2 þ y2 þ z2 ¼ 2 þ 2 þ 2 þ12 : ›t ›x ›y ›z a2 The basis system is formed by the products ð1Þ ð1Þ wð1Þ n1 ðxÞwn2 ðyÞwn3 ðzÞ; K ¼ 160 source points are distributed uniformly on the auxiliary sphere of the radius 0.95, Dt ¼ 0:05: The results are placed in Table 5. Example 5. Consider 2D time-dependent convection– diffusion equation:
are used in this section to approximate: (i) the space operator lðxÞ : 2cn ðxÞ and (ii) solutions of the corresponding Sturm – Liouville problems: 2wn ðxÞ: The approximation of lðxÞ pursues the following two goals: (1) to extend the initial PDE from region V to V0 : Note that in one-dimensional case this means a simple extension from the interval ½a; b to the bigger one [2 1, þ 1]. However, for a multi-dimensional PDE the initial irregular solution domain is replaced by a simple Cartesian one. (2) The second goal is to transform the initial PDE into the form convenient for applying the Galerkin projection procedure. The approximations are sought in the form of series: r ð0Þ ðxÞ ¼
›u ›2 u ›2 u ›u ›u ¼ 2 þ 2 2 wx 2 wy þ hðx; y; tÞ: ›t ›x ›y ›x ›y
N X
Ri ci ðxÞ;
i¼1
qð0Þ ðxÞ ¼
ð2Þ The basis system is wð2Þ n1 ðx; wx Þwn2 ðy; wy Þ: The solution domain in xy-plane is a disk: V ¼ {ðx; yÞ : x2 þ y2 # a2 }: The boundary condition of the Dirichlet type, the initial condition and the right hand side hðx; y; tÞ correspond to the exact solution: uex ðx; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ expð2yÞÞ cos t: Table 6 displays the maximal absolute error for the two cases: (1) wx ¼ wy ¼ 5; (2) wx ¼ 15; wy ¼ 215: The pairs ðM; xÞ are placed in parentheses. The time step Dt ¼ 0:01 and K ¼ 40 source points are used in both cases.
N X
pð0Þ ðxÞ ¼
N X
Pi ci ðxÞ;
i¼1
Qi ci ðxÞ:
i¼1
To get the unknown coefficients Ri ; Pi ; Qi we write the collocation conditions P for some inner points xj [ ½a; b; j ¼ 1; …; N $ N; e.g. Ni¼1 Ri ci ðxj Þ ¼ rðxj Þ and then solve the system by an appropriate numerical procedure. Note that the Atkinson’s technique mentioned above also can be used there. As a result one gets a new differential operator 1 › ð0Þ › ð0Þ l ðxÞ ¼ ð0Þ p ðxÞ 2 q ðxÞ ; ›x r ðxÞ ›x ð0Þ
Table 6 The maximal absolute error in solution of the 2D convection-diffusion equation: 2
2
›u › u › u ›u ›u ¼ 2 þ 2 2 wx 2 wy þ hðx; y; tÞ; ›t ›x ›y ›x ›y exact solution Dt ¼ 0:001 t
0.5 1.0 1.5 2.0
uex ðx; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ expð2yÞÞcos t;
wx ¼ 5; wy ¼ 5
wx ¼ 15; wy ¼ 215
(20,5)
(30,8)
(10,3)
(20,5)
(30,8)
3 £ 1022 2 £ 1022 3 £ 1023 1 £ 1022
3 £ 1023 2 £ 1023 3 £ 1024 1 £ 1023
5 £ 1026 3 £ 1026 6 £ 1027 2 £ 1026
1 £ 10þ3 7 £ 10þ2 1 £ 10þ2 5 £ 10þ2
3 £ 1023 2 £ 1023 3 £ 1024 1 £ 1023
3 £ 1024 1 £ 1024 2 £ 1025 1 £ 1024
ð33Þ
defined in [2 1, þ 1] which approximates lðxÞ when x [ ½a; b: We assume that lð0Þ ðxÞ continues to be a negative defined operator, i.e. the conditions (8) are true for r ð0Þ ðxÞ; pð0Þ ðxÞ; qð0Þ ðxÞ: Let us consider the Sturm – Liouville problem on [2 1, þ 1]: lð0Þ ðxÞ½v ¼ 2mv;
(10,3)
ð32Þ
vð21Þ ¼ vðþ1Þ ¼ 0:
ð34Þ
An approximate solution of the problem can be looked for in the form of Mth term expansion:
vðxÞ ¼
M X n¼1
W n wn ðxÞ:
ð35Þ
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Projecting with the weight r ð0Þ ðxÞ on the finite basis w1 ðxÞ; …; wM ðxÞ one gets the sequence of equations: klð0Þ ðxÞ½v þ mv; wn ðxÞl ¼ 0;
ð36Þ
where the scalar product is: ðþ1 r ð0Þ ðxÞuðxÞvðxÞdx: ku; vl ¼
ð37Þ
19
when lm 2 nl . N; i.e., when the number of harmonics M exceeds the number of terms in the expansion of the coefficients, then the matrices have a band-like form. It is well-known that under the assumptions mentioned m }; above problem (38) has exactly M solutions {mm ; w m ¼ 1; …; M such that the eigenvectors w m and corresponding functions
21 M X
Substituting Eq. (35) in Eq. (36) one gets the eigenvalue problem:
vm ðxÞ ¼
^w ^ w; S· ¼ mR·
satisfy the conditions [30]:
ð38Þ
kw k; w m lR ¼
w ¼ ½W ; …; W ; M T
ðþ1 21
r ð0Þ ðxÞwn ðxÞwm ðxÞdx:
ð40Þ
Substitution Eq. (32) in Eqs. (39) and (40) yields: Sn;m ¼
N X
½0:25nmp2 Pi Z1 ði 2 1; n; mÞ þ Qi Z2 ði 2 1; n; mÞ;
i¼1
N X
Ri;j Wki Wmj ¼ dkm ;
kvk ; vm l ¼ dkm :
ð46Þ
The first condition means orthogonality of the vectors w k with the weight matrix Ri;j and the second one means orthogonality of the functions vk ðxÞ in the sense of the scalar product (37). In particular orthogonality of the system vm ðxÞ; m ¼ 1; …;M means that a series f ðxÞ ¼ P M wn ðxÞ can be n¼1 Fn wn ðxÞ over the basis functions P v further re-expanded over vk ðxÞ : f ðxÞ ¼ M m¼1 Fm vm ðxÞ: The coefficients of both expansion Fn and Fmv are connected by the relations: Fmv ¼
M X
Tm;n Fn ;
n¼1
ð41Þ Rn;m ¼
M X i;j¼1
^ are positive defined matrices with the entries: S^ and R ðþ1 dw ðxÞ dwm ðxÞ Sn;m ¼Sm;n ¼ pð0Þ ðxÞ n dx dx 21 þ qð0Þ ðxÞwn ðxÞwm ðxÞ dx; ð39Þ Rn;m ¼ Rm;n ¼
ð45Þ
n¼1
where 1
Wmn wn ðxÞ;
Ri Z2 ði 2 1; n; mÞ;
ð42Þ
Tm;n ¼
M X
Fn ¼
M X
Wmn Fmv ;
m¼1
ð47Þ
Wmk Rn;k :
k¼1
i¼1
The following property of the v-basis plays a central role in the algorithm [24]:
where Z1 ði; n; mÞ ¼ ¼
ðþ1 21
ci ðxÞcn ðxÞcm ðxÞdx
klð0Þ ðxÞ½vm ðxÞ; vn ðxÞl ¼ 2mn dmn :
1 ½Dði þ n 2 mÞ þ Dði þ m 2 nÞ þ Dðn 2
þ m 2 iÞ þ Dði þ n þ mÞ; ðþ1 Z2 ði; n; mÞ ¼ ci ðxÞwn ðxÞwm ðxÞdx
ð43Þ
(
DðkÞ ¼
where Iðx; jÞ is given in Eq. (6). Re-expanding it over the v-basis and projecting to each vn ðxÞ one gets
1 ½Dði þ n 2 mÞ þ Dði þ m 2 nÞ 2 Dðn 2 þ m 2 iÞ 2 Dði þ n þ mÞ;
1
if k ¼ 0
0
otherwise
In particular it is possible to write an approximation of the trial function Fk ðxÞ using the v-basis. Let us consider the PDE for Fk ðxÞ : ðlð0Þ ðxÞ 2 pÞ½Fk ðxÞ ¼ Iðx; jk Þ;
21
¼
ð48Þ
ð44Þ
:
^ are It is important to note that entries of the matrices S^ and R calculated in an analytic form without any numerical integration due to the special approximation (32). In general case M 2 numerical integrations are needed to get the matrices ^ Moreover, it easy to prove that Sn;m ¼ 0; Rn;m ¼ 0 S^ and R:
Fk ðxÞ ¼
M X n¼1
dnv ðjk ;M; xÞvn ðxÞ;
dnv ðjk ;M; xÞ ¼ 2
cvn ðj;M; xÞ ; mn þ p
where cvn ðj;M; xÞ are obtained from cn ðj;M; xÞ using linear transformation (47). In the similar way in the 2D case one gets an analog of Eq. (18) where cn ðM; x; jk Þ is replaced by cvn ðM; x; jÞ ¼ cvn1 ðM; x; jÞcvn2 ðM; x; hÞ: In general, all the formulae derived in the previous section for differential operators with a known spectrum are valid when one takes the expansion over the v-basis.
20
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
Table 7 The maximal absolute error in solution of the equation: ›u › ›u ¼ ex cosðxÞ 2 x2 u þ hðx; tÞ; ›x ›t ›x
Table 8 The maximal absolute error in solution of the equation: ›u › ›u ›2 u ¼ ex cosðxÞ 2 x2 u þ 2 þ hðx; y; tÞ; ›x ›t ›x ›y
uex ðx; tÞ ¼ ½1 þ sinðxÞcosðtÞ; Dt ¼ 0:05
uex ðx; y; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ e2y Þcos t
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
x
M ¼ 10
M ¼ 20
M ¼ 30
M ¼ 40
0 1 2 3 4 5 6 7 8 9 10 11
Div Div Div 4.5 £ 1023 2.2 £ 1022 4.6 £ 1022 7.3 £ 1022 0.1 0.1 0.2 0.2 0.2
Div Div Div Div 3.3 £ 1022 2.2 £ 1025 5.5 £ 1025 4.4 £ 1024 1.4 £ 1023 3.0 £ 1023 5.4 £ 1023 8.6 £ 1023
Div Div Div Div Div 1.1 £ 1023 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025 1.9 £ 1025 4.5 £ 1025
Div Div Div Div Div Div Div 8.5 £ 1025 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025 1.4 £ 1025
N ¼ 5; Dt ¼ 0:05; K ¼ 10 4 5 6 7 8
0.1 0.2 0.3 0.4 0.5
Div 2.8 £ 1022 2.4 £ 1022 2.1 £ 1022 2.0 £ 1022
Div Div 4.4 £ 1022 4.1 £ 1022 3.9 £ 1022
Div Div 9.9 £ 1022 5.0 £ 1022 4.9 £ 1022
0.3 0.4 0.7 1.0 1.0 1.0
1.9 £ 1024 8.9 £ 1024 3.8 £ 1023 9.4 £ 1023 1.9 £ 1022 3.2 £ 1022
Div 3.6 £ 1023 6.7 £ 1024 2.1 £ 1024 2.0 £ 1024 2.1 £ 1024
Div Div 1.7 £ 1022 1.7 £ 1023 2.1 £ 1024 2.1 £ 1024
0.3 0.4 0.7 1.0 1.0 1.0
2.2 £ 1024 9.1 £ 1024 3.6 £ 1023 9.5 £ 1023 1.9 £ 1022 3.2 £ 1022
Div Div Div 1.6 £ 1024 1.8 £ 1024 2.6 £ 1024
Div Div Div Div 0.9 1.5 £ 1024
0.2 0.4 0.6 1.0 1.0 1.0
7.6 £ 1025 8.0 £ 1024 3.6 £ 1023 9.6 £ 1023 1.9 £ 1022 3.2 £ 1022
Div Div Div 3.3 £ 1026 2.6 £ 1025 1.3 £ 1024
Div Div Div Div 1.6 £ 1022 5.1 £ 1026
Example 6. Consider the following PDE: ›u › ›u ¼ ex cosðxÞ 2 x2 u þ hðx; tÞ; ›x ›t ›x with hðx; tÞ; boundary and initial data corresponding to the exact solution (26). Here the coefficients of the operator lðxÞ are: pðxÞ ¼ cosðxÞ;
qðxÞ ¼ x2 e2x ;
rðxÞ ¼ e2x :
To approximate these coefficients N ¼ 10 terms are used in expansion (32). This provides the maximal absolute errors in approximation of lðxÞ :8.1 £ 1026, 5.1 £ 1025 , 1.6 £ 1025 for pðxÞ; qðxÞ and rðxÞ correspondingly. The maximal absolute error in solution is shown in Table 7 and corresponds to the time step Dt ¼ 0:05: Example 7. As an example of 2D-problems consider the PDE:
›u › ›u ›2 u ¼ ex cosðxÞ 2 x2 e2x u þ 2 þ hðx; y; tÞ: ›x ›t ›x ›y Here the solution domain V is the ellipse bounded by the curve x ¼ 0:5 cos u; y ¼ 0:3 sin u; 0 # u # 2p: The function hðx; y; tÞ; initial and boundary conditions, are taken in accordance with the exact solution uex ðx; y; tÞ ¼ ð1 þ x2 þ 5x3 Þð1 þ y þ e2y Þcos t: The basis system is formed by the products vn ðxÞwð1Þ m ðxÞ: The source points ðjk ; hk Þ are placed on the equidistant ellipse x ¼ 0:95 cos u; y ¼ 0:75 sin u: To evaluate the maximal absolute error 100 checking points are taken inside V on each time-layer tj : In the upper part of Table 8 the results of the calculations with: N ¼ 5 2 number of terms in the approximation (35) of the coefficients of the differential operator; K ¼ 102number of terms in an approximation of the boundary data by the linear combinations like (4) and Dt ¼ 0:05 are presented. One can see that the increase of M
N ¼ 5; Dt ¼ 0:05; K ¼ 30 5 6 7 8 9 10 N ¼ 5; Dt ¼ 0:01; K ¼ 30 5 6 7 8 9 10 N ¼ 10; Dt ¼ 0:01; K ¼ 30 5 6 7 8 9 10
cannot reduce an error to a level below 4 4 5 £ 1022 : When K grows up to 30 the error decreases to a level 2 £ 1024. These data are placed in the second part of the table. An attempt to improve exactness of the calculations by diminishing Dt in five times from 0.05 to 0.01 is shown in the third part of the table. The error has not been reduced. It is due to low approximation of the coefficients of the differential operator. The increase of N from 5 to 10 allows to reduce the error ea approximately in 100 times. These data are presented in the last part of the table.
4. Conclusions A Trefftz type method has been proposed to solve timedependent PDEs. The method suggested follows the general scheme of the MFS with trial functions which satisfy the initial PDE only approximately. Using these trial functions we introduce an additional error in the approximate solution. But as the numerical experiments have shown,
S.Yu. Reutskiy / Engineering Analysis with Boundary Elements 28 (2004) 13–21
this error can be done of the same level as the one due to approximate satisfaction of the boundary condition.
References [1] Reutskiy S. A boundary method of Trefftz type with approximate trial functions. Engng Anal Bound Elem 2002;26(4):341–53. [2] Zhu S, Satravaha P, Lu X. Solving linear diffusion equations with the dual reciprocity method in Laplace space. Engng Anal Bound Elem 1995;13(1):1 –10. [3] Chen CS, Golberg MA, Hon YC. The method of fundamental solutions and quasi-Monte Carlo method for diffusion equations. Int J Numer Meth Engng 1998;43:1421–35. [4] Chen CS, Rashed YF, Golberg MA. A mesh free method for linear diffusion equations. Numer Heat Transfer, Part B 1998;33:469– 86. [5] Muleshkov AS, Golberg MA, Chen CS. Particular solutions of Helmholtz-type operators using higher order polyharmonic splines. Computat Mech 1999;23:411– 9. [6] Golberg MA, Chen CS, Muleshkov AS. The method of fundamental solutions for time-dependent problems. In: Chen CS, Brebbia CA, Pepper DW, editors. Boundary element technology, vol. XIII. Boston: WIT Press; 1999. p. 377– 86. [7] Golberg MA, Chen CS. The method of fundamental solutions for potential, Helmholtz and diffusion problems. In: Golberg MA, editor. Boundary inetgral methods—numerical and mathematical aspects. Southampton, UK: Computational Mechanics Publications; 1998. p. 103–76. [8] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comp Math 1998;9:69–95. [9] Atkinson KE. The numerical evaluation of particular solutions for Poisson’s equation. IMA J Numer Anal 1985;5:319–38. [10] Chen CS. The method of fundamental solutions and the quasi-Monte Carlo method for Poisson’s equation. In: Niederreiter H, Skiue P, editors. Lecture notes in statistics, vol. 106. Berlin: Springer; 1995. p. 158–67. [11] Chen CS, Golberg MA. A domain embedding method and the quasiMonte Carlo method for Poisson’s equation. In: Brebbia CA, Kim S, Osswald TA, Power H, editors. Boundary elements, vol. XVII. Southampton: Computational Mechanics Publications; 1995. p. 158–67. [12] Chen CS, Golberg MA, Hon YC. Numerical justification of the method of fundamental solutions and quasi-Monte Carlo method for Poisson-type equations. Engng Anal Bound Elem 1998;22:37–50.
21
[13] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary elements methods. Southampton, UK: Computational Mechanics Publications; 1992. [14] Chen CS, Brebbia CA, Power H. Dual reciprocity method using compactly supported radial basis functions. Commun Numer Meth Engng 1999;15:137– 50. [15] Chen CS, Golberg MA, Schaback RS. Recent developments of the dual reciprocity method using compactly supported radial basis functions. In: Rashed YF, editor. Transformation of domain effects to the boundary. Boston: WIT Press; 2002. [16] Chen CS, Rashed YF. Evaluation of thin plate spline based particular solutions of Helmholtz-type operators for the DRM. Mech Res Commun 1998;25:195 –201. [17] Herrera I. Trefftz method. In: Brebbia CA, editor. Progress in boundary element methods, vol. 3. New York: Wiley; 1983. [18] Zielinski AP, Zienkiewicz OC. Generalized finite element analysis with T-complete boundary solution functions. Int J Numer Meth Engng 1985;21:509– 28. [19] Cheung YK, Jin WG, Zienkiewicz OC. Direct solution procedure for solution of harmonic problems using complete, non-singular, Trefftz functions. Commun Appl Numer Meth 1989;5:159–69. [20] Plinter R. Recent development in the Trefftz method for finite element and boundary element applications. Adv Engng Software 1995;24: 107–15. [21] Portela A, Charafi A. Programming Trefftz boundary elements. Adv Engng Software 1997;28:509–23. [22] Doicu A, Eremin Y, Wriedt T. Acoustic and electromagnetic scattering analysis using discrete sources. San Diego: Academic Press; 2000. [23] Zielinski AP. On trial functions applied in the generalized Trefftz method. Adv Engng Software 1995;24:147 –55. [24] Reutskiy S, Tirozzi B. Quasi Trefftz spectral method for separable linear elliptic equations. Comp Math Appl 1999;37:47 –74. [25] Reutskiy S, Tirozzi B. A new boundary method for electromagnetic scattering from inhomogeneous bodies. J Quant Spectrosc Radiative Transfer 2002;72:837–52. [26] Lanczos C. Applied analysis. New York: Prentice Hall; 1956. [27] Bogomolny A. Fundamental solutions method for elliptic boundary value problems. SIAM J Numer Anal 1985;22:644 –69. [28] Morse PM, Feshbach H, Methods of theoretical physics, vol. 1 –2. New York: McGraw-Hill; 1953. [29] Carslaw HS, Jaeger JC. Conduction of heat in solids. London: Oxford University Press; 1995. [30] Strang G. Linear algebra and its applications. New York: Academic Press; 1976.