DISCONTINUOUS GALERKIN METHODS FOR
ORDINARY DIFFERENTIAL EQUATIONS
by
Russell E. Bauer
B.A., University of Northern Colorado, 1985
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
1995
This thesis for the Master of Science
degree by
Russell E. Bauer
has been approved
by
Weldon Lodwick
7/2 5/?S
Date
Bauer, Russell (M.S., Applied Mathematics)
Discontinuous Galerkin Methods for Ordinary Differential Equations
Thesis directed by Associate Professor Leo Franca
ABSTRACT
The discontinuous Galerkin method is applied to ordinary differential
equations and is shown to have the property of Astability. For polynomial
spaces of degree q 0, 1, 2 and 3, this method is applied to a model equation
and is shown to be stiff Astable and of order 2# + 1.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
in
To my parents
Contents
1 Introduction 1
1.1 AStable Methods...................................... 3
1.2 Stiff AStable Methods................................ 5
2 Discontinuous Galerkin Method 8
2.1 Astability for the DGM.............................. 11
3 Discontinuous Galerkin Method with piecewise polynomials
of degree q = 0,1,2,3 13
3.1 DGM(O)............................................... 15
3.2 DGM(l)............................................... 17
3.3 DGM(2)................................................20
3.4 DGM(3)............................................... 23
IV
4 Computational Results 29
5 Conclusion 39
V
List of Figures
2.1 Discretization strategy of the discontinuous Galerkin method . 9
4.1 error vs. number of intervals, N, for DGM(O) ...................32
4.2 error vs. number of intervals, N, for DGM(l) ...................33
4.3 error vs. number of intervals, IV, for DGM(2) .................34
4.4 error vs. number of intervals, N, for DGM(3) ..................35
4.5 DGM(O) for model equation with A = 30 and N = A...........36
4.6 DGM(2) for model equation with A = 30 and N = 4...........37
4.7 Trapezoidal method for model equation with A = 30 and
iV = 4........................................................... 38
vi
List of Tables
4.1 DGM(O) applied to model equation with A = 1 for various
number of intervals, N.....................................31
4.2 DGM(l) applied to model equation with A = 1 for various
number of intervals, N.....................................31
4.3 DGM(2) applied to model equation with A = 1 for various
number of intervals, N.....................................31
4.4 DGM(3) applied to model equation with A = 1 for various
number of intervals, N.................................... 32
4.5 Comparison of methods for A = 30 and n = 4 ..............33
4.6 Comparison of methods for A = 30 and n = 8 ..............34
4.7 Comparison of methods for A = 30 and N = 16 .............35
5.1 Comparison of discontinuous Galerkin methods .............39
Vll
ACKNOWLEDGEMENTS
I would like to thank my advisor, Professor Leopoldo Franca for his guid
ance and patience.
Vlll
Chapter 1
Introduction
In this work, we consider the approximate numerical solution of the ordinary
differential equation
y(t) = /(*,), t > 0, (11)
v( o) = yo
A numerical solution to this equation can be obtained by calculating the
approximate solution at a sequence of discrete points, say tn, n = 0,1,N.
Let yn represent this approximate solution at tn and y(tn) the exact solution
at tn. If these discrete points are equally spaced, then tn = nh, where the
step size h = tn tn~l.
Consider the model differential equation
y(t) = Ay, y(0) = 1, (1.2)
1
where A is a (complex) constant. The region of absolute stability for a nu
merical method is the set of (complex) values Ah for which all numerical
solutions of equation 1.2 will remain bounded as n oo [6j. A desirable
property of a numerical method is Astability. If the region of absolute sta
bility includes all values of Ah, where Re(Xh) < 0, then the method is said
to be Astable. In other words, a method is Astable if all numerical approx
imations of equation 1.2 tend to zero as n oo, when the method is applied
with fixed positive step size h to equation 1.2, where i?e(A) <0 [6].
If y(tn) yn = 0(hp+1), then the local truncation error is said to be of
order p + 1. The sum of the local errors is the global truncation error. If
the local error is 0(hp+1), then the global error is 0(hp), since the number
of local errors is proportional to h1.
Consider the power series Y,?La cizl representing a function f(z), so that
f(z) = Z!ci^
t=0
A Pade approximant, Ri>m for f{z) is a rational fraction [2, 5]
= Pi,m{z) = Qp + aiz + ... + aiz1
l,m Qi,m(z) bo + bxz + ... + bmzm
The Rn,n axe called the diagonal Pade approximants and Rni,n
(1.3)
are the subdi
2
agonal Pade approximation to ez [2]. Some properties of Pade approximants
are [2, 5, 9]
\Rl,m{z)\ < 1, whenever Re(z) < 0 for l = rti, m 1 or m 2 (1.4)
Rn,n > (1) as Re(z) oo (1.5)
Rl,m(z) 0 as Re(z) > oo for l = m lor m 2 (1.6)
W,/WI = 0(z'++1) (1.7)
1.1 AStable Methods
A general linear fcstep method for computing the approximate numerical
solution is defined by the formula
OfcJ/n+fc + akiyn+k 1 + ... + *0yn = h(/3kfn+k + + Po fn)i (18)
where a,, Pi are real constants, i = 0,1,k, 0, and h > 0 is the step
size, i.e. tm = mh, and fm = where ym is the approximation
to y(tm). If Pk = 0, then the method is explicit. The method is implicit
when Pk 7^ 0, since the right hand side of equation 1.8 would then contain a
function of yn+k.
3
Dahlquist [6] has shown that any explicit fcstep method is not Astable
and that the highest order of the global truncation error can not exceed 2 for
an Astable implicit linear multistep method. Therefore, there are no high
order Astable linear multistep methods. The smallest truncation error for
order p = 2 Astable methods is obtained for the trapezoidal rule; k = 1,
a0 = = 1 and j30 = f31 = 1/2 in equation 1.8. But, the use of linear
multistep methods with Richardson extrapolation can be used to obtain an
order of p = 4 Astable procedure [6].
The explicit RungeKutta methods are not Astable [10], but there are
implicit RungeKutta methods that are Astable [4, 8]. An rstage implicit
RungeKutta method has the form
kq = hf(yn + fiqjkj), q = 1,2,..., r j = 1,..., q 1
!/=!/ + T.H %
3=1
The implicit rstage RungeKutta method of order p = 2r applied to the
model equation 1.2 results in
2/n+l = Vn
where Rrtr(\h) is the rth Pade approximant to eXh [5, 3, 8, 10]. Therefore,
4
from equation 1.4
lR^rl < 1, if Re(Xh) < 0.
Hence, the methods are Astable. Also, rstage RungeKutta methods of
order p = 2r 1 and 2r 2 are Astable [3, 8].
Also, there is another class of Astable methods given by an extension of
Taylors series [10]. One case of this gives the trapezoidal rule. Another class
is a quadrature type [1]. For more discussion on Astability see [1, 3, 4, 5, 8,
9, 10, 13]
1.2 Stiff ^4Stable Methods
For a certain class of differential equations called stiff equations, another
desirable property of numerical methods is stiff Astability. Stiff equations
arise frequently in network analysis, simulation of chemical reactions and in
control systems. In these large systems, processes have significantly different
time scales, i.e. one part of the system may be decaying faster than another
part.
Accordingly, it is reasonable to consider the behavior for a numerical
method when Afe is large or when Re(Xh) > oo. In particular, if for
5
equation 1.2, yn+1/yn * 0 as Re(Xh) > oo, then the rapidly decaying
components would also decay rapidly in the numerical method. This is indeed
the definition for a method to be stiff Astable [1, 10].
The rstage implicit RungeKutta methods of order p = 2r are not stiff bi
stable, since for equation 1.2 one obtains yn+1 = Rr
the diagonal Pade approximant to e2, and Rrtr(p) * (l)r as Re(z) > oo
[3, 8, 9).
However, there are 9stage RungeKutta methods with order of 2r 1
that axe stiff Astable [10]. In this case
yn+1 = Rri,r(\h) y11.
Since
Rrif 0 as Re(z) > 00,
this class of methods is stiff Astable.
Also, the extension of Taylors series mentioned above can be made to
be stiff Astable [10]. The order is also reduced by one. The above methods
have the drawback of involving a lot of work.
For more discussion on these and on other methods that are stiff Astable
6
see [1, 3, 8, 9, 10].
In this work, we will show that the discontinuous Galerkin method [7,
11, 12] (DGM), exhibits the property of ^stability for any desired order of
global truncation error. We will also show that this method exhibits the
property of stiff Astability, for order of p = 1,3,5, and 7. We will then give
some examples where these properties are desirable.
7
Chapter 2
Discontinuous Galerkin
Method
To describe the DGM, consider the following notation:
In = (tn~\tn)
Pq{In) = { : V(t) = Y,VJ'
=0
vl = lim v(tn + a)
f s+0+
u" = lim v(tn + s)
s*0~
The discontinuous Galerkin method applied to the test equation
V = Ay, ?/(0) = 1
can be stated as: find yh P7(/n) for n = 1,2,N, such that
Bn{vh,yh) = Ln(vh), VUfe e Pq(In), (2.1)
8
(input)
yT1
(unknowns)
vl~l
yl
* discrete values
tnl tn
interval In
Figure 2.1: Discretization strategy of the discontinuous Galerkin method
where
Bn(v,y)= j {yhVXyhV^dt + yll'vll1
JIn
Ln(v) = ylZ^1
and where = 1, is the initial value. Note that yh and Vh Pq{In) for
n = 1,2,TV, that is yh and Vh are allowed to be discontinuous at the mesh
points, see figure 2.1.
To show how the DGM is related to the test equation 1.2, we define the
following space:
V = {w : v is continuous on the intervals In}
Let v(t) E V be an arbitrary function. Assume that y(t) is a solution to the
9
model equation.
We multiply the model equation by v to obtain
yv A yv = 0.
Then we integrate over the interval (0,1), with the continuity enforced weakly.
Thus, we have
E [ / (   y~l>T'
n=l UIn
= 0.
Since v varies independently on each interval this can be written as
[ (yv Ayv)dt + = y^v^1, n = 1,2,..., N.
J In
Hence, a variational problem for the model equation is: find y for n =
1,2,..., N, such that
Bn(v,y) = Ln, W G V.
With y replaced with yh Pq(In) and v replaced with Vh G Pg(In), we have
the DGM, equation 2.1.
The subscript h on y and v and the subscript n on B and L will be
dropped for brevity.
10
2.1 ^stability for the DGM
With v = y on the left hand side of equation 2.1, we have
B(y, y) = f (yy \y2)dt + (y+1)2. (2.2)
*'ir
By integration by parts,
Therefore,
JIn(vw=ife)2 for1)2!
Substituting this into equation 2.1 and taking the real part of the equation
gives
iie(A) Â£ y2dt + ifo!))2 + = VTV+1 < jfor1)2 + \(vVf
Then
Re(\)Jj2di + ifoll)2 < jfor1)2.
From this it follows that V Re(A) < 0, (yZ)2 iZ (yZ 1)2, or
lull < for11.
(2.3)
11
Therefore we have established that the method is Astable for any value of
q, (i.e., for an arbitrary high order accurate approximation), a result not
attainable for multistep methods, for example.
12
Chapter 3
Discontinuous Galerkin
Method with piecewise
polynomials of degree
9 = 0,1,2,3
Evaluating the integrals for B(v,y) can be much simplified by using the
transformation from the global domain / = (tn~1,tn) to a local domain.
This transformation Â£ : In (1,1) can be represented by
Â£(i) = ci + c2t,
where ci and c2 are constants. Thus,
et*"1) = i = d + c2tn_i
Â£(<) = 1 = Cl + C2tn
13
This yields
t =
21 t711 tn
h '
It follows from this that Â£* = 2/h. Solving for t gives
hÂ£ + fn_1 + tn
t =
Then = h/2 and dt = = (h/2)dÂ£. Thus
/ {yv ~ Ayv)dt = [ (y'iOtMO A
Jin J1
= Â£ A(W()d( (AA/2) Â£*({Wflif
Thus, the variational equation (2.1) can be stated as find y P9(1,1), for
n = 1,2,..., iV, such that
P(u,t/) = T(u), Vu Pg(1,1), (3.1)
where
b(v, y) = Â£  y Â£ vWMO# + vl^V <32)
and
i(u) =
1+
l
14
The shape functions are also transformed from the global to the local
coordinate system, so that y(Â£) and v(Â£) are expressed in terms of the shape
functions, $js, i = 1,..., q + 1. We will denote the discontinuous Galerkin
method with y(t) and v(t) as polynomials of degree q or less as DGM(g).
3.1 DGM(O)
In the case for q = 0, v(t) and y(t) are polynomials of degree zero. Hence,
y" = ?/"1 and v_ = u^_1. Denote these as yn and vn. The shape function
in the local system is $(Â£) = 1. So, $'(Â£) = 0. Thus, we have
y(0 = mvn = yn
y'{ 0 = nt)yn = o
(0 = = n
Then, equation 3.1 is
y ifVdJ + yV = y1v"
=S(W + S,>* = !,"V.
Since vn is arbitrary, we have
yn = M0yn~l =
1
.71 1
1A h
(3.3)
15
This is the implicit Euler method. The amplification factor, Mo, is the
Pade approximant i2o,i to eXh. From the properties of Pade approximants,
equations 1.4, 1.6 and 1.7, we have
\Mo(Xh)\ < 1, whenever Re(Xh) < 0
\Mo(Xh)\ > 0 as Re(Xh) > oo
M0(Xh) eXh
0{h2)
We have already shown that the DGM is Astable, but the second and third
expressions above show that the DGM(O) is stiff Astable and that the global
truncation error is of order p = 1.
Consider the model equation 1.2 for real negative values of A. The solution
is y(t) = eXt. Let 6t > 0 be some increment of time. Then,
y(t + St) < y(t).
In other words, the solution of equation 1.2 with real A is monotone. It would
be desirable to have the numerical method exhibit this property too.
For real values of A < 0,
> 
16
Recall that Mo < 1 for Re(A) < 0. Therefore,
yn < iT1.
That is, DGM(O) is a monotone method for real negative A.
3.2 DGM(l)
For q = 1, we have in the local coordinate system, the following shape
functions and their derivatives: $i(Â£) = 1/2(1 Â£), $2(0 = 1/2(1 + Â£),
$1(0 = 1/2 and $1(0 = 1/2. Also,
w'tO =
({)=mw1+
The following results will be needed.
fh $i$l# = 1/2 fh $2$1^ = 1/2
fh $i$l^ = 1/2 $1^ = 1/2
/21$K = l/3 fi $i$2^ = 1/3 A$lde = 2/3.
The first integral in equation 3.1 can be evaluated as follows:
/I v'(t)v(m=+*>:)
17
= (y+1 j\ i Vydi + yl Â£ vl1
+ (y+_1 j\ Wi# + V f *2&2dt) V
= (^r1+\y~>v+(vr1+\v)
For the second integral:
t L rtMtw=t /!1(*irI+^x^ior1+*2vi)t,(d(
\ h fl rl
= ^ter1 *?<*e+Â£ si^r1]
\A r1 /i
+Tfer *1
= yK^r1 + j^K'1 + (r + fjCKl
So the variational equation (3.1) is:
(^r+\y>r'
,AA (T!/+ 1 + ysCK
+ y+_1u+_1
+ (^r1 + \y)
1 + (fvr + y):
=* [(1/2 AtySJsrJ"1 + (1/2 AA/GJ^K1
+[(1/2 Xh/fyyl1 + (1/2 \hlZ)yn_]vn_ = y^'v
18
c +
Since v2 and 1 are arbitrary, we have two equations
/ 1 AA \
(2 3"^+ + ^2 ~ "6~ 2/_ ~~
1 Ah. , .1 Ah.
(~2 6~^+ + ^2 ~~ ~3~^_ ~ '
In matrix form, this can be written as
' 32A h 3A h ' ' vV ' ' Gyr1 '
3A h 32A h y 0
This matrix can then be row reduced to solve for y1 and y" in terms of
y"1. For the purpose of developing a numerical method only y" is needed.
It is
The amplification factor, M\, is the subdiagonal Pade approximant Ri<2(Xh)
to eAA. From the properties in equations 1.6 and 1.7, we have that
\Mi(Xh)\ 0 as Re(Xh) > oo
Mi eXh = 0(h2+1+1) = 0(h4).
Therefore, in addition to Astability, DGM(l) exhibits stiff Astability, and
has order of 3.
19
The numerator of Mi(Xh) has a zero at Ah = 3, and is negative for
Ah < 3. The denominator is positive for all real Ah < 0. Thus Mi(Xh) < 0
for real negative A whenever the step size h > 3/A. The sequence of
numerical solutions will then oscillate between negative and positive values.
Therefore, DGM(l) is not in general, monotone. But the method will exhibit
monotonicity for step size chosen small enough. For example, with A = 30,
choosing < 0.1 will result in this method being monotone.
3.3 DGM(2)
For q 2, the shape functions in the local coordinates are:
*i(0 = Â§(00
*2(0 = 0 +1
*3(0 = 1(0 + 0
The derivatives are:
m)=i 1
$'2) = ~2(
m)=+Â§
The test and weight functions expressed in terms of the shape functions are:
2/(0 = *!J/r1 + *2S/"1/2 + *32/
y'( 0 = *wrx + *2 yn~1/2 + *3 y
t7(0 = *1U+'1 + *2Unl/2 + *3^.
20
The following integrals will be needed:
fh *1*1% = 1/2 A = 2/3 A = 1/6
fli = 2/3 A $2*Â£dÂ£ = 0 A *2*^ = 2/3
A 3ide = i/6 A = 2/3 A wide = 1/2
A $2^ = 4/i5 A $2^ = 16/15 = 4/15
A*i
The first integral in equation 3.2 is evaluated using the above results as
follows:
J ^ '(Â£MÂ£)df
= f (^iy+_1 + ^y71172 + $32/")($i<1 + $2un_1/2 + $3tÂ£)#
= (ACMiyr1 + ^x^yn'1/2 + *i*sy)df) ^r1
+ (A^W1 + *2&2yn~1,i + $2$MdÂ£)
+ (A($3$iy+1 + 3&2yn~1/2 + $3$^)#)
= (Jyr1 + iy""1/2 Â£yK1
+ (y+1 + o t/n1/2 +
+ (iy+1fyn"1/2 + Â§y),;
For the second integral in equation 3.2,
y iMwm
21
= y Jj^y+~x + W~1/2 + + $2Vn~1/2 + $3 vn_)dÂ£
= f [/i^frr1 + i$2yn1/2 + $i*3Â£)<*f] ur1
+ f [/X1($1$2^1 + *lyn1/2 + *2*3y2K] u"1/2
+ f [/ic^i^r1 + ^2^3yn~1/2 + $hl)dt] vZ.
+
+
t (ir1 + 4n_1/ai)r1
fc^r1 + ify"172+^K~1/2
f (iyr1 + ^n"1/2 + r5y)v
So from the variational equation 3.1 or 3.2 and the fact that v" x, vn 1!2
and vZ are arbitrary, we have the three equations:
1 2Xh i 2 Xh 1/2 1 Aft = nl.
4 15 y+ ^ 3 15 ^ +tl 6 + 30 y~ V~
( Jy;1 y^1/2 + ( )j, = o
v 3 15 y+ 15 y v3 15 y~
.1 Ahs i 2 Xh < i /n /I 2A/i,
( + 80 W + ("3  ' + (2 TT* =
In matrix form, this becomes
' 15 4U 20 2Ah 5 + Ah ' V 1 ' 30y"1 '
10 Ah 8Ah 10 Ah = 0
5 + Ah 20 2Ah 15 4Ah I S  1 1 0 1
22
Upon solving this system, we obtain
yn = M2yT1 =
3(A 2h2 + 8A h + 20) B_x
X3h3 + 9X2h2 36Ah + 60 V~ '
In this case, M2(Xh) is the subdiagonal Pade approximant R2^(Xh) to
eXh. Again from equations 1.6 and 1.7,
M2(A*) 
M2(Xh) eXh
0 as Re(Xh) oo
= 0(h3+2+1) = 0(h6).
Therefore, DMG(2) is stiff Astable and has order p = 5.
The numerator of M2(Xh) has two complex roots and the denominator
has two complex roots and one positive real root. Then, since M2(0) = 1 > 0,
M2(Xh) > 0 for all negative real Ah. Therefore,
yn < yn 1
That is DGM(2) is monotone.
3.4 DGM(3)
For 9 = 3, the shape functions are
*i(0 = ^( + l)(3fl)l)
23
*a(f) = ^(f+l)(3Â£l)(Â£l)
*(Â£) = +1)(3Â£ +1)(Â£ 1)
*4(Â£) = ^(Â£ + l)(3Â£ + l)(3Â£l).
The derivatives are
m)
m)
m)
m)
27Â£2 18Â£ 1
16
9(9e2 2e 3)
16
9(9^2 2Â£ 3)
16
27Â£2 +18Â£ 1
16
And
y(0 = Ziyi'1 + *2 yn~2'z + ^yn~1/3 + $4y_
yW = ^r1 + yn~2/3 + % yn~1/3 + *iv"
v(() = + 2Vn2'3 + $3i0n"l/3 +
24
The following integrals will be needed:
/^l*!*!# = 1/2
/_i 1*1*4# = 7/80
1$2$3 # = 81/80
/_! 1*3*'2 # = 81/80
/_! 1*4*1# = 7/80
U l*4*i# = 1/2
/_! l*i*3# = 3//70
/_! 1*2*3# = 27/280
/_! 1*3*4# = 33/280
Then,
fh *i*'2# = 57/80
fh *2*1# = 57/80
fh *2*l# = 3/10
/i Widf = 0
/i *4*1# = 3/10
fh *?# = 16/105
fh *i*4# = 19/840
fh *2*4# = 3/70
fh *1# = 16/105
Â£ I^MO#
fh *1*3# = 3/10
fh *2*'2# = 0
/2a *3*1# = 3/10
/2a *3*;# = 57/80
fh *4*3# = 57/80
fh *1*2# = 33/280
fh *1# = 27/35
/_! *2# = 27/35
+ *i2)
= Â£ t*;^1 + *'2yn2/3 + *'yn1/3
t*!^1 + *2^n"2/3 + *3Un_1/3 + *4tÂ£.)#
= Â£(*a$iy1 + *! *lj/n2/3 + *i*^2/n1/3 + *i*l2/!!)#<1
+ + *2 *lj/n2/3 + *2*lj/n1/3 + *2*^)#wn'2/3
+ + *3*22/ 2/3 + *3*,3yn'1/3 + *3*ii/)#t;"1/3
+ Â£(*4*'12/r1 + *4*,2yTl_2/3 + *4 *3?/711/3 + *4*4lÂ£)#K
+ + o y"2/3 + y1,3s ^y>2/3
+ (^i1 S*'"2'3+0 y"~1,3y+I*)"1'3
25
+ yT1 + yn~2/3 yn~1/3y + yn)vn
K 1(T+ T i(T 8(T y T~} ~
Also,
Y Liy^v^d^
= Y Â£($ i^r1 + ^yn~2/z + $3 yn~1/3 + ^ynY
+ $2un_2/3 + 3vn~l/3 + $4 vn_.)d(
= % +i$2yn~2/3+$i$3yni/3+$i$4^k)]
+ f [A^SitfT1 + $&n213 + ^2$3y"1/3 + $2$42/!lMO] n_2/3
= f [/i^itfT1 + $3$2yn2/3 + $iy"1/3 + $3$4
= Â¥ [flA^Qiyr1 + $4$2yn~2/3 + $4$3yn1/3 + ^yH)#)] vi
= fiwsvr1 + mvn~m ~ r0yn1/3y + &K"1
+ f(~w0yV +1 yn~2/3 M>yn'1/3 ^K2/3
+ Â¥(hyl~x ~ myn~2'z +1yn~113 + iyX"1/3
+ f (ffivr1 IX2/3 + w0yn~1/3y + w5y)v
In a similar fashion as before, use the fact that u1, vn 2/3, vn 1^3 and
26
u" are arbitrary, to obtain the following matrix equation
840 128A/i 1197 99Ah 504 + 36A h 14719AA ' r vV 1
1197 99Ah 648Ah 170 + Xh 504 + 36A h y7l2/3
504 + 36A h 1701 +8lXh 648A h 1197 99A h y711/3
504 + 36AA 504 + 36A h 1197 99Ah 840 128Ah . . yl .
16807/r1
0
0
o
Upon solving this system we have for y
yn_ = M3{Xh)y
4(A 3h3 + 15A2h2 + 90 Xh + 210) n1
\4hi16\3h3 + m\2h2m\h + 840y ^
For q = 3, M3(Xh) is the Pade approximant R3i4(Xh) to eXh. therefore,
\M3(Xh)\
M3(Xh) 
v 0 as Re(Xh) > oo.
= 0(hA+3+1) = 0(h&)
Thus, DGM(3) is stiff Astable and has order p = 7.
The denominator of M3(Xh) has only complex zeros, therefore its sign
does not change. The numerator has two complex zeros, but one negative
real zero, x\ = (5\/65)1/3 (5S/6+5)1/3 5 which is approximately 5.65.
Since M3(0) = 1, M3(Xh) > 0 for Xh (a;i, 0], and M3 < 0 for Xh < 3.
Thus, DGM(3) is monotone for negative A only when the step size h < x\/X.
27
This is not a very strong restriction. For example, with A = 30, h can be
as large as 0.18.
28
Chapter 4
Computational Results
All calculations were done with 15 decimal places of accuracy. For A = 1,
tables 4.1, 4.2, 4.3 and 4.4 show the computational results for various n. The
last column of each table gives the absolute value of the error divided by
h raised to the power of the order of the method. Note that these values
are fairly constant. This verifies the results obtained earlier in this work for
the order of DGM. Figures 4.1, 4.2, 4.3 and 4.4 are loglog graphs of the
absolute value of the error versus the number of intervals, n. Since n = 1/h,
the negative of the slope on each graph gives approximately the order for the
method. Note that for DGM(3), the result is accurate to 15 decimal places
for n = 16. Therefore, since these calculations are accurate to 15 decimal
places, larger values of n do not increase the precision.
29
In tables 4.5, 4.6 and 4.7, the approximate value at each time step for four
methods is shown. The last column gives the true solution to four decimal
places. The four methods are: a method that is not Astable, Heuns method,
which is a twostage RungeKutta method; the order of p = 2 trapezoidal
method, which is Astable but is not stiff Astable; The DGM(O) method
or Eulers method with order only p = 1, which is both ^4stable and stiff
Astable; and the DGM(2) with order p = 5 and both stability properties
mentioned above. These methods were applied to the model equation (1.2),
but with A = 30. Figures 4.5, 4.6 and 4.7 show these results graphically.
Heims method gives absurd results, so it was not shown on a graph.
Notice that the DGM(O) actually gives better results than the trapezoidal
method, even though it is of lower order. Then, of course the DGM(2) comes
the closest to the true solution curve. The better results for the DGMs for
this test case is important in systems of equations that have some eigenvalues
with very large moduli. These methods will not restrict the step size as much
as other methods that are not stiff Astable.
30
N computed value of 2/(1) error error/h
1 5.000000000000000e01 1.321205588285577e01 1.321205588285577e01
2 4.444444444444444e01 7.656500327300209e02 1.531300065460042e01
4 4.096000000000000e01 4.172055882855769e02 1.668822353142307e01
8 3.897443431289458e01 2.186490195750346e02 1.749192156600277e01
16 3.790853319179361e01 1.120589074649381e02 1.792942519439009e01
32 3.735538614900618e01 5.674420318619422e03 1.815814501958215e01
64 3.707349329009732e01 2.855491729530868e03 1.827514706899755e01
128 3.693118106024922e01 1.432369431049851e03 1.833432871743810e01
Table 4.1: DGM(O) applied to model equation with A = 1 for various
number of intervals, N
N computed value of 2/(1) error error fk6
1 3.636363636363636e01 4.243077535078688e03 4.243077535078688e03
2 3.673094582185491e01 5.699829528932332e04 4.559863623145866e03
4 3.678043951904257e01 7.504598101665261e05 4.802942785065767e03
8 3.678697774589969e01 9.663712445429073e06 4.947820772059686e03
16 3.678782140046069e01 1.227166835482407e06 5.026475358135940e03
32 3.678792865263700e01 1.546450723588499e07 5.067409731054795e03
64 3.678794217611732e01 1.941026911422838e08 5.088285586680286e03
128 3.678794387401324e01 2.431309931161252e09 5.098826484754682e03
Table 4.2: DGM(l) applied to model equation with A = 1 for various
number of intervals, N
N computed value of y(l) error error/h?
1 3.679245283018868e01 4.508713044443224e05 4.508713044443224e05
2 3.678809236447542e01 1.482473311886778e06 4.743914598037691e05
4 3.678794891116255e01 4.794018315656956e08 4.909074755232723e05
8 3.678794426987461e01 1.527303783444012e09 5.004669037589338e05
16 3.678794412196591e01 4.821676391486562ell 5.055894143879414e05
32 3.678794411729567e01 1.514399716739945e12 5.081482231616974e05
64 3.678794411714896e01 4.729550084903167e14 5.078315734863281e05
128 3.678794411714437e01 1.332267629550188e15 4.577636718750000e05
Table 4.3: DGM(2) applied to model equation with A = 1 for various
number of intervals, N
31
N computed value of y(l) error error/h1
1 3.678792038435141e01 2.373279282541496e07 2.373279282541496e07
2 3.678794392443099e01 1.927132398105158e09 2.466729469574602e07
4 3.678794411559969e01 1.544542271858518ell 2.530578058212996e07
8 3.678794411713199e01 1.224575996161548e13 2.568121999502182e07
16 3.678794411714413e01 1.054711873393899e15 2.831220626831055e07
32 3.678794411714424e01 1.110223024625157e16 3.814697265625000e06
64 3.678794411714422e01 1.110223024625157e16 4.882812500000000e04
128 3.678794411714427e01 3.330669073875470e16 1.875000000000000e01
Table 4.4: DGM(3) applied to model equation with A = 1 for various
number of intervals, N
Figure 4.1: error vs. number of intervals, N, for DGM(O)
32
Figure 4.2: \error\ vs. number of intervals, TV, for DGM(l)
step DGM(O) DGM(2) trapezoidal Heuns 9 II C* >
0 l.OOOOe + 00 l.OOOOe+ 00 l.OOOOe+ 00 l.OOOOe+ 00 l.OOOOe+ 00
1 1.1765e 01 3.8748e 02 5.7895e 01 2.1625e + 01 5.5308e 04
2 1.3841e 02 1.5014e 03 3.3518e 01 4.6764e + 02 3.0590e 07
3 1.6283e 03 5.8177e 05 1.9405e 01 1.0113e + 04 1.6919e 10
4 1.9157e 04 2.2543e 06 1.1235e 01 2.1869e + 05 9.3576e 14
Table 4.5: Comparison of methods for A = 30 and n = 4
33
Figure 4.3: error vs. number of intervals, N, for DGM(2)
step DGM(O) DGM(2) trapezoidal Heuns (<") =
0 l.OOOOe + 00 l.OOOOe+ 00 l.OOOOe+ 00 l.OOOOe+ 00 l.OOOOe+ 00
1 2.1053e 01 3.2561e 02 3.0435e 01 4.2812e + 00 2.3518e 02
2 4.4321e 02 1.0602e 03 9.2628e 02 1.8329e + 01 5.5308e 04
3 9.3308e 03 3.4522e 05 2.8191e 02 7.8471e + 01 1.3007e 05
4 1.9644e 03 1.1241e 06 8.5799e 03 3.3596e + 02 3.0590e 07
5 4.1355e 04 3.6601e 08 2.6113e 03 1.4383e + 03 7.1941e 09
6 8.7064e 05 1.1918e 09 7.9473e 04 6.1578e + 03 1.6919e 10
7 1.8329e 05 3.8805e 11 2.4188e 04 2.6363e + 04 3.9790e 12
8 3.8588e 06 1.2635e 12 7.3614e 05 1.1287e + 05 9.3576e 14
Table 4.6: Comparison of methods for A = 30 and n = 8
34
Figure 4.4: \error\ vs. number of intervals, AT, for DGM(3)
step DGM(O) DGM(2) trapezoidal Heuns y(tn) = eM
0 1.0000e + 00 l.OOOOe + 00 l.OOOOe+ 00 l.OOOOe+ 00 l.OOOOe + 0
1 3.4783e 01 1.5415e 01 3.2258e 02 8.8281e 01 1.5335e 01
2 1.2098e 01 2.3761e 02 1.0406e 03 7.7936e 01 2.3518e 02
3 4.2081e 02 3.6626e 03 3.3567e 05 6.8803e 01 3.6066e 03
:
16 4.5898e 08 1.0160e 13 1.3747e 24 1.3611e 01 9.3576e 14
Table 4.7: Comparison of methods for A = 30 and N = 16
35
Figure 4.5: DGM(O) for model equation with A = 30 and N = 4
36
Figure 4.6: DGM(2) for model equation with A = 30 and N = 4
37
Figure 4.7: Trapezoidal method for model equation with A = 30 and N = 4
38
Chapter 5
Conclusion
In this work, we showed that the discontinuous Galerkin methods for ordi
nary differential equations are Astable. With piecewise polynomial spaces
of degree q = 0, 1, 2 and 3, we showed that the order is p = 2k + 1. Also,
in these cases, the methods are stiff Astable. Table 5.1 gives a summary of
thesis results.
Method Order Astable Stiffly Astable Monotone
DGM(O) 1 yes yes unconditionally
DGM(l) 3 yes yes for h < 3/A
DGM(2) 5 yes yes unconditionally
DGM(3) 7 yes yes for h < 5.65/A
Table 5.1: Comparison of discontinuous Galerkin methods
39
Bibliography
[1] 0. Axelsson. A class of Astable methods. BIT, 9:185199, 1969.
[2] G. A. Baker, Jr. and P. GravesMorris. Pade Approximants, Part I:
Basic Theory. AddisonWesley, 1981.
[3] J. C. Butcher. On Astable implicit RungeKutta methods. BIT, 17:375
378, 1977.
[4] J. C. Butcher. Implicit RungeKutta processes. Mathematics of Com
putations, 18:5064, 1978.
[5] J. C. Butcher. The Numerical Analysis of Ordinary Differential Equa
tions. John Wiley and Sons, 1987.
[6] G. Dahlquist. A special stability problem for linear multistep methods.
BIT, 3:2743, 1963.
[7] M. Delfour, W. Hager, and F. Trochu. Discontinuous galerkin meth
ods for ordinary differential equations. Mathematics of Computation,
36(154):455473, 1981.
[8] B. L. Ehle. High order Astable methods of the numerical solution of
systems of differential equations. BIT, 8:276278, 1968.
[9] B. L. Ehle. Astable methods and Pade approximations to the exponen
tial. SIAM Journal of Mathematical Analysis, pages 671680, 1973.
[10] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential
Equations. PrenticeHall, 1971.
40
[11] C. Johnson. Numerical Solution of Partial Differential Equations by the
Finite Element Method. Cambridge, 1987.
[12] P. Lesaint and P. A. Raviart. On a Finite Element Method for Solving
the Neutron Transport Equation, in Mathematical Aspects of Finite Ele
ments in Partial Differential Equations, pages 89123. Academic Press,
1974.
[13] H. A. Watts and L. F. Shampine. Astable block implicit onestep meth
ods. BIT, 12:252266, 1972.
41

PAGE 1
DISCONTINUOUS GALERKIN JVIETHODS FOR ORDINARY DIFFERENTIAL EQUATIONS by Russell E. Bauer B.A., University of Northern Colorado, A thesis submitted to the University of Colorado at. Denver in partial fulfillrnent of the requirements for t.he degree of of Science Applied Mathematics 1995
PAGE 2
This thesis for the Master of Science deg ree by Russell E. Bauer has been approved by Leopoldo Franca Thornas Rus sell vVeldon Lodwick Date
PAGE 3
Bauer, Russell (M.S., Applied Mathernatics) Discontinuous Galerkin for Ordinary Differential Equations Thesis directed by Associate Professor Leo Franca ABSTRACT The discontinuous Galerkin rnethod is applied to ordinary differential equations and is shown to have the property of Astability. For polynornial spaces of degree q = 0 1, 2 and :3, this rnethod is applied to a rnodel equation and is shown to be stiff Astable and of order 2q + 1. This abstract accurately represents the content of the candidate's thesis. I recommend its publication. Signed _________ Leopoldo Franca lll
PAGE 4
Contents 1 Introduction 1 1.1 AStable Methods :3 1.2 Stiff AStable Methods .J 2 Discontinuous Galerkin Method 8 2.1 Astability for the DGM . . . . . . . . 11 3 Discontinuous Galerkin Method with piecewise polynomials of degree q = 0 1 2,:3 13 :_u DC::l\'1(0) : _L2 DC::M(1) :3.:3 DC::M(2) : _L4 lV
PAGE 5
4 Computational Results 29 5 Conclusion 39 v
PAGE 6
List of Figures 2.1 Di scretization strateg y of the discontinuous Galerkin rnethod 4.1 !error! vs. nurnber of intervals, .1.V, for . . . . :32 4.2 !error! v s nurnber of interva ls, i"l, for DGM(1) . . . . :3:3 4.:3 !error! vs. nurnber of intervals, N, for DC::M(2) . . . . 4 .4 !error! v s nurnber of interval s N, for DGM(: n :3.5 4.5 DGM(O) for model equation with A= : w and N = 4 :36 4.6 DGM( 2) for model equation with A = : w and N = 4 :n 4. 7 Trapezoidal rnethod for model equation with A = :.H) and J.\' = 4 . . . . . . . . . . . . . . . :38 Vl
PAGE 7
List of Tables 4.1 DGM(O) applied to model equation with ). = 1 for various nmnber of intervals, N . . . . . . . . . . . :_n 4.2 DGM(1) applied to model equation with ). = 1 for various nmnber of intervals, N . . . . . . . . . . . :_n 4.:3 DGM(2) applied to model equation with ). = 1 for variou s nun1ber of intervals, N . . . . . . . . . . . :_n 4.4 D GM ( : n applied to model equation with ). = 1 for various nmnber of intervals, N . . . . . . . . . . . :32 4.5 Cornparison of rnethods for ). = : .HJ and n = 4 . . . . :3:3 4.6 Cornparison of rnethods for ). = :.H) and n = 8 . . . :34 4. 7 Cornparison of rnethods for ). = :.H) and N = 16 . . . :3.5 .5 .1 Cornparison of discontinuous Galer kin rnethods . . . . : .H) Vll
PAGE 8
ACKNOWLEDGEMENTS I would like to thank rny advisor, Professor Leopoldo Franca for his guidance and patience. Vlll
PAGE 9
Chapter 1 Introduction In this work, we consider the approxirnate nurnerical solution of the ordinary differential equation :iJ( t) = f(t,y), t 2 0, y(O) = :IJo. ( 1.1) A nurnerical solution to this equation can be obtained b y calculating the approxirnate solution at a sequence of discrete points s a y tn, n = 0, 1, ... ,N. Let y n represent this approximate solution at tn and :IJ( t n ) the exact solution at tn. If these discrete points are equally spaced, then tn = nh, where the step size h = t n tnl. Consider the rnodel differential equation :iJ( t) = A:IJ, y(O) = 1, 1 ( 1.2)
PAGE 10
where A. is a (complex) constant. The region of absolute stability for a numerical method is the set of ( cornplex) values A.h for which all nurnerical solutions of equation 1.2 will rernain bounded as n > oo [ G ] A desirable property of a nurnerical method is Astability. If the region of absolute stability includes all values of A.h, where Re(A.h) < 0, then the method is said to be Astable. In other words, a rnethod is Astable if all nurnerical approxirnations of equation 1.2 tend to zero as n > oo when the method is applied with fixed positive step size h to equation 1.2, where Re(A.) < 0 [ G ] If l y(f'1 ) ynl = O(h P +1), then the local truncation error is said to be of order p + 1. The sun1 of the local errors is the global truncation error. If the local error is O(hP+1), then the global error is O(hP), since the nurnber of local errors is proportional to h1 Consider the power series C i Z1 representing a function f(z), so that X l f(z ) = Lcizi. i=O A Pade approxirnant, R 1,m for f(z ) is a rational fraction [ 2, .5] Pt.m (z) R t,m = ____ c'QI,m(z ) I ao + a1z + ... + a1z bo + b1z + ... + bmzm (l.:.n The Rn,n are called the diagonal Pade approximants and RnI,n are the subdi2
PAGE 11
agonal Pade approxirnation to e z [2]. Son1e properties of Pade approxirnants are [ 2 .5, J R ,,m(z ) J < 1, whenever Rt:(z) < 0 for l = m, rn1 or rn2 (1.4) R n,n + ( lt as Re(z)+ oo (1..5) R i,m(z ) + 0 as Re(z) + oo for l = rnlor m 2 (l.G) 1.1 AStable Methods A general linear k:step method for computing the approximate nurnerical solution is defined by the fonnula n+k: n+k:l n } ( .D j' .D j' ) O'k:Y + O:k:1:lJ + .. + O:oY = l /Jk: n+k: + .. + l'o n ( 1.8) where o:i, 31 are real constants, i = 0, 1, ... ,k:,o:k: =f 0, and h > 0 i s the step s1ze, 1.e. tm = mh, and fm = f( tm, y m), where ym is the approximation to y ( tm). If #k: = 0, then the method is explicit. The method is implicit when #k: =f 0, since the right. hand side of equation 1.8 would then contain a function of yn+k.
PAGE 12
Dahlquist [6] has shown that any explicit k:step method is not A.stable and that the highest order of the global truncation error can not exceed 2 for an A.stable irnplicit linear nmltistep method. Therefore, there are no high order A.stable linear rnultistep rnethods. The srnallest truncation error for order p = 2 A.stable methods is obtained for the trapezoidal rule; k: = 1, n0 = n1 = 1 and ,30 = .31 = 1/2 in equation 1.8. But, the use of linear multistep methods with Richardson extrapolation can be used to obtain an order of p = 4 A.stable procedure [ 6 ] The explicit RungeKutta methods are not A.stable [10], but there are irnplicit RungeKutta rnethods that are A.stable [ 4 8 ] An rstage implicit RungeKutta rnethod has the form k:q = hj'(:ljn + ,3qjk:j), q = 1 2 ... r j = 1 ... q 1 'I' u. n+l u. n + """ 1 ... " L ;JI\.J j=l The irnplicit rstage RungeKutta rnethod of order p model equation 1.2 results in :IJn+ l = Rr,r(Ah) Yn 2r applied to the where Rr.r(>..h) is the rth Pade approximant to e>.h [.5, :3, 8 10] Therefore, 4
PAGE 13
fron1 equation 1. 4 J Rr,r J < 1, if Re(>.h) < 0. Hence, the methods are Astable. Also, rstage RungeKutta methods of order p = 2r1 and 2r2 are Astable [:3, 8 ] Also, there is another class of Astable rnethods given by an extension of Taylor's series [10] One case of this gives the trapezoidal rule. Another class is a quadrature type [ 1 ] For rnore discussion on Astability see [1, :3, 4, .5, 8, 1.2 Stiff AStable Methods For a certain class of differential equations called stiff equations, another desirable property of nmnerical methods is stiff Astability. Stiff equations arise frequently in network analysis, sirnulation of chernical reactions and in control systems. In these large systerns, processes have significantly different time scales, i.e. one part of the systen1 rnay be decaying faster than another part. Accordingly, it is reasonable to consider the behavior for a nmnerical method when J>.hJ is large or when Re(>.h) + CXJ. In particular, if for
PAGE 14
equation 1.2, yn+l jyn + 0 as Re( A. h) + oo, then the rapidly decaying cornponents would also decay rapidly in the nurnerical rnethod. This is indeed the definition for a rnethod to be stiff Astable [ 1, 10] The rstage irnplicit RungeKut ta rnethods of order p = 2r are not stiff Astable, since for equation 1.2 one obtains :1/'+1 = Rr.r(A.h);t/', where R r r ( z ) is the diagonal Pade approxirnant to ez, and Rr,.(fl) + ( 1t as Re(z) + oo However, there are qstage RungeKutta methods with order of 2r 1 that are stiff Astable [10] In this case n+l R (\/. ) n Y = rl,r A t Y Since Rrl,r+ 0 a.s Re(z)+ oo, this class of rnethods is stiff Astable. Also, the extension of Taylor's series mentioned above can be rnade to be stiff Astable [10] The order is also reduced by one. The above rnethods have the drawback of involving a lot of work. For rnore discussion on these and on other methods that are stiff Astable
PAGE 15
see [1, :3, 8, 9, 10]. In this work, we will show that the discontinuous Galer kin method [7, 11, 12] exhibits the property of Astability for any desired order of global truncation error. \Ve will also show that this method exhibits the property of stiff Astability, for order of p = 1, :3, 5, and 7. vVe will then give sorne exarnples where these properties are desirable. 7
PAGE 16
Chapter 2 Discontinuous Galerkin Method To describe the consider the following notation: q Pq(In) = {v : v( t) =I: viC i=O = lim v(C1 + .s) S>0+ = lim v(C1 + .s) s>OThe discontinuous Galerkin rnethod applied to the test equation ;/J = >.y, y(O) = 1 can be stated as: find Y h E Pq(In) for n = 1,2, ... ,N, such that. 8 (2.1)
PAGE 17
(input) (unknowns) '1.1 n.1 '1/n.1 ,':J,':J+ y"_'_ values interval I n Figure 2.1: Di scretization strategy of the discontinuous Galerkin rnethod where ( n.1 n.1 Ln. v) = ihv h + and where :IJt = 1, is the initial value. Note that Yh and 'Vh E Pq(In) for n = 1 2, ... N that is Yh and V h are allowed to be discontinuous at the rnesh points, see figure 2.1. To show how the is related to the test equation 1.2, we define the followin g s pace: V = {v : v i s continuous on the intervals I n} Let v( t) E V be an arbitrary function. Assmne that y ( t) is a solution to the 9
PAGE 18
model equation. vVe nmltiply the rnodel equation by v to obtain :iJv>.yv = 0. Then we integrate over the interval (0, 1), with the continuity enforced weakly. Thus, we have N """"' [ { ( \ l ( n1 n1 n1] L Jr :zrvAlJV )d + Yf)vf= 0. n=1 ln Since v varies independently on each interval I n, this can be written as Hence, a variational problen1 for the rnodel equation 1s: find y for n 1 2, ... }l, such that B n (v,y) = Ln, Vv E V \"lith :IJ replaced with Y h E P q(In) and v replaced with V h E Pq(In), we have the DGM, equation 2.1. The subscript h on y and v and the subscript n on B and L will be dropped for brevity. 10
PAGE 19
2.1 Astability for the DGM \Vith v = y on the left hand s ide of equation 2.1, we have ( 1(. \2{ ( n12 B y y) = lJ:lJ A:lJ )d + y + ) l n (2.2) B y integration b y p arts Therefore, Sub stituting this into equation 2.1 and taking the real p art of the equation giVes (,.1 2l 1 ( n 2 1( n12 n1 n1 1( n12 1( n12 R e A d = < ) l n .lJ + 2 ,lj) + 2 .lJ + ) .lJ,lj + 2 .lJ) + 2 .lJ + ) Then Fron1 this it follo ws tha t V Re( ,\) < 0 :=; (y:::1 ) 2 or I n I < I n11 .lJ_ .lJ_ (2.:.n 11
PAGE 20
Therefore we have established that the rnethod is A stable for any value of q (i.e. for an arbitrary high order accurate approximation), a res ult not attainable for rnultistep rnethods, for exarnple. 12
PAGE 21
Chapter 3 Discontinuous Galerkin Method with piecewise polynomials of degree q==0,1,2,3 Evaluating the integrals for B(v, :IJ) can be nmch sirnplified by usmg the transfonnation from the "global" dornain In = ( tnl, tn) to a "local" dornain. This transfonnation : In + ( 1, 1) can be represented by where c1 and c2 are constants. Thus,
PAGE 22
This yields 2t tn1 tn h It follows frorn this that (1 = 2 /h. Solving fort gives t = h( + tn1 + tn. 2 Then t.( = h / 2 and dt = t,(d( = (h/2)d(. Thus f .UrvA.yv)dt = 11 (:1/(()(tv(() A:IJ(()v(())t,(d( }In 1 Thus, the variational equation (2.1) can be stated as find y E P1 ( 1, 1), for n = 1,2, ... ,N, such that B(v, y) = L(v), Vv E P1 ( 1, 1), (:3.1) where B(v, :IJ) = 11 y' (()v(()d(11 y(()v(()d( + (:3.2) 1 1 and ( n1 n1 L v) = y_ v+ 14
PAGE 23
The shape functions are also transformed fron1 the global to the local coordinate systern, s o that :IJ( and v( are expressed in terms of the shape functions, i 's i = 1 ... q + 1. \Ve will denote the discontinuous Galerkin method with y(t) and v(t) as polynornials of degree q or less as DG.M(q). 3.1 DGM(O) In the case for q = 0, v( t) and :IJ( t) are polynornials of degree zero. Hence, :IJ: = yf,1 and vr:_ = vf,1 Denote these as yn and vn. The shape function in the local systen1 i s (() = 1. So, '(() = 0. Thus we have '( . '( . n :IJ () = ():IJ = 0 Then, equation :.u is ).h 11 n n z n n n1 n y v + y v = y v 2 1. .., Since vn is arbitrary, we have n n 1 1 n.1 :IJ = Afo y = 1 ).h y 1.5
PAGE 24
This is the irnplicit Euler rnethod. The amplification factor, A10 is the Pade approxirnant R0 1 to e>h. Fron1 the properties of P ade approxirnants, equations 1.4, 1.() and 1.7, we have IJfo(A.h)l S 1, whenever Re(A.h) S 0 IA1o( A.h) I + 0 as Re( A. h) + oo \Ve have already shown that the DGM is Astable, but the second and third expressions above show that the DGM(O) i s s tiff Astable and that the global truncation error is of order p = 1. Consider the rnodel equation 1.2 for real negative values of A.. The solution is y( t) = e>t Let 8t > () be some increment of tirne. Then, y(t + 8t) < y(t). In other words, the solution of equation 1.2 with real ,\ is rnonotone. It would be desirable to have the numerical method exhibit this property too. For real values of ,\ < 0, 1 Afo(A.h) = 1 ,\h > 0. 1()
PAGE 25
Recall that I .. M0I < 1 for Re(>..) < 0. Therefore, That is, DGM(O) is a rnonotone rnethod for real negative>... 3.2 DGM(l) For q = 1, we have in the local coordinate systern, the following shape functions and their derivatives: 1(0 = 1/2(1= 1/2(1 + 1/2 and 1/2. Also, '(t) = <[>' ( t) n1 + ' (t) n .lJ ., 1 ., .lJ+ 2 ., .lJ_ ( "" ( n1 "" ( n v = + The following results will be needed. = 1 / 2 = 1/2 1 = 1j:3 = The first integral in equation ;u can be evaluated as follow s : 17
PAGE 26
For the second integral: A.h jl ( ( l A.h jl (m. n1 n (m. n1 m. n l 2 _1 iJ ()v ()<{ = 2 _1 + <[>2;//_:) + A.h[ ( n ljl m.2z n j 1 A'. z n1] = ;2 Y+ "'1<( + Y1 1 1 A.h [ ( n1 jl A'. A'. z n jl m.2 z ') n l +;2 Y + + y_ v_ 1 1 So the variational equation 1) is: + n1 n1 n1 n1 Y+ v+ = y_ L'+ 18 n1 n1 = y_ v +
PAGE 27
Since v: and are arbitrary, we have two equations ( 1 A.h (1 ),h n 2 (j)Y+ + 2 = O. In rnatrix form, this can be writ ten as [ :3.2A.h .:3 ,>..h l [ :lJ.+:1 l = [ l .. >..h ,J 2A.h y 0 This matrix can then be row reduced to solve for and :lJ: in terms of For the purpose of developing a nmnerical method only y: is needed. It is n = Jf (A.h). = 2(:3 + A.h) 1 A.2h2 4A.h + 6 (:3. 4) The amplification factor, A11 is the subdiagonal Pade approximant R1,2 (>..h) to f/'h. From the properties in equations 1.6 and 1. 7, we have that Therefore, in addition to Astability, DGM(1) exhibits stiff Astability, and has order of :3. 19
PAGE 28
The numerator of 1\J1(>..h) has a zero at >..h = :3, and is negative for >..h < :.L The denorninator is positive for all real >..h < 0. Thus JJ1 (A. h) < 0 for real negative >.. whenever the step size h > :3/ >... The sequence of nurnerical solutions will then oscillate between negative and positive values. Therefore, is not in general, rnonotone. But the rnethod will exhibit monotonicity for step size chosen small enough. For exarnple, with>..= :w, choosing h < 0.1 will result in this method being monotone. 3.3 DGM(2) For q = 2, the shape functions in the local coordinates are: The derivatives are: 1(0 dO = e + 1 = + The test and weight functions expressed in tenns of the shape functions are: ,11(t). = unl +. ,11n1/2 + un ,'J ., 1.'1+ 2.'1 3.'120
PAGE 29
The following integrals will be needed: <1>1 = 1 / 2 = 2/:3 = 1 / 6 = 4/1.5 <1>1 = 2/:3 = 0 = 2 j: 3 = 16/1.5 = 2 /15 <1>1 = 1 /15 The first integral in equation : _L2 is evaluated u sing t.he above results as follows: [11 jl (;r;.l' nl + ;r;.l' n1/2 + ;r;.l' n) (;r;. nl + ;r;. nl/2 + ;r;. ',n) zs:. L'1'3L_ c._, 1 (J l (
PAGE 30
>.h [J (2.1jn1 + .. lj n1/2 + . y" )dt ] 'Vn1 2 1 1 + 1 2 1 3 .. .., + 4_./jn1 + 2 _./jn1/2 'Vn1 2 1 a + 1 a 1a + + >.h ( 2 .ljn1 + + 2 y )vn 1/2 2 1 o + 1o 1 o + + 2_/jn1/ 2 + 4 y ) v 2 1 o + 1 o 1a So from the variational equation ;u or :3.2 and the fact that v"112 and are arbitrary we have the three equations : ( 1 2A.h) n1 ( 2 A.h, n1/ 2 ( 1 A.h) n _, n1 2 15 Y+ + 3 15 )y + 6 + :.HJ y y ( 2 )._J.t). n 1 8A.. h n 1 ; 2 ( 2 A.h) n lj y + lj 0 ., 1:' + 1:' 'J 1:' '0 J J 0 J ( 1 A.h n 1 ( 2 A.h. n1 ; 2 ( 1 2A.h. n 6 + : .HJ )y+ + 3 15 )y + 2 15 )y = 0 In rnatrix fonn, this b e cor n e s [ 15 <1.\h. 2(.) :2A.. h 10A.h 8A.h .s + A.h 202A.h 5 + A.h l 10A.h 1.54A.h 22 + [ .ljn1 l :IJn 1/2 yfl_
PAGE 31
Upon solving this systern, we obtain In this case, 1\J2(A.h) is the subdiagonal Pade approxirnant R2 3(A.h) to e>.h. Again frorn equations 1.6 and 1. 7, Therefore, D:MG(2) is stiff Astable and has order p = .5. The nmnerator of 1H2 ( A.h) has two cornplex root s and the denorninator has two cornplex root s and one positive rea l root. Then, since JJ2(0) = 1 > 0, JJ2(A.h) > 0 for all negative real A.h. Therefore, That is DGM( 2) i s rnorwtone. 3.4 DGM(3) For q = :3, the shape functions are
PAGE 32
9 3(0 =(( + 1)(:3( + 1)((1) 4(() = 1 1 6 (( + 1)(:3( + 1). The derivatives are And ' ( ) = 21e 18( 1 1 16 9(9.C. 2 2.C. ?). ;(()=' .., 16.., ) 9(9.C. 2 2.C. ?). ;(()= '.., 16 .., ) ' ( .. ) = 21e + 18( 1. 4 16 u(t). = un1 +. un2/3 + u n1/3 + u n ,'j .., 1.'1+ 2 .'1 3.'1 4.'1u'(t). = '. n1 n 2/ 3 + n1/3 +'. n ,'j .., 1.11+ 2.11 3Y 4.1124
PAGE 33
The following integrals will be needed: I1 1 <1>1 = 1 / 2 = 7/80 J_11 = 81/80 = 81/80 J_ 1 1 <1> 4 = 7/80 J_ 1 1 <1> 4 = 1 I 2 .L11<1>1 :3//70 = 27/280 = :3:3/280 Then, <1>1 = 57 /80 =57 /80 = :3/10 = 0 = :3/10 = 16/10.5 = 19/ 840 = :3/70 = 16/10.5 <1>1 = :3/10 = 0 = :3/10 = .57/80 = .57/80 1 <1>1 = :3:3/280 1 .1 m.2 z '27/ '''" 1 = .Jn = 27 j:35 11 (' 'lln1 + 'lln2/3 + 'lln1/3 + ' 11n 1 + ' u n )dt)t'n 4 1 '1+ 4 2 '1 4 3.'1 4 4'1.., 1 ( ;3 n1 81 n2/3 O n1/3. .57. n ) n1/3 + 1QY+ 80 :lJ + :lf :lf + 80:lf_ v 2.5
PAGE 34
Also, A.h 11 2 1 :IJ( = A.h 11 ( n1 +., n2/3 + n1/3 + n). 2 1Y+ 2Y 3 Y 4Y1 ).h [1'1 (2,1jnl + . 'ljn2/3 + 'ljn1/3 + 'lJn)dt)] vn1 2 1 1 + 1 2 1 3.. 1 4 . .., + [1 ( 'ljn1 + . 'ljn2/3 + + 'ljn )dt)] Vn1/3 2 1 3 1. + 3 2. 3. 3 4. .., ).h [1'1 ( 'ljn1 + . 'ljn2/3 + 'ljn1/3 + 2yn )eft)'] vn 2 '1 4 1.. + 4 2 . 4 3 . 4 . .., 16_.1jnl + .33 ,ljn2/3 _}! .ljn1/3.lj + ltl yn)vn1 2 lOa + 280 tO 840 + + >.h(33 'ljn1 + 27,1jn2/3 27 'ljn1/3 JLyn) '[ln2/3 2 280 + 35 280 70 + >.h(_'L.ljn1_ 27.1jn2/3+ 27.1jnl/3+ 33.1Jn)vnl/3 2 70 + 280 35 280 + ( ltl .ljn1 }! .ljn2/3 + .33 ljn1/3.lj + 16_ljn )vn 2 840 + 10 280 lOa In a sirnilar fashion as before, use the fact that v'f1 vn2 / 3 vnl/3 and 26
PAGE 35
v: are arbitrary, to obtain the following rnatrix equation r 840 128Ah ;.197 o04 + .504 + :_WAh 1197 99Ah 648Ah 1701 + 81Ah .504 + :_WAh .504 + : wAh 170 + Ah 648Ah 119799Ah 1680:/j:1 1 0 0 0 Upon solving this system we have for y: 14719Ah 1 504 +:wAh 119799Ah 840 128Ah ,ljn1 1 + Y. n1/ 3 y :IJ: For q = :3, A13(Ah) is the Pade approxirnant R3 4(Ah) to e>h therefore, IA13(Ah)l+ 0 as Re(Ah)+ oo. Thus i s stiff A.stable and has order p = 7. The denorninator of .A13(Ah) has only complex zeros therefore its sign does not change. The nmnerator has two complex zeros but one negative real zero, x1 = (.5.J65)113(.5.J6+.5)1 135, which is approxirnately .5.65. Thus, DGM(:.n is rnonotone for negative A only when the step size h < x1 / A. 27
PAGE 36
This is not a very strong restriction. For example, with ,\ = :w, h can be as large as 0.18. 28
PAGE 37
Chapter 4 Computational Results All calculations were done with 15 decirnal places of accuracy. For >.. = 1, tables 4.1, 4.2, and 4.4 show the cornputational results for various n. The last colun1n of each table gives the absolute value of the error divided by h raised to the power of the order of the rnethod. Note that these values are fairly constant. This verifies the results obtained earlier in this work for the order of DGM. Figures 4.1, 4.2, 4.:3 and 4.4 are loglog graphs of the absolute value of the error versus the nurnber of intervals, n. Since n = 1/h, the negative of the slope on each graph gives approxirnately the order for the method. Note that for the result is accurate to 15 decirnal places for n = lG. Therefore, since these calculations are accurate to 1.5 decirnal places, larger values of r1. do not increase the precision. 29
PAGE 38
In tables 4.5, 4.6 and 4. 7, the approxirnate value at each tirne step for four methods is shown. The last colun1n gives the true solution to four decirnal places. The four rnethods are: a method that is not Astable, Heun's rnethod, which is a twostage RungeKutta rnethod; the order of p = 2 trapezoidal method, which is Astable but is not stiff Astable; The DGM(O) method or Euler's rnet.hod with order only p = 1, which is both Astable and stiff Astable; and the DGM(2) with order p = .5 and both stability properties mentioned above. These rnethods were applied to the model equation ( 1.2), but with A = : _H). Figures 4.5, 4.6 and 4. 7 show these results graphically. Heun's rnethod gives absurd results, so it was not shown on a graph. Notice that the DGM(O) actually gives better results than the trapezoidal method, even though it is of lower order. Then, of course the cornes the closest to the true solution curve. The better results for the DGM's for this test case is important in systems of equations that have sorne eigenvalues with very large rnoduli. These rnethods will not restrict the step size as rnuch as other rnethods that are not stiff Astable.
PAGE 39
N I cornputed value of y(1) I ctTor C'f'r'O'f'/h 1 5 .OOOOOOOOOOOOOOOe0 1 1.321205588285577e O 1 1.32120558828titi77eO 1 2 4.444444444444444eO 1 7.656500327300209e02 1.5313000654600eO 1 4 4.096000000000000e O 1 4.172055882855769e02 1.668822353142307e O 1 8 3.897 443431289458eO 1 2.186490 195750346e 02 1. 7 49192156600277eO 1 16 3. 790853319179361eO 1 1.1205890H6493 81e02 1. 7929 12519439009eO 1 32 3. 73553861 1900618eO 1 5 .67H20318619 122e 0 3 1.81581 19582ltie01 61 3. 7073 19329009732eO 1 2 .855 191729530868e 03 1.82751H0689975tieO 1 128 3 .69311 8 10602 1922eO 1 1 A32369431049851e03 1.8334328717 43810eO 1 Table 4.1: DC:M(O) applied to rnodel equation with ,\ = 1 for various nmnber of intervals, N N I cornputed value of y(l) ctTO'f' cnorjh3 1 3 .636363636363636eO 1 1.2 13077535078688e03 1.2 13077535078688e 03 2 3 .67309458218ti 191eO 1 5 .69982952 89 32332e04 4 .559863623145 866e03 4 3 .67804395190eO 1 7.504598101665261e 05 4.80294278506 5767 e03 8 3.67869777 4589969e O 1 9 .663712445429073e06 4.947820772059686e03 16 3 .67 8782 140046069eO 1 1.227166835482 4 07e 06 5 .026 4 75358135940e03 32 3 .678792865263700eO 1 1.546450723588499 e07 5 .067 409731054 795e 03 61 3 .67879 1217611732eO 1 1.911026911 122838e08 5 .088285586680286e 03 128 3 .67879 138H0132 1eO 1 2 ..1313099311612ti2e09 5 .09882M8 175 1682e 03 Table 4.2: applied to rnodel equation with ,\ = 1 for various nmnber of intervals, N N I cornputed value of y(l) I cnor cnorjh5 1 3.6792e01 1.5087130H13n1e0 5 1.5087130HH322 1e05 2 3 .678809236H7EieO 1 1 ..182 173311886778e 06 1. H391 1598037691e05 4 3 .67879489111625tieO 1 4. 7940 18315656956e08 4.90907 4 755232723e05 8 3 .67879442698H61eO 1 1.5273037834440 12e09 5 .004669037589338e 0 5 16 3 .678794412196ti91e O 1 4.821676391486562ell 5 .055894143879414 e05 32 3.678794 4 11729567e01 1.514399716739945e12 5.081482231616974e05 61 3.67879Hll7H8 96e01 1.72955008 1903167 e1 5.07831573 1863281e05 128 3.67879Hll71H37e01 1.332267629550188e15 1.577636718750000e05 Table 4.:3: DGM(2) applied to rn odel equation with ,\ = 1 for various nmnber of intervals, N :_n
PAGE 40
N I cornputed value of y(1) I crTor' crrorjl/ 1 3 .67879203843Ei1He O 1 2 .3732792825H496e 07 2 .373279282541496e 07 2 3 .67879eO 1 1. 92713239810Ei lti8e09 2 A66729f.1602e 07 3 .678791111559969e01 1 .5115122718Ei8Ei18e11 2 .530578058212996e 07 8 3 .678791111713199eO 1 1.22lti 18e13 2 .568121999502182e 07 16 3 .678794411714H3e O 1 1.054 711873393899e 15 2 .831220626831055e 07 32 3.678794411714424e 01 1.11022302e 16 3.814697265625000e 06 61 3 .678794411714e O 1 1.11022302462Ei1Ei7e 16 4.882812500000000e 04 128 3 .678794 4 11714 127e O 1 3.330669073875 4 70e 16 1.875000000000000e O 1 Table 4.4: DG applied to rnodel equation with >.. = 1 for various nurnber of intervals, N * n Figure 4.1: lerrorl vs. nurnber of intervals, N, for DGM(O)
PAGE 41
102 103 104 105 s106 107 108 109 10 101 n Figure 4.2: lerrorl vs. nmnber of intervals, N, for step DGM(O) DGM(2) trapezoidal Heun's y(t") = cAt 0 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 1 1.1765c01 3.8U8c02 t i.7895c01 2 .1625 c + 01 5.5308c0 4 2 1.3841c02 l.Ei014c 03 3.3518c 01 4.6764 c + 02 3.0590c07 3 1.6283c03 5.8177c 05 1.9 105c01 1.0113 c + 04 1.6919c1 0 4 1.9157c04 2.2Ei43c 06 1.12 3 5c 0 1 2 .1869 c + OEi 9.3576c14 Table 4 .. 5: Cornparison of methods for A = :_H) and n = 4
PAGE 42
* n Figure 4.:3: lerrorl vs. nmnber of intervals, N, for DGM(2) step DGM(O) DGM(2) trapezoidal Heun's y(t ) = c>..t 0 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 1 2.1053c01 3.2561c 02 3.0,135c 01 4.2812 c + 00 2.3518c02 2 4.4321c02 1.0602c 03 9.2628c 02 1.8329 c + 01 5.5308c04 3 9.3308c03 3Ati22c05 2.8191c02 7.8H1c + 01 1.3007c05 1 1.96Hc03 1.12Uc 06 8.Ei799c03 3.3596c + 02 3.0590c07 i ) 4.1355c04 3.6601c 08 2.6113c 03 1.4383 c + 03 7.1941c09 6 8.7064c05 1.1918c 09 7.9H3c04 6.1578 c + 03 1.6919c10 7 1.8329c05 3.880tic 11 2A188c 04 2.6363 c + 04 3.9790c12 8 3.8588c06 1.2635c 12 7.3614c 05 1.1287 c + 05 9.3576c14 Table 4.6: Cornparison of methods for A.= :.H) and n = 8
PAGE 43
* .:::* e * 10 10 1 102 103 n Figure 4.4: lerTorl vs. nurnber of intervals, N, for step DGM(O) DGM(2) trapezoidal Ileun's y(tn) = CAt 0 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 00 l.OOOOc + 0 1 3.4783c01 1.5415c01 3.2258c02 8.8281c 01 1.5335c 01 2 1.2098c01 2.3761c02 1.0 106c03 7.7936c 01 2.3Ei18c02 3 4.2081c02 3.6626c03 3.3567c05 6.8803c 01 3.6066c 03 16 .5898c 08 1.0160c13 1.37Hc2 1.36llc01 9.3ti76c U Table 4.7: Cornparison of rnethods for).= and N = 16
PAGE 44
0.9 0.8 0.7 0.6 >0.5 0.4 0.3 0.2 0.1 0.1 0.2 0.3 0.4 0.5 t 0.6 0.7 0.8 0.9 Figure 4 .. 5: DGM(O) for r nodel equation with ). = : _H) and N = 4
PAGE 45
0.9 0.8 0.7 0.6 >0.5 0.4 0.3 0.2 0.1 0.1 0.2 0.3 0.4 0.5 t 0.6 0.7 0.8 0.9 Figure 4 .6: DGM(2) for r nodel equation with ). = : _H) and N = 4
PAGE 46
0.8 0.6 0.4 >0.2 0 0.2 0.4 0.6 L__ ___.L_ __ __J__'.*"''''''_J._L____J 0 0.1 0.2 0.3 0.4 0.5 t 0.6 0.7 0.8 0.9 Figure 4.7: Tra pezoidal rnethod for rnodel equation with,\= JV = 4
PAGE 47
Chapter 5 Conclusion In this work, we showed that the discontinuous Galerkin rnethods for ordinary differential equations are Astable. \Vith piecewise polynomial spaces of degree q = 0, 1, 2 and ;3, we showed that the order is p = 2k: + 1. Also, in these cases, the rnethods are stiff Astable. Table .5.1 gives a sununary of thesi s results. hod Order Astable Stiffly A stable Monotone DGl\1(0) 1 yes yes unconditionally DGl\1(1) ;3 yes yes for h < :3/A. DGl\1(2) .J yes yes unconditionally DGM(:.n 7 yes ye s for h < 5.6.5/ A. Table .5.1: Cornparison of discontinuous Galerkin rnet.hods
PAGE 48
Bibliography [ 1 ] 0. Axelsson. A class of Astable methods. BIT, [ 2 ] G. A. Baker, .Jr. and P. Graves. Morris. Pade Approximanls1 Part I: Basic Theory. Addison\Vesley, 1981. [:3] .J. C. Butcher. On A stable implicit RungeKutta rnethods. BIT, 17::375:378, [4] .J. C. Butcher. Irnplicit Runge Kutta processes. J11alhenwlics of Comp'U.lations, 18:.5064, 1978. [.5] .J. C. Butcher. The Numerical Analysis of Ordinary Equations .John vViley and Sons, 1987. [G] G. Dahlquist. A special stability problen1 for linear multistep methods. BIT, :3:274:3, 196:3. [7] Delfour, vV. Hager, and F. Trochu. Discontinuous galerkin methods for ordinary differential equations. J11alhematics of Compu.lalion, :.W(154):4.5.54 7:3, 1981. [8] B. L. Ehle. High order Astable rnethods of the numerical solution of systen1s of differential equations. BIT, 8:276278, 1968. [9] B. L. Ehle. A.stable rnethods and Pade approxirnations to the exponen tial. 8IAA1 Journal of J11alhematical Analysis, pages 671680, 197:3. [10] C. \V. Gear. Nu.merical Initial VaZ.u.e Problems in Ordinary DUJ'erential Eq'U.alions. PrenticeHall, 1971. 40
PAGE 49
[11] C. Johnson. Numerical Sol'Uhon of Partial Dzfj'erential Equations by the Finite Element 111elhod. Carnbridge, 1987. [12] P. Lesaint and P. A. Raviart. On a Finite Element J11elhod for Solving the Neutron Transport Equation1 in klathematical Aspects of Finite Elements in Partial DZfJ'erential Equations, pages 12:. L Acadernic Press, [1:3] H. A. vVatts and L. F. Sharnpine. Astable block implicit onestep methods. BIT, 1972. 41
