An Improved Symmetric Positive-Definite Finite Element Method for the Complex Helmholtz Equation

0.54 MB·English

by  Russell B. Richins

AN IMPROVED SYMMETRIC POSITIVE-DEFINITE FINITE
ELEMENT METHOD FOR THE COMPLEX HELMHOLTZ
EQUATION
RUSSELL B. RICHINS
3
1
0
2
Abstract. Most discretizations of the Helmholtz equation result in a system of
n
linearequationsthathasanindefinitecoefficientmatrix. Muchefforthasbeenput
a
J into solving such systems of equations efficiently. In a previous work, the current
author and D.C. Dobson proposed a numerical method for solving the complex
8
1 Helmholtz equation based on the minimization variational principles developed by
Milton, Seppecher, and Bouchitt´e. This method results in a system of equations

] withasymmetricpositivedefinitecoefficientmatrix,butatthesametimerequires
A
solving simultaneously for the solution and its gradient, resulting in a much larger
N
system of equations. Herein is presented a method based on the saddle point vari-
h. ational principles of Milton, Seppecher, and Bouchitt´e, which produces symmetric
t positive definite systems of equations, but eliminates the necessity of solving for
a
m the gradient of the solution.
[

1 1. Introduction
v
5
The Helmholtz equation
3
4 ∇·L∇u = Mu,
4
isusefulinmodelingwavepropagationinproblemsarisingfrommanydifferentphysi-
.
1
cal situations. Suppose we wish to solve the Helmholtz equation in a domain Ω ⊂ Rd,
0
3 and assume that L and M are complex valued functions. A common source of nu-
1
merical methods for solving this equation is the variational principle
:
v (cid:90)
i (1) [−L∇u·∇v¯−Muv¯]dx = 0 ∀ v ∈ H1(Ω).
X 0
Ω
r
a Since this is a stationary principle, the resulting system of equations is indefinite,
and is generally more difficult to solve than a system of equations having a positive
definite coefficient matrix.
In their paper [7], Milton, Seppecher, and Bouchitt´e developed variational principles
that apply to the Helmholtz equation above, as well as the time-harmonic Maxwell
2010 Mathematics Subject Classification. Primary 65N30; Secondary 35A15.
1
2 RUSSELL B. RICHINS
equations and the equations of linear elasticity in lossy materials. To derive these
variational principles, we first define the dual variable
v = iL∇u.
Then
(cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
L 0 ∇u L∇u −iv
= = ,
0 M u Mu −i∇·v
or equivalently,
G = ZF,
where
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
∇u −iv L 0
F = , G = , Z = .
u −i∇·v 0 M
For a complex quantity z, we will write z(cid:48) = Re(z) and z(cid:48)(cid:48) = Im(z). Taking real and
imaginary parts, the constitutive relation becomes
G(cid:48) = Z(cid:48)F(cid:48) −Z(cid:48)(cid:48)F(cid:48)(cid:48) and G(cid:48)(cid:48) = Z(cid:48)F(cid:48)(cid:48) +Z(cid:48)(cid:48)F(cid:48),
which can be written in matrix form as
(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
G(cid:48)(cid:48) Z(cid:48)(cid:48) Z(cid:48) F(cid:48)
(2) = .
G(cid:48) Z(cid:48) −Z(cid:48)(cid:48) F(cid:48)(cid:48)
Solving this relation for the imaginary parts of F and G, we find that
(cid:18) (cid:19) (cid:18) (cid:19)
G(cid:48)(cid:48) F(cid:48)
(3) = L ,
F(cid:48)(cid:48) −G(cid:48)
where
(cid:18) (cid:19)
Z(cid:48)(cid:48) +Z(cid:48)(Z(cid:48)(cid:48))−1Z(cid:48) Z(cid:48)(Z(cid:48)(cid:48))−1
L = .
(Z(cid:48)(cid:48))−1Z(cid:48) (Z(cid:48)(cid:48))−1
The matrix L is positive definite as long as Z(cid:48)(cid:48) is positive definite (see [7]).
The previous approach in [8] was to use this constitutive relation and the correspond-
ing energy functional
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
F(cid:48) Z(cid:48)(cid:48) +Z(cid:48)(Z(cid:48)(cid:48))−1Z(cid:48) Z(cid:48)(Z(cid:48)(cid:48))−1 F(cid:48)
· dx
−G(cid:48) (Z(cid:48)(cid:48))−1Z(cid:48) (Z(cid:48)(cid:48))−1 −G(cid:48)
Ω
to formulate a numerical method. When this variational principle is discretized by
the finite element method, the result is a system of equations that has the block
form
    
A AT AT uˆ(cid:48) b
1 4 6 1
 A4 A2 AT5  vˆ1(cid:48)(cid:48)  =  b2 .
A A A vˆ(cid:48)(cid:48) b
6 5 3 2 3
AN IMPROVED METHOD FOR THE HELMHOLTZ EQUATION 3
A similar system of equations must be solved to find approximations for u(cid:48)(cid:48) and v(cid:48).
In all, to find u(cid:48) and u(cid:48)(cid:48), one must solve two positive definite systems of equations of
size 3N ×3N, where N is the number of nodes in the computational grid.
In this work, we develop a new method based on the saddle point variational prin-
ciples in [7] that does not require that v be solved for in order to find u. In this
formulation, we find u(cid:48) and u(cid:48)(cid:48) by solving two N × N symmetric positive definite
systems of equations. Since these are real matrices, the amount of work in solving
the two systems can be expected to be less than that required to solve the complex
N ×N system that results from (1) before taking positivity into account, since com-
plex multiplication is more than twice as expensive than real multiplication.
First, in Section 2, we will derive the saddle point variational principles from [7]
upon which our method is based. In Section 3, we will discuss the details of handling
Dirichlet, Neumann, and Robin boundary conditions with these variational princi-
ples. Section 4 contains the derivation of a standard bound on the error incurred
when the Helmholtz equation is solved using a finite element method that discretizes
the saddle-point variational principle. Section 5 outlines the numerical method and
discusses the conditioning of the system. In Section 6, we provide several numerical
examples, as well as numerical verification of the error bound from Section 4.
2. The Saddle Point Variational Principle
The derivation of the saddle point variational principle from [7] follows the same
steps presented in the introduction for the minimization principle. The difference is
that instead of continuing to the constitutive relation (3), we stop at equation (2).
As usual, we assume that Z(cid:48)(cid:48) is positive definite. From this constitutive relation, we
define the functional
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) F(cid:48)
Y(u(cid:48),u(cid:48)(cid:48)) = · dx.
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) F(cid:48)(cid:48)
Ω
Let u(cid:48),u(cid:48)(cid:48) ∈ H1(Ω) be the real and imaginary parts of a solution to the Helmholtz
equation. Let s ∈ H1(Ω) and define
0
(cid:18) (cid:19)
∇s
S = .
s
Then we have
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
F(cid:48) +S Z(cid:48)(cid:48) Z(cid:48) F(cid:48) +S
Y(u(cid:48) +s,u(cid:48)(cid:48)) = · dx
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) F(cid:48)(cid:48)
Ω
4 RUSSELL B. RICHINS
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) S
= Y(u(cid:48),u(cid:48)(cid:48))+2 · dx+Y(s,0).
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) 0
Ω
The integral in the line above can be rewritten as
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19) (cid:90) (cid:18) (cid:19) (cid:18) (cid:19)
G(cid:48)(cid:48) S −v(cid:48) ∇s
· dx = · dx
G(cid:48) 0 −∇·v(cid:48) s
Ω Ω
(cid:90) (cid:90) (cid:90)
= [−v(cid:48) ·∇s−∇·v(cid:48)s] dx = −∇·[v(cid:48)s] dx = −v(cid:48) ·ns dS = 0.
Ω Ω ∂Ω
Therefore,
(cid:90)
Y(u(cid:48) +s,u(cid:48)(cid:48)) = Y(u(cid:48),u(cid:48)(cid:48))+ S ·Z(cid:48)(cid:48)S dx,
Ω
and the last term must be nonnegative, since Z(cid:48)(cid:48) is assumed to be positive definite.
A similar calculation yields
(cid:90)
Y(u(cid:48),u(cid:48)(cid:48) +s) = Y(u(cid:48),u(cid:48)(cid:48))− S ·Z(cid:48)(cid:48)S dx.
Ω
This shows that (u(cid:48),u(cid:48)(cid:48)) is at a saddle point of the functional Y.
Suppose that (u(cid:48),u(cid:48)(cid:48)) is a saddle point of the functional Y. Then the functional
Q(s(cid:48),s(cid:48)(cid:48)) = Y(u(cid:48) + s(cid:48),u(cid:48)(cid:48) + s(cid:48)(cid:48)), defined for all s(cid:48),s(cid:48)(cid:48) ∈ H1(Ω), should have a saddle
0
point at s(cid:48) = s(cid:48)(cid:48) = 0. A necessary condition for this to happen is that the first
variation of Q should vanish. Let
(cid:18) (cid:19) (cid:18) (cid:19)
∇s(cid:48) ∇s(cid:48)(cid:48)
S(cid:48) = and S(cid:48)(cid:48) = .
s(cid:48) s(cid:48)(cid:48)
Then we must have
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) S(cid:48)
0 = · dx
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) S(cid:48)(cid:48)
Ω
(cid:90)
= [F(cid:48) ·Z(cid:48)(cid:48)S(cid:48) +F(cid:48) ·Z(cid:48)S(cid:48)(cid:48) +F(cid:48)(cid:48) ·Z(cid:48)S(cid:48) −F(cid:48)(cid:48) ·Z(cid:48)(cid:48)S(cid:48)(cid:48)] dx
Ω
(cid:90) (cid:20)(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
∇u(cid:48) L(cid:48)(cid:48) 0 ∇s(cid:48) ∇u(cid:48) L(cid:48) 0 ∇s(cid:48)(cid:48)
= · + ·
u(cid:48) 0 M(cid:48)(cid:48) s(cid:48) u(cid:48) 0 M(cid:48) s(cid:48)(cid:48)
Ω
(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)(cid:21)
∇u(cid:48)(cid:48) L(cid:48) 0 ∇s(cid:48) ∇u(cid:48)(cid:48) L(cid:48)(cid:48) 0 ∇s(cid:48)(cid:48)
+ · − · dx
u(cid:48)(cid:48) 0 M(cid:48) s(cid:48) u(cid:48)(cid:48) 0 M(cid:48)(cid:48) s(cid:48)(cid:48)
(cid:90)
= [∇u(cid:48) ·L(cid:48)(cid:48)∇s(cid:48) +u(cid:48)M(cid:48)(cid:48)s(cid:48) +∇u(cid:48) ·L(cid:48)∇s(cid:48)(cid:48) +u(cid:48)M(cid:48)s(cid:48)(cid:48) +∇u(cid:48)(cid:48)·
Ω
L(cid:48)∇s(cid:48) +u(cid:48)(cid:48)M(cid:48)s(cid:48) −∇u(cid:48)(cid:48) ·L(cid:48)(cid:48)∇s(cid:48)(cid:48) −u(cid:48)(cid:48) ·M(cid:48)(cid:48)s(cid:48)(cid:48)] dx
(cid:90)
= [(−∇·L(cid:48)(cid:48)∇u(cid:48) +M(cid:48)(cid:48)u(cid:48) −∇·L(cid:48)∇u(cid:48)(cid:48) +M(cid:48)u(cid:48)(cid:48))s(cid:48)
Ω
AN IMPROVED METHOD FOR THE HELMHOLTZ EQUATION 5
+(−∇·L(cid:48)∇u(cid:48) +M(cid:48)u(cid:48) +∇·L(cid:48)(cid:48)∇u(cid:48)(cid:48) −M(cid:48)(cid:48)u(cid:48)(cid:48))s(cid:48)(cid:48)] dx.
If we rewrite the equation ∇·L∇u = Mu in terms of real and imaginary parts, we
arrive at
∇·[L(cid:48)∇u(cid:48) −L(cid:48)(cid:48)∇u(cid:48)(cid:48) +i(L(cid:48)∇u(cid:48)(cid:48) +L(cid:48)(cid:48)∇u(cid:48))] = M(cid:48)u(cid:48) −M(cid:48)(cid:48)u(cid:48)(cid:48) +i(M(cid:48)u(cid:48)(cid:48) +M(cid:48)(cid:48)u(cid:48)),
or in other words,
∇·L(cid:48)∇u(cid:48) −∇·L(cid:48)(cid:48)∇u(cid:48)(cid:48) −M(cid:48)u(cid:48) +M(cid:48)(cid:48)u(cid:48)(cid:48) = 0
and
∇·L(cid:48)∇u(cid:48)(cid:48) +∇·L(cid:48)(cid:48)∇u(cid:48) −M(cid:48)u(cid:48)(cid:48) −M(cid:48)(cid:48)u(cid:48) = 0.
Notice that the left-hand sides of these equations are just the opposites of the ex-
pressions multiplying s(cid:48) and s(cid:48)(cid:48) in the integral above. Since the result of the integral
must be zero regardless of the choice of s(cid:48) and s(cid:48)(cid:48), the saddle point of Y must be a
solution to the Helmholtz equation.
3. Boundary Conditions
The calculations done above show that a saddle point of Y satisfying u(cid:48) = f(cid:48) and
u(cid:48)(cid:48) = f(cid:48)(cid:48) on ∂Ω is a solution of
(cid:26)
∇·L∇u = Mu in Ω
.
u = f(cid:48) +if(cid:48)(cid:48) on ∂Ω
We can also solve the Neumann problem
(cid:26)
∇·L∇u = Mu in Ω
.
v ·n = g(cid:48) +ig(cid:48)(cid:48) on ∂Ω
Let s(cid:48),s(cid:48)(cid:48) ∈ H1(Ω) be arbitrary test functions. Then we have
(cid:90)
0 = [(∇·v(cid:48) −∇·v(cid:48))s(cid:48) +(−∇·v(cid:48)(cid:48) +∇·v(cid:48)(cid:48))s(cid:48)(cid:48)] dx
Ω
(cid:90) (cid:90)
= [−v(cid:48) ·∇s(cid:48) −∇·v(cid:48)s(cid:48) +v(cid:48)(cid:48) ·∇s(cid:48)(cid:48) +∇·v(cid:48)(cid:48)s(cid:48)(cid:48)] dx+ [s(cid:48)v(cid:48) ·n−s(cid:48)(cid:48)v(cid:48)(cid:48) ·n] dx
Ω ∂Ω
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19) (cid:90)
G(cid:48)(cid:48) S(cid:48)
= · dx+ [s(cid:48)v(cid:48) ·n−s(cid:48)(cid:48)v(cid:48)(cid:48) ·n] dS
G(cid:48) S(cid:48)(cid:48)
Ω ∂Ω
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:90)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) S(cid:48)
= · dx+ [s(cid:48)v(cid:48) ·n−s(cid:48)(cid:48)v(cid:48)(cid:48) ·n] dS.
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) S(cid:48)(cid:48)
Ω ∂Ω
Therefore, in order to solve the Neumann problem, we solve the weak equation
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:90)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) S(cid:48)
· dx = [−s(cid:48)g(cid:48) +s(cid:48)(cid:48)g(cid:48)(cid:48)] dS ∀ s(cid:48),s(cid:48)(cid:48) ∈ H1(Ω).
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) S(cid:48)(cid:48)
Ω ∂Ω
6 RUSSELL B. RICHINS
To solve the Robin problem
(cid:26)
∇·L∇u = Mu in Ω
,
u+av ·n = g on ∂Ω
we begin with the weak form of the Neumann problem, which we will write as
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:90) (cid:18) (cid:19) (cid:18) (cid:19)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) S(cid:48) s(cid:48) v(cid:48) ·n
0 = · dx+ · dS
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) S(cid:48)(cid:48) s(cid:48)(cid:48) −v(cid:48)(cid:48) ·n
Ω ∂Ω
We split the boundary condition into its real and imaginary parts as
u(cid:48) +a(cid:48)v(cid:48) ·n−a(cid:48)(cid:48)v(cid:48)(cid:48) ·n = g(cid:48)
u(cid:48)(cid:48) +a(cid:48)v(cid:48)(cid:48) ·n+a(cid:48)(cid:48)v(cid:48) ·n = g(cid:48)(cid:48),
which we can write as
(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19)
u(cid:48) a(cid:48) a(cid:48)(cid:48) v(cid:48) ·n g(cid:48)
+ = .
u(cid:48)(cid:48) a(cid:48)(cid:48) −a(cid:48) −v(cid:48)(cid:48) ·n g(cid:48)(cid:48)
If the matrix in the equation above is called M, then
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
v(cid:48) ·n u(cid:48) g(cid:48)
= −M−1 +M−1 ,
−v(cid:48)(cid:48) ·n u(cid:48)(cid:48) g(cid:48)(cid:48)
so the weak form of the equation with Robin boundary conditions is
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:90) (cid:18) (cid:19) (cid:18) (cid:19)
F(cid:48) Z(cid:48)(cid:48) Z(cid:48) S(cid:48) u(cid:48) s(cid:48)
· dx− ·M−1 dS =
F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) S(cid:48)(cid:48) u(cid:48)(cid:48) s(cid:48)(cid:48)
Ω ∂Ω
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)
g(cid:48) s(cid:48)
− ·M−1 dS.
g(cid:48)(cid:48) s(cid:48)(cid:48)
∂Ω
The inverse of M is
(cid:18) (cid:19) (cid:18) (cid:19)
1 −a(cid:48) −a(cid:48)(cid:48) 1 a(cid:48) a(cid:48)(cid:48)
M−1 = = ,
−(a(cid:48))2 −(a(cid:48)(cid:48))2 −a(cid:48)(cid:48) a(cid:48) |a|2 a(cid:48)(cid:48) −a(cid:48)
so if we require that a(cid:48) be negative, the matrix that results from discretizing the left-
hand side will have the same block form as those that result from the other boundary
conditions.
4. Error Bound
We will make the following assumptions on Z(cid:48)(cid:48):
a. there is a constant γ such that |[Z(cid:48)(cid:48)] (x)| < γ for all i, j, and x ∈ Ω
(4) 1 ij 1
b. there is γ such that Z(cid:48)(cid:48)(x) > γ I for all x ∈ Ω.
2 2
AN IMPROVED METHOD FOR THE HELMHOLTZ EQUATION 7
Define the space V = [H1(Ω)]2, endowed with the norm
(cid:107)(u(cid:48),u(cid:48)(cid:48))(cid:107) = ((cid:107)u(cid:48)(cid:107)2 +(cid:107)u(cid:48)(cid:48)(cid:107)2 )1.
V H1(Ω) H1(Ω) 2
Also, we will assume that V and V are finite dimensional subspaces of H1(Ω),
N1 N2
and V = V ×V is the set in which we seek our numerical solution.
N N1 N2
Define an energy f(s(cid:48)) for s(cid:48) ∈ H1(Ω) as
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
1 1 S(cid:48) Z(cid:48)(cid:48) Z(cid:48) S(cid:48)
f(s(cid:48)) = Y(s(cid:48),u(cid:48)(cid:48))+Q(s(cid:48),u(cid:48)(cid:48)) = · dx+Q(s(cid:48),u(cid:48)(cid:48)),
2 2 F(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) F(cid:48)(cid:48)
Ω
where Q(s(cid:48),u(cid:48)(cid:48)) contains terms that arise from the enforcement of boundary condi-
tions and any inhomogeneous terms. We will write this as
1
f(s(cid:48)) = B(s(cid:48),s(cid:48))−F(s(cid:48),u(cid:48)(cid:48)),
2
where
(cid:90)
B(s(cid:48),s(cid:48)(cid:48)) = S(cid:48) ·Z(cid:48)(cid:48)S(cid:48)(cid:48) dx
Ω
and F(s(cid:48),u(cid:48)(cid:48)) contains the rest of the terms. If u(cid:48) is a minimizer of f(s(cid:48)), then
B(u(cid:48),s(cid:48)) = F(s(cid:48),u(cid:48)(cid:48)) ∀ s(cid:48) ∈ H1(Ω).
Therefore, we can write
1 1
f(s(cid:48)) = B(s(cid:48),s(cid:48))−F(s(cid:48),u(cid:48)(cid:48)) = B(u(cid:48),u(cid:48))−F(u(cid:48),u(cid:48)(cid:48))+ B(s(cid:48),s(cid:48))−F(s(cid:48),u(cid:48)(cid:48))
2 2
1
= B(u(cid:48),u(cid:48))−F(u(cid:48),u(cid:48)(cid:48))+ B(s(cid:48),s(cid:48))−B(u(cid:48),s(cid:48))
2
1 1 1
= B(u(cid:48),u(cid:48))−F(u(cid:48),u(cid:48)(cid:48))+ B(u(cid:48),u(cid:48))−B(u(cid:48),s(cid:48))+ B(s(cid:48),s(cid:48))
2 2 2
1 1
= B(u(cid:48),u(cid:48))−F(u(cid:48),u(cid:48)(cid:48))+ B(u(cid:48) −s(cid:48),u(cid:48) −s(cid:48))
2 2
Suppose that u(cid:48) ∈ V is such that
N N1
f(u(cid:48) ) = min f(s(cid:48)).
N
s∈VN1
Then
B(u(cid:48) −u(cid:48) ,u(cid:48) −u(cid:48) )1 = min B(u(cid:48) −s(cid:48),u(cid:48) −s(cid:48))1.
2 2
N N
s(cid:48)∈VN1
The inequalities above imply that
√ (cid:112) √
γ (cid:107)s(cid:48)(cid:107) ≤ B(s(cid:48),s(cid:48)) ≤ C γ (cid:107)s(cid:48)(cid:107) ∀ s(cid:48) ∈ H1(Ω),
2 H1(Ω) 1 H1(Ω)
8 RUSSELL B. RICHINS
and therefore,
√ √
γ (cid:107)u(cid:48) −u(cid:48) (cid:107) ≤ min C γ (cid:107)u(cid:48) −s(cid:48)(cid:107) .
2 N H1(Ω) 1 H1(Ω)
s(cid:48)∈VN1
Here and in what follows, C will represent any constant that does not depend on u(cid:48),
u(cid:48)(cid:48), or the grid spacing h.
LetF betheorthogonalprojectionfromH1(Ω)ontoV . Then(cid:107)F (cid:107) =
1 N1 1 B(H1(Ω),H1(Ω))
1, where B(H1(Ω),H1(Ω)) is the set of all bounded linear functions from H1(Ω) to
itself. We then take s(cid:48) = F u(cid:48) to obtain the inequality
1
(5) (cid:107)u(cid:48) −u(cid:48) (cid:107) ≤ C(cid:107)u(cid:48) −F u(cid:48)(cid:107) .
N H1(Ω) 1 H1(Ω)
If instead we use the energy
(cid:90) (cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
1 1 F(cid:48) Z(cid:48)(cid:48) Z(cid:48) F(cid:48)
f(s(cid:48)(cid:48)) = Y(u(cid:48),s(cid:48)(cid:48))+Q(u(cid:48),s(cid:48)(cid:48)) = · dx+Q(u(cid:48),s(cid:48)(cid:48)),
2 2 S(cid:48)(cid:48) Z(cid:48) −Z(cid:48)(cid:48) S(cid:48)(cid:48)
Ω
and perform calculations similar to those above, we obtain the bound
(6) (cid:107)u(cid:48)(cid:48) −u(cid:48)(cid:48) (cid:107) ≤ C(cid:107)u(cid:48)(cid:48) −F u(cid:48)(cid:48)(cid:107) ,
N H1(Ω) 2 H1(Ω)
where F is the orthogonal projection from H1(Ω) onto V .
2 N2
Combining inequalities (5) and (6), we find that
(cid:107)(u(cid:48) −u(cid:48) ,u(cid:48)(cid:48) −u(cid:48)(cid:48) )(cid:107)2 = (cid:107)u(cid:48) −u(cid:48) (cid:107)2 +(cid:107)u(cid:48)(cid:48) −u(cid:48)(cid:48) (cid:107)2
N N V N H1(Ω) N H1(Ω)
(cid:16) (cid:17)
≤ C (cid:107)u(cid:48) −F u(cid:48)(cid:107)2 +(cid:107)u(cid:48)(cid:48) −F u(cid:48)(cid:48)(cid:107)2 = C(cid:107)(u(cid:48) −F u(cid:48),u(cid:48)(cid:48) −F u(cid:48)(cid:48))(cid:107)2
1 H1(Ω) 2 H1(Ω) 1 2 V
and consequently,
(cid:107)(u(cid:48),u(cid:48)(cid:48))−(u(cid:48) −u(cid:48)(cid:48) )(cid:107) ≤ C(cid:107)(u(cid:48) −F u(cid:48),u(cid:48)(cid:48) −F u(cid:48)(cid:48))(cid:107) .
N N V 1 2 V
We partition Ω into subregions e , each of which can be viewed as a suitably shifted
l
and rotated version of a reference element eˆ, so that there exist affine changes of
variables F (x) = Bx+x such that F (eˆ) = e . In what follows, a hat over a function
l l l l
will denote the corresponding function defined over the reference element eˆobtained
by a change of variables.
We define the seminorm |·| by
s
|u|2 = [u,u] ,
s s
where
(cid:90)
(cid:88)
(7) [u,w] = Dαu·Dαw dx
s
eˆ
|α|=s
and α is a multiindex.
AN IMPROVED METHOD FOR THE HELMHOLTZ EQUATION 9
From [3] we get the inequality
(8) c−1hs−d|w| ≤ |wˆ| ≤ chs−d|w| ,
2 s,e s 2 s,e
l l
where c is a constant, w = wˆ◦F−1, and |·| denotes (7) with e in place of eˆ.
l s,el l
We now recall the following lemma from [1]:
Lemma 1 (Bramble-Hilbert Lemma). For some region Ω ⊂ Rd and some integer
k ≥ −1, let there be given a bounded linear functional
f : Hk+1(Ω) → R,
satisfying |f(u)| ≤ δ(cid:107)u(cid:107) for all u ∈ Hk+1(Ω) for some δ independent of u.
Hk+1(Ω)
¯
Suppose that f(u) = 0 for all u ∈ P (Ω). Then there exists a constant C, dependent
k
only on Ω such that
|f(u)| ≤ Cδ|u| , u ∈ Hk+1(Ω).
k+1
Let s ∈ {0,1} and fix w ∈ Hs(eˆ). Define the functionals
L (uˆ) = [uˆ−F uˆ,w] and L (uˆ) = [uˆ−F uˆ,w] .
1 1 s 2 2 s
Since
|L (uˆ)| ≤ |uˆ−F uˆ| |w| ≤ (|uˆ| +|F uˆ| )|w| ≤ ((cid:107)uˆ(cid:107) +(cid:107)F uˆ(cid:107) )|w|
j j s s s j s s H1(eˆ) j H1(eˆ) s
≤ 2(cid:107)uˆ(cid:107) |w| ≤ 2(cid:107)u(cid:107) |w| ,
H1(eˆ) s Hk+1(eˆ) s
and F u = u for polynomial functions u in V (j = 1,2), we see that the Bramble-
j Nj
Hilbert Lemma applies, and there exist constants such that
|L (uˆ(cid:48))| ≤ C|w| |uˆ(cid:48)| and |L (uˆ(cid:48)(cid:48))| ≤ C|w| |uˆ(cid:48)(cid:48)| ,
1 s k+1 2 s k+1
as long as k is small enough so that all polynomials of degree less than or equal to
k are contained in the span of the basis functions representing uˆ(cid:48) and uˆ(cid:48)(cid:48). Taking
w = uˆ(cid:48) −F uˆ(cid:48) in the first inequality and w = uˆ(cid:48)(cid:48) −F uˆ(cid:48)(cid:48) in the second yields
1 2
|uˆ(cid:48) −F uˆ(cid:48)| ≤ C|uˆ(cid:48)| and |uˆ(cid:48)(cid:48) −F uˆ(cid:48)(cid:48)| ≤ C|uˆ(cid:48)(cid:48)| .
1 s k+1 2 s k+1
Assuming that h ≤ 1 and using inequality (8), we see that
|u(cid:48) −F u(cid:48)| ≤ Chd−s|uˆ(cid:48) −F uˆ(cid:48)| ≤ Chd−s|uˆ(cid:48)| ≤ Chk−s+1|u(cid:48)|
1 s,e 2 1 s 2 k+1 k+1,e
l l
and
|u(cid:48)(cid:48) −F u(cid:48)(cid:48)| ≤ Chd−s|uˆ(cid:48)(cid:48) −F uˆ(cid:48)(cid:48)| ≤ Chd−s|uˆ(cid:48)(cid:48)| ≤ Chk−s+1|u(cid:48)(cid:48)|
2 s,e 2 2 s 2 k+1 k+1,e
l l
Consequently, the overall error satisfies
(cid:107)(u(cid:48),u(cid:48)(cid:48))−(u(cid:48) −u(cid:48)(cid:48) )(cid:107)2 ≤ C(cid:107)(u(cid:48) −F u(cid:48),u(cid:48)(cid:48) −F u(cid:48)(cid:48))(cid:107)2
N N V 1 2 V
10 RUSSELL B. RICHINS
(cid:88)(cid:2) (cid:3)
≤ C |u(cid:48) −F u(cid:48)|2 +|u(cid:48) −F u(cid:48)|2 +|u(cid:48)(cid:48) −F u(cid:48)(cid:48)|2 +|u(cid:48)(cid:48) −F u(cid:48)(cid:48)|2
1 0,e 1 1,e 2 0,e 2 1,e
l l l l
l
(cid:88)(cid:2) (cid:3)
≤ C h2k+2|u(cid:48)|2 +h2k|u(cid:48)|2 +h2k+2|u(cid:48)(cid:48)|2 +h2k|u(cid:48)(cid:48)|2
k+1,e k+1,e k+1,e k+1,e
l l l l
l
≤ Ch2k(|u(cid:48)|2 +|u(cid:48)(cid:48)|2 ).
k+1,Ω k+1,Ω
We have now proved
Theorem 1. Under the assumptions (4) on Z(cid:48)(cid:48), if the solution (u(cid:48),u(cid:48)(cid:48)) ∈ [Hk+1(Ω)]2
and the finite element subspace used in the numerical method contains [P (Ω¯)]2, then
k
there exists a constant C such that the error satisfies
(cid:107)(u(cid:48),u(cid:48)(cid:48))−(u(cid:48) −u(cid:48)(cid:48) )(cid:107)2 ≤ Ch2k(|u(cid:48)|2 +|u(cid:48)(cid:48)|2 ),
N N V k+1,Ω k+1,Ω
where h ≤ 1 is the grid spacing.
4.1. Regularity. For the error bound above to be useful, we must have k ≥ 1,
which means that u(cid:48),u(cid:48)(cid:48) ∈ H2(Ω), at the least. Classical regularity theory, such as
found in [5] tells us that if L is positive definite, bounded, and has C1 regularity,
then u(cid:48),u(cid:48)(cid:48) ∈ H2(Ω).
5. The Numerical Method
We will examine the numerical solution of the model problem
(cid:26)
∇·L∇u = Mu in Ω
u = f(cid:48) +if(cid:48)(cid:48) on Ω
The first step in solving the problem is to select a set of finite element basis functions.
The numerical examples presented here will use a rectangular grid with bilinear basis
functions of the form
(cid:40) (cid:16) (cid:17)(cid:16) (cid:17)
1− |x−xj| 1− |y−yk| if |x−x | ≤ h, |y −y | ≤ h
h h i k ,
0 otherwise
where h is the grid spacing.
Regardless of how the basis is chosen, we will assume that the basis functions are
labeled as {ψ }K1 and we assume that the solution has the form
k k=1
(cid:18) u(cid:48) (cid:19) (cid:18) ψ(cid:48) +(cid:80)α(cid:48)ψ (cid:19)
u(cid:48)(cid:48) = ψ0(cid:48)(cid:48) +(cid:80)αk(cid:48)(cid:48)ψk ,
0 k k