NASA Contractor Report 187584

ICASE Report No. 91-47 iiiii

ICASE _-T_

_!-7- _

EXPONENTIALLY STABLE _ROXIMATIONS

OF WEAKLY DAMPED WA_EQUATIONS

H. T. Banks

K. Ito

C. Wang

x,,x,

iiii

Contract No. NAS1-18605

June 1991 ..........

Institute for Computer Applications in _nce and Engineering

NASALangleyResearchCenter

Hampton, Virginia 23665-5225

Operated by the Universities Space Re_h Association

:iiiili_i!

:ii_i121

National Aeronautics and

Space Administration

langley Research Center

Hampton, Virginia 23665-5225

.,,Dr',:.... ,v i' ;_Ttf"w c "i _ ;,'"LY -'A"_T _' WAV_

_ ,_',TY _",!:_ t-i_=:t _ 7i ;clrt (lC.t,_,=) #+1

{.";eL

"> t - 7 "
_r
EXPONENTIALLY STABLE APPROXIMATIONS OF WEAKLY

DAMPED WAVE EQUATIONS 1

H.T. Banks, K. Ito, C. Wang

Center of Applied Mathematical Sciences

Department of Mathematics

University of Southern California

Los Angeles, CA 90089

ABSTRACT

We consider wave equations with damping in the boundary conditions.

Techniques to ascertain the uniform preservation under approximation of exponential stability are presented. Several schemes for which preservation can

be guaranteed are analyzed. Numerical results that demonstrate the lack of

stability under approximation for several popular schemes (including standard

finite difference and finite element schemes) are given.

1This research was supported in part by the Air Force Office of Scientific Research under

grant AFOSR-90-0091, in part by the National Science Foundation under grant NSFDMS

8818530, and in part by the National Aeronautics and Space Administration under NASA

grant NAG-l-Ill6. Parts of this research were carried out while the authors were visitors

at the Institute for Computer Applications in Science and Engineering (ICASE), which is

supported under NASA Contract No. NAS1-18605.

1 Introduction

In this paper, we consider approximation methods for the boundary damped normalized (i.e., wave speed c = 1) wave equation system

(1.1) 02,u(t,x) = A,u(t,x), t > 0, x e f/: R",

(1.2) u(t,x) =0, t > O, x E F0 c ¢9Y_,

(1.3) a,u(t,x) + a,u(t,z) = 0, t > 0, • E rl : on,

where the domain 12 is an open bounded subset of R ", F0 is a relatively closed subset

of the boundary aft of the domain fl, and F1 is the complementary subset of F0 in

the boundary c912. The symbol a, represents the directional derivative operator in

the outward normal direction of the boundary. Since the boundary condition (1.2)

is often referred to as a reflecting boundary condition, F0 is called the reflecting

boundary. Similarly, F1 is called a partially absorbing boundary.

The system (1.1)-(1.3) arises in many important models for distributed parameter control problems. In particular, in the model of a vibrating flexible membrane, the solution u(t, x) represents the transverse displacement of the membrane,

and in models for acoustic pressure fields, the solution u(t, x) represents the fluid

pressure (see, [BKS], [BKSW], [BPS], ILl], for more detailed examples). It can

be shown that for a given initial state u(O,x) = uo(x) and u,(O,x) = Vo(X), the

solution of (1.1)-(1.3) decays exponentially in time and the decay rate is uniform

for all initial states (u0, v0) in a certain function space (see [C], [L]). The stability

of the solutions of (1.1)-(1.3) plays an essential role in several control theoretical

issues ([BK],[BW]). In this paper, we are interested in approximation methods that

uniformly preserve the exponential stability of the solutions of (1.1)-(1.3) for the

approximate solutions.

The equations (1.1)-(1.3) can be written abstractly as a differential equation

d

-_w(t) = Aw(t), t > O,w(t) e ;t,

in an infinite dimensional function space j7. In this context, we consider a linear

control system given by

(1.4)

d

_w(t) = Aw(t) + Bh(t), h(t) e R "_,

where h is a control input and B is a linear operator from R "_ into _/. The most

common approach for the approximation of a control problem involving (1.4) is to
formulate a sequenceof finite dimensionalcontrol systemsof the form

(1.5) _N(t) = a"w_Ct) + BNhCt), t > O, w"(t) e _N,

where the dimension of the space )/u increases toward infinity as N tends to infinity. In general, equation (1.5) is derived from (1.4) using space discretization

techniques such as finite difference, finite elements or spectral methods developed

for the approximation of the solutions of (1.1)-(1.3). A control strategy is then

designed for the finite dimensional control problem involving (1.5). This control

is used as an approximation to the desired control function for the control problem involving (1.4) (for example, see [G], [BK], [BW]). One of the most practical

conditions to assure the well-posedness of the finite dimensional control problems,

as well as the convergence of the approximate controls to the desired control for

the infinite-dimensional system is that the solutions of (1.5) for h -- 0 preserve the

exponential decay of the solutions of (1.1)-(1.3).

From a stability analysis of the solutions of (1.1)-(1.3), it is easy to see that

the energy dissipation in (1.1)-(1.3) comes exclusively from the boundary condition

(1.3). Since no medium damping exists in (1.1), we refer to (1.1)-(1.3) as a weakly

damped wave equation. Although many approximation techniques can provide convergent approximations for the solutions of (1.1)-(1.3) by the solutions of (1.5) with

h _- 0, the nature of the dissipation in (1.1)-(1.3) makes it very difficult to preserve

the uniform exponential decay of the solutions of (1.1)-(1.3). Numerical results

from our investigations reveal that most of the popular discretization techniques

can not maintain a uniform decay rate in the solutions of (1.5) with h - 0 as the

dimension of the approximate system (1.5) increases.

Although the preservation of the stability or the stabilizability of hyperbolic

type control systems is well-known to be a delicate approximation problem, there

exist, to our knowledge, only a few analytical results and these are for the approximation of "hyperbolic" delay equations (see [IK]). In this paper, we present a

general approach for the analysis of uniform preservation of exponential stability of

approximation systems for weakly damped wave equations. This approach is later

used to show that a particular mixed finite element method and polynomial based

Galerkin methods preserve a uniform exponential decay rate in the approximate

solutions as the dimension of the approximate system increases to infinity.

An outline of the paper is as follows. In Section 2, we specify the class of

approximation methods for the weakly damped wave equation considered in this

paper. Then we present a general approach for the analysis of decay rate in the

solutions of these approximate wave equations. In Section 3, the approach developed

in the proceeding section is used to prove that a mixed finite element method for
the 1-dimensional problem preserves uniformly an exponential decay rate in the

approximate solutions. In Section 4, polynomial based Galerkin approximations of

weakly damped wave equations in hypercubical domains are analyzed. These are

also shown to preserve uniformly an exponential decay rate. In Section 5, several

popular discretization methods are investigated numerically; a sharp distinction

between the methods analyzed in the earlier sections and the other popular methods

emerges from our numerical findings. In Section 6, our concluding remarks offer a

perspective of the results presented here.

We finish this introduction by observing that interest in the weakly damped

wave equation is also strongly motivated by the question of exact controllability of

the wave equation via partial boundary control. As it was pointed out by Datko (see

[D]), exact controllability or stabilizability via boundary control may be extremely

sensitive to perturbations. However, we view equations (1.1)-(1.3) as the model

of a physical system, where the control input is a nonhomogeneous term in (1.1)

as it is formulated in ([BFS], [BKSW], [BKS], etc.). Thus, our control input does

not introduce changes in the boundary condition (1.3). More discussion on the

significance of our results for the preservation of boundary stabilizability of the

wave equation is given in Section 6.
2 Estimation of decay rate for approximate solutions of the

wave equation: A general approach

In many analytical studies of control problems, the system (1.1)-(1.3) is taken in

the sense of distributions and one seeks mild or weak solutions. As a second order

equation, the state space for the mild solution (u(t),ut(t)) of (1.1)-(1.3) is taken to

be _ - H_o(fl ) × n2(fl), where the Hilbert space H_.o(fl ) is defined by

gl{o(fl ) = {u(.) • Hl(fl)luCz) = O,z • r0}.

For simplicity of analysis, we assume that F0 has a non empty relative interior; thus

an inner product in H_o (fl) that is equivalent to the usual H 1(fl) inner product can

be defined by

< u,v >1= fo v_(_). Vv(_)d_, u,. H_o(fi).

The infinitesimal generator for the mild solution semigroup of (1.1)-(1.3) is defined

by

P(A)={ (u'v) a,_(_) E_llvEH_'°(fl)'AuEL'(fl)' + _vC_)= 0,_ } • r,

A: D(A) C )4 ,.--, )4, A(u,v) = (v,A=u).

It is easy to verify that A generates a C0-contraction semigroup S(t) in _. It is

shown in [C] and ILl that there exist constants M _> 1,w > 0 such that

IIS(t)llLc_l_ M e-_'', for t > 0.

Alternatively and equivalently, for the initial state (uo, vo) • P(A), the solution

(u(t), v(t)) = S(t)(uo, vo) satisfies the following variational equality,

/

(2.1) d d_ < u(t),f >1 = < v(t),f >1 t

< vCt),9>,,cnl = - < _(t),9 >,- I_ _v(t,_)g(_)d_

,tl 1

for all (f,g) • H_o(fl ) × H_o(n ). Equation (2.1) is often referred to as a weak

formulation of (1.1)-(1.3) and it is a natural formulation in which to consider approximation methods.

Consider two finite dimensional spaces H_ and H_ of functions defined on

O given by

H N = span{¢N}g:,, H( = span{¢_}f:,.
We consider approximations to the solutions of (1.1)-(1.3) in the form

N N

k=l k=l

In most cases, Hff and H_ are chosen to be subspaces of H_o(fl ) and L2Cfl),

respectively.However, thiscondition isnot essentialinthe analysis presented below.

Let us denote the column vectors of the coefficients u N(t), v_ (t) by

u_Ct))

: ,0_Ct) =

v_ (t))

:

v_Ct)

Then system (1.1)-(1.3) is approximated by an ordinary differential equation of the

form

[,0 M01(:., N b_v(t ))= [0. -K g -B g oN(t) "

The matices K N and M N are symmetric and positive definite by construction.

The matrix B _ is symmetric nonnegative definite. It arises from the boundary

integral term in (2.1). In general, mappings qN: H_ _-* H_ are chosen [IK]. The

approximate solution (u_(t),vgct)) is required to satisfy the following variational

system which is analoguous to (2.1):

d

d'-_

d

< v_(t),¢_ >L_C_I= aq N (vN (t))qN (¢:) dz.

In particular, if the mappings qN axe defined as qN(¢ff) __ ¢_,j = 1,...,N, then,

the matrices K N, M N and B N are defined by

K,_ = /o v_,_(_),v_'(_)a_,

M,_ - f._,_(_)_(z)dz,

B'N - fr, acNcx)¢N(x)dz"

In the case of Galerkin type of methods, i.e., finite elements, spline based Galerkin

methods, polynomial based Galerkin methods, H g = H N, and therefore, qN __ I

the identity operator.
We are interested in the exponential decay rate of the approximate solutions

(uN,vlV). First, we need to assume that the appropriate norm for the approximate solutions is used. In general, the convergence properties of the approximation

schemes can be stated as follows: Let i jv be a mapping from _N = H_ x H_ into

_. Then (u N, v N) is said to converge to (u, v), if

v")- (u,v)II ' o,

as N tends to infinity.

Definition 2.1 (Uniform exponentially stable approximation) A given approximation method is said to preserve uniformly the exponential stability of the solutions

of (1.1)-(1.3), if there exist constants M and > o independent of N such that for

any initial state (Uo, Vo) the corresponding aprozimate solutions satisfy

])iN(uNCt),vg(t))])_ <<_Me- 'll(uo, vo)Hx,N = 1,2, ....

It follows from the positivity of the matrices K N and M _ that a norm is

defined on H_ x H_ by

))(uN,vN)IIn_xH_ =< KNftN,_t N >RN + < MN_N,# N >R_,

where fin and VN are vector representations of u N and v N with respect to the bases

of H_ and H N, respectively. We require the following condition to hold.

Condition 2.1 (Norm compatibility) There exist constants cl, c2 independent of N

such that for all (u N, v N) E H_ x H_,

We note that in the case of Galerkin type of methods mentioned previously,

the above condition holds with cl = c_ = 1. Under the above condition, in order

to establish preservation of exponential stability, it is sufficient to show that the

approximate solutions (uN(t), vN(t)) have a uniform exponential decay rate under

the norm II " IIH_×R_.

In the analysis of the stability of solutions of (1.1)-(1.3), a generalized Lyaponov

function Q(t) is used in ([C], [L)). It is defined by

Q(t) = -_([((u(t),v(t))[l_) + 2v(t)(i-Vu(t)) + _,,- 1 u(t)v(t) dx
where l is a vector field defined on fl that represents certain characteristics of the

domain fl and its boundary, _,i(x) = Oz,_(x), and the term E(t) - I[(u(t), v(t)}[_

represents the energy associated with a solution (u(t),v(t)) at time t. The main

steps of the analysis consist in establishing the following results:

(L1) There exists a constant 7"1 > 0 such that for all t >/'1,

Q(t) >_ t E(t).

(L2) There exists a constant T2, such that for all t > T2, (_(t) <_ 0.

(L3) For any given T > 0, there exists a constant M such that for all t < T,

Iq(t)l _< ME(t).

From (L1)-(L3), we can find a constant C such that for all initial conditions and

for all t > T = max{T1, T2} we have

E(t) < 4ME(T) < cE(O).

t - t

Using the semigroup property of S(t), e.g., see [P, p.l161, we conclude that S(t) is

exponentially stable as t tends to infinity.

Since (2.2) is an approximation of (1.1)-(1.3), it is natural to attempt to

follow a similar idea to analyze the stability of the solutions of (2.2). Thus, we

introduce a function QN(t) defined by

QN(t) t KNfz N , MN_N(t),_N(t ) >RN) 5(< (t) aN(t) >,N + <

+ < wN_N(t),fzN(t) >RN

where W N is an N × N matrix. Obviously, W Iv plays a role analoguous to that of

the integral term in Q(t) involving the vector field £(.). Then if we take Eg(t) =

II(u N (t), v N(t))II_×H_, by establishing (L1)-(L3) for some constants 7'1, T2, and M,

independent of N, we can obtain the uniform exponential stability of the approximate solutions following an argument similar to that outlined previously for the

solutions of (1.1)-(1.3). The following condition is required if (L1) and (L3) hold

with constants/'1 and M independent of N.

Condition 2.2 (Boundedness o]W N) There exist constants 81 and fl: independent

of N such that that for all fN,gN E R N,

I "4 wNfN,g N >R N I <- 81 < KNgN,g N >R N q-_2 < MNfN, f N >R r: •
Lemma 2.1 Under Condition 2.17, there exists a constant 7"1 independent of N such

that for all t >_ 7"1,

ON(t) > 1⁄4ENCt).

Moreover, for any given T > O, there exists a constant M independent of N such

that for all t < T,

IQ_'(t)l < MEN(t).

Proof: It is sufficient to take 7"1 = 4 max{f_l,fl2} and M > T/2 + max{fll,fl2}.

[]

As it is in the case of the analysis of the decay rate for solutions of (1.1)-(1.3),

condition (L2) is the most difficult to establish (see [C], [L]). In fact, by using (2.2),

we can compute the derivative of QN as follows

QNCt) 1

5(< KNuN,Ct N >RN + < MNoN, o N >RN)

--t < BNo N,O N >R N + < WNv N,O N >RN

- < + B" N),a "

The above derivation involves use of the expressions

u=N (t) __ oN (t), vN(t) = --(MN)-'CKNftN(t) + BNoN(t)).

obtained from (2.2). By adding and subtracting

we obtain

(2.3) 4"(t)

(< KN_tN,ft _v >R_ + < MN_N,_ N >RN)/2,

_ 1 (< KNft_r(t),_jvCt) >R N + < MN_VCt),O_V(t) >RN)

2

+ < gN_v(t),aN(t) >R,_ -- < KN_tN(t),cNfzN(t) >RN

+ < MNON(t),oN(t) >_,, + < MNOg(t),C_VJv(t) >R"

-t < BNoN(t),ON(t) >R" - < BNON(t),C_raN(t) >m _,

where the matrix C _r is defined by (cN) r = W jv(MN) -x.

The following conditions can be of practical use in establishing (L2).

Condition 2.3 There exists constants ax,a2,a3,#, and 6, with (# + 6) < 1⁄2, each

independent of N such that
(,91) For any f,v 6 R N,

< MNfN,f N >RN + < MN/N,cNf N >nN (a 1 q- (22) < BNfN,f N >R t_

#

+it < MN f N,fN >nN;

(Se} For any fN gN 6 R N

< KNgN, g N >RN -- < KNg N, cNg N >RN -- < BNf N, CNg N ?>RN

<_ (6+#) < KNg N,9 N >aN +as < BNfN ,fN >a N .

lz

We note that a necessary (but not sufficient) condition that ($2) hold (take

(f_¢ -- 0) is

(S3) For anygNER N,

< gNgN,g N >RN -- < KNgN, CNg N >RN_ 6 < KNgN, 9 N >RN .

If a matrix W N can be found such that Condition 2.3 holds, we can show

that (L2) holds immediately.

Lemma 2.2 Under Condition 2.3 there exists a constant T_ independent of N such

that QN(t) < 0 for all t > T2.

Proof: Let # be chosen such that 6+# _<1/2, then by taking T2 = al +a2//z+as/Iz,

using (2.3), we have QN(t) _< 0 for all t _> T2.

[]

In order to obtain a uniform exponential decay rate for the approximate

solution semigroups sN(.), where (uN(t),vN(t)) = sN(t)(uN(O),vN(O)), we also

need the following Lax-stability condition for the approximation scheme.

Condition 2.4 There exist constants C > 1 and w > 0 independent of N such that

IIsNCt)IlzcH,"×H )_ce w'

for all t > O.

The above condition is automatically satisfied if the approximation is convergent. Finally, we state the main result of this section.

9
Theorem 2.1 Suppose that a matrix W N can be found for each N such that Conditions 2.£_-_.4 hold. Then, there exists constants L > 1 and a > 0 such that

for all t > 0 and for all N.

Proof: By Lemma 2.1 and 2.2, we can find constants T and M independent of N

such that for all t > T, the following statements hold.

(a) QN(t) > tEN(t)/4.

(b) 0NCt)< o.

(c) QNCT) < MEN(T).

By Condition 2.4 and (c), we obtain

QN (T) < ME N (T) < C2Me 2_TB y (0).

Combining the above inequality with (a) and (b), we obtain

E N (t) < 4C2ie:WTE N (0).

-t

If we let ,k2 = 4C:Me 2wT, the above inequality can be rewritten as

,k2

liS N (t)(u N (0), v N (0)) [l_,_ ×H_ < Tfl( uN (0), vN (0)) rl_,_ ×,,e,

for all (u N(0), v N(0)) E g_ × H_ and t :> T. As a consequence, for r = 4._ _, we

have

1

IISN( )ll ("ffx" ) < 5'

for all N. Using the semigroup property of sN(.) and Condition 2.4, we obtain the

inequality

IISN(t)IIL(H_×H¢) < L. e -'t

by taking a = In 2/4,_ 2 and L is a constant independent of N.

[]

It is clear now that the difficult part of the analysis for a given approximation

is to estabish Conditions 2.1-2.4. In the subsequence sections, we will specify the

required matrix W N for two different types of approximation methods and verify

that Conditions 2.1-2.4 hold.

10
3 Mixed finite element methods

In the most common implementations of finite element methods, the two components of the mild solutions of (1.1)-(1.3) are approximated by functions uN(t,x)

and vN(t, x) with the same smoothness in spatial variable x, or more precisely, the

approximation spaces H N and H_ are often chosen to be identical. However, both

the analysis for the existence of the solutions of (1.1)-(1.3) and nature of the wave

propagation suggest that u and v have different smoothness in x. The methods

analyzed here use different approximation spaces H_ and H_ which give different smoothness in x to the functions u N and v N. The term "mixed finite element

method" is used here to indicate that two different types of approximation elements

are involved in these schemes.

We first consider the case of a one-dimensional wave equation. Let fl =

(0, 1), so that the equations (1.1)-(1.3) can be written for a typical set of boundary

conditions(i.e.r0 = {0},rl = {1})

0 2 O:

(3.1) _-su(t,x) = _x2U(t,z), t > O,x e (0,1),

(3.2) = o, t > o,

(3.3) _zu(t,1) = -aOu(t,1), t >0.

The following approximation method has been proposed by Ito and Kappel

in [IK]. The basis elements ¢_r and cN for spaces g N and H N, respectively are

defined by

cN (x) = 0, 1-NIx-xkl, otherwise, x• [xk-l,Xk+,]l"l[0,1],

and

cN(x) = 1, 0, xotherwise, E [zk-l,xk.l]Fl[0,1J,

where k = 1, "'" ,N and xk = k/N. It is easy to see that {¢kN }*=1 N forms a basis

for the set of continuous piecewise linear polynomials (i.e. linear splines) on [0, 1]

corresponding to the uniform mesh {k/N}_= o with function value equal to zero

at x = O. The ¢_'s are piecewise constant functions with the same support as

the corresponding CN's. Therefore, H_ c H_(fl) = {¢ E H1(0,1)1¢(0) = 0} and

g N = span{¢g}N= 1 C L2(fl). The matrices MN, K N, and B g of the previous

section are defined by

1

11
BiN = a¢_(1)¢_(1).

Before choosing the matrix W N, we first establish the following useful equalities.

Lemma 3.1 For any fN G R N, the following equalities hold.

(i) Let gN = MN IN, then the component g_ is given by

k=l

(2f_ + f_)/N, N N

gkN = (f:-I + 2f_ + fk+')/ ' k -- 2,...,Y- I,

(fNN + f__l)/N, k = N.

(ii) Let gN = KNf N, then the component g_ is given by

N(2fg_ fN),

gg = N(-f__, + 2f_- f_,),

N(fNN -- fN_I),

k=l

k =2,...,N- 1,

k=N.

(iii) For any fN E R N, we have

1 1 N

< M NfN,fN >RN= _(fN)2 -t- _ )--_ (fN + f:_l)2.

k=2

(iv) For any fN E R N, we have

N

< K NIN,I N >R,,= N(I_) 2+ N F_,(I_- fff-_)_"

k=2

Proof: Equalities (i) and (ii) follow directly from the definitions of the matrices M N

and K _¢, respectively. The following equalities are used frequently in this section,

and in particular, are useful in establishing (iii) and (iv). Let ak, bk be given numbers

for k --- 1,...,N. Then

N-1

(3.4) a2bl + aNbN + __, (ak + ak+l)bk

k=2

n-1

(3.s) - a_b,+ aNbN+ _(ak --_k+,)b_

k=2

N

_k(bk + bk_,),

k=2

N

a_(bk -- bk-1).

k=2

12
From (i), we have

< MN/N, fN >R N

= _ 1 N N N

1 (2f_ + fN)fN + ._(f_ + f_-l)f_ N

1 N-1

+"_ _ ((fff + fff-1) + (fff+l + fff))fff"

k=2

By taking ak = fN + f___,bk = f_, and applying (3.4), we obtain

, __1(fN_ 2 1

< MNfN fN >R N= NWl J +

N

E(fN + f__,)2.

k=2

From (ii), we have

< KNfN, fN >R _ N(2f: - f2N)f: + N(f_r _ f_,__)f_NN

N-I

+u E ((:Z- -/:))s:.

k=2

By taking ak = fN _ fff-1, and bk = fN, and applying (3.5), we obtain

N

< KNfN,f u >Ru= N(f:) 2 + N_-_(f N- f:_l) 2.

k---2

[]

Now we can define a matrix C N (and hence the matrix W N = (CN)TM lv)

as follows. For any vector fN E R N, the vector gN = cNfN is given by

A", gN= f:_,),

2N(f r -

k = 2,...,N- 1,

k=N.

The definition of C N can be seen as a discrete approximation of the term in the

multiplier Q(t) involving the vector field L In fact, in our case, we can take

_(z) = 2x. In fact, for any function fN(x) = _k=1N f_¢_(x) e H_, we approximate 2£(x)df N(x)/dx by

eN(x) _- cN(x,)T . cN fN"

where cg (x)r is a row vector valued function given by cJV(x)r = (¢N(x),..., CNN(X)).

One can verify from the definition of C N, cN and cN, that on each sub-interval

[xk-l,xk}, for k = 3,..., N - 1, we have

flV(Xk) - flV(Xk__) eN(x) = 2 Zk_ 1

Xk -- Xk-2

+ Xk

fN(xk+l)- fN(Xk-l)) ,

_k+l -- 2::k-1

13
which can be regarded as a discrete approximation to the first integral in the definition of Q(t).

Thus, for any v N E H2n with

N

,,-"(=)= :E ff _' (=),

k=l

and eN as defined above, the following equality holds

fo_eN(=)v_C(_,)d= >R,, •

WNvN,f n

The Lax-stability for this approximation scheme is given by the convergence

analysis in [IK], and the norm compatibility condition holds with cl = c_ = 1. We

need to establish Conditions 2.2 and 2.3.

Lernrna 3.2 For any vector fN, gN E R N, the following inequality holds

[ < MNf N, CNg N >nN I <---2 < MNf N, fN >n N +2 < KNg N, gN >it N .

Proof: By definition of the matrices M n and C N, we have

1 (2fN + f_)gn + 2(fg + fNN_,)(gNN --gN_,)

< MNf N, CNg N >RN =

1 N-I

k=2

Using the inequality I_bl_<(_/N + lVb2)/2,weobtain

1

(3.6) I(a, "-Fa_)(bl + b2)l <_ --.._(a_ + a_) + N(b_ + b_).

By (3.6), we obtain

1<: MN f N, cNg N _>R N I < :((f,,)2 +(f:+f(),+ (:_+fL,)=) -E

1 N-x

+_ Z; ((sg,+ :_)' + of; + s__,)_)

k=2

+N ((gf)' + (gf- g_)'+ (g_- gL,):)

N-1

+_' Z_((gf+,- g_)'+ (gf - gf-,)')

k=2

14
< _ (:P)'+ (::+/:_,)'

+2N ((g_)' + _(gf - gN_,),)

k=2

[]

Lemma 3.3 For any fN E R N, the following equalities hold.

(3.7)

(3.8)

< MNf N,CNf N >R N

< KNf N, cNf N >RN

= _ < MNfN, fN >RN

+4(f_+ f#-_)fZ, - (/Z,+ :L,)_,

= < KNfN,fN >R N +N2(fN_ fff_,)2.

Proof: First consider (3.7) and observe, as in the previous proof, that we have

N N (2:_+ .:_).:g + 2(:g+ :;,__)(f;, - :___) < MNf N,CN/N >RN=

i N-1

+-__ k(_ +_) +¢:_+___))(¢:_ +s;l -(_ +___)).

k=2

N-1

The sum _k=_ on the right hand side can be written as

1 N-1

k=2

2 1 N

= (:Z+ :L:)'- -_(::+/_)'--__,(I_+ :_'-,)'.

k=3

Therefore, we obtain

1 i¢N_21

< MNf N,CNf N >RN=--_Wl J N

Thus using Lemma 3.1 (iii), we find

N

_(f: +fL,) N 2 +4 (::_N +f;-,)/;-(/; N N N +:;-_) N 2.

k=2

N N

< MNfN,cNfN >n_= --< MNIN, f N >RN +4(fff + f_-l)f_ -(fff + fff-1)2"

Next considering (3.8), we have

< K NI N 'C NfN >R_v __ g(2f_- :_)f2N + 2g'(f N - fff i)2

N-I

+_v5:k((s:-::__)-(:_,-s:))((s:-s__,) +(:_-s:)).

k=2

15
N-1

In a manner similar to that above, the sum _k=2 can be rewritten as

N-1

N :E[(k-1)(::-::_,)'- k(:_,-:_)'+(:,"-:,"_,)']

k--2

N-I

•- -N(N- 11(/_- f__,)' + N(f:- f:)' + N _ (f:- f:_,)'

k=2

Using the above equality, we have

< KN/N,CNf N >RN= g (f_)' + _(f_- I_-,)' + N'(I_- N_,) _.

k=2

Using Lemma 3.1, we obtain (3.8).

By the definition of the matrix B N we have for any fN gN E R N,

< BNIN, f N >n_= a(fNN) 2, < BNfN, CNg N >RN= 20tN/NN(g N --gN_l).

Combining the results in Lemma 3.3, we can obtain Condition 2.3.

Lemma 3.4 Condition 2.3 holds.

Proof: Consider

< MNfN, fN >RN + < MNf N, cNf N >RN

Also, we have

[]

= -(:_ + :L,)'+ 4:$(f,_ + :L,)

_< 4(fff) 2 =-- 4 < BN f _, fN >R_ •

Ot

< KNgN,g N >RN -- < KNgN, cNg N >R,¢ -- < BNfN,CNg N ?>R,_¢

= -N'(o_ - gN_1)' --2aNf_(g_ - g_-1)

< _(/_)2= _ < BN/N,/N >R" •

[]

We can thus invoke Theorem 2.1 to conclude that the method described

above preserves uniformly exponential stability of the solutions of equations (3.1) -

(3.3).

16
In the case of a two-dimensional wave equation on a rectangular domain

_l = {(x, y) : 0 < x < 1,0 < y < 1}, we may consider the following special case of

(1.1)-(1.3):

0 _ O: 0 2

(3.9) _-_u(t,z, y) - 0x, u(t'x'Y)+ b-_zy_u(t,_,y),(z,y) e n,t > 0,

(3.10) u(t,x,y) = O, (z,y) e F0, t > 0,

(3.11) Onu(t,x,y) = -a_u(t,x,y), (x,y) e Fl,t >0,

where F0 = {(x,y) E _,x = O, or y = 0} and F, = 13/F0. The mixed finite element

method is defined as follows. For each given integer N, let (xi, Yi) be the grid points

defined by x, = i/g,y i = j/N, for i,j = 0,1,...,g. For each pair of indices (i,j)

with 1 < i,j < N, a neighborhood of the grid point (xi,yj) is defined as

D,_= {(_,y)• _: Ix- _,J< 1 i ]_ _,+ y;I< 1 _ _,l,-y;I < _, -y- __).

Two families of functions are defined by

cg(x' Y) = 0, 1-Nmax{lx-xil, ly-yi],lx-y-xi-yj[}, (x,y) otherwise, eDi N,

Y I _'i, (x,y) • D,_, ¢ii (x,y) = , 0, otherwise,

where the constants "Y_i are chosen such that we have the equality

fo =fo

We then choose H N = span(¢,y},.j=,, N _ H g = span{_b0 N }i,j=l N We prefer,

for simplicity of presentation, to denote an element of an N 2 × N 2 matrix A _ by

A(ii),(k 0N . Using this notation, the matrices M N, K N, and B N are defined by

Mi'i)(kO N = fn ¢_in(x'Y)¢g(x'y)dxdy'

Kc,_).c_, ) : v,,;(_,_), v,_c_._)e_ey,

B(ii)'(_0 = fr .¢_(_,u)¢_Cx,y)ar,.

Although our numerical results indicate that this method preserves uniformly exponential stability of the solutions of (3.9)-(3.11), we are, to date, unable to confirm it

analytically. We have attempted to use arguments analogous to the one-dimensional

case, but the number of terms to compute and the associated tedium quickly become

overwhelming.

17
4 Polynomial based Galerkin approximation

In this section, we present a stability analysis for a general class of Galerkin type

approximations of the wave equation in R n. Several technical assumptions on the

domain 12 and the approximation spaces H_, H_ are needed for exponential stability of the approximate solutions. In particular, hypercubic domains in R" with all

reflecting sides sharing at least one common corner satisfy these conditions. The requirements on the approximation spaces can be satisfied if we take H_ = H_ = H N,

and H N is a carefully defined subspace spanned by the tensor products of polynomial functions.

Since the arguments developed here can also be useful for the analysis of

other approximation methods, the assumptions are made in a general form. We

will indicate in each case how polynomial based Galerkin methods satisfy these

assumptions.

Now consider the wave equation given by (1.1)-(1.3) where we make the

following assumptions on the domain 12.

Assumption 4.1 There exists a constant r > 0 such that for any e > O, there

exists a vector field _(.) E C'(fi; R"), with the properties:

(Yl) The matrix L(x) defined by

Lii(x) lf_.(___x) _gti(x)/= 1

= 2' ax; + _ ' _(e,,;(_)+,;,,(_)]

satisfies L(x) - I > 0 on the domain 12;

(V2) For x E Of 2, the outward normal unit vector rl(x) satisfies

tCz). _(x) < 0, /or z e r0, t(_).,(_) _>r, /or _ e r,;

[V3} The following inequalities hold

Tt I'1

12: t,.,jl_<,, 12: t,.,_;l<_,.

i,i=1 1,i=1

where

a 2 a s

18
We note that (V1)-(V3) are exactly the same conditions as these required

in [C] where the author treats the exponential decay of the solution of the wave

equation. For a given domain fl, if there exists a point x0 E R" such that (x -

x0). r/(x) < 0 on F0 and (x- x0)-r/(x) >_ r, on rl, by taking l(x) = x- x0, the

conditions (V1)-(V3) are clearly satisfied. In particular, consider the hypercube I2

given by fl = {x E R",0 < x0 < 1}; if F0 is the union of a collection of sides that

share at least one common point x0, then one can verify that t(x) = x - z0 satisfies

(V1)-(V3). In fact, in the case of a hypercube, if F0 is the union of a collection of

sides without any common point, it can be shown that the solutions of the partial

differential equation (1.1)-(1.3) are not exponentially stable.

Now consider subspaces H g = H_ = H N C H2(fl)A H_o(fl ) given by

H N = span{¢N}g=l, where the ¢_v are a general family of approximation elements.

The matrices M N, K N, B N of Section 2 are defined by

M_ = /,_,_,"(x)_(x)dx,

Ki_ = fn VCCx). v_cx)_x

B,_ = fr _C(xl_f(x)ex.

Before we introduce the multiplier matrix C _', we present a few useful special

properties of the matrices M N, K N, B N.

Lemma 4.1 Consider a sequence {0N(-)}N=I of functions in H N and let V N be the

matriz defined by

Then, the matrix D N - (vN)I"(MN)-aBN can also be given by

D,'_= fr, "C(x)_7(x)dx"

Proof: Since the 8g(.)'s are elements of H N, there exist constants 0_ such that

N

C(x) = _ _TCx)O_.

.i=1

Let 0 N be the matrix with entries 0;7; then we have V N = MNo N. In fact,

v N = ¢_(x)¢N(xldx)0_ = __, MhNo,j.N

= l=l

19
Therefore, we obtain D _' = (oN)rB N and hence

Di_ = N N n fr

l=l 1 = l

_O,"(_)_'(x)d_.

[]

Lemma 4.2 Let V N be the same matrix as defined above. Then, the matrix F n =

(VN)r(Mn)-IK N is also given by

Jrf 0in (z)___¢i° N

n = (x)dx - [ CCz)ACfCz)dzJIl

Proof: We only need to consider F N = (oN)TK N where the matrix O N is the same

as before. We have

N

/=1

x l=l l=l

Therefore, we obtain

[]

Now, by imitating the multiplier used in Section 2 for Q(t), we define the

matrix W N = (cN)TM N as

/1

w,_=/°{2(:. vC)_[(x)+(E _,_,- _)C(_)_[(_)}d_ k=l

where £N(.) is a vector field which satisfies Assumption 4.1. In order to be able to

use the previous matrix equalities, we need further technical assumptions on /_N.

Assumption 4.2 There exist constants r > 0 and M > 0 independent of N such

that there exists a vector field lN(-) 6 C4(fl; R") which satisfies Assumption 4.1 and

furthermore

2¢_• vC +_ d_i,i¢ N Nk 6 H N, with II:(_)IIR° <-- M, i_eNCx)l 11/It <-- M,

i=1 i=1

for all x C (2 and k = 1,...,N.

2O
The aboveassiamptionis satisfied if there exists x0 E R n such that _N (x) =

x - Xo satisfies Assumption 4.1 and if H N is defined by

H" = f(_,,...,_n) ; p_(xk)e H_0(a), where

pk(.) are polynomials of degree less than or equal to m_.

In the case of a hypercube, one can easily construct H N by taking p_ E P_k (0, 1),

where P,,_ (0, 1) are polynomials of degree less than equal to mk which satisfy the

appropriate Dirichlet boundary condition.

The following lemma will establish Condition 2.2.

Lemma 4.3 Under Assumption _.1, there exists constants cl, c2 depending only on

the domain t2 such that

J< MU f _,cNg N >R_ [

for all fN, gN E R N.

= 1< WNyN,gN >R" I

< c, < KNgN,g N >RN +C2 < MNfN,f N >RN,

Proof: Let us associate each vector fN E R N with an element fN(.) E H N given by

fN(x ) = _N=, fff CN(x). It iS easy to verify that for any vectors fg gg E R y

< KNgN,g N >RN

< MNIN, fN >RN

Similarly, we have

= f. IlVg"Cx)ll_.d_= 119_(-)lll,

= II/N(-)II_:(n)

< MNfN,CNgN >R_=< WNfN,g _r >RN

= Sla[2fN(x)(,N(x).Vg/v(x))-J - (i__i<'{x)- 1)fN(x)gN(x)d_:.

From standard results, there exists a constant e depending only on fl such that for

all u E H_o(f2), ]lull_(,) < cllu]l _. By Assumption 4.2, we obtain

I< MNf_,C_9 N>_ I < _-II/N(")

+ £_(.) - 1 + nIIF'(')IIL-I.I 119"(')ll_

L_[n)

-- c2 < MN fN, f N >tlN q-C, < KNgN,g N >R_¢,

wherec_= a/2,cl = ,,lle_(.)ll,.-C,-,) + _/211E,"_-, e_(.) - lllL-c,).

The following two lemmas are useful in establishing Condition 2.3.

[]

21
Lemma 4.4 Let Assumption 4.2 hold. For any fN E R N, let IN(.) be the Junction

in n N given by/N(x) = 2_=, fff¢_(x). Then

< MNf N,CNf N >RN = < wNfN,fN >R _

= _ < MNfN, fN >n" + Jr[, f" C,.)'e"(_).o(_)d_.

Proof." Let G N = W N + MN; then

k=l

Then, we have

< GNfN, fN >R N

= /n,__l(2(tN'V¢_f_)¢N(z)f_

- fn{2(eN.Vf_(x))fN(z)

= fodiv(fN(x)=£N(z)) dz

= Jr[ fN(_)'e"(x).,7(_)d_,

-J-_'_£_k (X)fJV(X)'} dx

k=l

dx

because fN(x) = 0 on I'0. []

Let us define the matrix L_V(z) by (recall condition (V1) of Assumption 4.1)

oq(..))

L_(_)= _ _, o_; + _ ]"

Lemma 4.5

H N given by oN(x) = _kN=lgkCN(x). Then we have

< KNgN,CJVg _v >RN -- < KNoN,g N >R_

-- 2 fn(VgN)T(LN(x) -- I)VgN(x)dx-- _l f. (Z '_ £i,,j,)g N N (x)'dx

i,3'= I

_Ffr,_ (I_ lNy_Tj)gN(x)_d x + fr,IVg_l:e_ •,Td_.

i,j = 1

Under Assumption 4.2, for any gN E R N, let gN (.) be the Junction in

22
Proof: We have < KNgN,CNg N >RN=< WN(MN)-IKNgN,p N >RN. By taking

0_(,::r) -- 2l N- V_(x) -1- t_A:(X ) -- 1 ¢_(X)

we have

where V N is as defined in Lemmas 4.1 and 4.2. By Lemma 4.2, we have

< KNgN,CNg N >aN=< FNgN,g N _>RN,

where the matrix F N is given by

As a consequence, we obtain

(4.1) •_ KNg N, cNg N >R N

n

= - foc2t" •vg_(.) + (_ C,_- 1)9_(_))Ag_(_)e_

k=l

l k=l

Consider first the term

k=l

For all functions ¢ E H2(f_)

= div 2V¢(t N • V¢) -IV¢l't N + ¢(_t_)V¢

i=1

_2(v¢)rLNV¢ _ iv 0¢

- t,,,j_-_C i,j= 1

23
Therefore, we obtain

: =

+i., [(.<.. +

k=, 0." o'}r/

- ,rio IvgNI=tN "ndz,

LNVgN(x) + ,,i=1 klN _"# -NOgN_ + gN AgN(z) } dz

- IvgNI2eN •.] dx

where we have used Assumption 4.2 to observe that 2(_ N- Vg N) q- _=1 lk,k9 N is in

H N and hence vanishes on Fo. Since gN vanishes on ro, we have that also 2(gN'Vg N)

must vanish on F0. Moreover, since gN(x) = 0 on to, VgN(x) = +]VgN[_7 on F0

and hence [vgN[2_ N "_7 = q-[vgNI _N" _gN = 0 on r 0. As a consequence, we obtain

_ - _i,ii-_xig + gN/hgN dx

i,.i= l

k=l OT/

Integrating the terms

k _N i)g N N

id=l ;'O'-ff_-z/g and gNAgN(z)

by parts once, we obtain

I N = - ,_{2 < :v:, v: >_o-Iv:l' - (._..e,,,;;) }dx

id=l

, • t_.;,i) TM _ dz. - IvgNI2eN " + (..= 2

Finally combining the above result with (4.1), we obtain

< KNg N, cNg N _>RN

, - _e_,,ij(g_) = dz

id=l

+ Iv¢_l=l" .,7+ _,,,:j)(g_)= d...

1 =1

24
This least equality is equivalent to

< KNgN,CNg N >aN -- < KNgN,9 N >aN

= fn2 (< LNVg N, Vg _v >,. -IVgN] ') dx

--2 (_ Nl_,Oi)(glv ,)dx+ [VgN[_£ N.rldx

i]=1 1

+.2if ° t,,,.n,)(9"

I j=l

[]

By using Lemma 4.1 and the same notation as for Lemma 4.4, 4.5, we find

< B_cfN' cNgN _'RN -_- fr Ot (2( tN " Vg Iv + (_ ntNk,k -- 1)g N )f Ndz'

1 k:l

< BNfN'IN >RN = fr a(fN)2dx"

1

Lemma 4.6 Under Assumption 4.2, Condition 2.8 holds.

Proof: (i) By the boundedness of t_y and Lemma 4.4, (S1) holds. (ii) By taking e

sufficiently small in Assumption 4.1 and 4.2, we have

IVgN(x)I2tN(X) _(X) + (, =1

and

- 9 (x)'<_6 /oIvgN(x)12dz

2 ij=l

for 5 arbitarily small. Therefore ($3) following Condition 2.3 holds.

Schwartz inequality, < BNf N, cN9 N >RN can be bounded by

(iii) By the

I< BNIN, CNg N >R,, ] < Cl fr l]N(z)l'dz+c_U fr ]VgN(z)l 2dr

12 1 1

where the constants cl, c_, and u are independent of N and u can be arbitrarily

chosen. Therefore, ($2) of Condition 2.3 holds.

As a consequence of Lemma 4.3 and 4.6, on the domain fl and the subspace H N, we conclude that under Assumption 4.2 these approximation methods

uniformly preserve the exponential stability of the weakly damped wave equation.

25
5 Numerical studies of approximation methods

As has been demonstrated in the previous three sections, establishing a uniform

decay rate estimate for a given approximation method can be complex and tedious.

In fact the general approach outlined in Section 2 offers only a broad direction

toward the selection of the appropriate multiplier QN. On the other hand, for any

given approximation method, by computing the eigenvalues of the matrices A N in

(1.5), one can observe the trends in the location of the eigenvalues as N increasess.

In several cases presented below, the numerical results clearly indicate that some of

the eigenvalues of A N are approaching the imaginary axis as N increases, and it is

therefore unlikely that a uniform decay rate can be preserved.

The example used for demonstration in this section is the two-dimensional

wave equation (3.9)-(3.11) where the parameter a is taken to be 1. The construction

of the matrix A N is presented for each approximation method. The eigenvalues are

then computed using the subroutine F02AFF of the FORTRAN software library

NAG. The results were compared with the results using the IMSL library, and no

visible difference can be observed. All computations were carried out on a IBM

3081 computer.

5.1 Polynomial based Galerkin approximation methods

Let {Qk('))k=0g be the Legendre polynomials of degree less than or equal to N on

the interval [-1,1]. We define a sequence of polynomials Pk by

= +

fork= 1,...,N.

It is easy to see that pg __ span{pk}N=l is the set of polynomials of degree

less than or equal to N vanishing at s = -1. Let us define matrices M N, K S as

follows

The matrices M, _, K, N can be easily computed using the recursive properties and

the orthogonality of the Legendre polynomials. Now let ¢_y(x,y)g = d2_iN(x,y) =

P,(x)Ps(y ). We assume a scaling of the domain fl onto [-1,1] × [-1,1] is made.

26
Figure 5.1: Locations of the eigenvalues of the matrix A N for the polynomial based

Galerkin method.

N=3

k

i

t

!

I

N=4

N=6

°_

o.

N=9

Q•

.-

! !

i _

!

27
The matrices M N, K N, and B N defined in the previous section can be computed as

follows:

[M"v](iS),(kt)

[K'']Cii),(ki)

[BN] C_i),C_)

= [MN],,,[M;N],:,,,

= rMN1 rKN1 [K_]ikIM_]j' +t 0 j,kt ° Jj',

= a[MN],,.Pi(1)P,(1) + aPi(1)Pk(1)[M_]it.

The matrix AN is given by

_( MN)-0I KN _( MN)-"-

'

B .'v )•

The dimension of A N is N 2. For N = 3,4,6,,9, the locations of eigenvalues are

displayed in Figure 5.1. We note that as would be predicted from the analysis of

Section 4, a uniform margin between the eigenvalues of A N and the imaginary axis

is maintained for all N. In the Table 5.1, we list the margin of stability for each N.

Table 5.1:

axis.

Margin between the eigenvalues of the matrix A N

N

4

5

6

7

8

9

max{Re,k, A e a(AN)}

-0.6230

-0.6215

-0.6232

-0.6097

-0.5815

-0.5628

and the imaginary

5.2 Polynomial spline based Galerkin approximation methods

For N a given integer, the interval [0, 1] is divided into equal sized subintervals

[xl,-1,xk] with xl, = k/N for k = 1,...,N. Let B _r''' be the set of polynomial

splines of order m corresponding to the grid {xk} that vanish at x = 0. Consider

a basis "t-'-"kr'N''"_N+"-llk=lfor S N'''' on fl = [0,1] × [0, 1]. The matrices M N.,, K No are

defined as

[MN],,:= f01 B,N.,,,(x)B i.',,.,,, (x)dx, i,j=l,...,N+m-1,

28
Figure 5.2: Locations of the eigenvalues of the matrix A N for the linear spline based

Galerkin method.

N=4

I

I

I

Ii t

"I'ti

N=5

• ?

° I _ 1

N--7

w

N=lO

aim

o °

I

i

29
Figure 5.3: Locations of the eigenvalues of the matrix A N for the cubic spline based

Galerkin method.

N=3

• • 1

• e °

• T

N=5

•o

e• •

• e e_e

IP

I

N=6

ee o

N=9

••o

°e•

o•

30
[K_],i = Io'2 , B '_(x B '"_(x)dx, i,j = 1,...,N + rn- 1.

The matrices M, N_ ,K, N can be easily evaluated. We define ¢,N(x,y) = ¢,;(_,u) N =

B_'"_(x)BT"n(y). The the matrices MN,K N and S N are defined as

[MN](iD'(k') = fn CN(x'Y)¢_(x'y)dxdy'

[KN](il)'(k') = fn V¢i_(x'Y)" V¢_(x,y)dxdy

[BN]('i)(k') = fr, a¢,j(x,y)¢g(z,y)dr,.

Again, the matrices M N, K N, and B N can be computed using M N and K W. Since

only the first derivative of the spline function is required, either linear spline or

cubic spline functions can be used. In Figures 5.2 and 5.3, the locations of the

eigenvalues of the matrix A N using linear and cubic spline functions, respectively,

are depicted. We note the similarity between the locations of the eigenvalues with

small imaginary part and the locations of the eigenvalues of the matrix A N using

polynomials depicted in Figure 5.1. However, the eigenvalues with large imaginary

part tend toward the imaginary axis. Another interesting observation is that unlike the eigenvalues of A N using polynomials and cubic splines, the eigenvalues of

A N using linear splines do not have exceedingly large negative real parts; further

discussion concerning this observation will be given in the final section.

5.3 Finite element method

The classical finite element method requires that we subdivide the domain 12 into

triangles where functions linear on each triangle are used for the approximation.

More precisely, for any 1 < i < N and 1 _< j < N, a function ¢_ is defined as

¢;iN { 1-Nmax{[x-x,[,[y-Yil, lx-y-x,+Y1[}, (x,y) eD, N, (x, y) = 0, otherwise,

where the support D_ is defined by

D,_ = {(x,y) e (1: ]x- xi] < 1 1 1 _ _,lu-u_l _<_,lx-u-x, + u;I_<7)

with xi -_ i/N, yi = j/N. We choose the functions ¢,i(x,y) N = CN(x,y) and the

matrices M N, K N, and B N are defined as in the other Galerkin methods.

31
Figure 5.4: Loaction of the eigenvalues of the matrix A N for the finite element

method.

N-4

:...:

N-6

!

!

N--7

I °

N--lO

!

.g !

e •

• |

I

P

P

32
The computation of the matrices M N, K n, and B n uses finite element techniques based on elementary matrices.

The location of the eigenvalues are depicted in Figure 5.4. We note that

some eigenvalues with large imaginary part tend toward the imaginary axis which

suggest that no uniform decay rate of solutions is preserved. We also observe that

only a few eigenvalues have large negative real part.

5.4 Mixed finite element method.

The mixed finite element method presented in Section 3 was also implemented. In

Figure 5.5, the location of eigenvalues are depicted.

Figure 5.5: Location of the eigenvalues of the matrix A N for the mixed finite element

method.

N=5

"# -,

• Q

N=6

• I e • 1

| e e

I .t

N=8

• • • ee

• II •

$

. '6, j

e

'l

'I!

Although a uniform margin

N=10

I Q

• • e I' i

! '"

_p

• !

|

• I •

between the eigenvalues of the matrix A n and

33
the imaginary axis is maintained, the locations of the eigenvalues are dramatically

different from the previous approximation schemes at lower values of N. This might

suggest a slower convergence rate for the approximation scheme.

In Table 5.2, the size of the margin between the eigenvalues of A N and the

imaginary axis is reported.

Table 5.2: Margin between the eigenvalues of the matrix A N and the imaginary axis

for the mixed finite-element method.

N max{ReA, A E a(AN)}

4 -0.4689

5 -0.4596

6 -0.4552

7 -0.4527

8 -0.4512

9 -0.4969

10 -0.4496

15 -0.4632

5.5 Finite difference approximation method

Unlike the other approximation methods presented so far, in the finite difference

method one attempts to approximate the solution of the wave equation at the grid

points directly. In fact, for any given integer N, let U,i(t),v,y(t N N ) be the approximation of the solution values u(t,z,,yi) , v(t,x, yi) -- ut(t,x,,yj) of (3.9)-(3.11). The

wave equation is approximated by

d N

(5.1) = v Ct),

(5.2) d n N 2

fori=2,...,N-1 andj=2,...,N-1.

At the boundary x = 0 or y = 0, by taking u0Ni = u m = 0, the equation (5.2) i,0

still holds. However, at the boundary x = 1 or y = 1, the second derivatives are

34
Figure 5.6: Locations of the eigenvalues of the matrix A N for the finite-difference

method.

N=4

I,q

!

!

I

N=7

!

!

!

,...j'

|

!

N=8

.e

i

N=10

F

,....(

-i!

35
approximated by one sided finite differences. In particular, we take

_,,,, = (t)+ .__,N(t)+

+N (avi_,N(t) - NCu_N(t) - u _v,,N__Ct)))

dN = N _

-'_v_v,i (-2u_,iCt) + u_,j-1(t) + u_v,j+a(t))

+N (,_v'L(t)-NC,,L(t)- ,,__,,_(t)))

for i = 1,...,N - 1,j = 1,...,N- 1. At the corner (1,1), we take

d N

-_V_C,N = N (av_,N(t) - NCu_c.N(t) - u_,N_l(t)))

+N (avN,g(t) -- NCu_,NCt ) -- u_c_l,N(t))) .

The locations of the eigenvalues of the matrix A _¢ are depicted in Figure 5.6.

It was a surprise to us that the eigenvalues of A N for this finite difference

approximation are so different from those of the other schemes. Many variations of

the finite difference methods exist; we believe that an adjustment of the method used

here may well produce different locations for eigenvalues. However, our experience

with the above standard finite difference approximation method suggests that no

uniform margin between the eigenvalues of A N and imaginary axis is maintained

for all N.

36
6 Concluding remarks

Even in view of the analysis for the methods given in Sections 3 and 4, one might

ask why these two methods succeed while several other popular approximation

schemes fail. One possible explanation is the following. As approximation methods

of a second order system, many approximation methods approximate the second

component v which has the same smoothness as Vu theoretically by functions with

equal order of smoothness as the approximation for u. In another words, if u0N 6 H N

and v0N 6 H_ are taken to be the initial values of the equations (1.1)-(1.3), the

solution v is much less smooth than elements in H_. This is not the case in the

mixed finite element method and the polynomial based Galerkin methods.

We should also comment on the implications of the present results with respect to the boundary control problem. If the boundary condition (1.3) is considered

as a boundary feedback control, with an appropriate approximation of the input

operator, the results reported here can be used to construct finite dimensional approximate control systems that are uniformly exponentially stabilizable. However,

this does not guarantee that the infinite dimensional system under the approximate control is exponentially stable. The inherent sensitivity of the conservative

wave equation with respect to perturbations of the boundary conditions remains

one major drawback of this type of model.

Finally, we observe that in both polynomial based Galerkin methods and

the mixed finite element method, a large number of extraneous eigenvalues with

large negative real parts are introduced. Therefore, the resulting finite dimensional

ordinary differential equation generated by these methods can be very stiff. This can

cause problems in integration algorithms. Since our main concern is to solve linear

quadratic control problems using these approximation schemes, this observation

should not substantially diminish the attractiveness of these two methods.

37
References

[BFS]

[BK]

[BKS]

[BKSW]

[BPS]

[BW]

[c}

ID}

[C}

[IK]

[L]

ILl}

[P]

H.T. Banks, W. Fang and R.J. Silcox, Active control of noise: A time

domain approach, in Proc. Conf. on Recent Advances in Active Control of

Sound and Vibration (VPISU, Blacksburg) April, 1991, pp. $43-$47.

H.T. Banks and K. Kunisch, The linear regulatorproblem for parabolic

systems, SIAM J. Control and Opt. 22(1984), pp.684-698.

H.T. Banks, S.L. Keeling and R.J. Silcox,Optimal controltechniques for

active noise suppression, Proc. 27th IEEE Conf. on Dec. and Control,

Austin, TX, 1988, pp. 2006-2011.

H.T. Banks, S.L. Keeling, R.J. Silcoxand C. Wang, Linear quadratic tracking problem in Hilbert space: Application to optimal active noise suppression, Proc. 5th IFAC Syrup. on Control of DPS (A. EI-Jai, M. Amouroux,

eds.), Perpignan, France, June, 1989, pp. 1%22.

H.T. Banks, G. Propst and R.J. Silcox, A comparison of time domain

boundary conditions for acoustic waves in wave guides, ICASE Rep No.

91-27,1991, NASA Langly Res. Ctr., Hampton, VA.

H.T. Banks and C. Wang, Optimal feedback control of infinite-dimensional

parabolic evolution systems: approximation techniques, SIAM J. Control

and Optim., 27(1989), pp. 1182-1219.

G. Chen, Energy decay estimates and exact boundary value controllability

for wave equation in a bounded domain, J. Math. Pures Appl., 5(1979),

pp. 249-274.

R. Datko, Not all feedback stabilized hyperbolic systems are robust with

respect to small time delays in their feedback, SIAM J. Control and Opt.

26 (1988), pp.69%713.

J.S. Gibson, The Riecati integral equations for optimal control problems on

Hilbert spaces_ SIAM J. Control Optim., 17(1979), pp. 537-565.

K. Ito and F. Kappel, On variational formulations of the Trotter-Kato

theorem, CAMS Rep. _91-7, Univ. of Southern California, April, 1991.

J. Lagnese, Decay of solutions of wave equations in a bounded region with

boundary dissipation, J. Diff. Eqns., 50(1983), pp. 163-182.

J.L. Lions, Exact controllability, stabilization and perturbation for distributed systems, SIAM Review, 30(1988), pp. 1-68.

A. Pazy, Semigroups of Linear Operator and Applications to Partial Differential Equations, Springer, 1983.

38

4

Report Documentation Page

1 Report NO 2 Government Accession No 3 RecipienT's Catalog No

NASA CR-187584

ICASE Report No. 91-47

t

4. Title and Subtme

EXPONENTIALLY STABLE APPROXIMATIONS OF WEAKLY

DAMPED WAVE EQUATIONS

7 Authorlst

H. T. Banks

K. Ito

C. Wang

9 Performing Organization Name and Address

Institute for Computer Applications in Science

and Engineering

Mail Stop 132C, NASA Langley Research Center

Hampton, VA 23665-5225

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Langley Research Center

Hampton, VA 23665-5225

5 Report Date

June 1991

6 Performing Organizatton Code

8 Performing Organization Report No

91-47

10 Work Unit No

505-90-52-01

11 Contract or Grant No

NASI-18605

13. Type of Report and Period Covered

Contractor Report

14 Sponsoring Agency Code

15. Supplementaw Notes

Langley Technical Monitor:

Michael F. Card

Final Report

To appear in Proc. Intl Conference

on Control of DPS, Voran (July

1990), Birkhauser

16. Abstract

We consider wave equations with damping in the boundary conditions. Techniques

to ascertain the uniform preservation under approximation of exponential stability

are presented. Several schemes for which preservation can be guaranteed are analyzed. Numerical results that demonstrate the lack of stability under approximation

for several popular schemes (including standard finite difference and finite element

schemes) are given.

17. Key Words (Sug_sted by Author(s))

wave equations, approximations,

exponentially stable

19. SecuriW Cla_if. (of this report)

Unclassified

18. Distribution Statement

64 - Numerical Analysis

66 - Systems Analysis

Unclassified - Unlimited

20. SecuriW Cla_if. Iof this pa_) 21. No. of _s

Unclassified 40

22. Price

A03

NASA FORM IiI_ OCT _6

NASA-Langley, I991