Citation
Discontinuous Galerkin methods for ordinary differential equations

Material Information

Title:
Discontinuous Galerkin methods for ordinary differential equations
Creator:
Bauer, Russell E
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
viii, 41 leaves : illustrations ; 29 cm

Subjects

Subjects / Keywords:
Galerkin methods ( lcsh )
Differential equations ( lcsh )
Differential equations ( fast )
Galerkin methods ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 40-41).
Thesis:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Applied Mathematics
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Russell E. Bauer.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
34018869 ( OCLC )
ocm34018869
Classification:
LD1190.L622 1995m .B38 ( lcc )

Downloads

This item has the following downloads:


Full Text
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 A-stability. 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 A-stable 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 A-Stable Methods...................................... 3
1.2 Stiff A-Stable Methods................................ 5
2 Discontinuous Galerkin Method 8
2.1 A-stability 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, (1-1)
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 A-stability. If the region of absolute sta-
bility includes all values of Ah, where Re(Xh) < 0, then the method is said
to be A-stable. In other words, a method is A-stable 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 h-1.
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 Rn-i,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 A-Stable Methods
A general linear fc-step method for computing the approximate numerical
solution is defined by the formula
OfcJ/n+fc + ak-iyn+k 1 + ... + *0yn = h(/3kfn+k + + Po fn)i (1-8)
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 fc-step method is not A-stable
and that the highest order of the global truncation error can not exceed 2 for
an A-stable implicit linear multistep method. Therefore, there are no high
order A-stable linear multistep methods. The smallest truncation error for
order p = 2 A-stable 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
multi-step methods with Richardson extrapolation can be used to obtain an
order of p = 4 A-stable procedure [6].
The explicit Runge-Kutta methods are not A-stable [10], but there are
implicit Runge-Kutta methods that are A-stable [4, 8]. An r-stage implicit
Runge-Kutta method has the form
kq = hf(yn + fiqjkj), q = 1,2,..., r j = 1,..., q 1
!/=!/ + T.H %
3=1
The implicit r-stage Runge-Kutta 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
l-R^rl < 1, if Re(Xh) < 0.
Hence, the methods are A-stable. Also, r-stage Runge-Kutta methods of
order p = 2r 1 and 2r 2 are A-stable [3, 8].
Also, there is another class of A-stable 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 A-stability see [1, 3, 4, 5, 8,
9, 10, 13]
1.2 Stiff ^4-Stable Methods
For a certain class of differential equations called stiff equations, another
desirable property of numerical methods is stiff A-stability. 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 A-stable [1, 10].
The r-stage implicit Runge-Kutta 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 9-stage Runge-Kutta methods with order of 2r 1
that axe stiff A-stable [10]. In this case
yn+1 = Rr-i,r(\h) y11.
Since
Rr-if 0 as Re(z) > 00,
this class of methods is stiff A-stable.
Also, the extension of Taylors series mentioned above can be made to
be stiff A-stable [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 A-stable
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 A-stability, 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
tn-l tn
interval In
Figure 2.1: Discretization strategy of the discontinuous Galerkin method
where
Bn(v,y)= j {yhV-XyhV^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 A-stable 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 t71-1 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 <3-2)
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 = y-1v"
=S-(-W + S,>* = !,"-V.
Since vn is arbitrary, we have
yn = M0yn~l =
1
.71 1
1-A 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 A-stable, but the second and third
expressions above show that the DGM(O) is stiff A-stable 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 f-i $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 £ vl-1
+ (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^K-1
+[(-1/2 Xh/fyyl-1 + (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
' 3-2A h 3-A h ' ' vV ' ' Gyr1 '
-3-A h 3-2A h y- 0
This matrix can then be row reduced to solve for y-1 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 A-stability, DGM(l) exhibits stiff A-stability, 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 = §(0-0
*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 + ^y71-172 + $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 £y-K-1
+ (-|y+-1 + o t/n-1/2 +
+ (iy+1-fyn"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$2yn-1/2 + $i*3|£)<*f] ur1
+ f [/X1($1$2^-1 + *lyn-1/2 + *2*3y2K] u"-1/2
+ f [/ic^i^r1 + ^2^3yn~1/2 + $hl)dt] vZ.
+
+
t (ir1 + 4n_1/a-i-)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 A-stable 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)(3f-l)-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(() = + 2Vn-2'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 + *'2yn-2/3 + *'yn-1/3
t*!^-1 + *2^n"2/3 + *3Un_1/3 + *4t£.)#
= £(*a$iy-1 + *! *lj/n-2/3 + *i*^2/n-1/3 + *i*l2/!!)#<-1
+ + *2 *lj/n-2/3 + *2*lj/n-1/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 + |y-1,3s ^y>-2/3
+ (^i-1 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$3yn-i/3+$i$4^k)]
+ f [A^SitfT1 + $&n-213 + ^2$3y"-1/3 + $2$42/!lMO] n_2/3
= f [/i^itfT1 + $3$2yn-2/3 + $iy"-1/3 + $3$4 = ¥ [flA^Qiyr1 + $4$2yn~2/3 + $4$3yn-1/3 + ^yH)#)] vi
= fiwsvr1 + mvn~m ~ r0yn-1/3y + &-K"1
+ f(~w0yV +1 yn~2/3 M>yn'1/3 ^-K2/3
+ ¥(-hyl~x ~ myn~2'z +1yn~113 + i|y-X"1/3
+ f (ffivr1 IX2/3 + w0yn~1/3y + w5y-)v-
In a similar fashion as before, use the fact that u-1, 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 147-19AA ' r vV 1
-1197 -99Ah 648Ah 170 + Xh -504 + 36A h y7l2/3
504 + 36A h -1701 +8lXh -648A h 1197 99A h y71-1/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) n-1
\4hi-16\3h3 + m\2h2-m\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 A-stable 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 (5-S/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 log-log 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 A-stable, Heuns method,
which is a two-stage Runge-Kutta method; the order of p = 2 trapezoidal
method, which is A-stable but is not stiff A-stable; The DGM(O) method
or Eulers method with order only p = 1, which is both ^4-stable and stiff
A-stable; 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 A-stable.
30


N computed value of 2/(1) error error/h
1 5.000000000000000e-01 1.321205588285577e-01 1.321205588285577e-01
2 4.444444444444444e-01 7.656500327300209e-02 1.531300065460042e-01
4 4.096000000000000e-01 4.172055882855769e-02 1.668822353142307e-01
8 3.897443431289458e-01 2.186490195750346e-02 1.749192156600277e-01
16 3.790853319179361e-01 1.120589074649381e-02 1.792942519439009e-01
32 3.735538614900618e-01 5.674420318619422e-03 1.815814501958215e-01
64 3.707349329009732e-01 2.855491729530868e-03 1.827514706899755e-01
128 3.693118106024922e-01 1.432369431049851e-03 1.833432871743810e-01
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.636363636363636e-01 -4.243077535078688e-03 -4.243077535078688e-03
2 3.673094582185491e-01 -5.699829528932332e-04 -4.559863623145866e-03
4 3.678043951904257e-01 -7.504598101665261e-05 -4.802942785065767e-03
8 3.678697774589969e-01 -9.663712445429073e-06 -4.947820772059686e-03
16 3.678782140046069e-01 -1.227166835482407e-06 -5.026475358135940e-03
32 3.678792865263700e-01 -1.546450723588499e-07 -5.067409731054795e-03
64 3.678794217611732e-01 -1.941026911422838e-08 -5.088285586680286e-03
128 3.678794387401324e-01 -2.431309931161252e-09 -5.098826484754682e-03
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.679245283018868e-01 4.508713044443224e-05 4.508713044443224e-05
2 3.678809236447542e-01 1.482473311886778e-06 4.743914598037691e-05
4 3.678794891116255e-01 4.794018315656956e-08 4.909074755232723e-05
8 3.678794426987461e-01 1.527303783444012e-09 5.004669037589338e-05
16 3.678794412196591e-01 4.821676391486562e-ll 5.055894143879414e-05
32 3.678794411729567e-01 1.514399716739945e-12 5.081482231616974e-05
64 3.678794411714896e-01 4.729550084903167e-14 5.078315734863281e-05
128 3.678794411714437e-01 1.332267629550188e-15 4.577636718750000e-05
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.678792038435141e-01 -2.373279282541496e-07 -2.373279282541496e-07
2 3.678794392443099e-01 -1.927132398105158e-09 -2.466729469574602e-07
4 3.678794411559969e-01 -1.544542271858518e-ll -2.530578058212996e-07
8 3.678794411713199e-01 -1.224575996161548e-13 -2.568121999502182e-07
16 3.678794411714413e-01 -1.054711873393899e-15 -2.831220626831055e-07
32 3.678794411714424e-01 1.110223024625157e-16 3.814697265625000e-06
64 3.678794411714422e-01 -1.110223024625157e-16 -4.882812500000000e-04
128 3.678794411714427e-01 3.330669073875470e-16 1.875000000000000e-01
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 A-stable. 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 A-stable. Table 5.1 gives a summary of
thesis results.
Method Order A-stable Stiffly A-stable 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 A-stable methods. BIT, 9:185-199, 1969.
[2] G. A. Baker, Jr. and P. Graves-Morris. Pade Approximants, Part I:
Basic Theory. Addison-Wesley, 1981.
[3] J. C. Butcher. On A-stable implicit Runge-Kutta methods. BIT, 17:375-
378, 1977.
[4] J. C. Butcher. Implicit Runge-Kutta processes. Mathematics of Com-
putations, 18:50-64, 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:27-43, 1963.
[7] M. Delfour, W. Hager, and F. Trochu. Discontinuous galerkin meth-
ods for ordinary differential equations. Mathematics of Computation,
36(154):455-473, 1981.
[8] B. L. Ehle. High order A-stable methods of the numerical solution of
systems of differential equations. BIT, 8:276-278, 1968.
[9] B. L. Ehle. A-stable methods and Pade approximations to the exponen-
tial. SIAM Journal of Mathematical Analysis, pages 671-680, 1973.
[10] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential
Equations. Prentice-Hall, 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 89-123. Academic Press,
1974.
[13] H. A. Watts and L. F. Shampine. A-stable block implicit one-step meth-
ods. BIT, 12:252-266, 1972.
41


Full Text

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 A-stability. 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 A-stable 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 A-Stable Methods :3 1.2 Stiff A-Stable Methods .J 2 Discontinuous Galerkin Method 8 2.1 A-stability 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 guid-ance 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 A-stability. If the region of absolute stability includes all values of A.h, where Re(A.h) < 0, then the method is said to be A-stable. In other words, a rnethod is A-stable 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 h-1 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 Rn-I,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 A-Stable 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 multi-step methods with Richardson extrapolation can be used to obtain an order of p = 4 A.-stable procedure [ 6 ] The explicit Runge-Kutta methods are not A.-stable [10], but there are irnplicit Runge-Kutta rnethods that are A.-stable [ 4 8 ] An r-stage implicit Runge-Kutta 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 r-stage Runge-Kutta 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 A-stable. Also, r-stage Runge-Kutta methods of order p = 2r-1 and 2r2 are A-stable [:3, 8 ] Also, there is another class of A-stable 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 A-stability see [1, :3, 4, .5, 8, 1.2 Stiff A-Stable Methods For a certain class of differential equations called stiff equations, another desirable property of nmnerical methods is stiff A-stability. 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 A-stable [ 1, 10] The r-stage irnplicit Runge-Kut 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 q-stage Runge-Kutta methods with order of 2r 1 that are stiff A-stable [10] In this case n+l R (\/. ) n Y = r-l,r A t Y Since Rr-l,r--+ 0 a.s Re(z)--+ -oo, this class of rnethods is stiff A-stable. Also, the extension of Taylor's series mentioned above can be rnade to be stiff A-stable [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 A-stable

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 A-stability for any desired order of global truncation error. \Ve will also show that this method exhibits the property of stiff A-stability, 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) = ih-v 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 ( n-1 n-1 n-1] L Jr :zrv-AlJV )d + Y-f)v-f= 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 A-stability 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( n-12 n-1 n-1 1( n-12 1( n-12 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 n-11 .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 = ( tn-l, 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( + tn-1 + 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 ( n-1 n-1 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 n-1 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 A-stable, but the second and third expressions above show that the DGM(O) i s s tiff A-stable 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) n-1 + ' (t) n .lJ ., 1 ., .lJ+ 2 ., .lJ_ ( "" ( n-1 "" ( 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. n-1 n (m. n-1 m. n l 2 _1 iJ ()v ()<{ = 2 _1 + <[>2;//_:) + A.h[ ( n -ljl m.2z n j 1 A'. z n-1] = --;-2 Y+ "'1<( + Y-1 -1 -1 A.h [ ( n-1 jl A'. A'. z n jl m.2 z ') n l +--;-2 Y + + y_ v_ -1 -1 So the variational equation 1) is: + n-1 n-1 n-1 n-1 Y+ v+ = y_ L'+ 18 n-1 n-1 = 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 A-stability, DGM(1) exhibits stiff A-stability, 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). = un-l +. ,11n-1/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' n-l + ;r;.l' n-1/2 + ;r;.l' n) (;r;. n-l + ;r;. n-l/2 + ;r;. ',n) zs:. L-'1'3L-_ c._, -1 (J l (
PAGE 30

>.h [J (2.1jn-1 + .. lj n-1/2 + . y" )dt ] 'Vn-1 2 1 1 + 1 2 1 3 .. .., + 4_./jn-1 + 2 _./jn-1/2 'Vn-1 2 1 a + 1 a 1a -+ + >.h ( 2 .ljn-1 + + 2 y )vn -1/2 2 1 o + 1o 1 o -+ + 2_/jn-1/ 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) n-1 ( 2 A.h, n-1/ 2 ( 1 A.h) n _, n-1 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. n-1 ; 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 -10-A.h 8A.h .s + A.h -20-2A.h 5 + A.h l 10-A.h 1.54A.h 22 + [ .ljn-1 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 A-stable 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). = un-1 +. un-2/3 + u n-1/3 + u n ,'j .., 1.'1+ 2 .'1 3.'1 4.'1u'(t). = '. n-1 n -2/ 3 + n-1/3 +'. n ,'j .., 1.11+ 2.11 3Y 4.1124

PAGE 33

The following integrals will be needed: I--1 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 (' 'lln--1 + 'lln--2/3 + 'lln--1/3 + ' 11n -1 + ' u n )dt)t'n 4 1 '1+ 4 2 '1 4 3.'1 4 4'1--.., ----1 ( ;3 n--1 81 n--2/3 O n--1/3. .57. n ) n--1/3 + 1QY+ 80 :lJ + :lf :lf + 80:lf_ v 2.5

PAGE 34

Also, A.h 11 2 -1 :IJ( = A.h 11 ( n-1 +., n-2/3 + n-1/3 + n). 2 1Y+ 2Y 3 Y 4Y--1 ).h [1'1 (2,1jn-l + . 'ljn-2/3 + 'ljn-1/3 + 'lJn)dt)] vn-1 2 -1 1 + 1 2 1 3.. 1 4 . .., + [1 ( 'ljn-1 + . 'ljn-2/3 + + 'ljn )dt)] Vn-1/3 2 -1 3 1. + 3 2. 3. 3 4. .., ).h [1'1 ( 'ljn-1 + . 'ljn-2/3 + 'ljn-1/3 + 2yn )eft)'] vn 2 '-1 4 1.. + 4 2 . 4 3 . 4 . .., 16_.1jn-l + .33 ,ljn-2/3 _}! .ljn-1/3.lj + ltl yn)vn-1 2 lOa + 280 tO 840 -+ + >.h(33 'ljn-1 + 27,1jn-2/3 27 'ljn-1/3 JLyn) '[ln-2/3 2 280 + 35 280 70 + >.h(_--'L.ljn-1_ 27.1jn-2/3+ 27.1jn-l/3+ 33.1Jn)vn-l/3 2 70 + 280 35 280 -+ ( ltl .ljn-1 }! .ljn-2/3 + .33 ljn-1/3.lj + 16_ljn )vn 2 840 + 10 280 lOa --In a sirnilar fashion as before, use the fact that v'f--1 vn-2 / 3 vn-l/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 1197-99Ah 840 -128Ah ,ljn-1 1 + Y. n-1/ 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.J6-5)113-(.5.J6+.5)1 13-5, 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 log-log 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 A-stable, Heun's rnethod, which is a two-stage Runge-Kutta rnethod; the order of p = 2 trapezoidal method, which is A-stable but is not stiff A-stable; The DGM(O) method or Euler's rnet.hod with order only p = 1, which is both A-stable and stiff A-stable; 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 A-stable.

PAGE 39

N I cornputed value of y(1) I ctTor C'f'r'O'f'/h 1 5 .OOOOOOOOOOOOOOOe-0 1 1.321205588285577e O 1 1.32120558828titi77e-O 1 2 4.444444444444444e-O 1 7.656500327300209e-02 1.5313000654600e-O 1 4 4.096000000000000e O 1 4.172055882855769e-02 1.668822353142307e O 1 8 3.897 443431289458e-O 1 2.186490 195750346e -02 1. 7 49192156600277e-O 1 16 3. 790853319179361e-O 1 1.1205890H6493 81e-02 1. 7929 12519439009e-O 1 32 3. 73553861 1900618e-O 1 5 .67H20318619 122e -0 3 1.81581 19582ltie-01 6-1 3. 7073 19329009732e-O 1 2 .855 191729530868e -03 1.82751H0689975tie-O 1 128 3 .69311 8 10602 1922e-O 1 1 A32369431049851e-03 1.8334328717 43810e-O 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 .636363636363636e-O 1 1.2 13077535078688e-03 -1.2 13077535078688e -03 2 3 .67309458218ti 191e-O 1 -5 .69982952 89 32332e-04 -4 .559863623145 866e-03 4 3 .67804395190e-O 1 7.504598101665261e -05 -4.80294278506 5767 e-03 8 3.67869777 4589969e O 1 -9 .663712445429073e-06 -4.947820772059686e-03 16 3 .67 8782 140046069e-O 1 1.227166835482 4 07e 06 -5 .026 4 75358135940e-03 32 3 .678792865263700e-O 1 1.546450723588499 e-07 -5 .067 409731054 795e -03 6-1 3 .67879 1217611732e-O 1 -1.9-11026911 122838e-08 5 .088285586680286e -03 128 3 .67879 138H0132 1e-O 1 -2 ..1313099311612ti2e-09 -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.6792e-01 1.5087130H13n1e-0 5 1.5087130HH322 1e-05 2 3 .678809236H7Eie-O 1 1 ..182 173311886778e -06 1. H391 1598037691e-05 4 3 .67879489111625tie-O 1 4. 7940 18315656956e-08 4.90907 4 755232723e-05 8 3 .67879442698H61e-O 1 1.5273037834440 12e-09 5 .004669037589338e -0 5 16 3 .678794412196ti91e O 1 4.821676391486562e-ll 5 .055894143879414 e-05 32 3.678794 4 11729567e-01 1.514399716739945e-12 5.081482231616974e-05 6-1 3.67879Hll7H8 96e-01 1.72955008 1903167 e-1 5.07831573 1863281e-05 128 3.67879Hll71H37e-01 1.332267629550188e-15 1.577636718750000e-05 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 .67879e-O 1 -1. 92713239810Ei lti8e-09 2 A66729f.1602e 07 3 .67879-1-111559969e-01 -1 .5-1-15122718Ei8Ei18e-11 2 .530578058212996e 07 8 3 .67879-1-111713199e-O 1 -1.22lti 18e-13 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 6-1 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

10-2 10-3 104 10-5 -s10-6 10-7 108 10-9 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.7895c-01 2 .1625 c + 01 5.5308c-0 4 2 1.3841c02 l.Ei014c 03 3.3518c 01 4.6764 c + 02 3.0590c-07 3 1.6283c03 5.8177c 05 1.9 105c-01 1.0113 c + 04 1.6919c1 0 4 1.9157c04 2.2Ei43c 06 1.12 3 5c 0 1 2 .1869 c + OEi 9.3576c-14 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 3Ati22c-05 -2.8191c02 7.8H1c + 01 1.3007c-05 1 1.96Hc-03 1.12-Uc 06 8.Ei799c-03 3.3596c + 02 3.0590c-07 i ) 4.1355c04 3.6601c -08 2.6113c 03 1.4383 c + 03 7.1941c09 6 8.7064c05 1.1918c -09 7.9H3c-04 6.1578 c + 03 1.6919c10 7 1.8329c05 3.880tic 11 2A188c -04 2.6363 c + 04 3.9790c-12 8 3.8588c06 1.2635c -12 7.3614c 05 1.1287 c + 05 9.3576c-14 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.2098c-01 2.3761c02 1.0 106c-03 7.7936c 01 2.3Ei18c02 3 4.2081c-02 3.6626c03 3.3567c05 6.8803c -01 3.6066c 03 16 .5898c 08 1.0160c13 1.37Hc-2 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 A-stable. \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 A-stable. Table .5.1 gives a sununary of thesi s results. hod Order A-stable 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 A-stable 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 Runge-Kutta rnethods. BIT, 17::375-:378, [4] .J. C. Butcher. Irnplicit Runge Kutta processes. J11alhenwlics of Comp'U.lations, 18:.50-64, 1978. [.5] .J. C. Butcher. The Numerical Analysis of Ordinary Eq-uations .John vViley and Sons, 1987. [G] G. Dahlquist. A special stability problen1 for linear multistep methods. BIT, :3:27-4:3, 196:3. [7] Delfour, vV. Hager, and F. Trochu. Discontinuous galerkin methods for ordinary differential equations. J11alhematics of Comp-u.lalion, :.W(154):4.5.5-4 7:3, 1981. [8] B. L. Ehle. High order A-stable rnethods of the numerical solution of systen1s of differential equations. BIT, 8:276-278, 1968. [9] B. L. Ehle. A.-stable rnethods and Pade approxirnations to the exponen tial. 8IAA1 Jo-urnal of J11alhematical Analysis, pages 671-680, 197:3. [10] C. \V. Gear. Nu.merical Initial VaZ.u.e Problems in Ordinary DUJ'erential Eq'U.alions. Prentice-Hall, 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. A-stable block implicit one-step methods. BIT, 1972. 41