A finite element of Galerkin type semidiscrete method is proposed for numerical solution of a linear hyperbolic partial differential equation. The question of stability is reduced to the stability of a system of ordinary differential equations for which Dahlquist theory applied. Results of separating the part of numerical solution which causes the spurious oscillation near shock-like response of semidiscrete scheme to a step function initial condition are presented. In general all methods produce such oscillatory overshoots on either side of shocks. This overshoot pathology, which displays a behavior similar to Gibb's phenomena of Fourier series, is explained on the basis of dispersion of separated Fourier components which relies on linearized theory to be satisfactory. Expository results represented

Durgun, K.

English

NASA Technical Reports Server

N86-14083
:\
'_ FINITE ELEMENT OR GALERKIN TYPE SBMIDISCRETE SCHE_IES_ri
_: Kanat Durgun
r
, ABSTRACT
'A finite element or Galerkin type semidiscrete method is proposed for numerical
solution of a linear hyperbolic partial differential equation. The question of
stability is reduced to the stability of a system of ordinary differential
_ equations for which Dahlquist theory applies.
_-. We also present some results of separating the part of numerical solution which
- causes the spurious oscillation near shock-like response of semidiscrete scheme
to a step function initial condition. In general all methods produce such
Z
oscillatory overshoots on either side of shocks. This overshoot pathology,
• !
=
which displays a behaviour similar to Gibb's ?henomena of Fourier series, is !
!
explained on the basis of dispersion of separated Fourier components which
I
relies on linearized theory to be satisfactory. We present expository results,
polished formal proofs will appear elsewhere.
i
INTRODUCTION _' 4
Our model of one and two dimensional linear hyperbolic equations are
•(I) _-_+c_--v-_=o
_, (2) _--9-u+ a%V.u+ _-u: o
,_ _" NASA Summer Faculty Fellew
_ Associate Professor of _.lathematics,University of Arkansas at Little Rock.
5-1
1986004609-043
https://ntrs.nasa.gov/search.jsp?R=19860004614 2020-03-20T16:34:55+00:00Z
2 I.
Introducing c= oz_-_&=,C-Cf.QsoC and Cd=cS,_oL equation (2) can be written as
(2') _U_t c_ _aU c_-_ o--;-+G ._ + =
Galerkin or finite element semidiscretization [1], [2] , [3] , [4] ,seeks an
approximate solution for equation (2') in the form
"" where
(4) _)_= {0 0therwise.
4
We obtain a system of ordinary differential equations by requiring that the
-. residual _,,_--_tcC_._._.cS.,_]__,beorthogonal to the basis functions _O i.e
_? <_L_ , R,>=O. Candidates for L_m are too many producing algorithms with
increasing complexity proportional with their smoothness. We only present bilinear
finite elements on squares. The orthogonality requirement y;elds, say in one
dimensional case.
where _¢and L kare discrete Toeplitz operators with eigenvectors _).
If _is an identity operator then scheme is explicit, otherwise implicit. If
the real part of the corresponding eigenvalue _0_) is zero then the scheme is _
conservative [6] , [7] . The quantity
( _) _ (.c_ :. _ ,1._':),__
u)
is the velocity of propagation of numerical solutions in comparison with exact
propagation velocity C in (l). The quotient _'_);c_or difference C_}-C in
an appropriate northis the measure of spurious oscillations and dispersions in
numerical solutions. Purely mathematical treatment without the effects of
discretization i.e nonnumerical can be fecundin [8].
In the next sections, to study the response of semidiscrete scher_eto sharp
¢
5-2
1986004609-044
;;gradientchangeswe simulatea shockby a step functioninitialconditionin (I)
herewe presentan heuristicargumentfor the cause of parasiticoscillations
A
' arounda pointof discontinuity.
_, Considerthe weightedGalerkinsemidiscretization
_ (7) ._. 4.(,- +_ =
i of our model equation(l), where _e[o, 1]is a parameter.Note thatoc-ocorresponds
i :to the equation
(8) _w---_=centereddifferenceapproximationto C-C'_LL_• _L- . _ ,,' '
il Since for any n , in equation(7), indicestake three successiveinteger
4
-_ valueswe may relabelthem for n even as LL_and for n odd as t_ we then
-
_ obtainrespectivelythe followingsystems
(9)
d,,, c u.,2) i.
;'t;for=.o,and ))
•_ zh.. -
,- (lO) _T +  :-
'* (_-_') +'_" + _t; / u.,,,,
; These equationsare consistentapproximationsfor the followingsystems
1
" I
I
., I,. _ "a,(
Eliminatingu or I/In (Illwe obtainrespectively
(13) _'L" "__._
; .--- C7-r_'z'w,
_t_ _
5-3
] 986004609-045
Showing that in a doubly spaced grid wave equation is satisfied. This indicates
that finite differencing is consistent with (13) rather than (1). Also adding
' the equations in (ll}
; we see that discretization is consistent for the avarage of the solutions at
two successive grid points• However _ubtracting equations in (ll) we obtain
This shows that due to discretlzation difference, however small, of two successive
C solutions propagates as an error wave in tilediscrete medium in the opposite
direction.
" For (12), adding we obtain
and subtracting we find
which is the cause of oscillations in general.
t
GALERKIN S_IDISCRETIZATION FOR EQUATION (2') !
On the square with vertices (_,_.,__ n_,), (_r_+l,_1)' (_m+l,%n-i) and
(_-l_ _-_ ) we take basis functions to be {
, I
%-
i 0 otherwise
] 986004609-046
i. These are pyramids whose basc.is a square with vertices are given above, centered
: at (_, u ) with unit height.Forming the inner products with the residual we
obtainh
,4
Only nonvanishing terms come for the values of indices I-_=n_..t_m.,fr_4and L, n-_,_.,_'_t.
" l',husequation (19) reduces to
_.._-_ --_ .'_'_,,_ _%.,.._>:o
Computation of inner products as double integrals are straightforward but
" tedious. Replacing the values of various integrals in equation (20), we obtain
_. (21) _-__ u_p.,-l-_m.,,n-_ . .
where B=c_-'T-_G_h'=C_-c_'_G '3 _ and e_3= :ZC_._.
i: oystem of equations (21) can be written in matrix notation on a rectangle [ O,(M+l)h_
T
_. X[O,(N+l)h in various ways. Let LJl:_,lal6 .... _I] ' _.:_,:L,..N , IN be the_
NxN identity matrix, and LN=[_U_Nx N where
_.i_{I j-L-t0 otherwise,
_4.Q
superscript T indicates transposition. Then
_ _E.t4 (.L,_+ _i _.LI) U, + ( L,_+ 4::i., + L_) U_1 , ,n',
(,c*_-,+_-_.I,  oLI)U,1- k t'_°*'_'+_'-('_L, $ˆ].,,,-eLi)_Jo_-
4-
.\
] 986004609-047
3_ _t
wherewe introducedvectorsin IR,14
4
: " "_T •
:e _%-- - .lu'_O_'"., ' t , _= %_Z,...,
_k= [_a , ,.... aE,_+_ , ...,
. "T
. _=[ o .. o _.. 1 ,_-o,_...,_*_
-_.
We 1et
"4 l o _
_, 4., 0 -_o 0
"i-N=LH+4Z,+L,= , _= ',.
1,I (:_-- N
--(') "'., 0 -% .V
I 4-. _ -,_o
5-6
9 o .. "®
1986004609-048
+\
" ,¢. Thenthe system (22)in vectorand block tridiaganalmatrixnotationbecomes
• w .
36 at _-_ I
• 1
"-" % " 0 .e. i '
. t.,_.., - 4T_I.U. .{_,,% i .+;,. I I1i ,
(23) . .
-e_ "Uo: '.,o+,,+_<. _-r. u<>i
o 0, t + o °: t" -_ i-t- I1 I._ • '
. i,I 3G 44: " i+J, . i
" O " " ] 0 .i"d 0 o o ;
1
1
where -_ i
" = I , _r t "i "L'A,*_,":++.>_, <_-"_" ' _>""_" %_+"_ iiwl
Hereentriesof matricesare NXN matricesand entriesof vectorsare N vectors.
Note thatfor time independentboundaryconditionsthe last term in this equation
vanishes,Furthersimplificationis obtainedby introducingNM dimensionalvectors
or M dimensionalcompoundvectorsi.e wctors whose componentsare N dimensional
vectors,
u= , ,_}= ,_: J,N=. ,_=• k[z :+x]. . , Z:- .44,
and the squarematricesof order NM,_ for the matrixon the left and(_ for
_+the matrixon the right hand side of equation(23),
_ The linearsystem (22)or equivalently(23)can be writtenas
;_, (24) _ :-]=_ "i-
#,#
+ _ 5-7
1986004609-049
._ withthe initialcondition_L:_-_o whent:0.
' Thissystemhasa uniquesolutionforJ.--_. Letting_= A_+_K:,wlth
-TN ckl, .o T.
-r.0 T.o 0
,b
0 o ,".
. -r. T.o.
directmultiplication_howsthat A1 and A2 co,,,,,ute,thisis a directcon-
.., sequenceof bothbeingT_eplltzmatrices,thereforetheyhavethe sameelgenvectors. '
Elgenvaluesof Al, as easilyverified,are
_._.=_Cz_-_ _.L__ , k. _+z,...,_.NM+rl ..
with correspondingeigenvectors
k
_+ -- Ib • " It _ '
': Nhk+X N_rX
"_ J t
T- f
+ Let(_kbeanelgenvalueof A2 associat+d.it, t,e e,genvector_<k'sinceA,: A_ )
.+
we have
!
"T
t
_kL)O , hence_ Is nonslngular[g].It Is knownthatJ)kis unltarlly
similarto a dlagoml matrix_) wlthelgenvaluesof_ , whicharethesumof the
etgenvalues of A1 and A2, are thedlagona_entries.This slmllarltjtransfor-
•"_.._,,_ ] l.e columns of S arematlonis performedby taking_ .... "_N_
-I
. ,5"_elgenvectorsof,_ LettingS I_L=_.)andmultiplying(24) by theInitial )
value problemreducesto +
(2s)
5" IJL.-L1
dl *oNotethatonedoesnotne:_to compute , since S -._
Forthesolutionof (25)one stepmethodssuchas Runge-Kuttamethodcanbe used.
+ Alsoa largenumberof multlstepmethods,implicitor explicitin time(predictor- ',
J 5-8 , '.
1986004609-050
ilit:
corrector),o.ce a starting procedure is realized by a one step method, are
available,and their stability theory is well understood and detailed treatment
ca_ be found in [lO] , [Ill .
._ To show that the finite differencing scheme is conservative, we must show that
the eigenvalues_(_J,_)ofGalerkin difference operators in (21) belonging to
eigenvectors expI.[_c_t_'ru_<l_rL_,are purely imaginary where _-u._c_, and
cz ,_c_ L
_<_=u_S;_cx.Substituting _L _CE_=L(t)cX_._ x. _u_,_n_ in (21) after some
manipulationyields for the left hand side
_ and for the right hand side
_-_ R.H.S=-_-_._¢ " " L _ --
I
"_. ,.:
Hence !
. where
which is imag,nary. Setting u_r(_,_L_,_J_n._L(.%oc)wefind the numerical solution I_
The discrepancy between C(_,d,} ar,d C or more precisely the order of zero of
C_c_,o(_-C about cu,_-.C)is the the order of accuracy of the semidiscrete method.
To show that this method is of order four, K'eexpand _(_o,eL) in a Taylor series
L and a straightforwardcomputation shows that
;_ 5-9
1986004609-051
C ERROR ESTIMATES FOR PRE A_IDPOST 03CILLATIONS ABOUT DISCONTINUITIES
To justify the heuristic argument _reser.tedearlier we may assume that this
spurious oscillations are rapresented by small perturbations in _0 , and rep]ace
_) by u_ Œtrial solutions
(28) _L¢o%£)= 0.__ ,
Expanding ¢'_u_-r_) in a Taylor series about g_:O and retaining only the
linear terms we obtain
L Sincecu>>F.., terms of order L_can be neglected and introducing group velocity
". (29)_c_): _ C__c_)
= (28) can be written as
Straightforward computation shows that eigenvalues of (7) corresponding to
eigenvectors { _'L_A_ , are
%
ar.dtherefore
Using (29),_u_) is easil_ computed as
',_l)_= c _+(_-,_),Co_o_
Due to the discretization of the domain of the equation, the group velocity
corresponding to 2h wavelenght, from equation (31) is
which is the same as depicted in (17) and for _=_ in (15).
To obtain estimates on loca_ and global error of numerical solution we recall
the definitions. [5] p.43, [13_. We say an infinite series __._k_is (C,1)
5-I0
• ®
-- _-_._._T._-"
i986004609-052
11plp_j iI_% . • "t •
:S,. E
su_.anableif lim O'_-_lim£°+S'+-'_'=S existswhere Sk_'KT-z_u_. in this case
C
.i we write ____t_--.S (C,l)sense.
An infiniteseries ZtL_. is said robesummableby Abel'smethod ( some say,.
. Poisson's)or simplyA-summableto s, if _-_tEr_ is convergentfor_r_<:'land
lim_-tt F_ . k
=11,_(l-r)___skr=Swhere Skis definedabove.
We need two results,the firstis that the series
I
is (C,l)and also A-summableto zero.
It is known that, [5] p.20
L and from the t.-igonometricdentity
it followsthat
2_
: Thus _ C
Stn.Ck+ -.r..= .... ac z ,vr
le.--_
and lim(l'_O_.-{oshowA-summability,recallthe Poisson'sformula[5], p.6l;
_r_ -_(xa
- e_
I r_ _- r_
Lettingr._,(we see that the assertionis true.
The secondresultis that
is (C,l)and also A-summableto 4"(aa_
It is knownthat [5], p.21;
! O.A;_.x - _
R.
Using the trigonometricidentity
-z
m 2.
we find
_ = ..... = _co_c_C ....
and the result followsby letting r_.-_
: 5-11
986004609-053
To showA-summabilitywe use Poisson'sformula
F
• n=_ 4 - _LFF'-cL_JC+r Z
and let r-_l.
We now estimatetheJ_znormof the globalerror.As a directconsequenceof
Parseval'sidentity,it is knownthat Fouriertransformis an isometricisomorphism
betweenthe Hilbertspacesinvolved[16] p.51°52,[15] p.25.Thereforeit
sufficesto co:,_putehe &norm of the Fouriertransformof the error.To simulate
the shock,we let the initialconditionto be the step function
: Withoutlossof generalitywe may assumethat the discreteFouriertransformof
E
-; the net initialcondition_L_(o)is equal to the Fouriertransformof (35) , and
: )
i
we obtain
(36) U(.uJo_--_r_C_u,o_-- ,,
.eo 13.=0 t t '
• It followsfrom the proofsof statementsconcerning equations(33)and (34)that I
f
l
serieson the Hght of equation(36) is (C,l)and henceA-sunTnableto
: (37)U(_o)= _(_,o)=
_L 5L_
Fromequation(1),Fouriertransform,of the exact solutionis easilycomputed
^ ', -_.c_l:
uc_o,_)=U(_,o)e "4
LL_(_)is the solutionof semidiscreteequation,for simplicitywe assume Kh
to be the identityoperator,takingthe discreteFouriertransformof the
semidiscreteequationand solvingthe resultingdiffer._ntialequationone
obtains
^ ^ %£_")_:
For conservativeschemes _k(uo3=-_._C'{._), thereforethe_. norm of the global
error is __
I 1
'-, tlEIi -
'Ell _'
_.. 5-12
1986004609-054
• • ' i'_I " _. . . + . . - . .4
• 4
i'_ Introducing dimensionless variables _._tc and _=_)_ , a straightforward
computation shows that
--  Ce+( -c)
_ • r_ _._2__ ay.
.'_ o 2.II
: CONCLUSION
The semidiscrete method proposed here has a reasonable Courant number and a fourth
t
: order accuracy. Results are theoretically conclusive. Computational evidence for
detailed comparison of this method with conventional methods will await our
: numerical experiments.
The measure of oscillations in the numerical solution, in a neighborhood of sharp
+
"S changes is the pointwise error. We were able to show with a lenghty argument,
• although there are some gaps in details of proofs, that maxima of the difference
between the exact and the numerical solutions continually diminish and minima
]
: continually increase in an interval of lenght 4h on each side of the sharp i
_ gradient change. Numerical solution is approximately 0.28h larger in the
I
upstream direction.
ACKNOWLEDGEMENT
+
I wish to express my sincere appreciation to Dr. W. Goodrich for providing me
the encouragement and this gratifying experience, i would also like to acknowledge
the pleasant working environment provided for NASA at Johnson Space Center.
5-13
1986004609-055
!REFERENCES
i. Babuska, I. and Aziz, A. K.: Survey Lectures on the Mathematical Foundatior_
of the Finite Elements Method with Applications to Partial Differential
Equations Academic Press, NY 1972.
2. Mitchell, A. R. and Wait, R.: The Finite Element Method in Partial
Differential Equations, Wiley, 1977.
3. Thomee, V.: Convergence Estimates for Semidiscrete Galerkin Methods for
. Variable Coefficient Initial Value Problems. Lecture notes in Math., 333,
_ Springer, Berlin, 1973, pp. 243-262. *
4. Layton, W. J. : Stable Galerkin Methods for Hyperbolic Systems. SIkM J.
Numerical Analysis, Vol. 20, April 1983, pp. 221-233.
" 5. Zygmund, A.: Trigonometrical Series. Dover, 1955.
6 Hyman, J M.: A Method of Lines Approach to the Numerical Solution of !
Conservation Laws. Advances in Computer Methods for Partial Differential
-_ Equations. IMACS, New Brunswick, NJ, 1979.
/
7. Lax, P. D.: Hyperbolic Systems of Conservation Laws II. Comm. on Pure and
Aoplied Math., Vol. 7, 1959, pp. 159-193.
P.. Morawetz, C. S. : Notes on Time Decay and Scattering for _J_me Hypersolic
- Problems. SI_ Regional Conference Series in Applied Mathematics.
9. Varga, R. S.: Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs,
NJ, 1962.
i0. Dahlquist, G.: Stability and Error Bounds in the Numerical Integration of
Ordinary Differential Equations. Inaugural Dissertation. Appsala, Sweden,
1958.
Ii. Henrici, P.: Error Propagation for Difference Methods. SIAM Series in l
Applied Mathematics. John Wiley, 1963. _ _
12. Vichnevetsky, R. and Bowles, J. B. : Fouries Analysis of Numerical
Approximations of Hyperbolic Equations. SIAM Studies in Applied "
Mathematics. 1982.
13. Hardy, G. H. and Rogoslnski, W. W.: Fourier Series. Cambridge. 1968.
14. Whitham, G. D. : Linear and Nonlinear Waves, Wiley. 1974.
15. Sneddon, I. N.: Fourier Transforms. McGraw Hill. 1951.
16. Richtmyer, R. D.: Difference Methods for Initial Value Problems.
Interscience. 1962.
!
_: 5-14
1986004609-056


Finite element or Galerkin type semidiscrete schemes

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server