- Email: [email protected]
A least-squares method for the Helmholtz equation
Computer methods in applied mechanics and engineering Comput. Methods Appl. Mech. Engrg. 175 (1999) 12 I - 136
EJ ,SEVIER
www.elsevier.com/locate/cma
A least-squares method for the Helmholtz equation P. M o n k * , D a - Q i n g W a n g Department qf Mathematical Sciences, University of Delaware, Newark, DE 19716, USA Received 9 March 1998
Abstract
We investigate the use of least-squares methods to approximate the Helmholtz equation. The basis used in the discrete method consists of st lutions of the Helmholtz equation (either consisting of plane waves or Bessel functions) on each element of a finite element grid. Unlike p~evious methods of this type, we do not use polynomial based finite elements. The use of small elements (and relatively few basis functions per element) allows us to prove convergence theorems for the method and, to some extent, control the conditioning of the resulting linear s} ~tem. Numerical results show the efficiency of the new method and suggest that it may be possible to obtain accurate results with a coarser grid than is usual for standard finite element methods. © 1999 Elsevier Science S.A. All rights reserved.
1. Introduction
When computing the solution of the Helmholtz equation or Maxwell's equations using finite element or finite difference methods we are faced with the problem that there must be a sufficient number of discretization points per wavelength to resolve the solution. A standard rule of thumb is to use at least ten grid points per wavelength, but even this is not sufficient for large problems in which the computational domain is many wavelengths in diameter [13]. This has prompted attempts to use other methods such as corrected least-squares methods, special basis functions or higher-order methods (see e.g. [14,1,12]). In another direction of particular interest to us, Cessenat and Despres [5] suggest using exact solutions of the partial differential equation on each element and choose the expansion coefficients to obtain approximate continuity of the function and its flux across element boundaries. This is done using an approach similar to a mixed finite element method call the 'ultra-weak' wtriational method. The authors report very exciting results on the efficiency and convergence rate of the method. Nevertheless, the method is rather complicated. In this paper we also want to investigate the use of local complete families of solutions of the Helmholtz equation to approximate a scattering problem. But instead of using the 'ultra-weak' variational method of C,zssenat and Despres [5] we will enforce continuity across element boundaries by minimizing a simple quadratic functional. This approach leads to solving a linear system with a Hermitian matrix to obtain the discrete solution, and allows us to use the conjugate gradient scheme to solve the discrete problem. To fix ideas, let /2 C ~2 be a bounded smooth domain with boundary F consisting of two disconnected components F~ and F,. We wish to approximate the solution u of the model scattering problem
Au+k2u=O
in ~ ,
(1.1)
u =g c3u
on F,,
(1.2)
on I_~
(1.3)
;Ju
iku = 0
~' Corresponding author. 00~5-7825/99/$ - see front matter PII: S 0 0 4 5 - 7 8 2 5 ( 9 8 ) 0 0 3 2 6 - 0
© 1999 Elsevier Science S.A. All rights reserved.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 17.5 (1999) 121-136
122
Here, g is a given smooth function and k is a fixed, real positive parameter called the wave-number. The boundary condition (1.3) is a simple absorbing boundary condition. Although not the best absorbing boundary condition available, this condition is simple to implement and we can derive exact solutions on some domains for comparison with numerical results. We shall test two methods for approximating (1.1)-(1.3). In both methods we first mesh the region ,(2 using triangles (we could use more general meshes since there is no need to use a standard finite element grid). The two methods differ according to the expansion functions used on each element. In one case we use a basis of plane waves as in [5] and in the other case we use a Bessel function basis. Using a duality estimate we can estimate the error in the solution from a knowledge of the approximation properties of the appropriate basis. We shall prove the following: (1) For the plane wave basis, convergence under mesh refinement follows from a theorem of [5]. Convergence for a fixed mesh but increasing the number of basis functions is proved using a theorem from 1181. (2) For Bessel functions we prove convergence both under mesh refinement and use a theorem from [18] to prove convergence as the number of basis functions increases on a fixed mesh. Using a fixed grid, our Bessel function method is similar in spirit to the scheme for computing interior eigenvalues presented by Driscoll [9]. If only one element is used the method becomes the Trefftz method analyzed by Eisenstat [10] (see also l11]). Similar methods include the weighted residual method [6] and the boundary approximation method [15-17]. As in [9], the flexibility offered by using a grid can be used to improve the convergence rate of the method. Our estimates, and experience, shows that both the methods we propose converge at high order to smooth solutions of (1.1)-(1.3) if sufficient basis functions are used on each element. A similar approach to ours has been presented in [20]. The least-squares functional is different to ours (in essence we also fit tangential derivatives along element boundaries) and no proof of convergence is given. However, [20] contains several useful extensions of least-squares methods including a discussion of problems with corners (we restrict ourselves to smooth domains). A related method is the partition of unity finite element method of Melenk and Babu~ka [18,2]. We shall use theoretical results from these papers in our work. The method differs from ours in that the partition of unity is used to create a conforming subspace for the standard weak form of (1.1)-(1.3). Hence, a standard Galerkin scheme can be used and the method fits in to the theoretical framework of standard finite element methods. However, the Galerkin formulation does not produce a Hermitian positive definite linear system. Finite element least-squares methods for this problem (see e.g. [14]) and for general elliptic and first-order systems are very well developed (see e.g. [3,4]). The difference here is that we do not use finite elements. We end the paper with a series of numerical which show that the scheme can be efficient. We also comment on the conditioning of the linear system for the problem. We complete the introduction by remarking about the notation for function spaces and norms in the paper. Let S be a Lipschitz domain in the plane with boundary OS. We denote by HP(S) the usual Sobolev space of functions with p derivatives in L2(S). The norm on this space is Similarly, there is also a space of functions H~(OS) on the boundary, and we denote the norm I " I~,.~,s. Sometimes, we will need to use semi-norms. We denote by I/" Ilp,~.s the pth order semi-norm (i.e. the square root of the sum of the squares of the norms of the pth order derivatives).
II-II,, .
2. The method Suppose we cover 12 with a finite element mesh % of regular and quasi-uniform elements of maximum diameter h. In fact we do not need to be so restrictive, since the grid could contain a mixture of element shapes and there is no need for the usual finite element restrictions on vertex placement (see [7]). Here, we shall present results for triangular meshes. The edges in the mesh are divided into three disjoint sets
e
[ Eo 12' ~
if e is in the interior of ~,
boundary
~
if e is on the boundary ~ .
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
123
Each edge e has an associated normal t,. For a sufficiently smooth function f defined in the neighborhood of e we define jump of f across e, for given x E e, by [f](x) =lim (f(x + eu ) - f(x - ev~)) e ~l)
On each element K E 7,, we assume that we have N smooth basis functions {tpeK}p,v= which are solutions of the Helmholtz equation on K, and vanish outside K (so they are discontinuous). Then, we approximate u by u, defined by N
E E a,fC,,f(x), KGvt~
p -
I
where the coefficients aeK need to be determined. We do this by a simple least-squares technique. Let a be the vector of length M = N card(~,) listing all the unknown coefficients. Then, we define J(a)= ~
"~,
9
2
{l[ u,,l[o,,+A~,l[u,,][o,~} + ~
e E=-£11
"~
eEE
2
A~.[(u,,-g)lo,,,+ ~
eEE_~
~Ua
-~-iku,,
I2 O,e
,
(2.1)
where A is a positive coefficient on each edge which might depend on N, k, etc. In this paper we shall choose A,=k. We have tested other choices computationally. For example, we tried ,~ = 1/h~ where h, is the length of edge e but did not find a choice that was a marked improvement over h~ = k. The correct weighting of terms (perhaps involving different norms) is likely to be important to improve conditioning and accuracy of the method. This aspect needs to be investigated further. The solution of the least squares problem a = arg rain J(b)
(2.2)
gives the coefficients of u. and hence the numerical solution of our problem. As discussed in the Introduction, candidates for the basis functions are as follows:
Plane waves Here, we choose N directions ~,, 1 ~
~bj, (x) = exp(ikdp • (x - xK)),
(2.3)
where x s. is the centroid of K. In this paper we shall choose d = (cos(0t,), sin(0p)) r ,
(2.4)
where Op = 2"rr(p - 1)/(N - 1). It is possible that a non-uniform choice of the angle parameter Or (perhaps even limited to a suitable sector) might help satisfy the outgoing radiation condition but we will not investigate this aspect here.
~essel functiorzs Here, we choose N = 2P + 1 then K
~p(x)=,l~p.e_l~(kr)exp(i(p-P)O),
l<~p<~N,
where J, is the cylindrical Bessel function of first kind and order p - P - 1 , cylindrical coordinates with their origin at the centroid x K.
(2.5) and r = I x - x K [ a n d 0 are
:,. E r r o r e s t i m a t e s
We start by proving a simple duality result that provides the basic error estimate for the method. Next, we apply the duality result together with approximation properties of the basis functions to prove the desired error estimates.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
124
THEOREM 3.1. Let u be a smooth solution of (1.1)-(1.3) and let u a be a discrete function computed using either the plane wave or Bessel function basis with coefficients given by a then
Ch "20(a))"2
II. - uoll,,.~ ~
PROOF. Let v satisfy the dual problem At) + kZu = 05
in ~(2,
u=0
on C ,
Ov - - +ikv = 0
on C ,
0v
where 05 is an arbitrary function in Lz(O). Using integration by parts and the fact that on each element u, is a solution of the Helmholtz equation, we can write (U -- U,,, 05) = (U -- a ,,, AV + k2v)
= ~
(u - u,,, Av + k2v)~
K G "i"h
(
io
= ~ (a(u-.,,)+k~(u-uo),v)+ u-uo, ~-~ ~K- --~(U-Uo),V K E,rh
)
OK
Ov K E'r h
Here, (., -)s is the L2(S) inner product (or duality pairing as appropriate) and vx is the unit outward normal to K. Now, rewriting these sums as sums over edges we obtain
e~b"
The true solution, and its gradient, are continuous across interior edges. In addition v = 0 and u - - g on ~. Furthermore, we can rearrange the sum on F= using the absorbing boundary condition to obtain
eEL'()
+
~ (
e
e
~',Ov)
-
E
(( O~_ik)(u_u,),v) i
•
5
e
e~E~
e
Hence, using the Cauchy-Schwarz inequality we obtain
Ov
eEEv
By the Cauchy-Schwarz inequality again
O,e
eEE~
,
P. Monk, D.-Q. Wang / Comput. Methods" AppL Mech. Engrg. 175 (1999) 121-136
)(u - uo, ?))1 ~
ir + 11 I( ) I
~l[u,,[lo~ +
u,,
e - u,I;,,, ~ + E A,,Ig
+ E
125
e~l'-,
~O p -ik
u,~
e~E_~
-Z
O.e
e
I°°1
2
+ J~l 2
O,e
)]
1/2
0,e
Hence, we obtain the bound (3.1)
[(u - u., 6)1 <~(J(a))~'21110111 where
We estimate Jl[vl[I using the trace estimate that for any suitably smooth function w lw]o.~ ~< C x (h s. Vw ~ ~¢ + h~'l}w o x ) ,
e
'
~
t
(3.2)
'
where K~ is a triangle with e as an edge, hx," is the diameter of K, and CK~ is given by
CK'
hK. Px
where PK,, is the radius of the largest inscribed circle in K e (see e.g. [19]). Then, we have that
IIio111 ~ where
c~,
1 (hK, Ilotl2,,K 2 +
h-'l[Vvll~,) + [hK W ~ + h~,,1 Ivlt'~) ~
ll'lJz.~,~,, is the 2-seminorm on K,. Thus, if a,. = k
-112
We obtain the desired estimate by using this estimate in (3.1).
D
3./. Plane waves 7n this section we analyze the use of the plane wave basis given in (2.3) with uniformly distributed propagation angles given by (2.4).
3.1.1. Convergence under mesh refinement ttere, we assume that a fixed N has been chosen, and we study convergence under classical finite element grid refinement (i.e. h--) 0). The error estimates of Cessenat and Despres [5} can be used to prove the following convergence estimate. THEOREM 3.2. Suppose u ~ C ~u ~J2-~(~) then there exists a constant C independent o f h such that
II. - .,Ak,~ ~< c h'~' ,,,2,
(3.3)
REMARK. Suppose N = 2n + I so N is odd. If N = 5 we are guaranteed O(h) convergence, if N = 7 we get O(h:~), and if N = 1 1 we get O(h4). PROOF. In the proof of Theorem 3.7 of Cessenat and Despres [5] it is shown that on there is a function u,, such that if u satisfies (1.1) then for each x E / 2 ,
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
126
lu(x) - u,,(x)[ <~ Ch "+'llullC2,,- ,,,, IV(u(x) - u.,(x))] ~<
Ch"llu[Ic2,, ~,,,,,
Hence, using this estimate, on the edges of the mesh o
Ch~,,+,
eE~ h
ull 2
.
e~Ttl
Using the fact that in a quasi-uniform mesh there are at m o s t C / h 2 edges we obtain "~
E
e ~ rzj
lu -- u ,}~,.~+ IV(u - u .)1o.,, ~< Ch
2n
1
2
}[Ullc2,,+,,~,).
Using this estimate, we can easily see that if A = k and N = 2n + 1 then (min J(a)) 1/2 <~
Ch"-"211ullc2.....
[]
3.2. Convergence tinder increasing order In this case we assume that a fixed grid is chosen, and that convergence is to be obtained by increasing the number of directions per cell N. We can use a Theorem of Babu~ka and Melenk [2] to prove convergence. In particular (using our notation), they state the following result: THEOREM 3.3 (Babu.(ka and Melenk, Theorem 5). Let K be an element in the mesh and let the exterior angles of K be bounded below by ATr, 0 < , t < 2 . Assuming uEH'~(K), s > l then there is a function u x = N K K Y>-~ I al, tgl, (x) (see (2.3)) such that [ In 2 N'~ ~'~ )> [lu - u K,Ili.~ <~ ~ - / Ilull,.K
j = 0 . . . . . s - 1.
Using this result and our duality estimate we can easily estimate the error in the least-squares method: THEOREM 3.4. Assume that u E H"(f2) f o r some s ~ 2. Let u,, be the approximate solution computed using plane wave functions on a .fixed mesh. We have the following error estimate: f ln2 N'~ a,.io,.,. -2~
Ill,
-
u, IJ,,.,, ~< c
~ - )
llull, ,,
where Ami.'rr is a lower bound for the smallest exterior angle of any triangle in the mesh (see Theorem 3.3). PROOF. We use Theorem 3.1. Let a ' denote the coefficients of the plane wave function expansion guaranteed by Theorem 3.4 on each triangle. Then
}lu - u I},,,. <~ Ch '/2(j(a'))'/2 To estimate
J(a') note that it suffices to estimate
I[Vu,,,]h2,.~ + Ih,,[u,,,]], 2,....
(3.4)
Suppose e is common to elements K~ and K 2, and denote by u,,,~ and u , 2 the Bessel function approximations on K~ and K 2, respectively
I~v,,,,llo,,, ~< Iv(u - u,, )l,,.,, + IV(u - u., :)L. But by standard trace results for Sobolev spaces (recalling that the grid is fixed) the two terms on the right-hand side can be estimated as follows (we only give details of the first estimate): Iv(u -- u,,., )1~,.~ ~< Cllu - u,,. ,11~.,~ .
But using Theorem 3.3 this is just
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-1~6
127
{ In2 N'~2ax, ~'-e~
IV(u - u ,.?, )l,, ,.
ct
Ilull:.,,,
where A~, is the angle parameter from Theorem 3.3 for triangle K~. Hence, we conclude (using similar estimates for the term involving [u - u,,.] in (3.4)) that
I
j
,,,. + l
<-
where A.,i. = min a x . Summing over all edges (and using similar estimates for the boundary terms) gives the estimate //In 2 N'~ a,..,,,~.- 2)
__
}
which completes the proof.
Ilull ,,, []
3.3. Bessel functions In this section we prove convergence of the method using the Bessel function basis given in (2.5). We assume that the solution u of the boundary value problem can be extended to satisfy the Helmholtz equation on a domain /2 containing the closure o f / 2 in its interior. This rules out domains with corners, but we shall discuss this case later. We also assume that the triangulation is sufficiently fine that, for each element, the smallest ball containing K is also contained in jrl. 3.3.1. Convergence under mesh refinement In this section we shall assume that a fixed number N of Bessel functions is used on each triangle. Convergence is obtained by decreasing the mesh size h. We start by proving an estimate of the error in approximating a solution of Laplace's equation by a harmonic polynomial. Results of this type are well known (e.g. [10]). This version has explicitly the dependence on element size and number of basis functions. THEOREM 3.5. Suppose N = 2n + 1 > 1 and w ~ H " + I ( I ~ ) satisfies Aw = 0 in 1). Let K be an element in the mesh of radius h x and let w^,(x, y) = ~ m
~
a~r I'< exp(im0) -
n
where (r, 0) are radial coordinates relative to the centroid of the element K, and r~a,,,~. . . . ,, are arbitrary constants for each triangle K. We have the following approximation result K~n
h ~L____- 0" +__~1 inf WN [{w - WNII,,.K <~ C (n + 1)! llOr n + ] tlo.~ for O<~p <~2 where B K is the ball centered at the centroid o f K and having radius h K. PROOF. Since K C / 2 there is an R > 0 such that the ball Bn of radius R centered at the centroid of K is such that/~ C Bn C / ) . We use radial coordinates centered at the centroid of K and denote by h K the distance from the centroid to the furthest vertex of the triangle. On B R we can expand u as follows: w(r, O) =
k
frhl., Iexp(im0),
w',,k-R}
where the series converges at r = R in the L 2 sense. Now we choose
128
P. Monk, D.-Q. Wang I Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
wx(r, O) =
m = -- n
wm
exp(imO).
Hence
IIw - w~llkK =
=•,_ m I
tl I
( r'~",
w,.\~-)
<~ 2 ~
exp(im0) + m~ , = --~
(Iw,,l" + Iw m] ) -~
~2{n4 1) <~ z~rn x
fh'k"
~
(,
) re-'n+ I
ZTi'F"2i"+l) /K' ~ . . ~'"K . . . ((n + ])!)2 ao m 2"rrhx~''+ 11
.
.
((m
=11+ 1
3" + 1
II .... + , ,
- U, 7
n! (n S-m)!
~
or
dp
[ [ I )(;n-n--l) "~ 2 P w m ~ + w ,, --~
t'X <~ 21rh2^~fl~" ao ...... +~
~<
o
{ r'~rml " exp(im0) ,o.x Wmk, R }
)
2
2 dp
2 .~ p (Iwm] + }W-ml') "~
m, 1)1)9- (Iw,.I 2 + --n-*
(n - m ) !
ni
)
2 dp
~ Iw_,.l 2 , ( P m.-,,-i)2(k (m - , /~-l'-
l)r
)
do
2 IIO,B K
where B,~ is the ball of radius h K centered at the centroid of K. For m = 1 or 2, we note that the intermediate value theorem shows that for any w u (since w u is a solution of Laplace's equation) there is a constant C independent of K, u, N and h~ such that C
II(w - w~),llo.~ ~
IIw- .,&K
and similarly for the y derivative (see e.g. [8 p. 274]).
[]
This theorem has an immediate corollary: C O R O L L A R Y 3.6. Under the assumptions o f Theorem 3.5, f o r 0 ~ p <~ 2 h ,,+j e
0., +
inf IIw - wNII,, <- c - w,v (n + l)r
~ u
0.fi
P R O O F . We square and add the estimates in Theorem 3.5 over all triangles. The overlap of the balls for different K does not change the order of convergence since the grid is regular and quasi-uniform (so there is a minimum element internal diameter and hence a maximum finite number of triangles that can be inside B~ and outside K independent of h). However, the constant in the estimate is increased to allow for this overlap. []
Now, we shall relate the error estimates for Laplace's equation to similar results for the Helmholtz equation. The key tool is the following Vekua transformation [211 which transforms a solution of the Helmholtz equation to the solution of Laplace's equation. For a solution (in radial coordinates) u(r, O) of the Helmholtz equation, we define uo(r, O) = u(r, 0) -
r r 0 , u(p, O) P ~ p Jo(k~fp(p - r)) dp
f{
(3.5)
where Jo is the Bessel function of first kind and order zero. The function u 0 satisfies Au 0 = 0. Vekua shows that this integral identity can be inverted as follows: 0 Jo(k~/r(r , - p)) d O . u(r, O) = uo(r, O) - f ( uo( p, O) -~p More interestingly, if u o is a harmonic monomial so that
(3.6)
P. Monk. D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
uo(r, 0) = r" exp(in0),
129
n t>0
the corresponding solution of the Helmholtz equation is u(r, O) = 2"n !J,(kr) exp(in0 ) .
Thus, the Bessel function basis for the Helmholtz equation analyzed next is equivalent to using a harmonic polynomial basis for approximating Laplace's equation, and we shall use this correspondence to prove error estimates. Using our results for harmonic polynomials and the transformation of Vekua (3.5) we can use arguments like those in [10] to prove the following approximation theorem for solutions of the Helmholtz equation. T H E O R E M 3. Z Suppose n is fixed and n > O. Let u ~ H ~+ t(g)) satisfy Au + k2u = 0 in ~ . Let K be an element it+ the mesh o f radius h K and let amJm(kr) exp(imO)
ux(x, y) = m ~ -- n K
N
wlTere {a. },,=o are arbitrary constants f o r each triangle K. We have the following approximation result inf Ilu - uuIl ,,.K <~ C h ~+ ' - P[lull n + ,.BK , uN
f o r 0 <<-p <~ 2 where B x is the ball centered at the centroid o f K and having radius h K. PROOF. Let us consider the case p = 0 first. We denote by w the solution of Laplace's equation on B e corresponding to u via the Vekua transformation (3.5). Now, we choose w N to be the harmonic polynomial approximation of w given in the proof of Theorem 3.5. Finally, we denote by u~v the Bessel function expansion obtained from w N via the inverse transformation (3.6). Next, we use the boundedness of transformations (3.5) and (3,6) in the Sobolev norms for L2(K) and H"+~(Bx). This follows from the representation results once we assume that all elements have diameter less than a fixed h o < 1. Then we can obtain a uniform bound for the appropriate (analytic) kernel that is independent of the diameter of the element. Using the boundedness, we have
JJu - uNIIo.,,
Crfw - wNl/0,,
f-~Jr ~ |
. . . .
on+
I
w
O'BK
n+l
A .,;imilar argument holds for p = 1 and p = 2.
D
Finally, we can state and prove the main theorem of this section: TH,~OREM 3.8. Assume that u has the necessary smoothness and extension properties given at the beginning o f this section. Let u, be the approximate solution o f ( 1 . 1 ) - ( 1 . 3 ) computed using (2.2) with Bessel functions. We have the following error estimate."
Ilu - u,,llo,,,
Ch°- 'llull,,+ ,
where[2©~. REMARK. It is astonishing that the order of convergence is exactly the same as for plane waves (see (3.3)). The optimality of either result is not known. PROOF. The proof is very similar to that of Theorem 3.4 with some modifications to deal with the fact that the elements are changing size. Let a ' denote the coefficients of the Bessel function expansion used in the proof of Theorem 3.7 on each triangle. Then by Theorem 3.1 it suffices to estimate J(a'). As before we note that it is enough to estimate
P. Monk. D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
130
I[Vu.,l[~,.. + ILlu,,,ll~,. •
(3.7)
Suppose e is common to elements Kj and K e, and denote by uo,~ and u,,,2 the Bessel function approximations on K~ and K 2, respectively
I[Vu,,4,,,. ~< [V(u - u . , , )10., + IV(u - u,, 2)1o~. But by the trace estimate (3.2) the two terms on the right-hand side can be estimated as follows (we only give details of the first estimate): [V(u - u,, , )l,,.~ <~ C(h ~x/,2]l(u - u,,,.~)]lz.,.,x, + h K,' /21](u -- U,'.,r[,.~.K,) But using Theorem 3.7 this can be estimated by I/2
n-I
h
IV(u - u ,, ),,,~ <- C(hK, hK, Ilull,, +, ,¢, + ..,~,
1221/11
-,~, Ilull,,+, ,¢,
Hence, we conclude (using similar estimates for the term involving [u - u,,] in (3.7)) that
I[Vu..l[~.~ + la,.lu,,,ll~.,." <~Ch 2" 'IluIIL, ,~,uK~. Summing over all edges (and using similar estimates for the boundary terms) gives the estimate J(a') <~ Ch ~'- ' l l u l i ~ t , . .
Use of Theorem 3.1 completes the proof.
[]
3.4. Convergence under increasing order
In this section we suppose that the mesh is fixed, and the order of the Bessel function basis (N = 2n + 1) is increased to provide convergence. Again, we start with a result of Babu~ka and Melenk [2] to prove convergence. In particular (using our notation), they state the following result: T H E O R E M 3.9 (Babushka and Melenk, Theorem 5). Let K be an element in the mesh and let the exterior angles o f K be bounded below by ATr, O < a < 2 . Assuming u ~ H ' ( K ) , s > 1 then there is a function u~ = ~N K K p : j a~ @, (x) (see (2.5)) such that
l[u _ u Kll,.x ~< Cik~-)/'ln N\
fluJl.~,K, j
= 0
. . . . . s - 1.
Using this result and our duality estimate exactly as in the proof of Theorem 3.4 we can easily estimate the error in the least squares method: T H E O R E M 3.10. Assume that u E H'(~(2) f o r some s >! 2. Let u, be the approximate solution computed using plane wave ,functions on a fixed mesh. We have the following error estimate: {lnN\,%~,,c, -2~
Ilu - u,,ll0 ,, ~< ck-fi-
)
tlutl,.,,
where Ami,.rr is a lower bound f o r the smallest exterior angle o f any triangle in the mesh (see Theorem 3.9).
4. Numerical results In this section, we describe the implementation of the least-squares method for the Helmholtz equation and discuss some numerical results obtained for a model problem with a known solution.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech, Engrg. 175 (1999) 121-136
131
4.1. M o d e l p r o b l e m
We solve (1.1)-(1.3) in an annular domain 12 C R 2 with the circular boundaries ~ and F~. The radius of F, and F~ are r m and roo ,, respectively. The incident wave is given by g(x) = - u i ( x )
= -exp(-ikx~),
x = (x,,x2).
The computational domain ~(2 is triangulated by a triangular mesh (see Fig. 1). In each triangular cell K, we K approximate the solution u of (1.1)-(1.3) as described in the introduction where the basis functions ~b,,(x) are either plane waves or Bessel functions as described in Section 2. The parameter a,. = k is chosen to obtain a better numerical stability. The other choice mentioned in the Introduction, Ae = 1/h~ does not yield good results. The best choice for A, to obtain good numerical stability is still unknown and will be one of our future research topics. On both boundaries F and I~ curved elements are used. This is easy using the method outlined here since it is only necessary to perform line integrals along the curved boundary. Indeed, the error estimates proved earlier in the paper also hold for curvilinear elements with an appropriate redefinition of the mesh parameter h, etc. The minimization of J(a) results in a positive definite Hermitian linear system in the unknown vector a which is solved by the standard conjugate gradient method. Next, we briefly outline this method for the complexHermitian linear system. Let A be an N × N positive definite Hermitian matrix and x, b ~ C N, then the solution of A x = b is the same as the minimizer of the following quadratic functional 1 x*Ax
-
Re(b'x)
where the asterisk represents the transposed conjugate. The pseudo code for the complex conjugate gradient method is as follows: Initialize
X o a n d Po ~& 0
f o r k = O, I ..... N
do
r k =- b - A x k if
Ilrkl[ ~ l 0
,5 o r pk = = 0, b r e a k
Re(r* Apk ) Ol--
p~Ap~
Pk+l = O t p k - - r k
C
-0.5
-1
-1.fi
-2 -2
-1.5
-1
4.5
0
05
I
1.5
2
Fig. 1. A typical mesh used in the computations presented here.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
132
Re(pk*+ , r , ) fl-
* A P*+! P,+l
xk+l =Xk + fiPk+l end do.
4.2. Results using plane waves In this subsection, we present the numerical results using plane waves. We shall choose N propagation directions dp as follows: dp = (cos(0p), sin(tgp)) t ,
where
0p = 27r(p - 1)/N,
for 1 <~p<~N.
4.2.1. Mesh refinement We compute the least-squares solution u, in a domain /2 with rin = 1 and roo t = 2. For comparison, the exact solution u~ solution is computed through a series expansions in terms of cylindrical Bessel functions. We use three meshes with h ..... = 0.412581 for the coarsest mesh and h m a x = 0.105434 for the finest mesh. The meshes are quite high quality, with maximum aspect ratio of the elements between 4 to 5. We present results for N = 5, 7 and 11 basis functions per element and for wave numbers k = 1, 10. We shall present tables of the relative error (denoted Err) defined by IIU! - - U J[L2{n )
Err =
The relative errors and the average rates of convergence are shown in Tables 1 and 2. The average rates in these experiments almost verify or exceed the estimate of Theorem 3.2 but are rather variable. One possible cause of this might be the rather large condition number of the linear system (particularly for N = 11). Good preconditioners are needed to obtain the desired results when the size of the system increases, but are currently not available for this problem.
4.2.2. Increasing order Tables 1 and 2 also show that the errors go down when the number of the basis functions per element N increase. Table 3 and Fig. 2(a) give the relative errors and the condition numbers of the resulting linear system from the least-squares method when r~n = 0.3, rou , = 1.0 with a wave number of k = 50. We use the a fairly Table I Relative error Err a n d average c o n v e r g e n c e rate when k = 1. The predicted rate is f r o m T h e o r e m 3.2
N = 5 N = 7
0.4125
0.2069
0.1054
Ave. rate
Pred. rate
8.531 e-3 3.747e-3
1.997e-3 9.092e-4
4.923e-4 2.357e-4
2.09 2.07
1 2
Table 2 Relative error Err and average c o n v e r g e n c e rate when k = 10. The predicted rate is f r o m T h e o r e m 3.2 hmax
N = 5
N = 7 N = 11
0.4125
0.2069
0.1054
Ave. rate
Pred. rate
0.795 0.3707 2.351e-3
0.6596 5.027e-2 9.701e-5
0.4246 2.859e-3 2.04e-5
0.461 3.57312 3.46
2 3
l
P. Monk, D.-Q. Wang
/
Comput. Methods Appl. Mech. Engrg. 17,5 (1999) 121-136
133
Tzble 3 Relative error Err and condition number of the matrix for k = 50 with h .... = 0.2246 N C~md. no
9 0.88 1.9e3
11 0.74 2.8e3
13 0.37 4.3e4
15 7.7e-2 1.2e6
N E~r Cond. no
19 1. I e-3 4.4e8
21 9.3e-4 2.5e9
23 9.927e-4 3.3e9
25 9.24e-4 8.1 e9
Err
IJml , ~
L,eml $qu~'e 1MI~ I r ~ e , ~ n g Orders I
0.~
W ~ In~u~in~ ~Ordm'o
7
l
6
!°
0.7
~
17 6.6e-3 2.4e7
0.6
1,
\
• 0.5
/
0,4
!
0.3 0.2 1 01 0
1~0 The nurNoer ot line wave basra
~
26
(a) L i n e a r - l i n e a r .
0
IN, "The num~,r 01 p l m e wave ~ , .
(b) L o g - l i n e a r .
Fig. 2. Plots of the relative error (err) against number of plane waves in the basis N. On the left is a linear-linear plot and on the right a log-linear plot.
e l e m e n t a l basis from N = 9 to 25 and the directions are u n i f o r m l y chosen unit vectors. Despite such a coarse mesh, the error w h e n N = 17 is less than 1%. This discretization corresponds to a p p r o x i m a t e l y 2 degrees of f r e e d o m per wavelength. F i g , 2fb~ is the p l o t o f -In(Err) with resoect o f N. F r o m about N = 13 to N = I8 the g r a o h is nearly linear, which w o u l d be e x p e c t e d if the c o n v e r g e n c e rate is faster than polynomial. As N increases the conditioning of the system deteriorates and this is likely to be the cause o f the s l o w d o w n in c o n v e r g e n c e as N increases. Nevertheless, the results show that the error decreases very rapidly w h e n N increases.
4!.3: l~esut~s ust'ng B'esset'T~mcti'ons We also tested the feast squares s c h e m e s w h e n Besse[ functions o f first ~ ' n d are u s e d as basis timctions on
~b~(x)= J ( p _ e _ ] ) ( k r ) e x p ( i ( p - P)O),
l <~p<<-g,
The s a m e incident w a v e as in the plane w a v e case is used.
,!.~..f. l~es~ refinement W e used the s a m e d o m a i n (}'in = [ and rou t -~ 2) and three m e s h e s as appeared in the plane w a v e case. We restricted ourselves to k = 10 and investigated the use o f N = 7, 9 and 11 basis functions per element. The relative errors and the a v e r a g e rate o f c o n v e r g e n c e are s h o w n in T a b l e 4, For N = 9 we have close to the predicted rate but for N = 7 we o b s e r v e a better rate and for N = 11 a poorer
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
134
Table 4 Relative error Err and average convergence rate when k = 10. The predicted rate of convergence is from Theorem 3.8 h,,,,x
N= 7 N = 9 N = 11
0.4125
0.2069
0.1054
Ave. rate
Pred. rate
0.219 1.83e-2 5.75e-3
2.229e-2 4.504e-4 2.28e-4
1.225e-3 7.01e-5 6.95e-5
3.806 4.063 3.219
2 3 4
Table 5 Relative error Err and condition number of the matrix for k = 50 with h ...... = 0.2246 N
Err Cond. no N
Err Cond. no
II 0.63 3.4e3
13 0.22 5.3e4
15 3.4e-2 1.7e6
17 3.3e-3 4.2e7
19 1.53e-3 3.1e9
21 1.45e-3 1.3ell
23 1.44e-3 9.8el2
25 1.43e-3 9.2e14
Leasl So,au~s vWlhInaeasing the numbw basks (Be~e~ lunction)
~01 01 ~Og( E~rot~ vv~tl ~
~
~
bamil
0.$ ~-
/ 0,5
/ /
l
l
J
/
i ~ 0.3
/ / /
i
/ /
02.
i /
/
(11 ./
15
20 Nunt~er ol Basis
(a) Linear-linear
25
15 "T'heNumt~ oltoeses
i 2O
2G
(b) Log-linear
Fig. 3. Plots of the relative error (err) against number of Bessel functions in the basis N. On the left is a linear-linear plot and on the right a log-linear plot.
rate of convergence than predicted. The latter is most likely due to ill-conditioning of the linear system. Nevertheless, the results s h o w an accuracy of better than 2% for N = 9 and 11 even on the coarsest mesh.
4.3.2. Increasing order
We use the Bessel functions basis to compute the problem with k = 50, r~, = 0.3 and rout = 1 (the same problem in Section 4.2.2). The order of the basis is varied from N = 1 1 to 25. The results are shown in Table 5. As in Section 4.2.2 useful accuracy is attained when N = 17 at a low density of degrees of freedom. Fig. 3(a) is the plot of the relative error for Table 5 with respect to the number of basis functions per element while Fig. 3(b) is the plot of - l o g ( E r r ) with N. Again, there is a linear regime in Fig. 3(b)) suggesting the faster than polynomial convergence predicted by Theorem 3.10. We believe that if a good preconditioner were used, the plot in Fig. 3(b)) should be y = C In(N) which verifies the estimate in Theorem 3.8.
P. Monk. D.-Q. Wang / Comput. Methods" Appl. Mech. Engrg. 175 (1999) 121-136
135
5. Extension to domains with corners For polygonal scatterers, we could try the approach o f Driscoll [9] which takes into account the corner singularities. Briefly, the domain exterior to the polygon is meshed using elements such that each vertex of the polygon is a vertex of the mesh. The element containing a vertex of the polygon must be entirely contained in a disc centered at vertex which does not contain any other vertices. In this element we could use a fractional order Bessel function expansion. This expansion is centered about the domain vertex, takes into account the boundary data exactly, and has the correct singularity. In the remaining elements, we can use the plane wave or Bessei function expansions of this paper. At present it is not known if this approach of using only fractional order Bessel functions at each corner is sufficient to approximate the problem (this is known for the eigenvalue problem studied by Driscoll [9]). This approach is currently under study. If curvilinear edges meet at a corner, the approach of Driscoll may no longer work without modification. We expect to need an augmented fractional order basis, but this needs to be studied further.
6. Conclusion On the theoretical level we have provided an error analysis for the least-squares method in a number of cases. The error analysis suggests a high-order approximation as the dimension of the local basis is increased, and this is observed in practice. There are still many gaps however. In particular, a better understanding of how to precondition the least-squares problem is still needed. On the practical level, we see that the method can suffer from a large condition number but the actual condition number depends on a variety of factors. While both the methods converge rapidly, they are not as robust as the standard finite element method in that the predicted order of convergence is not always seen. It is clear from our results for both bases that for smooth solutions the most efficient way to obtain good accuracy is to increase the dimension of the local basis rather than refine the mesh. Our results suggest that the plane wave basis is more practically useful than the Bessel function basis. This is because various integrals can be evaluated in closed form and plane wave function evaluations are rapid compared to Bessel function evaluations. Yet the error in the solution using either method is comparable. Clearly, an interesting case is the full 3D Maxwell system. In this case a Bessel function basis is easy to construct and will be tested.
Acknowledgment and disclaimer Effort of Peter Monk and Da-Quing Wang sponsored by the Air Force Office of Scientific Research, Air Force Materials Command, USAF, under grant number F49620-96-1-0039. The US Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Air Force Office of Scientific Research or the US Government.
References [1] 1. Babu~ka, F. Ihlenburg, E.T. Paik and S.A. Sauter, A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128 (1995) 325-359. [21 I. Babu~ka and J.M. Melenk, The partition of unity method, Int. J. Numer. Methods Engrg. 40 (1997) 727-758. [3] J. Bramble, R. Lazarov and J. Pasciak, A least-squares approach based on a discrete minus one inner product for first order systems, Math. Comput. 66 (1997) 935-955. [41 J.H. Bramble, R.D. Lazarov and J.E. Pasciak, Least-squares for second-order elliptic problems, Comput. Methods Appl. Mech. Engrg. 152 (1998) 195-210. [51 O. Cessenat and B. Despres, Application of the ultra-weak variational formulation of elliptic PDEs to the 2-dimensional Helmholtz problem, SIAM J. Numer. Anal. 35 (1998) 255-299.
136
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
[6] Y.K. Cheung, W.G. Jin and O.C. Zienkiewicz, Solution of the Helmholtz equation by Trefftz method, Int. J. Numer. Methods Engrg. 32 (1991) 63-78. [7] P.G. Ciarlet and P.-A. Raviart, A mixed finite element method for the biharmonic equation, in: C. de Boor, ed., Mathematical Aspects of Finite Elements in Partial Differential Equations (Academic Press, New York, 1974) 125-145. [8] R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. II (John Wiley, 1966). [9] T. Driscoll, Eigenmodes of isospectral drums, SIAM Rev. 39 (1997) 1-17. [10] S.C. Eisentat, On the rate of convergence of the Bergman-Vekua method for the numerical solution of elliptic boundary value problems, SIAM J. Numer. Anal. 11 (1974) 654-680. [11] I. Herrera, Boundary Methods: An Algebraic Theory (Pitman, 1984). [12] F. lhlenburg and I. Babu~ka, Finite element solution of the Helmholtz equation with high wave number, Part II: The h - p version of the FEM, SIAM J. Numer. Anal. 34 (1997) 315-358. [13] H. Ihlenburg and I. Babu~ka, Finite element solution of the Helmholtz equation with high wave number, Part I: The h-version of the FEM, Comput. Math. Applic. 30 (1995) 9-37. [14] B.-N. Jiang, J. Wu and L.A. Povinelli, The origin of spurious solutions in computational electromagnetism, J. Comput. Phys. 125 (1996) 104-123. [15] Z.-C. Li, R. Mathon and P. Sermer, Boundary methods for solving elliptic problems with singularities and interfaces, SIAM J. Numer. Anal. 24 (1987) 487-498. 116] Zi-Cai Li, Numerical Methods for Elliptic Problems with Singularities (World Scientific, 1993). [17] Zi-Cai Li, Boundary approximation methods for the Helmholtz equation with degeneracy, ZAMM (1996) 453-454. [18] J.M. Melenk and I. Babu~ka, The partition of unity finite element method: Basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289-314. [19] P. Monk and E. Siili, The adaptive computation of far field patterns by a posteriori error estimation of linear functionals, SIAM J. Numer. Anal., to appear. [20] M. Stojek, Least-squares Trefftz-type elements for file Helmholtz equation, Int. J. Numer. Methods Engrg. 41 (1998) 831-849. [21] I.N. Vekua, New Methods for Solving Elliptic Equations (North-Holland, 1967).
EJ ,SEVIER
www.elsevier.com/locate/cma
A least-squares method for the Helmholtz equation P. M o n k * , D a - Q i n g W a n g Department qf Mathematical Sciences, University of Delaware, Newark, DE 19716, USA Received 9 March 1998
Abstract
We investigate the use of least-squares methods to approximate the Helmholtz equation. The basis used in the discrete method consists of st lutions of the Helmholtz equation (either consisting of plane waves or Bessel functions) on each element of a finite element grid. Unlike p~evious methods of this type, we do not use polynomial based finite elements. The use of small elements (and relatively few basis functions per element) allows us to prove convergence theorems for the method and, to some extent, control the conditioning of the resulting linear s} ~tem. Numerical results show the efficiency of the new method and suggest that it may be possible to obtain accurate results with a coarser grid than is usual for standard finite element methods. © 1999 Elsevier Science S.A. All rights reserved.
1. Introduction
When computing the solution of the Helmholtz equation or Maxwell's equations using finite element or finite difference methods we are faced with the problem that there must be a sufficient number of discretization points per wavelength to resolve the solution. A standard rule of thumb is to use at least ten grid points per wavelength, but even this is not sufficient for large problems in which the computational domain is many wavelengths in diameter [13]. This has prompted attempts to use other methods such as corrected least-squares methods, special basis functions or higher-order methods (see e.g. [14,1,12]). In another direction of particular interest to us, Cessenat and Despres [5] suggest using exact solutions of the partial differential equation on each element and choose the expansion coefficients to obtain approximate continuity of the function and its flux across element boundaries. This is done using an approach similar to a mixed finite element method call the 'ultra-weak' wtriational method. The authors report very exciting results on the efficiency and convergence rate of the method. Nevertheless, the method is rather complicated. In this paper we also want to investigate the use of local complete families of solutions of the Helmholtz equation to approximate a scattering problem. But instead of using the 'ultra-weak' variational method of C,zssenat and Despres [5] we will enforce continuity across element boundaries by minimizing a simple quadratic functional. This approach leads to solving a linear system with a Hermitian matrix to obtain the discrete solution, and allows us to use the conjugate gradient scheme to solve the discrete problem. To fix ideas, let /2 C ~2 be a bounded smooth domain with boundary F consisting of two disconnected components F~ and F,. We wish to approximate the solution u of the model scattering problem
Au+k2u=O
in ~ ,
(1.1)
u =g c3u
on F,,
(1.2)
on I_~
(1.3)
;Ju
iku = 0
~' Corresponding author. 00~5-7825/99/$ - see front matter PII: S 0 0 4 5 - 7 8 2 5 ( 9 8 ) 0 0 3 2 6 - 0
© 1999 Elsevier Science S.A. All rights reserved.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 17.5 (1999) 121-136
122
Here, g is a given smooth function and k is a fixed, real positive parameter called the wave-number. The boundary condition (1.3) is a simple absorbing boundary condition. Although not the best absorbing boundary condition available, this condition is simple to implement and we can derive exact solutions on some domains for comparison with numerical results. We shall test two methods for approximating (1.1)-(1.3). In both methods we first mesh the region ,(2 using triangles (we could use more general meshes since there is no need to use a standard finite element grid). The two methods differ according to the expansion functions used on each element. In one case we use a basis of plane waves as in [5] and in the other case we use a Bessel function basis. Using a duality estimate we can estimate the error in the solution from a knowledge of the approximation properties of the appropriate basis. We shall prove the following: (1) For the plane wave basis, convergence under mesh refinement follows from a theorem of [5]. Convergence for a fixed mesh but increasing the number of basis functions is proved using a theorem from 1181. (2) For Bessel functions we prove convergence both under mesh refinement and use a theorem from [18] to prove convergence as the number of basis functions increases on a fixed mesh. Using a fixed grid, our Bessel function method is similar in spirit to the scheme for computing interior eigenvalues presented by Driscoll [9]. If only one element is used the method becomes the Trefftz method analyzed by Eisenstat [10] (see also l11]). Similar methods include the weighted residual method [6] and the boundary approximation method [15-17]. As in [9], the flexibility offered by using a grid can be used to improve the convergence rate of the method. Our estimates, and experience, shows that both the methods we propose converge at high order to smooth solutions of (1.1)-(1.3) if sufficient basis functions are used on each element. A similar approach to ours has been presented in [20]. The least-squares functional is different to ours (in essence we also fit tangential derivatives along element boundaries) and no proof of convergence is given. However, [20] contains several useful extensions of least-squares methods including a discussion of problems with corners (we restrict ourselves to smooth domains). A related method is the partition of unity finite element method of Melenk and Babu~ka [18,2]. We shall use theoretical results from these papers in our work. The method differs from ours in that the partition of unity is used to create a conforming subspace for the standard weak form of (1.1)-(1.3). Hence, a standard Galerkin scheme can be used and the method fits in to the theoretical framework of standard finite element methods. However, the Galerkin formulation does not produce a Hermitian positive definite linear system. Finite element least-squares methods for this problem (see e.g. [14]) and for general elliptic and first-order systems are very well developed (see e.g. [3,4]). The difference here is that we do not use finite elements. We end the paper with a series of numerical which show that the scheme can be efficient. We also comment on the conditioning of the linear system for the problem. We complete the introduction by remarking about the notation for function spaces and norms in the paper. Let S be a Lipschitz domain in the plane with boundary OS. We denote by HP(S) the usual Sobolev space of functions with p derivatives in L2(S). The norm on this space is Similarly, there is also a space of functions H~(OS) on the boundary, and we denote the norm I " I~,.~,s. Sometimes, we will need to use semi-norms. We denote by I/" Ilp,~.s the pth order semi-norm (i.e. the square root of the sum of the squares of the norms of the pth order derivatives).
II-II,, .
2. The method Suppose we cover 12 with a finite element mesh % of regular and quasi-uniform elements of maximum diameter h. In fact we do not need to be so restrictive, since the grid could contain a mixture of element shapes and there is no need for the usual finite element restrictions on vertex placement (see [7]). Here, we shall present results for triangular meshes. The edges in the mesh are divided into three disjoint sets
e
[ Eo 12' ~
if e is in the interior of ~,
boundary
~
if e is on the boundary ~ .
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
123
Each edge e has an associated normal t,. For a sufficiently smooth function f defined in the neighborhood of e we define jump of f across e, for given x E e, by [f](x) =lim (f(x + eu ) - f(x - ev~)) e ~l)
On each element K E 7,, we assume that we have N smooth basis functions {tpeK}p,v= which are solutions of the Helmholtz equation on K, and vanish outside K (so they are discontinuous). Then, we approximate u by u, defined by N
E E a,fC,,f(x), KGvt~
p -
I
where the coefficients aeK need to be determined. We do this by a simple least-squares technique. Let a be the vector of length M = N card(~,) listing all the unknown coefficients. Then, we define J(a)= ~
"~,
9
2
{l[ u,,l[o,,+A~,l[u,,][o,~} + ~
e E=-£11
"~
eEE
2
A~.[(u,,-g)lo,,,+ ~
eEE_~
~Ua
-~-iku,,
I2 O,e
,
(2.1)
where A is a positive coefficient on each edge which might depend on N, k, etc. In this paper we shall choose A,=k. We have tested other choices computationally. For example, we tried ,~ = 1/h~ where h, is the length of edge e but did not find a choice that was a marked improvement over h~ = k. The correct weighting of terms (perhaps involving different norms) is likely to be important to improve conditioning and accuracy of the method. This aspect needs to be investigated further. The solution of the least squares problem a = arg rain J(b)
(2.2)
gives the coefficients of u. and hence the numerical solution of our problem. As discussed in the Introduction, candidates for the basis functions are as follows:
Plane waves Here, we choose N directions ~,, 1 ~
~bj, (x) = exp(ikdp • (x - xK)),
(2.3)
where x s. is the centroid of K. In this paper we shall choose d = (cos(0t,), sin(0p)) r ,
(2.4)
where Op = 2"rr(p - 1)/(N - 1). It is possible that a non-uniform choice of the angle parameter Or (perhaps even limited to a suitable sector) might help satisfy the outgoing radiation condition but we will not investigate this aspect here.
~essel functiorzs Here, we choose N = 2P + 1 then K
~p(x)=,l~p.e_l~(kr)exp(i(p-P)O),
l<~p<~N,
where J, is the cylindrical Bessel function of first kind and order p - P - 1 , cylindrical coordinates with their origin at the centroid x K.
(2.5) and r = I x - x K [ a n d 0 are
:,. E r r o r e s t i m a t e s
We start by proving a simple duality result that provides the basic error estimate for the method. Next, we apply the duality result together with approximation properties of the basis functions to prove the desired error estimates.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
124
THEOREM 3.1. Let u be a smooth solution of (1.1)-(1.3) and let u a be a discrete function computed using either the plane wave or Bessel function basis with coefficients given by a then
Ch "20(a))"2
II. - uoll,,.~ ~
PROOF. Let v satisfy the dual problem At) + kZu = 05
in ~(2,
u=0
on C ,
Ov - - +ikv = 0
on C ,
0v
where 05 is an arbitrary function in Lz(O). Using integration by parts and the fact that on each element u, is a solution of the Helmholtz equation, we can write (U -- U,,, 05) = (U -- a ,,, AV + k2v)
= ~
(u - u,,, Av + k2v)~
K G "i"h
(
io
= ~ (a(u-.,,)+k~(u-uo),v)+ u-uo, ~-~ ~K- --~(U-Uo),V K E,rh
)
OK
Ov K E'r h
Here, (., -)s is the L2(S) inner product (or duality pairing as appropriate) and vx is the unit outward normal to K. Now, rewriting these sums as sums over edges we obtain
e~b"
The true solution, and its gradient, are continuous across interior edges. In addition v = 0 and u - - g on ~. Furthermore, we can rearrange the sum on F= using the absorbing boundary condition to obtain
eEL'()
+
~ (
e
e
~',Ov)
-
E
(( O~_ik)(u_u,),v) i
•
5
e
e~E~
e
Hence, using the Cauchy-Schwarz inequality we obtain
Ov
eEEv
By the Cauchy-Schwarz inequality again
O,e
eEE~
,
P. Monk, D.-Q. Wang / Comput. Methods" AppL Mech. Engrg. 175 (1999) 121-136
)(u - uo, ?))1 ~
ir + 11 I( ) I
~l[u,,[lo~ +
u,,
e - u,I;,,, ~ + E A,,Ig
+ E
125
e~l'-,
~O p -ik
u,~
e~E_~
-Z
O.e
e
I°°1
2
+ J~l 2
O,e
)]
1/2
0,e
Hence, we obtain the bound (3.1)
[(u - u., 6)1 <~(J(a))~'21110111 where
We estimate Jl[vl[I using the trace estimate that for any suitably smooth function w lw]o.~ ~< C x (h s. Vw ~ ~¢ + h~'l}w o x ) ,
e
'
~
t
(3.2)
'
where K~ is a triangle with e as an edge, hx," is the diameter of K, and CK~ is given by
CK'
hK. Px
where PK,, is the radius of the largest inscribed circle in K e (see e.g. [19]). Then, we have that
IIio111 ~ where
c~,
1 (hK, Ilotl2,,K 2 +
h-'l[Vvll~,) + [hK W ~ + h~,,1 Ivlt'~) ~
ll'lJz.~,~,, is the 2-seminorm on K,. Thus, if a,. = k
-112
We obtain the desired estimate by using this estimate in (3.1).
D
3./. Plane waves 7n this section we analyze the use of the plane wave basis given in (2.3) with uniformly distributed propagation angles given by (2.4).
3.1.1. Convergence under mesh refinement ttere, we assume that a fixed N has been chosen, and we study convergence under classical finite element grid refinement (i.e. h--) 0). The error estimates of Cessenat and Despres [5} can be used to prove the following convergence estimate. THEOREM 3.2. Suppose u ~ C ~u ~J2-~(~) then there exists a constant C independent o f h such that
II. - .,Ak,~ ~< c h'~' ,,,2,
(3.3)
REMARK. Suppose N = 2n + I so N is odd. If N = 5 we are guaranteed O(h) convergence, if N = 7 we get O(h:~), and if N = 1 1 we get O(h4). PROOF. In the proof of Theorem 3.7 of Cessenat and Despres [5] it is shown that on there is a function u,, such that if u satisfies (1.1) then for each x E / 2 ,
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
126
lu(x) - u,,(x)[ <~ Ch "+'llullC2,,- ,,,, IV(u(x) - u.,(x))] ~<
Ch"llu[Ic2,, ~,,,,,
Hence, using this estimate, on the edges of the mesh o
Ch~,,+,
eE~ h
ull 2
.
e~Ttl
Using the fact that in a quasi-uniform mesh there are at m o s t C / h 2 edges we obtain "~
E
e ~ rzj
lu -- u ,}~,.~+ IV(u - u .)1o.,, ~< Ch
2n
1
2
}[Ullc2,,+,,~,).
Using this estimate, we can easily see that if A = k and N = 2n + 1 then (min J(a)) 1/2 <~
Ch"-"211ullc2.....
[]
3.2. Convergence tinder increasing order In this case we assume that a fixed grid is chosen, and that convergence is to be obtained by increasing the number of directions per cell N. We can use a Theorem of Babu~ka and Melenk [2] to prove convergence. In particular (using our notation), they state the following result: THEOREM 3.3 (Babu.(ka and Melenk, Theorem 5). Let K be an element in the mesh and let the exterior angles of K be bounded below by ATr, 0 < , t < 2 . Assuming uEH'~(K), s > l then there is a function u x = N K K Y>-~ I al, tgl, (x) (see (2.3)) such that [ In 2 N'~ ~'~ )> [lu - u K,Ili.~ <~ ~ - / Ilull,.K
j = 0 . . . . . s - 1.
Using this result and our duality estimate we can easily estimate the error in the least-squares method: THEOREM 3.4. Assume that u E H"(f2) f o r some s ~ 2. Let u,, be the approximate solution computed using plane wave functions on a .fixed mesh. We have the following error estimate: f ln2 N'~ a,.io,.,. -2~
Ill,
-
u, IJ,,.,, ~< c
~ - )
llull, ,,
where Ami.'rr is a lower bound for the smallest exterior angle of any triangle in the mesh (see Theorem 3.3). PROOF. We use Theorem 3.1. Let a ' denote the coefficients of the plane wave function expansion guaranteed by Theorem 3.4 on each triangle. Then
}lu - u I},,,. <~ Ch '/2(j(a'))'/2 To estimate
J(a') note that it suffices to estimate
I[Vu,,,]h2,.~ + Ih,,[u,,,]], 2,....
(3.4)
Suppose e is common to elements K~ and K 2, and denote by u,,,~ and u , 2 the Bessel function approximations on K~ and K 2, respectively
I~v,,,,llo,,, ~< Iv(u - u,, )l,,.,, + IV(u - u., :)L. But by standard trace results for Sobolev spaces (recalling that the grid is fixed) the two terms on the right-hand side can be estimated as follows (we only give details of the first estimate): Iv(u -- u,,., )1~,.~ ~< Cllu - u,,. ,11~.,~ .
But using Theorem 3.3 this is just
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-1~6
127
{ In2 N'~2ax, ~'-e~
IV(u - u ,.?, )l,, ,.
ct
Ilull:.,,,
where A~, is the angle parameter from Theorem 3.3 for triangle K~. Hence, we conclude (using similar estimates for the term involving [u - u,,.] in (3.4)) that
I
j
,,,. + l
<-
where A.,i. = min a x . Summing over all edges (and using similar estimates for the boundary terms) gives the estimate //In 2 N'~ a,..,,,~.- 2)
__
}
which completes the proof.
Ilull ,,, []
3.3. Bessel functions In this section we prove convergence of the method using the Bessel function basis given in (2.5). We assume that the solution u of the boundary value problem can be extended to satisfy the Helmholtz equation on a domain /2 containing the closure o f / 2 in its interior. This rules out domains with corners, but we shall discuss this case later. We also assume that the triangulation is sufficiently fine that, for each element, the smallest ball containing K is also contained in jrl. 3.3.1. Convergence under mesh refinement In this section we shall assume that a fixed number N of Bessel functions is used on each triangle. Convergence is obtained by decreasing the mesh size h. We start by proving an estimate of the error in approximating a solution of Laplace's equation by a harmonic polynomial. Results of this type are well known (e.g. [10]). This version has explicitly the dependence on element size and number of basis functions. THEOREM 3.5. Suppose N = 2n + 1 > 1 and w ~ H " + I ( I ~ ) satisfies Aw = 0 in 1). Let K be an element in the mesh of radius h x and let w^,(x, y) = ~ m
~
a~r I'< exp(im0) -
n
where (r, 0) are radial coordinates relative to the centroid of the element K, and r~a,,,~. . . . ,, are arbitrary constants for each triangle K. We have the following approximation result K~n
h ~L____- 0" +__~1 inf WN [{w - WNII,,.K <~ C (n + 1)! llOr n + ] tlo.~ for O<~p <~2 where B K is the ball centered at the centroid o f K and having radius h K. PROOF. Since K C / 2 there is an R > 0 such that the ball Bn of radius R centered at the centroid of K is such that/~ C Bn C / ) . We use radial coordinates centered at the centroid of K and denote by h K the distance from the centroid to the furthest vertex of the triangle. On B R we can expand u as follows: w(r, O) =
k
frhl., Iexp(im0),
w',,k-R}
where the series converges at r = R in the L 2 sense. Now we choose
128
P. Monk, D.-Q. Wang I Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
wx(r, O) =
m = -- n
wm
exp(imO).
Hence
IIw - w~llkK =
=•,_ m I
tl I
( r'~",
w,.\~-)
<~ 2 ~
exp(im0) + m~ , = --~
(Iw,,l" + Iw m] ) -~
~2{n4 1) <~ z~rn x
fh'k"
~
(,
) re-'n+ I
ZTi'F"2i"+l) /K' ~ . . ~'"K . . . ((n + ])!)2 ao m 2"rrhx~''+ 11
.
.
((m
=11+ 1
3" + 1
II .... + , ,
- U, 7
n! (n S-m)!
~
or
dp
[ [ I )(;n-n--l) "~ 2 P w m ~ + w ,, --~
t'X <~ 21rh2^~fl~" ao ...... +~
~<
o
{ r'~rml " exp(im0) ,o.x Wmk, R }
)
2
2 dp
2 .~ p (Iwm] + }W-ml') "~
m, 1)1)9- (Iw,.I 2 + --n-*
(n - m ) !
ni
)
2 dp
~ Iw_,.l 2 , ( P m.-,,-i)2(k (m - , /~-l'-
l)r
)
do
2 IIO,B K
where B,~ is the ball of radius h K centered at the centroid of K. For m = 1 or 2, we note that the intermediate value theorem shows that for any w u (since w u is a solution of Laplace's equation) there is a constant C independent of K, u, N and h~ such that C
II(w - w~),llo.~ ~
IIw- .,&K
and similarly for the y derivative (see e.g. [8 p. 274]).
[]
This theorem has an immediate corollary: C O R O L L A R Y 3.6. Under the assumptions o f Theorem 3.5, f o r 0 ~ p <~ 2 h ,,+j e
0., +
inf IIw - wNII,, <- c - w,v (n + l)r
~ u
0.fi
P R O O F . We square and add the estimates in Theorem 3.5 over all triangles. The overlap of the balls for different K does not change the order of convergence since the grid is regular and quasi-uniform (so there is a minimum element internal diameter and hence a maximum finite number of triangles that can be inside B~ and outside K independent of h). However, the constant in the estimate is increased to allow for this overlap. []
Now, we shall relate the error estimates for Laplace's equation to similar results for the Helmholtz equation. The key tool is the following Vekua transformation [211 which transforms a solution of the Helmholtz equation to the solution of Laplace's equation. For a solution (in radial coordinates) u(r, O) of the Helmholtz equation, we define uo(r, O) = u(r, 0) -
r r 0 , u(p, O) P ~ p Jo(k~fp(p - r)) dp
f{
(3.5)
where Jo is the Bessel function of first kind and order zero. The function u 0 satisfies Au 0 = 0. Vekua shows that this integral identity can be inverted as follows: 0 Jo(k~/r(r , - p)) d O . u(r, O) = uo(r, O) - f ( uo( p, O) -~p More interestingly, if u o is a harmonic monomial so that
(3.6)
P. Monk. D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
uo(r, 0) = r" exp(in0),
129
n t>0
the corresponding solution of the Helmholtz equation is u(r, O) = 2"n !J,(kr) exp(in0 ) .
Thus, the Bessel function basis for the Helmholtz equation analyzed next is equivalent to using a harmonic polynomial basis for approximating Laplace's equation, and we shall use this correspondence to prove error estimates. Using our results for harmonic polynomials and the transformation of Vekua (3.5) we can use arguments like those in [10] to prove the following approximation theorem for solutions of the Helmholtz equation. T H E O R E M 3. Z Suppose n is fixed and n > O. Let u ~ H ~+ t(g)) satisfy Au + k2u = 0 in ~ . Let K be an element it+ the mesh o f radius h K and let amJm(kr) exp(imO)
ux(x, y) = m ~ -- n K
N
wlTere {a. },,=o are arbitrary constants f o r each triangle K. We have the following approximation result inf Ilu - uuIl ,,.K <~ C h ~+ ' - P[lull n + ,.BK , uN
f o r 0 <<-p <~ 2 where B x is the ball centered at the centroid o f K and having radius h K. PROOF. Let us consider the case p = 0 first. We denote by w the solution of Laplace's equation on B e corresponding to u via the Vekua transformation (3.5). Now, we choose w N to be the harmonic polynomial approximation of w given in the proof of Theorem 3.5. Finally, we denote by u~v the Bessel function expansion obtained from w N via the inverse transformation (3.6). Next, we use the boundedness of transformations (3.5) and (3,6) in the Sobolev norms for L2(K) and H"+~(Bx). This follows from the representation results once we assume that all elements have diameter less than a fixed h o < 1. Then we can obtain a uniform bound for the appropriate (analytic) kernel that is independent of the diameter of the element. Using the boundedness, we have
JJu - uNIIo.,,
Crfw - wNl/0,,
f-~Jr ~ |
. . . .
on+
I
w
O'BK
n+l
A .,;imilar argument holds for p = 1 and p = 2.
D
Finally, we can state and prove the main theorem of this section: TH,~OREM 3.8. Assume that u has the necessary smoothness and extension properties given at the beginning o f this section. Let u, be the approximate solution o f ( 1 . 1 ) - ( 1 . 3 ) computed using (2.2) with Bessel functions. We have the following error estimate."
Ilu - u,,llo,,,
Ch°- 'llull,,+ ,
where[2©~. REMARK. It is astonishing that the order of convergence is exactly the same as for plane waves (see (3.3)). The optimality of either result is not known. PROOF. The proof is very similar to that of Theorem 3.4 with some modifications to deal with the fact that the elements are changing size. Let a ' denote the coefficients of the Bessel function expansion used in the proof of Theorem 3.7 on each triangle. Then by Theorem 3.1 it suffices to estimate J(a'). As before we note that it is enough to estimate
P. Monk. D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
130
I[Vu.,l[~,.. + ILlu,,,ll~,. •
(3.7)
Suppose e is common to elements Kj and K e, and denote by uo,~ and u,,,2 the Bessel function approximations on K~ and K 2, respectively
I[Vu,,4,,,. ~< [V(u - u . , , )10., + IV(u - u,, 2)1o~. But by the trace estimate (3.2) the two terms on the right-hand side can be estimated as follows (we only give details of the first estimate): [V(u - u,, , )l,,.~ <~ C(h ~x/,2]l(u - u,,,.~)]lz.,.,x, + h K,' /21](u -- U,'.,r[,.~.K,) But using Theorem 3.7 this can be estimated by I/2
n-I
h
IV(u - u ,, ),,,~ <- C(hK, hK, Ilull,, +, ,¢, + ..,~,
1221/11
-,~, Ilull,,+, ,¢,
Hence, we conclude (using similar estimates for the term involving [u - u,,] in (3.7)) that
I[Vu..l[~.~ + la,.lu,,,ll~.,." <~Ch 2" 'IluIIL, ,~,uK~. Summing over all edges (and using similar estimates for the boundary terms) gives the estimate J(a') <~ Ch ~'- ' l l u l i ~ t , . .
Use of Theorem 3.1 completes the proof.
[]
3.4. Convergence under increasing order
In this section we suppose that the mesh is fixed, and the order of the Bessel function basis (N = 2n + 1) is increased to provide convergence. Again, we start with a result of Babu~ka and Melenk [2] to prove convergence. In particular (using our notation), they state the following result: T H E O R E M 3.9 (Babushka and Melenk, Theorem 5). Let K be an element in the mesh and let the exterior angles o f K be bounded below by ATr, O < a < 2 . Assuming u ~ H ' ( K ) , s > 1 then there is a function u~ = ~N K K p : j a~ @, (x) (see (2.5)) such that
l[u _ u Kll,.x ~< Cik~-)/'ln N\
fluJl.~,K, j
= 0
. . . . . s - 1.
Using this result and our duality estimate exactly as in the proof of Theorem 3.4 we can easily estimate the error in the least squares method: T H E O R E M 3.10. Assume that u E H'(~(2) f o r some s >! 2. Let u, be the approximate solution computed using plane wave ,functions on a fixed mesh. We have the following error estimate: {lnN\,%~,,c, -2~
Ilu - u,,ll0 ,, ~< ck-fi-
)
tlutl,.,,
where Ami,.rr is a lower bound f o r the smallest exterior angle o f any triangle in the mesh (see Theorem 3.9).
4. Numerical results In this section, we describe the implementation of the least-squares method for the Helmholtz equation and discuss some numerical results obtained for a model problem with a known solution.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech, Engrg. 175 (1999) 121-136
131
4.1. M o d e l p r o b l e m
We solve (1.1)-(1.3) in an annular domain 12 C R 2 with the circular boundaries ~ and F~. The radius of F, and F~ are r m and roo ,, respectively. The incident wave is given by g(x) = - u i ( x )
= -exp(-ikx~),
x = (x,,x2).
The computational domain ~(2 is triangulated by a triangular mesh (see Fig. 1). In each triangular cell K, we K approximate the solution u of (1.1)-(1.3) as described in the introduction where the basis functions ~b,,(x) are either plane waves or Bessel functions as described in Section 2. The parameter a,. = k is chosen to obtain a better numerical stability. The other choice mentioned in the Introduction, Ae = 1/h~ does not yield good results. The best choice for A, to obtain good numerical stability is still unknown and will be one of our future research topics. On both boundaries F and I~ curved elements are used. This is easy using the method outlined here since it is only necessary to perform line integrals along the curved boundary. Indeed, the error estimates proved earlier in the paper also hold for curvilinear elements with an appropriate redefinition of the mesh parameter h, etc. The minimization of J(a) results in a positive definite Hermitian linear system in the unknown vector a which is solved by the standard conjugate gradient method. Next, we briefly outline this method for the complexHermitian linear system. Let A be an N × N positive definite Hermitian matrix and x, b ~ C N, then the solution of A x = b is the same as the minimizer of the following quadratic functional 1 x*Ax
-
Re(b'x)
where the asterisk represents the transposed conjugate. The pseudo code for the complex conjugate gradient method is as follows: Initialize
X o a n d Po ~& 0
f o r k = O, I ..... N
do
r k =- b - A x k if
Ilrkl[ ~ l 0
,5 o r pk = = 0, b r e a k
Re(r* Apk ) Ol--
p~Ap~
Pk+l = O t p k - - r k
C
-0.5
-1
-1.fi
-2 -2
-1.5
-1
4.5
0
05
I
1.5
2
Fig. 1. A typical mesh used in the computations presented here.
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
132
Re(pk*+ , r , ) fl-
* A P*+! P,+l
xk+l =Xk + fiPk+l end do.
4.2. Results using plane waves In this subsection, we present the numerical results using plane waves. We shall choose N propagation directions dp as follows: dp = (cos(0p), sin(tgp)) t ,
where
0p = 27r(p - 1)/N,
for 1 <~p<~N.
4.2.1. Mesh refinement We compute the least-squares solution u, in a domain /2 with rin = 1 and roo t = 2. For comparison, the exact solution u~ solution is computed through a series expansions in terms of cylindrical Bessel functions. We use three meshes with h ..... = 0.412581 for the coarsest mesh and h m a x = 0.105434 for the finest mesh. The meshes are quite high quality, with maximum aspect ratio of the elements between 4 to 5. We present results for N = 5, 7 and 11 basis functions per element and for wave numbers k = 1, 10. We shall present tables of the relative error (denoted Err) defined by IIU! - - U J[L2{n )
Err =
The relative errors and the average rates of convergence are shown in Tables 1 and 2. The average rates in these experiments almost verify or exceed the estimate of Theorem 3.2 but are rather variable. One possible cause of this might be the rather large condition number of the linear system (particularly for N = 11). Good preconditioners are needed to obtain the desired results when the size of the system increases, but are currently not available for this problem.
4.2.2. Increasing order Tables 1 and 2 also show that the errors go down when the number of the basis functions per element N increase. Table 3 and Fig. 2(a) give the relative errors and the condition numbers of the resulting linear system from the least-squares method when r~n = 0.3, rou , = 1.0 with a wave number of k = 50. We use the a fairly Table I Relative error Err a n d average c o n v e r g e n c e rate when k = 1. The predicted rate is f r o m T h e o r e m 3.2
N = 5 N = 7
0.4125
0.2069
0.1054
Ave. rate
Pred. rate
8.531 e-3 3.747e-3
1.997e-3 9.092e-4
4.923e-4 2.357e-4
2.09 2.07
1 2
Table 2 Relative error Err and average c o n v e r g e n c e rate when k = 10. The predicted rate is f r o m T h e o r e m 3.2 hmax
N = 5
N = 7 N = 11
0.4125
0.2069
0.1054
Ave. rate
Pred. rate
0.795 0.3707 2.351e-3
0.6596 5.027e-2 9.701e-5
0.4246 2.859e-3 2.04e-5
0.461 3.57312 3.46
2 3
l
P. Monk, D.-Q. Wang
/
Comput. Methods Appl. Mech. Engrg. 17,5 (1999) 121-136
133
Tzble 3 Relative error Err and condition number of the matrix for k = 50 with h .... = 0.2246 N C~md. no
9 0.88 1.9e3
11 0.74 2.8e3
13 0.37 4.3e4
15 7.7e-2 1.2e6
N E~r Cond. no
19 1. I e-3 4.4e8
21 9.3e-4 2.5e9
23 9.927e-4 3.3e9
25 9.24e-4 8.1 e9
Err
IJml , ~
L,eml $qu~'e 1MI~ I r ~ e , ~ n g Orders I
0.~
W ~ In~u~in~ ~Ordm'o
7
l
6
!°
0.7
~
17 6.6e-3 2.4e7
0.6
1,
\
• 0.5
/
0,4
!
0.3 0.2 1 01 0
1~0 The nurNoer ot line wave basra
~
26
(a) L i n e a r - l i n e a r .
0
IN, "The num~,r 01 p l m e wave ~ , .
(b) L o g - l i n e a r .
Fig. 2. Plots of the relative error (err) against number of plane waves in the basis N. On the left is a linear-linear plot and on the right a log-linear plot.
e l e m e n t a l basis from N = 9 to 25 and the directions are u n i f o r m l y chosen unit vectors. Despite such a coarse mesh, the error w h e n N = 17 is less than 1%. This discretization corresponds to a p p r o x i m a t e l y 2 degrees of f r e e d o m per wavelength. F i g , 2fb~ is the p l o t o f -In(Err) with resoect o f N. F r o m about N = 13 to N = I8 the g r a o h is nearly linear, which w o u l d be e x p e c t e d if the c o n v e r g e n c e rate is faster than polynomial. As N increases the conditioning of the system deteriorates and this is likely to be the cause o f the s l o w d o w n in c o n v e r g e n c e as N increases. Nevertheless, the results show that the error decreases very rapidly w h e n N increases.
4!.3: l~esut~s ust'ng B'esset'T~mcti'ons We also tested the feast squares s c h e m e s w h e n Besse[ functions o f first ~ ' n d are u s e d as basis timctions on
~b~(x)= J ( p _ e _ ] ) ( k r ) e x p ( i ( p - P)O),
l <~p<<-g,
The s a m e incident w a v e as in the plane w a v e case is used.
,!.~..f. l~es~ refinement W e used the s a m e d o m a i n (}'in = [ and rou t -~ 2) and three m e s h e s as appeared in the plane w a v e case. We restricted ourselves to k = 10 and investigated the use o f N = 7, 9 and 11 basis functions per element. The relative errors and the a v e r a g e rate o f c o n v e r g e n c e are s h o w n in T a b l e 4, For N = 9 we have close to the predicted rate but for N = 7 we o b s e r v e a better rate and for N = 11 a poorer
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
134
Table 4 Relative error Err and average convergence rate when k = 10. The predicted rate of convergence is from Theorem 3.8 h,,,,x
N= 7 N = 9 N = 11
0.4125
0.2069
0.1054
Ave. rate
Pred. rate
0.219 1.83e-2 5.75e-3
2.229e-2 4.504e-4 2.28e-4
1.225e-3 7.01e-5 6.95e-5
3.806 4.063 3.219
2 3 4
Table 5 Relative error Err and condition number of the matrix for k = 50 with h ...... = 0.2246 N
Err Cond. no N
Err Cond. no
II 0.63 3.4e3
13 0.22 5.3e4
15 3.4e-2 1.7e6
17 3.3e-3 4.2e7
19 1.53e-3 3.1e9
21 1.45e-3 1.3ell
23 1.44e-3 9.8el2
25 1.43e-3 9.2e14
Leasl So,au~s vWlhInaeasing the numbw basks (Be~e~ lunction)
~01 01 ~Og( E~rot~ vv~tl ~
~
~
bamil
0.$ ~-
/ 0,5
/ /
l
l
J
/
i ~ 0.3
/ / /
i
/ /
02.
i /
/
(11 ./
15
20 Nunt~er ol Basis
(a) Linear-linear
25
15 "T'heNumt~ oltoeses
i 2O
2G
(b) Log-linear
Fig. 3. Plots of the relative error (err) against number of Bessel functions in the basis N. On the left is a linear-linear plot and on the right a log-linear plot.
rate of convergence than predicted. The latter is most likely due to ill-conditioning of the linear system. Nevertheless, the results s h o w an accuracy of better than 2% for N = 9 and 11 even on the coarsest mesh.
4.3.2. Increasing order
We use the Bessel functions basis to compute the problem with k = 50, r~, = 0.3 and rout = 1 (the same problem in Section 4.2.2). The order of the basis is varied from N = 1 1 to 25. The results are shown in Table 5. As in Section 4.2.2 useful accuracy is attained when N = 17 at a low density of degrees of freedom. Fig. 3(a) is the plot of the relative error for Table 5 with respect to the number of basis functions per element while Fig. 3(b) is the plot of - l o g ( E r r ) with N. Again, there is a linear regime in Fig. 3(b)) suggesting the faster than polynomial convergence predicted by Theorem 3.10. We believe that if a good preconditioner were used, the plot in Fig. 3(b)) should be y = C In(N) which verifies the estimate in Theorem 3.8.
P. Monk. D.-Q. Wang / Comput. Methods" Appl. Mech. Engrg. 175 (1999) 121-136
135
5. Extension to domains with corners For polygonal scatterers, we could try the approach o f Driscoll [9] which takes into account the corner singularities. Briefly, the domain exterior to the polygon is meshed using elements such that each vertex of the polygon is a vertex of the mesh. The element containing a vertex of the polygon must be entirely contained in a disc centered at vertex which does not contain any other vertices. In this element we could use a fractional order Bessel function expansion. This expansion is centered about the domain vertex, takes into account the boundary data exactly, and has the correct singularity. In the remaining elements, we can use the plane wave or Bessei function expansions of this paper. At present it is not known if this approach of using only fractional order Bessel functions at each corner is sufficient to approximate the problem (this is known for the eigenvalue problem studied by Driscoll [9]). This approach is currently under study. If curvilinear edges meet at a corner, the approach of Driscoll may no longer work without modification. We expect to need an augmented fractional order basis, but this needs to be studied further.
6. Conclusion On the theoretical level we have provided an error analysis for the least-squares method in a number of cases. The error analysis suggests a high-order approximation as the dimension of the local basis is increased, and this is observed in practice. There are still many gaps however. In particular, a better understanding of how to precondition the least-squares problem is still needed. On the practical level, we see that the method can suffer from a large condition number but the actual condition number depends on a variety of factors. While both the methods converge rapidly, they are not as robust as the standard finite element method in that the predicted order of convergence is not always seen. It is clear from our results for both bases that for smooth solutions the most efficient way to obtain good accuracy is to increase the dimension of the local basis rather than refine the mesh. Our results suggest that the plane wave basis is more practically useful than the Bessel function basis. This is because various integrals can be evaluated in closed form and plane wave function evaluations are rapid compared to Bessel function evaluations. Yet the error in the solution using either method is comparable. Clearly, an interesting case is the full 3D Maxwell system. In this case a Bessel function basis is easy to construct and will be tested.
Acknowledgment and disclaimer Effort of Peter Monk and Da-Quing Wang sponsored by the Air Force Office of Scientific Research, Air Force Materials Command, USAF, under grant number F49620-96-1-0039. The US Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Air Force Office of Scientific Research or the US Government.
References [1] 1. Babu~ka, F. Ihlenburg, E.T. Paik and S.A. Sauter, A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128 (1995) 325-359. [21 I. Babu~ka and J.M. Melenk, The partition of unity method, Int. J. Numer. Methods Engrg. 40 (1997) 727-758. [3] J. Bramble, R. Lazarov and J. Pasciak, A least-squares approach based on a discrete minus one inner product for first order systems, Math. Comput. 66 (1997) 935-955. [41 J.H. Bramble, R.D. Lazarov and J.E. Pasciak, Least-squares for second-order elliptic problems, Comput. Methods Appl. Mech. Engrg. 152 (1998) 195-210. [51 O. Cessenat and B. Despres, Application of the ultra-weak variational formulation of elliptic PDEs to the 2-dimensional Helmholtz problem, SIAM J. Numer. Anal. 35 (1998) 255-299.
136
P. Monk, D.-Q. Wang / Comput. Methods Appl. Mech. Engrg. 175 (1999) 121-136
[6] Y.K. Cheung, W.G. Jin and O.C. Zienkiewicz, Solution of the Helmholtz equation by Trefftz method, Int. J. Numer. Methods Engrg. 32 (1991) 63-78. [7] P.G. Ciarlet and P.-A. Raviart, A mixed finite element method for the biharmonic equation, in: C. de Boor, ed., Mathematical Aspects of Finite Elements in Partial Differential Equations (Academic Press, New York, 1974) 125-145. [8] R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. II (John Wiley, 1966). [9] T. Driscoll, Eigenmodes of isospectral drums, SIAM Rev. 39 (1997) 1-17. [10] S.C. Eisentat, On the rate of convergence of the Bergman-Vekua method for the numerical solution of elliptic boundary value problems, SIAM J. Numer. Anal. 11 (1974) 654-680. [11] I. Herrera, Boundary Methods: An Algebraic Theory (Pitman, 1984). [12] F. lhlenburg and I. Babu~ka, Finite element solution of the Helmholtz equation with high wave number, Part II: The h - p version of the FEM, SIAM J. Numer. Anal. 34 (1997) 315-358. [13] H. Ihlenburg and I. Babu~ka, Finite element solution of the Helmholtz equation with high wave number, Part I: The h-version of the FEM, Comput. Math. Applic. 30 (1995) 9-37. [14] B.-N. Jiang, J. Wu and L.A. Povinelli, The origin of spurious solutions in computational electromagnetism, J. Comput. Phys. 125 (1996) 104-123. [15] Z.-C. Li, R. Mathon and P. Sermer, Boundary methods for solving elliptic problems with singularities and interfaces, SIAM J. Numer. Anal. 24 (1987) 487-498. 116] Zi-Cai Li, Numerical Methods for Elliptic Problems with Singularities (World Scientific, 1993). [17] Zi-Cai Li, Boundary approximation methods for the Helmholtz equation with degeneracy, ZAMM (1996) 453-454. [18] J.M. Melenk and I. Babu~ka, The partition of unity finite element method: Basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289-314. [19] P. Monk and E. Siili, The adaptive computation of far field patterns by a posteriori error estimation of linear functionals, SIAM J. Numer. Anal., to appear. [20] M. Stojek, Least-squares Trefftz-type elements for file Helmholtz equation, Int. J. Numer. Methods Engrg. 41 (1998) 831-849. [21] I.N. Vekua, New Methods for Solving Elliptic Equations (North-Holland, 1967).