- Email: [email protected]
Continuous and discontinuous finite element methods for burgers' equation
COMPUTER METHODS IN APPLIED @ NORTH-HOLLAND PUBLISHING
MECHANICS COMPANY
AND ENGINEERING
25 (1981) 65-84
CONTINUOUS AND DISCONTINUOUS FINITE ELEMENT METHODS FOR BURGERS’ EQUATION*
Paul ARMINJON Dipartement
and Claude BEAUCHAMP
de Mathe’matiques et de Statistique, Universite’ de Monkal, Canada
Revised
Received 12 December 1979 manuscript received 11 February
B.P. 6128, Montre’al H3C 3J7,
1980
We apply two discontinuous finite element methods to the inviscid Burgers’ equation and to the full equation with viscosity. In both cases we compare with a continuous space-time finite element method previously studied. For v = 0 discontinuous methods give better results, while the reverse prevails for the viscous equation.
1. Introduction
In previous reports [l], [2], a space-time continuous finite element procedure derived from Bonnerot and Jamet [4] was applied to the numerical solution of Burgers’ equation in one space dimension: 2
O
g+ug=v$$,
with initial and boundary
t>O
(la)
u(O,t)=u(l,t)=O,
(lb)
conditions and
u (x, 0) = uo(x) and in two space dimensions: vdu, v, + uv, + vu, = vdv, ut + uux + vu, =
(2)
with corresponding initial and boundary conditions; here A stands for the 2-dimensional Laplace operator. The solutions were regular, and the method had an (observed) order of accuracy equal to 2. In the present paper we shall apply two finite element methods of the discontinuous type to the inviscid Burgers’ equation in one space dimension (3)
ur + uu, = 0 with both regular and irregular initial functions * This research has been supported the province of Quebec.
by the National
t&(x).
Research
Council
of Canada
and the Ministry
of Education
of
66
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
We shall also apply to (3) the space-time continuous finite element method described in [l], for comparison. Finally we shall apply one of the discontinuous finite element methods presented here to the full equation (1) with viscosity, and compare the results with those obtained in [l]. Discontinuous finite elements were introduced earlier for the neutron transport equation by Reed and Hill [15], Lesaint and Raviart [ll]; for problems in nonlinear elasticity by Wellford and Oden [16], [17], [18]; for the equations of hydrodynamics by Fortin [8], [9]; and for ordinary differential equations by Rachford and Wheeler [14], Delfour, Hager and Trochu [6], [7], and others. Section 2 gives some properties of the inviscid Burgers’ equation which are of concern for our study, sections 3 and 4 describe the two discontinuous finite element methods considered here; these methods are called method F and method L-R here since they are obtained from two methods for ordinary differential equations by Fortin and Lesaint-Raviart, respectively. Section 5 describes some of the numerical experiments which have been performed. It is well known that solutions of (3) may present discontinuities even for regular initial data, while solutions of (1) tend to be more regular than the corresponding initial data, due to the smoothing effect of viscosity. The present paper aims at showing the advantage of discontinuous finite elements to eventually get a better fit of possible discontinuities or singularities such as those currently met in gas dynamics and shock tubes, while preparing for the use of these methods to study the numerical solution of the system of equations governing the atmospheric flow in numerical weather forecasting, in a manner designed to give a good account of the irregularities (“quasi-discontinuities”) arising in the vertical coordinate, particularly for the lowest kilometer of the atmosphere (boundary layer effect). Our results indicate that (i) for the inviscid equation (3) the discontinuous finite element methods considered here are more efficient and accurate than the continuous space-time finite element method of reference [l] (which will be called method C here). We also obtained good results for a discontinuous initial function. (ii) for the full Burgers’ equation with viscosity, the discontinuous methods F and L-R are very clearly inferior to the continuous method C.
2. Some properties of the inviscid Burgers’ equation We consider the initial-boundary ut + uu, = 0,
t >o,
value problem x > 0,
(3)
u (x70) = uo(x ),
(3a)
u(0, t)=
W)
0
and assume that this problem has a regular solution u(x, t). The characteristic curves of equation (3) are defined as the integral curves of the ordinary differential equation dx - = u(x, t).
(4)
dt
On such a characteristic
curve we have x = x(t), and the solution u(x, t) of problem
(3) is a
P. Anninjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
67
function of t only, for which du au au dx --+--P-+u~~O, dt- at ax dt
au
by (3). The solution u(x, t) is thus a constant along each characteristic curve. But since u = a constant on each such integral curve of dx/dt = u(x, t), we have dx/dt = u (x0, to) = uO(xO)so that the characteristic curves are the lines x(t) = Uo(z$ + x0,
(5)
or t = [uo(xO)]-‘(x - x0). If we consider for example the case X
z&J(x)=
l-x i 0
x 5 l/2, 1/25x 5 1, elsewhere,
0 s
(6)
we have the following situation (see fig. 1).
Fig. 1. Characteristics
for solution of eq. (3).
68
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
(i) If 0 < x0 -=cl/2, the characteristic lines t = and do not intersect each other. For all P = P = (x, t) with 0 < x < 1 and t > 1 there is only the value of the solution u(x, t) is completely constant value): u(x, t) = u(x0, 0) =
have slopes varying from +m to +2 (x, t) in the trapezoidal region OMBA or all one characteristic line passing through P, and determined by this line (along which u has a (x - x0)/x0
u&J = x0
(since 0 5 x0 I l/2) with t = (x - XO)/XO, and therefore ~0 = x/(1 + t) = u (x, t). (ii) For P = (x, t) in the triangular region MBC there is also exactly one characteristic going through P originating at (x0, 0) such that t = (x - x,)/x, or x0 = (x - t)/(l - t) with l/2 I x0 I 1, and therefore u(x, t) = 1 - x0 = (1 - x)/(1 - t). Notice that each such characteristic goes through B = (1,l). Summing up these two cases, we have &
for
O< x < 1
and
t >
2x - 1 (region OMBA or above),
U(X, t) = ’
(7)
and
t <
2x - 1 (triangular region MBC).
(iii) For x > 1 and t < 2x - 1 there is no characteristic going through (x, t) and originating at a point (x0, 0) where uO(xO)would be #O. In fact, the only characteristic going through such a point Q is the parallel to the t-axis, x = x Q. Since uO(xQ)= 0, we have u(Q) = 0. (iv) For x > 1 and t > 2x - 1 (triangular region T) there are two characteristic lines through P = (x, t) originating in (xl, 0) and (x2, 0), respectively, with 0 < x1 < l/2 and l/2 < x2 < 1. This corresponds to a jump discontinuity at time t with the values
Thus even for the regular initial function uO(x) defined by (6) Burgers’ inviscid equation can generate singular solutions. This example suggests the use of a discontinuous finite element method even if we limit ourselves, in this paper, to the values of u(x, t) for which 0 I x I 1, 0 5 t I 1, thus leaving the first singular point x = 1, t = 1 on the boundary of the domain of computation.
3. Method F. Discontinuous
finite elements derived from a classical variational
equivalent
We consider eq. (3) with initial and boundary conditions (3a)--(3b), and first perform a discretization in time by approximating the time derivative uI by a difference quotient, thus
P. Arminjon, C. Beau&,
Continuous and discontinuous finite element methods for Burgers’ equation
69
reducing the partial differential equation (3) to a system of ordinary differential equations to which we shall apply a method described by Fortin in [8], [9]. Let 0 = t” < t’ < . . . < t” . . . be an ordered sequence of points on the t-axis with At = t”+’- t” (not necessarily a constant). Setting U”(X) = U(X, t”), we use the difference quotient approximation
au
in
_z
-d-
(9)
At
at
to obtain the following approximate n
’
ad
in
system
~"-1
(10)
ax +dt= At’
uO(x> = uo(x), u”(0) = 0,
(11) 12
=0,1,2,...
(12)
To obtain our variational equivalent of this system, we take an arbitrary function !P = V(x) in CZ[O, l] = {P E C’[O, 11, ?-P(O)= V(1) = 0}, multiply (10) by !P and integrate by parts to obtain
-
=[,‘$+‘dx.
(13)
Let O=xO
=
OsjsM,
Sj,
F(Xj)=O,
[ Vi,(Xj)
OsjsM, =
OsjsM,
0,
$$ (Xj) =
6,,
OsjsM
(144
W)
70
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
?P$‘, restricted
for i = A4 we take the function i = 1,2, . . . . A4 - 1 by
(y>‘[
3-&(x
-xi4
Xi-1 I
2(?)1_3(?_&%)2+1,
Vi(x) =
0
q;(x) =
(x
0
0
-p
-xi)[q.f-y+
Xi,
xi 6 x 5 xi+1,
1
1 )
Xi-]
11,
Xi 5
I
X I Xi, X s
Xi+19
(W
elsewhere,
(x -;;-d2
I
X I
are given, for
elsewhere,
(x--;-d’ [”
?Py”(x>=
to [xIMPI,~~1. These functions
[x --y-1_
1],
elsewhere.
The trial functions will generally have two different values ui,_ = u(x;) and Ui.+= u(x+) at the knots x1, x2, . . . , xMpl. At the last knot X~ we shall only consider the left-hand value Us,_ as a degree of freedom, and at x0 we take uo.- = u o,+= uo(0), where the last term denotes the given value of the initial function uo(x) at x = 0. This gives a total of 2(M - l)+ 1 = 2M - 1 degrees of freedom, and the dimension of the space Ps[O, l] of test functions is also equal to 2M - 1. If we now take ly = !I”‘;, P = ?Pyi,in (13) and use appropriate quadrature formulae, we obtain a system of 2M - 1 nonlinear equations in the 2M - 1 unknown nodal values uy_ and u;, (1 I i I M - 1) and uL,_ (for n = 1,2,. . .). Since the functions we must integrate in (13) are piecewise polynomials, for each integral we choose a Newton-Cotes formula which gives its exact value; this yields the following system: u;,_=uI;,+=u(O,t”)=O
and
I
by(3b),
(IQ
=&[3u:_;‘++7u:~‘+7u:;‘+3u~~~_],
i=l,2
,...,
M-l,
I (u:_,,)‘+;u;_I.+u:--~(u:-)2+~(u:+)2+~u~+u~+~,-+(u:,,y] +s u:,l,-+3(u” 2 I,+ [
-u”_)-uy_l I,
I
(17’4
.+
, I =%Iu:;:._+q(u;;‘+u;_L)-u:_;I+ 1
i = 1,2, . . . , M - 1,
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
71
and for i = A4 we have
(17c) Eq. (17a) is obtained by taking P = !Pj in (13) while (17b) results from the choice ly = Vi. We note that if we assume that uy_ = uy+ (continuity of the trial functions), (17b) becomes ,&
I’ du;.
n u;+‘,l+2 ui+12U:-1+6
uy-‘_l .
u;+,2
UY-1
1 (18)
which is a consistent difference scheme for the equation uu, + uI = 0. If one wants to approximate a continuous solution, we shall see in section 5 that we get a better approximation of the nodal values ul = u”(xi) if we choose the approximate value
u; =&4;_+u;+).
(19)
System (17a)-(17c) consists of 2M - 1 nonlinear algebraic equations in 2M - 1 unknowns, which we solve with a simplified form of Newton’s method for systems. Specifically, to solve a system F(X) = 0, where F: Iw” +w” and X = (Xi)E Iw”, we start from a first approximation X0 = (Xp) and use the following iteration: X!+l= X: -A(Fi(Xk)/$$
(Xc)),
(20)
I
where A is a constant chosen to accelerate convergence; we found out that A = l/2 gives faster convergence than the usual value A = 1. When trying to approximate X = [u;,_, UT,+,. . . , uL-~,+, uh,-1’ we use the first approximation X0 = [uY,:*, u;;‘, . . . , uEll,+, ub:?]‘.
4. Method GR.
Discontinuous elements derived from a modified variational equivalent
We first perform the same time discretization as in method F and apply to the corresponding system of ordinary differential equations a modified variational formulation derived from the one appearing in Lesaint-Raviart [ll] by Delfour and Trochu (see [6], [7]). This leads to a family of numerical schemes depending on some real parameters ai. For Burgers’ equation
72
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
we use the values ai = 0 for all concerned value problem
g 1
=f(x, 4
m
u(0) = u”
indices i, or ai = 1 for all
i.
We consider the initial
x E (0, T), given
and a partition 0 = x0 < x1 < . . . < x,+, = T of the interval [0, T]. We set Ji = [xi-l, xi] and write Ui(X) for the restriction of the unknown function u = u(n) to the interval Ji (i = 1,2, . . . , M). We do not assume continuity of u at the nodes, so that Ui(Xi) may differ from Ui+l(xi). It can then be shown (see Delfour and Trochu [6]) that problem (P) is equivalent to the M + 1 variational problems i
(
u. + (1 - ao)[ul(0)-
-
(1 -(Yi)fJi(Xi)[&+l(Xi)-
uo] = u”
given, i.e.
ui(xi)l
=
O,
aOuo + (1 - cyo)ul(0) =
lsi
d’,
(224
Wb)
for all functions Ui
EHl(Ji)={EL2(Ji):$EL2(Jj)}. U
The numbers uo, Ui(Xi-I), Ui(Xi) (i = 1,2, . . . , M) and uM+, are the unknowns; u. may be considered as the left-hand value of our approximate solution u at x = 0, and u~+~ as the right-hand value at x = xM. At the node xi the left-hand value is Ui(Xi)and the right-hand value is Ui+l(Xi) (see fig. 2). Observe that if we assume the solution u to be continuous on [0, T], the last two terms in (22b) vanish, and the equivalence of (P) with (22) is obvious. But the interest of the method lies in its ability to handle discontinuities at the nodes; terms like Ui+I(Xi)- Ui(Xi)represent the jump at Xi. Now to find an approximation of the solution of (22) we seek functions u! E P“(A) C H’(4) for which (22) is satisfied for all ViE P”(A). Pk(Ji) is the space of polynomials of degree ok defined on J, = [xi-l, xi]; its dimension is k + 1. In practice we use here the case k = 1 corresponding to piecewise linear functions: we approximate the solution u of (22) by a piecewise linear function without continuity conditions at the nodes; this choice enables us to compare methods F and L-R. A basis for P’(J) is then given, for 1~ i 5 M by the functions (fig. 3):
u!(x) = h;‘(x u’(x) = -h;‘(x with hi =
X; -
xi~l.
-xi-,), -x,),
Xi-1 I
X
I
xi,
xi_, I x 5 xi,
(234 GW
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
8
I
I
0 =
x0
9
Ji
x. l-l
xi
I
x. 1+1
S-1
I I
73
,x
34
Fig. 2. Solution to problem (F’).
Every function u E P’(4) can be written U(X)’ U(Xi-l)Vp(X)+ U(Xi)Ui’(X).
(24)
For each point x1, x2, . . . , xM there are two basis functions for a total of 2M basis functions. The M unknown functions ul, . . , , uM define 2M degrees of freedom, and there are one degree of freedom u. at x0 = 0 and one degree uM+l at x,+, for a total of 2M + 2. Using (22a) and substituting the 2M basis functions o !, UT for Vi in (22b) lead to a system of 2M + 1 equations. The number of unknowns ~0, ui(xi-l), Ui(xi), uM+Ithat are really present in this system
X.
c-1
hi
Fig. 3. A basis for P’(4).
yi
74
P. Arminjon, C. Beauchamp, Continuous and discontinuous finite element methods for Burgers’ equation
depends on the choice of the parameters ai. We observe that u. and Us+] only appear in the first and last equation, respectively. We distinguish four cases. Case 1.
(Y~=O and
(Ye = 1.
The unknowns u. and uM+l disappear from our system, leaving 2M + 1 equations in 2M unknowns. In this case we shall omit the basis function ZJh and the corresponding equation. This amounts to considering the trial function uM and test functions ?_&in the subspace of H1(Jw) consisting of functions which vanish at x~. Case2
ao=O,
aMfl.
Case 3.
a0 # 0,
(Ye = 1.
In these two cases one unknown disappears: with a square system of order 2M + 1. Cuse4.
ao#O
and
u. in case 2 and uM+l in case 3, and we are left
~y~=l.
There are now 2M + 2 unknowns and only 2M + 1 equations. additional condition on u, for instance u. = 0 or &+l = 0.
We simply
take one
4.1. Computation of the integralsin (22) Using (23) and dui/dx = [ui(xi) - ui(xi_l)]/(Xi - Xi-l), we obtain
I dx dui
v{(x)dx
=
ui(xi)-2ui(%l),
j
=
1,2.
(25)
Ji
To compute _l’,,vi(x)f (x, ui) dx, we use the three-point hi [ ij-g 5
I
J, g(X)dX=T
5
+gg
(
l+A66_+1-X6% 2 1 1
( l-VSx,_ 2
11
+1+V6% 2
2
Gauss-Legendre xi)
Xi
+$
formula
g(x’-l;X.)
>I *
(26)
This leads to ~ ui(X)f(X, ui)dX -3 I, + + ui(Ci)f
(Ci,
Ui (Ci))]
[$ vC(ui)f(uiyui(ui))+t uj(bi)f(bi, ui(bi))
3
(27)
where ai, bi and ci are the three arguments of the function g appearing in the right member of
P. Arminjon,
C. Beauchamp,
(26) respectively.
we
get
75
Setting
l-A&% 2
y1=
Continuous and discontinuous finite element methods for Burgers’ equation
!
y3=
Y2 =;,
’
1+V% 2
(28)
,
IJ !(a,>
=
yl,
21 !(bi)
=
72,
21 I!(Ci>
=
737
2, ?(ai)
=
737
vf(bi)
=
y29
uf(Ci)
=
Yl,
(29)
and from (24) N(R)
=
y3ui(Xi-1)
+
ylUi(Xi),
&(bi)
=
y2”i(xi-l)
+
y2&(Xi),
h(G)
=
ylui(Xi-1)
+
y3&(Xi)*
(30)
Letting Ui(Xi)= u(x;) = Ui,- and Ui+l(xi) = obtain the following system of equations: CYRUS.-+
(1 - (Y~)u~,+ = u”
-+
+
+ $
Yf(ci,
y3Ui,-
-w
$
[$
+
5
YJ(Ci3
ylUi.-
+
yJ(ai7
YlUi-I,+)]
y3Ui-I,+)
= Ui,+
and introducing
(25) to (30) into (22) we
(given),
+ 1 [i yf(ai,
+
u(x+)
ylUi,-
-
(1
(3la) +
Y3Ui-I,+) + $ Yf(bi,
-
cYi)(Ui.+ - Ui,-) =
-
(Yi-l(Ui-I,+
- Ui-I,_) =
(31’4
0,
YlUi,- + Y3Ui-_1,+) + $ YJ(b,
1
Y2Ui,- + Y2+1,+)
YzUi,_ + YzUi-*,+)
0.
In our numerical experiments with method L-R we use the above procedure (lo)-(12) rewritten in the form du” _&c2w_1, dx
At u”(x)
u”(0) = 0.
(3Ic) to solve problem
(32)
We only consider the cases ai = 0 (i = 0,1, . . . , M) and ai = 1 (i = 0,1, . . . , M) here. They are particular cases of cases 2 and 3, respectively, discussed earlier. Remark. In the numerical experiments the approximation of u”-‘(x) on the interval Ji caused some difficulty. We first chose to use the natural approximation
U”-lI,,(X)= u~r~+v~(x)+u~~‘u~(x),
(33)
76
P. Arminjon. C. Beauchamp.
Continuous and discontinuous finite element methods for Burgers’ equation
but this led to an unstable scheme. To improve our scheme, we tentatively replaced the values uY:,‘.+and UT:’ in (33) by the values
and (34) which are precisely the numerical approximations of the solution u at t = t” at points Xi-1 and x,. With this modification our scheme became stable. 4.2. Resolution of the nonlinear system (31) Case (i).
Qyi= 0,
i=O,l,...,
M.
In this case, since (Y()= 0, (31a) immediately gives uo,+ = U” (given), and uo,_ does not appear. The last term of (31~) also vanishes, and we use the known value of uo,+ to determine ul,_ by solving (31~) with i = 1, which is an equation in one unknown ul,_. We can then find ul.+ by solving (31b). At the ith step we determine Ui,- from (31~) and then Ui,+ from (31b) (i = 1,2, . . . , M). Case (ii).
Cui=
i=O,l,...,
1,
M.
From (31a) we have uo,_ = u” (given); since (Y;= 1 (1 I i I M), the last term of (31b) vanishes. For each i = 1,2, . . . , M, (31b) and (31~) now form a system of 2 equations in 2 unknowns u,_~,+, Ui,- which is solved by Newton’s method: we first determine uo,+ and ul,_, then ul,+ and u2,_, and so on. (Notice that in this case the variational formulation (22) of the L-R method really coincides with that of Lesaint and Raviart [ll]. Otherwise (22) is a generalization due to Delfour and Trochu [6].)
5. Numerical experiments
We have performed several types of experiments to test the efficiency of methods F and L-R presented in sections 3 and 4 in the frame of problem (3) for Burgers’ inviscid equation as well as for the viscous equation (1). We considered three different initial functions for the inviscid equation (3): x
(11)
OSxSi, (34)
i
(12)
for
uo(x) = l-x
u&)=4x(1-x)
for
$SxSl, OlXll,
(35)
P. Arminjon,
C. Beauchamp,
U”(X)=
(13)
Continuous and discontinuous finite element methods for Burgers’ equation
2x
05x<;,
1-X
+x51.
I
77
(36)
This last initial function has a jump discontinuity at x = l/2. To find the analytic solutions for 0 < x < 1 and t > 0, we used the property that u (x, t) is constant along each characteristic line of equation x = u,(x,)t + x0 in the (x, t)-plane, as we did in section 2 for the initial function (11); in this case the solution is given by (7). For (12) we have
I-
u(x, t) = 1 + 4t - ((1 + 4t)‘- 16xt)l”][ 1 _ 1 + 4t - ((1 + 4t)* - 16xt)1’2 2t 8t
(37)
We shall deal with the case (13) in a separate subsection. 5.1. Discontinuous elements for the inviscid equation For methods F and L-R and initial functions (11) and (12) we have computed the mean relative error at all grid points (xi, t”) and the maximum relative error at the nodes (xi, tN = 1) for different values of Ax and At for 0 5 t 5 1 and 0 I x I 1. The results appear in tables l-5, with the exception of those for the case (i) above, i.e. when (Y~= 0 (i = 0,1, . . . , M), for the L-R method; in that case the L-R method was rather unstable. For most values of Ax and At it diverged; for the initial function (11) we obtained convergence only if the particular values of x given by x = 1/2(t” + 1) were nodal points. These values are precisely those for which the derivative du(x, t)/ax is discontinuous. The number (t” + 1)/(2Ax) is then a positive integer i (M/2 < i CM) for each n = 0,1, . . . , N. If Ax has the form (l/m) (m E IV), this implies that At/Ax must be an integer since (t” + l)/Ax = nAt/Ax + l/Ax must then equal 2i (n = 0, 1, . . . , N). Even when the L-R method does converge, the results are not very good. For instance, with Ax = 10e2, At = 2 x lo-* the mean relative error at the nodes is 4% and the maximum relative error is 25% ; these mediocre results are the best obtained by method L-R for the case ai = 0 (i = 0,1, . . . , M). The large value of the maximum error is probably due to the discontinuity of u(x, t) at x = 1, t = 1 since the maximum error occurs near this point. Tables 2 and 4 show the results obtained for the initial function (11) and the inviscid equation (3) by methods F and L-R, respectively. In this case method F gives better results; Table 1. Discontinuous finite elements (method F. initial function (11)) Mean relative error
Fl
Ax
At
l/20 l/20 l/40
l/loo l/500 l/200
left 10.3% 7.8% 6.56%
right 12.4% 9.4% 7.76%
Computing time
26 set 65 set 93 set
78
P. Arminjon,
C. Beauchamp,
Table 2. Discontinuous
Continuous and discontinuous finite element methods
finite elements Mear relative
Ax = l/10
Ax = l/20
Ax= l/40
Ax = l/80
(method
F, initial function
Maximum
relative
for Burgers’ equation
11, cr = l/2)
Computing
time
At
error
error at t = 1
l/50
l/100 l/150 l/200
2.99% 3.39% 3.64% 3.79%
27% 34% 36% 37.8%
7.87 set 12.23 set 16.54 set 20.8 set
l/50
1.49%
13.3%
l/100 l/200
1.23% 1.35%
24.6% 31.12%
17.13 set 26.6 set 40.6 set
l/50 l/100 l/150 l/200
0.76% 0.61% 0.56%
13% 20.2% 24%
71.77 set 81.83 set 94.95 set
l/200
0.401%
13.72%
295
diverge
set
here we have used the approximate value defined by (19) for method F, while we only retained the left and right estimates uy_ and uy+ in table 1. Tables 3 and 5 show the corresponding results for the initial function (12). In this case method L-R is slightly superior to method F and appears to be more stable since we did not find a case where it would diverge. The accuracy is Table
3. Discontinuous
finite
elements
(method
F, initial
function
12, (Y =
l/2)
At
Mean relative error
Maximum relative error at t = 1
Computing time
Ax = l/10
l/50 l/100 l/150 l/200
1.96% 1.54% 1.75% 2.07%
3.23% 1.66% 2.79% 4.29%
8.8 14.7 19.7 24.6
Ax = l/20
l/50 l/100 l/l50 l/200
0.86% 0.67% 0.59%
diverge 1.41% 0.94% 0.70%
36.6 set 46 set 55 set
l/50 l/100 l/150 l/200
0.43%
diverge diverge diverge 0.67%
180 set
l/400
0.217%
0.328%
885 set
Ax = l/40
Ax = l/80
set set set set
P. Arminjon, C. Beauchamp, Continuous and discontinuous finite element methods for Burgers’ equation Table 4. Discontinuous finite elements (method L-R, OStIl,cwi=l,i=O,l,..., M)
At
Mean relative error right left
79
initial function (11)
Maximum relative error left right
Computing time
Ax = l/20
l/50 l/70 l/80 l/90 l/100 l/120 l/200
1.67% 1.93% 1.59% 1.97% 2.00% 2.04% 2.20%
1.61% 1.83% 1.52% 1.87% 1.90% 1.97% 2.15%
26.5% 25% 28% 26.4% 27.3% 28% 31.3%
25.2% 20.4% 27% 23% 24% 25% 29.7%
10.5 12.4 13.5 17.4 18.7 21.5 32.8
set set set set set set set
Ax = l/40
l/10 l/20 l/40 l/80 l/100 l/150 11200
4.8% 2.77% 1.65% 1.17% 1.11% 1.09% 1.09%
4.88% 2.81% 1.64% 1.14% 1.08% 1.05% 1.08%
28.63% 23.30% 21.24% 23.50% 25% 28.10% 30.20%
27.6% 21.3% 17.27% 18.10% 20% 24.15% 27.20%
7.88 set 11.45 set 17.80 set 28.30 set 33.20 set 47.82 set 62 set
Ax = l/50
l/150 l/200
0.90% 0.89%
0.88% 0.87%
26.9% 29%
22.30% 25%
58.7 set 76 set
Ax = l/80
l/160 l/320
0.65% 0.59%
0.64% 0.58%
24% 30.4%
19.6% 26.8%
96.9 set 187 set
comparable, but the computing times are better; by contrast, method F was often divergent for (12) except if At/Ax was chosen small. As mentioned earlier, the results for method L-R with cyi= 1 (tables 4-5) were obtained by solving system (31b), (31~) with Newton’s method, where we computed the Jacobian and its inverse at each iteration. To reduce the computing time, we applied Newton’s method with the same Jacobian matrix at each iteration. This led to the same results and reduced the computing time by approximately 35%. The results appear in table 6. 5.2. Continuous space-time finite elements for the inviscid equation For a comparison, we also applied to (3) with initial function (11) a continuous space-time finite element method presented in [l] for the viscous equation (1). Since v = 0 here we drop the boundary condition ~(1, t) = 0 and introduce a degree of freedom at x = 1 by considering a basis function qM = (xM - xM_$‘(x - xMP1)for xMel I x 5 xM. For Ax = l/20 and At = l/100 the mean relative error was 56%) for a computing time of 113 set, as compared to (1.23%) 26 set) and (2%) 18 set) for methods F and L-R respectively. The maximum errors were respectively 63%) 24% and 27%) showing the clear superiority of the discontinuous methods for the inviscid equation.
P. Arminjon,
80 Table
C. Beauchamp,
5. Discontinuous
finite
Continuous and discontinuous finite element methods for Burgers’ equation
elements
(method
L-R,
initial
function
(12),
OltSl,(Y,=l) Mean relative error left right
At
Maximum relative error at t = 1 left right
Computing time
Ax = l/10
l/20 l/40 l/50 l/80 l/120 l/160 l/200 l/1000
2.67% 1.80% 1.95% 2.63% 3.14% 3.44% 3.64% 4.37%
2.52% 1.42% 1.52% 2.24% 2.83% 3.18% 3.41% 4.3%
3.96% 1.27% 1.7% 2.49% 3.04% 3.36% 3.57% 4.4%
4.0% 1.19% 1.63% 2.42% 2.97% 3.29% 3.5% 4.3%
5.9 set 7.8 set 8.7 set 10.9 set 13.47 set 16.35 set 19.1 set 77.6 set
Ax = l/20
l/20 l/40 l/70 l/80 l/100 l/500
3.20% 1.53% 0.98% 0.96% 1 .O% 2.0%
3.24% 1.51% 0.88% 0.84% 0.86% 1.94%
5.17% 2.2% 0.90% 0.68% 0.76% 1.82%
5.19% 2.21% 0.91% 0.69% 0.75% 1.81%
7.84 set 11.6 set 16.35 set 17.6 set 20 set 77 set
Ax = l/40
l/70 l/80 l/l40 l/l60 l/200
0.987% 0.82% 0.52% 0.51% 0.51%
0.88% 0.83% 0.49% 0.47% 0.47%
0.9% 1.11% 0.49% 0.37% 0.35%
0.91% 1.15% 0.49% 0.37% 0.35%
13.9sec 30.6sec 47.16 set 52.4 set 62.5 set
Ax = l/80
l/160
0.43%
0.43%
0.59%
0.59%
101 set
S.3. Discontinuous elements for the viscous equation
For the full Burgers’ equation (1) with boundary conditions ~(0, t) = ~(1, t) = 0 for t 10 we cannot apply the L-R method, but we did apply method F with the same regular initial Table 6. Discontinuous Jacobian matrix)
Ax
I1 11 11 11 I1 12
l/20 l/40 l/80 l/80 l/160 l/40
At
l/80 l/160 l/160 l/200 l/200 l/160
finite
elements
Mean relative error
(method
L-R,
(Y,= 1 with
Maximum relative error at t = 1
left
right
left
right
1.95% 1.09% 0.66% 0.62% 0.45% 0.51%
1.85% 1.06% 0.65% 0.61% 0.45% 0.47%
25.8% 28.6% 24.6% 26% 22% 0.37%
21.8% 24.9% 19.6% 21.5% 18.7 % 0.37%
constant
Computing time
12.8 set 36.7 set 69.4 set 84 set 162 set 39 set
P. Aminjon,
C. Beauchamp,
Continuous
and discontinuous
finite element methods for Burgers’ equation
81
Table 7. Errors Discontinuous (method At
a
F
l/100 l/200 l/500 l/1000
1.5 0.3 0.3 0.5
39% 19% 6.2% 3.015%
elements F) Emax
Continuous space-time elements Computing
59.4% 20% 9.7% 4.95%
26 46 70 107
time
set set set set
c 0.78% 0.49% 0.84% Note: reached 0.085
Computing
for the
time
2 set 3.5 set 5 set At = l/140, E minimum value
function as in [l]. With v = 1 and Ax = l/20 we obtain the results in table 7 for the mean and maximum relative errors C, E,,,. Here we chose to approximate ul by
my- + (1 - CX))u~+.
(38)
These results show the superiority of continuous space-time elements for the viscous equation. This is in agreement with the intuitive observation that regular functions should be better approximable by regular approximating functions.
Fig. 4. Characteristics
for problem
with discontinuous
initial function.
82
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
5.4. Discontinuous elements for the inviscid equation with the discontinuous initial function (13)
The initial function (13) (see (36)) has a jump at x = l/2, and we can determine the solution only in some regions of the (x, t)-plane. There is only one characteristic line going through each point of regions R1 and R3 (fig. 4) so that the solution u(x, t) can be uniquely determined at these points. But for region RZ each point (x, t) is the intersection of two characteristic lines, one originating in the interval 0 < x < l/2, the other in l/2 < x < 1. We find
j-$$ for uW)=
(x,t)ER1,
l_x
l_t
for
(x,t)ER3.
(39)
:
We computed the solution with method F and with method L-R with ai = 1. In both cases the numerical solution is well-behaved in regions R1 and R3 and converges to the value given by (39). Figs. 5 and 6 show the graphs of the numerical results obtained for u(x, t) at time t = 0.3 set
J(X,O.3)
Fig. 5. Behaviour
of u(x, 0.3) (method F, Ax = l/40, At = l/200).
P. Arminjon, C. Beauchamp,
83
Continuous and discontinuous finite element methods for Burgers’ equation
0.9
0.5
x 0.5
Fig. 6. Behaviour
1
of u(x, 0.3) (method L-R, Ax = l/40, At = l/160).
by methods F and L-R, respectively, and for Ax = l/40, with At = l/200 (method F), At = l/160 (method L-R). The mean relative error at points where the solution is known is 0.28% for method F and 0.35% for method L-R; the corresponding computing times were 47 set and 13 sec. 6. Concluding remarks Our results indicate that although the two discontinuous finite element methods presented here do not have the same accuracy and efficiency as other space-time finite element methods for continuous solutions of the viscous equation (see [l], [2]), they seem to be superior in the case of the ill-behaved solutions of the inviscid equation. This should make them well suited to problems in fluid mechanics where the viscosity can be neglected, and problems associated with the system of primitive equations of dynamic meteorology. We are presently testing these methods for some boundary layer problems where the viscosity is very small. References [l] P. Arminjon and C. Beauchamp, A finite element method for Burgers’ equation in hydrodynamics, Numer. Meths. Eng. 12 (1978) 415-428. [2] P. Arminjon and C. Beauchamp, Numerical solution of Burgers’ equation in two space dimensions, Meths. Appl. Mech. Eng. 19 (1979) 351-365.
Int. J. Comp.
84
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
[3] P. Anninjon, Finite element analysis of moving loads on a floating ice shell, Topics in Numerical Analysis, Proceedings 3rd Conference Royal Irish Academy (Academic Press, London, 1977) 299-314. [4] R. Bonnerot and P. Jamet, A second order finite element method for the one-dimensional Stefan problem, Int. J. Numer. Meths. Eng. 8 (1974) 811-820. [5] A. Chalifour, Analyse dune methode d’elements finis discontinus pour les equations differentielles lineaires, Report no. CRM-643 (Centre de Recherches Mathematiques, Univ. Montreal, 1976). [6] M. Delfour and F. Trochu, Discontinuous approximation of ordinary differential equations and applications to optimal control problems, Report no. CRM-751 (Centre de Recherches Mathematiques, Univ. Montreal, 1977). [7] M. Delfour, W.W. Hager and F. Trochu, Discontinuous piecewise polynomial Galerkin methods for ordinary differential equations, Report no. CRM-830 (Centre de Recherches Mathematiques, Univ. Montreal, 1978). [8] M. Fortin, Resolution numerique des equations de Navier-Stokes par des elements finis de type mixte, Rapport IRIA-LABORIA 76-184 (1976). [9] M. Fortin, Mixed finite element methods for incompressible flow problems, J. Comp. Phys. 31 (1979). [lo] P. Jamet and R. Bonnerot, Numerical solution of the Eulerian equations of compressible flow by a finite element method which follows the free boundary and the interfaces, J. Comp. Phys. 18 (1975) 21-45. [ll] P. Lesaint and P.A. Raviart, On a finite element method for solving the neutron transport equation, Proceedings of the Symposium Mathematical Aspects of Finite Elements in Partial Differential Equations, Madison, Wisconsin (Academic Press, New York, 1974). [12] P. Lesaint, Continuous and discontinuous finite element methods for solving the transport equation, The mathematics of finite elements and applic. II, Proc. 2nd Brunel Univ. Conf. Inst. Math. Appl., Uxbridge, 1975 (Academic Press, London, 1976) 151-161. [13] J.T. Oden and J.N. Reddy, An introduction to the mathematical theory of finite elements (Wiley, New York, 1976). [14] H.H Rachford and M.F. Wheeler, An H-’ Galerkin procedure for the two-point boundary value problem, in: C. de Boor (ed.), Mathematical aspects of finite element methods in partial differential equations (Academic Press, New York, 1975) 353-382. [15] W.H. Reed and T.R. Hill, Triangular mesh methods for the neutron transport equation, LA UR 73479 (Los Alamos Lab., 1973). [16] L.C. Wellford and J.T. Oden, Discontinuous finite element Galerkin approximations of shock waves in nonlinear elastic solids, TICOM Report 74-6 (The Texas Institute for Computational Mechanics, Univ. Texas at Austin, 1974). [17] L.C. Wellford and J.T. Oden, Discontinuous finite element approximations for the analysis of shock waves in nonlinear elastic materials, J. Comp. Phys. 19 (1975). [18] L.C. Wellford and J.T. Oden, A theory of discontinuous finite element Galerkin approximations of shock waves in nonlinear elastic solids -Part 1. Variational theory, Comp. Meths. Appl. Mech. Eng. 8 (1976) 1-16. [19] O.C. Zienkiewicz, The finite element method in engineering science (McGraw-Hill, London, 1971). [20] R. Courant and D. Hilbert, Methods of mathematical physics (Interscience, Wiley, New York, 1962). [21] P. Arminjon, Numerical solution of a boundary layer problem with continuous and discontinuous finite elements. To appear.
MECHANICS COMPANY
AND ENGINEERING
25 (1981) 65-84
CONTINUOUS AND DISCONTINUOUS FINITE ELEMENT METHODS FOR BURGERS’ EQUATION*
Paul ARMINJON Dipartement
and Claude BEAUCHAMP
de Mathe’matiques et de Statistique, Universite’ de Monkal, Canada
Revised
Received 12 December 1979 manuscript received 11 February
B.P. 6128, Montre’al H3C 3J7,
1980
We apply two discontinuous finite element methods to the inviscid Burgers’ equation and to the full equation with viscosity. In both cases we compare with a continuous space-time finite element method previously studied. For v = 0 discontinuous methods give better results, while the reverse prevails for the viscous equation.
1. Introduction
In previous reports [l], [2], a space-time continuous finite element procedure derived from Bonnerot and Jamet [4] was applied to the numerical solution of Burgers’ equation in one space dimension: 2
O
g+ug=v$$,
with initial and boundary
t>O
(la)
u(O,t)=u(l,t)=O,
(lb)
conditions and
u (x, 0) = uo(x) and in two space dimensions: vdu, v, + uv, + vu, = vdv, ut + uux + vu, =
(2)
with corresponding initial and boundary conditions; here A stands for the 2-dimensional Laplace operator. The solutions were regular, and the method had an (observed) order of accuracy equal to 2. In the present paper we shall apply two finite element methods of the discontinuous type to the inviscid Burgers’ equation in one space dimension (3)
ur + uu, = 0 with both regular and irregular initial functions * This research has been supported the province of Quebec.
by the National
t&(x).
Research
Council
of Canada
and the Ministry
of Education
of
66
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
We shall also apply to (3) the space-time continuous finite element method described in [l], for comparison. Finally we shall apply one of the discontinuous finite element methods presented here to the full equation (1) with viscosity, and compare the results with those obtained in [l]. Discontinuous finite elements were introduced earlier for the neutron transport equation by Reed and Hill [15], Lesaint and Raviart [ll]; for problems in nonlinear elasticity by Wellford and Oden [16], [17], [18]; for the equations of hydrodynamics by Fortin [8], [9]; and for ordinary differential equations by Rachford and Wheeler [14], Delfour, Hager and Trochu [6], [7], and others. Section 2 gives some properties of the inviscid Burgers’ equation which are of concern for our study, sections 3 and 4 describe the two discontinuous finite element methods considered here; these methods are called method F and method L-R here since they are obtained from two methods for ordinary differential equations by Fortin and Lesaint-Raviart, respectively. Section 5 describes some of the numerical experiments which have been performed. It is well known that solutions of (3) may present discontinuities even for regular initial data, while solutions of (1) tend to be more regular than the corresponding initial data, due to the smoothing effect of viscosity. The present paper aims at showing the advantage of discontinuous finite elements to eventually get a better fit of possible discontinuities or singularities such as those currently met in gas dynamics and shock tubes, while preparing for the use of these methods to study the numerical solution of the system of equations governing the atmospheric flow in numerical weather forecasting, in a manner designed to give a good account of the irregularities (“quasi-discontinuities”) arising in the vertical coordinate, particularly for the lowest kilometer of the atmosphere (boundary layer effect). Our results indicate that (i) for the inviscid equation (3) the discontinuous finite element methods considered here are more efficient and accurate than the continuous space-time finite element method of reference [l] (which will be called method C here). We also obtained good results for a discontinuous initial function. (ii) for the full Burgers’ equation with viscosity, the discontinuous methods F and L-R are very clearly inferior to the continuous method C.
2. Some properties of the inviscid Burgers’ equation We consider the initial-boundary ut + uu, = 0,
t >o,
value problem x > 0,
(3)
u (x70) = uo(x ),
(3a)
u(0, t)=
W)
0
and assume that this problem has a regular solution u(x, t). The characteristic curves of equation (3) are defined as the integral curves of the ordinary differential equation dx - = u(x, t).
(4)
dt
On such a characteristic
curve we have x = x(t), and the solution u(x, t) of problem
(3) is a
P. Anninjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
67
function of t only, for which du au au dx --+--P-+u~~O, dt- at ax dt
au
by (3). The solution u(x, t) is thus a constant along each characteristic curve. But since u = a constant on each such integral curve of dx/dt = u(x, t), we have dx/dt = u (x0, to) = uO(xO)so that the characteristic curves are the lines x(t) = Uo(z$ + x0,
(5)
or t = [uo(xO)]-‘(x - x0). If we consider for example the case X
z&J(x)=
l-x i 0
x 5 l/2, 1/25x 5 1, elsewhere,
0 s
(6)
we have the following situation (see fig. 1).
Fig. 1. Characteristics
for solution of eq. (3).
68
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
(i) If 0 < x0 -=cl/2, the characteristic lines t = and do not intersect each other. For all P = P = (x, t) with 0 < x < 1 and t > 1 there is only the value of the solution u(x, t) is completely constant value): u(x, t) = u(x0, 0) =
have slopes varying from +m to +2 (x, t) in the trapezoidal region OMBA or all one characteristic line passing through P, and determined by this line (along which u has a (x - x0)/x0
u&J = x0
(since 0 5 x0 I l/2) with t = (x - XO)/XO, and therefore ~0 = x/(1 + t) = u (x, t). (ii) For P = (x, t) in the triangular region MBC there is also exactly one characteristic going through P originating at (x0, 0) such that t = (x - x,)/x, or x0 = (x - t)/(l - t) with l/2 I x0 I 1, and therefore u(x, t) = 1 - x0 = (1 - x)/(1 - t). Notice that each such characteristic goes through B = (1,l). Summing up these two cases, we have &
for
O< x < 1
and
t >
2x - 1 (region OMBA or above),
U(X, t) = ’
(7)
and
t <
2x - 1 (triangular region MBC).
(iii) For x > 1 and t < 2x - 1 there is no characteristic going through (x, t) and originating at a point (x0, 0) where uO(xO)would be #O. In fact, the only characteristic going through such a point Q is the parallel to the t-axis, x = x Q. Since uO(xQ)= 0, we have u(Q) = 0. (iv) For x > 1 and t > 2x - 1 (triangular region T) there are two characteristic lines through P = (x, t) originating in (xl, 0) and (x2, 0), respectively, with 0 < x1 < l/2 and l/2 < x2 < 1. This corresponds to a jump discontinuity at time t with the values
Thus even for the regular initial function uO(x) defined by (6) Burgers’ inviscid equation can generate singular solutions. This example suggests the use of a discontinuous finite element method even if we limit ourselves, in this paper, to the values of u(x, t) for which 0 I x I 1, 0 5 t I 1, thus leaving the first singular point x = 1, t = 1 on the boundary of the domain of computation.
3. Method F. Discontinuous
finite elements derived from a classical variational
equivalent
We consider eq. (3) with initial and boundary conditions (3a)--(3b), and first perform a discretization in time by approximating the time derivative uI by a difference quotient, thus
P. Arminjon, C. Beau&,
Continuous and discontinuous finite element methods for Burgers’ equation
69
reducing the partial differential equation (3) to a system of ordinary differential equations to which we shall apply a method described by Fortin in [8], [9]. Let 0 = t” < t’ < . . . < t” . . . be an ordered sequence of points on the t-axis with At = t”+’- t” (not necessarily a constant). Setting U”(X) = U(X, t”), we use the difference quotient approximation
au
in
_z
-d-
(9)
At
at
to obtain the following approximate n
’
ad
in
system
~"-1
(10)
ax +dt= At’
uO(x> = uo(x), u”(0) = 0,
(11) 12
=0,1,2,...
(12)
To obtain our variational equivalent of this system, we take an arbitrary function !P = V(x) in CZ[O, l] = {P E C’[O, 11, ?-P(O)= V(1) = 0}, multiply (10) by !P and integrate by parts to obtain
-
=[,‘$+‘dx.
(13)
Let O=xO
=
OsjsM,
Sj,
F(Xj)=O,
[ Vi,(Xj)
OsjsM, =
OsjsM,
0,
$$ (Xj) =
6,,
OsjsM
(144
W)
70
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
?P$‘, restricted
for i = A4 we take the function i = 1,2, . . . . A4 - 1 by
(y>‘[
3-&(x
-xi4
Xi-1 I
2(?)1_3(?_&%)2+1,
Vi(x) =
0
q;(x) =
(x
0
0
-p
-xi)[q.f-y+
Xi,
xi 6 x 5 xi+1,
1
1 )
Xi-]
11,
Xi 5
I
X I Xi, X s
Xi+19
(W
elsewhere,
(x -;;-d2
I
X I
are given, for
elsewhere,
(x--;-d’ [”
?Py”(x>=
to [xIMPI,~~1. These functions
[x --y-1_
1],
elsewhere.
The trial functions will generally have two different values ui,_ = u(x;) and Ui.+= u(x+) at the knots x1, x2, . . . , xMpl. At the last knot X~ we shall only consider the left-hand value Us,_ as a degree of freedom, and at x0 we take uo.- = u o,+= uo(0), where the last term denotes the given value of the initial function uo(x) at x = 0. This gives a total of 2(M - l)+ 1 = 2M - 1 degrees of freedom, and the dimension of the space Ps[O, l] of test functions is also equal to 2M - 1. If we now take ly = !I”‘;, P = ?Pyi,in (13) and use appropriate quadrature formulae, we obtain a system of 2M - 1 nonlinear equations in the 2M - 1 unknown nodal values uy_ and u;, (1 I i I M - 1) and uL,_ (for n = 1,2,. . .). Since the functions we must integrate in (13) are piecewise polynomials, for each integral we choose a Newton-Cotes formula which gives its exact value; this yields the following system: u;,_=uI;,+=u(O,t”)=O
and
I
by(3b),
(IQ
=&[3u:_;‘++7u:~‘+7u:;‘+3u~~~_],
i=l,2
,...,
M-l,
I (u:_,,)‘+;u;_I.+u:--~(u:-)2+~(u:+)2+~u~+u~+~,-+(u:,,y] +s u:,l,-+3(u” 2 I,+ [
-u”_)-uy_l I,
I
(17’4
.+
, I =%Iu:;:._+q(u;;‘+u;_L)-u:_;I+ 1
i = 1,2, . . . , M - 1,
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
71
and for i = A4 we have
(17c) Eq. (17a) is obtained by taking P = !Pj in (13) while (17b) results from the choice ly = Vi. We note that if we assume that uy_ = uy+ (continuity of the trial functions), (17b) becomes ,&
I’ du;.
n u;+‘,l+2 ui+12U:-1+6
uy-‘_l .
u;+,2
UY-1
1 (18)
which is a consistent difference scheme for the equation uu, + uI = 0. If one wants to approximate a continuous solution, we shall see in section 5 that we get a better approximation of the nodal values ul = u”(xi) if we choose the approximate value
u; =&4;_+u;+).
(19)
System (17a)-(17c) consists of 2M - 1 nonlinear algebraic equations in 2M - 1 unknowns, which we solve with a simplified form of Newton’s method for systems. Specifically, to solve a system F(X) = 0, where F: Iw” +w” and X = (Xi)E Iw”, we start from a first approximation X0 = (Xp) and use the following iteration: X!+l= X: -A(Fi(Xk)/$$
(Xc)),
(20)
I
where A is a constant chosen to accelerate convergence; we found out that A = l/2 gives faster convergence than the usual value A = 1. When trying to approximate X = [u;,_, UT,+,. . . , uL-~,+, uh,-1’ we use the first approximation X0 = [uY,:*, u;;‘, . . . , uEll,+, ub:?]‘.
4. Method GR.
Discontinuous elements derived from a modified variational equivalent
We first perform the same time discretization as in method F and apply to the corresponding system of ordinary differential equations a modified variational formulation derived from the one appearing in Lesaint-Raviart [ll] by Delfour and Trochu (see [6], [7]). This leads to a family of numerical schemes depending on some real parameters ai. For Burgers’ equation
72
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
we use the values ai = 0 for all concerned value problem
g 1
=f(x, 4
m
u(0) = u”
indices i, or ai = 1 for all
i.
We consider the initial
x E (0, T), given
and a partition 0 = x0 < x1 < . . . < x,+, = T of the interval [0, T]. We set Ji = [xi-l, xi] and write Ui(X) for the restriction of the unknown function u = u(n) to the interval Ji (i = 1,2, . . . , M). We do not assume continuity of u at the nodes, so that Ui(Xi) may differ from Ui+l(xi). It can then be shown (see Delfour and Trochu [6]) that problem (P) is equivalent to the M + 1 variational problems i
(
u. + (1 - ao)[ul(0)-
-
(1 -(Yi)fJi(Xi)[&+l(Xi)-
uo] = u”
given, i.e.
ui(xi)l
=
O,
aOuo + (1 - cyo)ul(0) =
lsi
d’,
(224
Wb)
for all functions Ui
EHl(Ji)={EL2(Ji):$EL2(Jj)}. U
The numbers uo, Ui(Xi-I), Ui(Xi) (i = 1,2, . . . , M) and uM+, are the unknowns; u. may be considered as the left-hand value of our approximate solution u at x = 0, and u~+~ as the right-hand value at x = xM. At the node xi the left-hand value is Ui(Xi)and the right-hand value is Ui+l(Xi) (see fig. 2). Observe that if we assume the solution u to be continuous on [0, T], the last two terms in (22b) vanish, and the equivalence of (P) with (22) is obvious. But the interest of the method lies in its ability to handle discontinuities at the nodes; terms like Ui+I(Xi)- Ui(Xi)represent the jump at Xi. Now to find an approximation of the solution of (22) we seek functions u! E P“(A) C H’(4) for which (22) is satisfied for all ViE P”(A). Pk(Ji) is the space of polynomials of degree ok defined on J, = [xi-l, xi]; its dimension is k + 1. In practice we use here the case k = 1 corresponding to piecewise linear functions: we approximate the solution u of (22) by a piecewise linear function without continuity conditions at the nodes; this choice enables us to compare methods F and L-R. A basis for P’(J) is then given, for 1~ i 5 M by the functions (fig. 3):
u!(x) = h;‘(x u’(x) = -h;‘(x with hi =
X; -
xi~l.
-xi-,), -x,),
Xi-1 I
X
I
xi,
xi_, I x 5 xi,
(234 GW
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
8
I
I
0 =
x0
9
Ji
x. l-l
xi
I
x. 1+1
S-1
I I
73
,x
34
Fig. 2. Solution to problem (F’).
Every function u E P’(4) can be written U(X)’ U(Xi-l)Vp(X)+ U(Xi)Ui’(X).
(24)
For each point x1, x2, . . . , xM there are two basis functions for a total of 2M basis functions. The M unknown functions ul, . . , , uM define 2M degrees of freedom, and there are one degree of freedom u. at x0 = 0 and one degree uM+l at x,+, for a total of 2M + 2. Using (22a) and substituting the 2M basis functions o !, UT for Vi in (22b) lead to a system of 2M + 1 equations. The number of unknowns ~0, ui(xi-l), Ui(xi), uM+Ithat are really present in this system
X.
c-1
hi
Fig. 3. A basis for P’(4).
yi
74
P. Arminjon, C. Beauchamp, Continuous and discontinuous finite element methods for Burgers’ equation
depends on the choice of the parameters ai. We observe that u. and Us+] only appear in the first and last equation, respectively. We distinguish four cases. Case 1.
(Y~=O and
(Ye = 1.
The unknowns u. and uM+l disappear from our system, leaving 2M + 1 equations in 2M unknowns. In this case we shall omit the basis function ZJh and the corresponding equation. This amounts to considering the trial function uM and test functions ?_&in the subspace of H1(Jw) consisting of functions which vanish at x~. Case2
ao=O,
aMfl.
Case 3.
a0 # 0,
(Ye = 1.
In these two cases one unknown disappears: with a square system of order 2M + 1. Cuse4.
ao#O
and
u. in case 2 and uM+l in case 3, and we are left
~y~=l.
There are now 2M + 2 unknowns and only 2M + 1 equations. additional condition on u, for instance u. = 0 or &+l = 0.
We simply
take one
4.1. Computation of the integralsin (22) Using (23) and dui/dx = [ui(xi) - ui(xi_l)]/(Xi - Xi-l), we obtain
I dx dui
v{(x)dx
=
ui(xi)-2ui(%l),
j
=
1,2.
(25)
Ji
To compute _l’,,vi(x)f (x, ui) dx, we use the three-point hi [ ij-g 5
I
J, g(X)dX=T
5
+gg
(
l+A66_+1-X6% 2 1 1
( l-VSx,_ 2
11
+1+V6% 2
2
Gauss-Legendre xi)
Xi
+$
formula
g(x’-l;X.)
>I *
(26)
This leads to ~ ui(X)f(X, ui)dX -3 I, + + ui(Ci)f
(Ci,
Ui (Ci))]
[$ vC(ui)f(uiyui(ui))+t uj(bi)f(bi, ui(bi))
3
(27)
where ai, bi and ci are the three arguments of the function g appearing in the right member of
P. Arminjon,
C. Beauchamp,
(26) respectively.
we
get
75
Setting
l-A&% 2
y1=
Continuous and discontinuous finite element methods for Burgers’ equation
!
y3=
Y2 =;,
’
1+V% 2
(28)
,
IJ !(a,>
=
yl,
21 !(bi)
=
72,
21 I!(Ci>
=
737
2, ?(ai)
=
737
vf(bi)
=
y29
uf(Ci)
=
Yl,
(29)
and from (24) N(R)
=
y3ui(Xi-1)
+
ylUi(Xi),
&(bi)
=
y2”i(xi-l)
+
y2&(Xi),
h(G)
=
ylui(Xi-1)
+
y3&(Xi)*
(30)
Letting Ui(Xi)= u(x;) = Ui,- and Ui+l(xi) = obtain the following system of equations: CYRUS.-+
(1 - (Y~)u~,+ = u”
-+
+
+ $
Yf(ci,
y3Ui,-
-w
$
[$
+
5
YJ(Ci3
ylUi.-
+
yJ(ai7
YlUi-I,+)]
y3Ui-I,+)
= Ui,+
and introducing
(25) to (30) into (22) we
(given),
+ 1 [i yf(ai,
+
u(x+)
ylUi,-
-
(1
(3la) +
Y3Ui-I,+) + $ Yf(bi,
-
cYi)(Ui.+ - Ui,-) =
-
(Yi-l(Ui-I,+
- Ui-I,_) =
(31’4
0,
YlUi,- + Y3Ui-_1,+) + $ YJ(b,
1
Y2Ui,- + Y2+1,+)
YzUi,_ + YzUi-*,+)
0.
In our numerical experiments with method L-R we use the above procedure (lo)-(12) rewritten in the form du” _&c2w_1, dx
At u”(x)
u”(0) = 0.
(3Ic) to solve problem
(32)
We only consider the cases ai = 0 (i = 0,1, . . . , M) and ai = 1 (i = 0,1, . . . , M) here. They are particular cases of cases 2 and 3, respectively, discussed earlier. Remark. In the numerical experiments the approximation of u”-‘(x) on the interval Ji caused some difficulty. We first chose to use the natural approximation
U”-lI,,(X)= u~r~+v~(x)+u~~‘u~(x),
(33)
76
P. Arminjon. C. Beauchamp.
Continuous and discontinuous finite element methods for Burgers’ equation
but this led to an unstable scheme. To improve our scheme, we tentatively replaced the values uY:,‘.+and UT:’ in (33) by the values
and (34) which are precisely the numerical approximations of the solution u at t = t” at points Xi-1 and x,. With this modification our scheme became stable. 4.2. Resolution of the nonlinear system (31) Case (i).
Qyi= 0,
i=O,l,...,
M.
In this case, since (Y()= 0, (31a) immediately gives uo,+ = U” (given), and uo,_ does not appear. The last term of (31~) also vanishes, and we use the known value of uo,+ to determine ul,_ by solving (31~) with i = 1, which is an equation in one unknown ul,_. We can then find ul.+ by solving (31b). At the ith step we determine Ui,- from (31~) and then Ui,+ from (31b) (i = 1,2, . . . , M). Case (ii).
Cui=
i=O,l,...,
1,
M.
From (31a) we have uo,_ = u” (given); since (Y;= 1 (1 I i I M), the last term of (31b) vanishes. For each i = 1,2, . . . , M, (31b) and (31~) now form a system of 2 equations in 2 unknowns u,_~,+, Ui,- which is solved by Newton’s method: we first determine uo,+ and ul,_, then ul,+ and u2,_, and so on. (Notice that in this case the variational formulation (22) of the L-R method really coincides with that of Lesaint and Raviart [ll]. Otherwise (22) is a generalization due to Delfour and Trochu [6].)
5. Numerical experiments
We have performed several types of experiments to test the efficiency of methods F and L-R presented in sections 3 and 4 in the frame of problem (3) for Burgers’ inviscid equation as well as for the viscous equation (1). We considered three different initial functions for the inviscid equation (3): x
(11)
OSxSi, (34)
i
(12)
for
uo(x) = l-x
u&)=4x(1-x)
for
$SxSl, OlXll,
(35)
P. Arminjon,
C. Beauchamp,
U”(X)=
(13)
Continuous and discontinuous finite element methods for Burgers’ equation
2x
05x<;,
1-X
+x51.
I
77
(36)
This last initial function has a jump discontinuity at x = l/2. To find the analytic solutions for 0 < x < 1 and t > 0, we used the property that u (x, t) is constant along each characteristic line of equation x = u,(x,)t + x0 in the (x, t)-plane, as we did in section 2 for the initial function (11); in this case the solution is given by (7). For (12) we have
I-
u(x, t) = 1 + 4t - ((1 + 4t)‘- 16xt)l”][ 1 _ 1 + 4t - ((1 + 4t)* - 16xt)1’2 2t 8t
(37)
We shall deal with the case (13) in a separate subsection. 5.1. Discontinuous elements for the inviscid equation For methods F and L-R and initial functions (11) and (12) we have computed the mean relative error at all grid points (xi, t”) and the maximum relative error at the nodes (xi, tN = 1) for different values of Ax and At for 0 5 t 5 1 and 0 I x I 1. The results appear in tables l-5, with the exception of those for the case (i) above, i.e. when (Y~= 0 (i = 0,1, . . . , M), for the L-R method; in that case the L-R method was rather unstable. For most values of Ax and At it diverged; for the initial function (11) we obtained convergence only if the particular values of x given by x = 1/2(t” + 1) were nodal points. These values are precisely those for which the derivative du(x, t)/ax is discontinuous. The number (t” + 1)/(2Ax) is then a positive integer i (M/2 < i CM) for each n = 0,1, . . . , N. If Ax has the form (l/m) (m E IV), this implies that At/Ax must be an integer since (t” + l)/Ax = nAt/Ax + l/Ax must then equal 2i (n = 0, 1, . . . , N). Even when the L-R method does converge, the results are not very good. For instance, with Ax = 10e2, At = 2 x lo-* the mean relative error at the nodes is 4% and the maximum relative error is 25% ; these mediocre results are the best obtained by method L-R for the case ai = 0 (i = 0,1, . . . , M). The large value of the maximum error is probably due to the discontinuity of u(x, t) at x = 1, t = 1 since the maximum error occurs near this point. Tables 2 and 4 show the results obtained for the initial function (11) and the inviscid equation (3) by methods F and L-R, respectively. In this case method F gives better results; Table 1. Discontinuous finite elements (method F. initial function (11)) Mean relative error
Fl
Ax
At
l/20 l/20 l/40
l/loo l/500 l/200
left 10.3% 7.8% 6.56%
right 12.4% 9.4% 7.76%
Computing time
26 set 65 set 93 set
78
P. Arminjon,
C. Beauchamp,
Table 2. Discontinuous
Continuous and discontinuous finite element methods
finite elements Mear relative
Ax = l/10
Ax = l/20
Ax= l/40
Ax = l/80
(method
F, initial function
Maximum
relative
for Burgers’ equation
11, cr = l/2)
Computing
time
At
error
error at t = 1
l/50
l/100 l/150 l/200
2.99% 3.39% 3.64% 3.79%
27% 34% 36% 37.8%
7.87 set 12.23 set 16.54 set 20.8 set
l/50
1.49%
13.3%
l/100 l/200
1.23% 1.35%
24.6% 31.12%
17.13 set 26.6 set 40.6 set
l/50 l/100 l/150 l/200
0.76% 0.61% 0.56%
13% 20.2% 24%
71.77 set 81.83 set 94.95 set
l/200
0.401%
13.72%
295
diverge
set
here we have used the approximate value defined by (19) for method F, while we only retained the left and right estimates uy_ and uy+ in table 1. Tables 3 and 5 show the corresponding results for the initial function (12). In this case method L-R is slightly superior to method F and appears to be more stable since we did not find a case where it would diverge. The accuracy is Table
3. Discontinuous
finite
elements
(method
F, initial
function
12, (Y =
l/2)
At
Mean relative error
Maximum relative error at t = 1
Computing time
Ax = l/10
l/50 l/100 l/150 l/200
1.96% 1.54% 1.75% 2.07%
3.23% 1.66% 2.79% 4.29%
8.8 14.7 19.7 24.6
Ax = l/20
l/50 l/100 l/l50 l/200
0.86% 0.67% 0.59%
diverge 1.41% 0.94% 0.70%
36.6 set 46 set 55 set
l/50 l/100 l/150 l/200
0.43%
diverge diverge diverge 0.67%
180 set
l/400
0.217%
0.328%
885 set
Ax = l/40
Ax = l/80
set set set set
P. Arminjon, C. Beauchamp, Continuous and discontinuous finite element methods for Burgers’ equation Table 4. Discontinuous finite elements (method L-R, OStIl,cwi=l,i=O,l,..., M)
At
Mean relative error right left
79
initial function (11)
Maximum relative error left right
Computing time
Ax = l/20
l/50 l/70 l/80 l/90 l/100 l/120 l/200
1.67% 1.93% 1.59% 1.97% 2.00% 2.04% 2.20%
1.61% 1.83% 1.52% 1.87% 1.90% 1.97% 2.15%
26.5% 25% 28% 26.4% 27.3% 28% 31.3%
25.2% 20.4% 27% 23% 24% 25% 29.7%
10.5 12.4 13.5 17.4 18.7 21.5 32.8
set set set set set set set
Ax = l/40
l/10 l/20 l/40 l/80 l/100 l/150 11200
4.8% 2.77% 1.65% 1.17% 1.11% 1.09% 1.09%
4.88% 2.81% 1.64% 1.14% 1.08% 1.05% 1.08%
28.63% 23.30% 21.24% 23.50% 25% 28.10% 30.20%
27.6% 21.3% 17.27% 18.10% 20% 24.15% 27.20%
7.88 set 11.45 set 17.80 set 28.30 set 33.20 set 47.82 set 62 set
Ax = l/50
l/150 l/200
0.90% 0.89%
0.88% 0.87%
26.9% 29%
22.30% 25%
58.7 set 76 set
Ax = l/80
l/160 l/320
0.65% 0.59%
0.64% 0.58%
24% 30.4%
19.6% 26.8%
96.9 set 187 set
comparable, but the computing times are better; by contrast, method F was often divergent for (12) except if At/Ax was chosen small. As mentioned earlier, the results for method L-R with cyi= 1 (tables 4-5) were obtained by solving system (31b), (31~) with Newton’s method, where we computed the Jacobian and its inverse at each iteration. To reduce the computing time, we applied Newton’s method with the same Jacobian matrix at each iteration. This led to the same results and reduced the computing time by approximately 35%. The results appear in table 6. 5.2. Continuous space-time finite elements for the inviscid equation For a comparison, we also applied to (3) with initial function (11) a continuous space-time finite element method presented in [l] for the viscous equation (1). Since v = 0 here we drop the boundary condition ~(1, t) = 0 and introduce a degree of freedom at x = 1 by considering a basis function qM = (xM - xM_$‘(x - xMP1)for xMel I x 5 xM. For Ax = l/20 and At = l/100 the mean relative error was 56%) for a computing time of 113 set, as compared to (1.23%) 26 set) and (2%) 18 set) for methods F and L-R respectively. The maximum errors were respectively 63%) 24% and 27%) showing the clear superiority of the discontinuous methods for the inviscid equation.
P. Arminjon,
80 Table
C. Beauchamp,
5. Discontinuous
finite
Continuous and discontinuous finite element methods for Burgers’ equation
elements
(method
L-R,
initial
function
(12),
OltSl,(Y,=l) Mean relative error left right
At
Maximum relative error at t = 1 left right
Computing time
Ax = l/10
l/20 l/40 l/50 l/80 l/120 l/160 l/200 l/1000
2.67% 1.80% 1.95% 2.63% 3.14% 3.44% 3.64% 4.37%
2.52% 1.42% 1.52% 2.24% 2.83% 3.18% 3.41% 4.3%
3.96% 1.27% 1.7% 2.49% 3.04% 3.36% 3.57% 4.4%
4.0% 1.19% 1.63% 2.42% 2.97% 3.29% 3.5% 4.3%
5.9 set 7.8 set 8.7 set 10.9 set 13.47 set 16.35 set 19.1 set 77.6 set
Ax = l/20
l/20 l/40 l/70 l/80 l/100 l/500
3.20% 1.53% 0.98% 0.96% 1 .O% 2.0%
3.24% 1.51% 0.88% 0.84% 0.86% 1.94%
5.17% 2.2% 0.90% 0.68% 0.76% 1.82%
5.19% 2.21% 0.91% 0.69% 0.75% 1.81%
7.84 set 11.6 set 16.35 set 17.6 set 20 set 77 set
Ax = l/40
l/70 l/80 l/l40 l/l60 l/200
0.987% 0.82% 0.52% 0.51% 0.51%
0.88% 0.83% 0.49% 0.47% 0.47%
0.9% 1.11% 0.49% 0.37% 0.35%
0.91% 1.15% 0.49% 0.37% 0.35%
13.9sec 30.6sec 47.16 set 52.4 set 62.5 set
Ax = l/80
l/160
0.43%
0.43%
0.59%
0.59%
101 set
S.3. Discontinuous elements for the viscous equation
For the full Burgers’ equation (1) with boundary conditions ~(0, t) = ~(1, t) = 0 for t 10 we cannot apply the L-R method, but we did apply method F with the same regular initial Table 6. Discontinuous Jacobian matrix)
Ax
I1 11 11 11 I1 12
l/20 l/40 l/80 l/80 l/160 l/40
At
l/80 l/160 l/160 l/200 l/200 l/160
finite
elements
Mean relative error
(method
L-R,
(Y,= 1 with
Maximum relative error at t = 1
left
right
left
right
1.95% 1.09% 0.66% 0.62% 0.45% 0.51%
1.85% 1.06% 0.65% 0.61% 0.45% 0.47%
25.8% 28.6% 24.6% 26% 22% 0.37%
21.8% 24.9% 19.6% 21.5% 18.7 % 0.37%
constant
Computing time
12.8 set 36.7 set 69.4 set 84 set 162 set 39 set
P. Aminjon,
C. Beauchamp,
Continuous
and discontinuous
finite element methods for Burgers’ equation
81
Table 7. Errors Discontinuous (method At
a
F
l/100 l/200 l/500 l/1000
1.5 0.3 0.3 0.5
39% 19% 6.2% 3.015%
elements F) Emax
Continuous space-time elements Computing
59.4% 20% 9.7% 4.95%
26 46 70 107
time
set set set set
c 0.78% 0.49% 0.84% Note: reached 0.085
Computing
for the
time
2 set 3.5 set 5 set At = l/140, E minimum value
function as in [l]. With v = 1 and Ax = l/20 we obtain the results in table 7 for the mean and maximum relative errors C, E,,,. Here we chose to approximate ul by
my- + (1 - CX))u~+.
(38)
These results show the superiority of continuous space-time elements for the viscous equation. This is in agreement with the intuitive observation that regular functions should be better approximable by regular approximating functions.
Fig. 4. Characteristics
for problem
with discontinuous
initial function.
82
P. Arminjon,
C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
5.4. Discontinuous elements for the inviscid equation with the discontinuous initial function (13)
The initial function (13) (see (36)) has a jump at x = l/2, and we can determine the solution only in some regions of the (x, t)-plane. There is only one characteristic line going through each point of regions R1 and R3 (fig. 4) so that the solution u(x, t) can be uniquely determined at these points. But for region RZ each point (x, t) is the intersection of two characteristic lines, one originating in the interval 0 < x < l/2, the other in l/2 < x < 1. We find
j-$$ for uW)=
(x,t)ER1,
l_x
l_t
for
(x,t)ER3.
(39)
:
We computed the solution with method F and with method L-R with ai = 1. In both cases the numerical solution is well-behaved in regions R1 and R3 and converges to the value given by (39). Figs. 5 and 6 show the graphs of the numerical results obtained for u(x, t) at time t = 0.3 set
J(X,O.3)
Fig. 5. Behaviour
of u(x, 0.3) (method F, Ax = l/40, At = l/200).
P. Arminjon, C. Beauchamp,
83
Continuous and discontinuous finite element methods for Burgers’ equation
0.9
0.5
x 0.5
Fig. 6. Behaviour
1
of u(x, 0.3) (method L-R, Ax = l/40, At = l/160).
by methods F and L-R, respectively, and for Ax = l/40, with At = l/200 (method F), At = l/160 (method L-R). The mean relative error at points where the solution is known is 0.28% for method F and 0.35% for method L-R; the corresponding computing times were 47 set and 13 sec. 6. Concluding remarks Our results indicate that although the two discontinuous finite element methods presented here do not have the same accuracy and efficiency as other space-time finite element methods for continuous solutions of the viscous equation (see [l], [2]), they seem to be superior in the case of the ill-behaved solutions of the inviscid equation. This should make them well suited to problems in fluid mechanics where the viscosity can be neglected, and problems associated with the system of primitive equations of dynamic meteorology. We are presently testing these methods for some boundary layer problems where the viscosity is very small. References [l] P. Arminjon and C. Beauchamp, A finite element method for Burgers’ equation in hydrodynamics, Numer. Meths. Eng. 12 (1978) 415-428. [2] P. Arminjon and C. Beauchamp, Numerical solution of Burgers’ equation in two space dimensions, Meths. Appl. Mech. Eng. 19 (1979) 351-365.
Int. J. Comp.
84
P. Arminjon, C. Beauchamp,
Continuous and discontinuous finite element methods for Burgers’ equation
[3] P. Anninjon, Finite element analysis of moving loads on a floating ice shell, Topics in Numerical Analysis, Proceedings 3rd Conference Royal Irish Academy (Academic Press, London, 1977) 299-314. [4] R. Bonnerot and P. Jamet, A second order finite element method for the one-dimensional Stefan problem, Int. J. Numer. Meths. Eng. 8 (1974) 811-820. [5] A. Chalifour, Analyse dune methode d’elements finis discontinus pour les equations differentielles lineaires, Report no. CRM-643 (Centre de Recherches Mathematiques, Univ. Montreal, 1976). [6] M. Delfour and F. Trochu, Discontinuous approximation of ordinary differential equations and applications to optimal control problems, Report no. CRM-751 (Centre de Recherches Mathematiques, Univ. Montreal, 1977). [7] M. Delfour, W.W. Hager and F. Trochu, Discontinuous piecewise polynomial Galerkin methods for ordinary differential equations, Report no. CRM-830 (Centre de Recherches Mathematiques, Univ. Montreal, 1978). [8] M. Fortin, Resolution numerique des equations de Navier-Stokes par des elements finis de type mixte, Rapport IRIA-LABORIA 76-184 (1976). [9] M. Fortin, Mixed finite element methods for incompressible flow problems, J. Comp. Phys. 31 (1979). [lo] P. Jamet and R. Bonnerot, Numerical solution of the Eulerian equations of compressible flow by a finite element method which follows the free boundary and the interfaces, J. Comp. Phys. 18 (1975) 21-45. [ll] P. Lesaint and P.A. Raviart, On a finite element method for solving the neutron transport equation, Proceedings of the Symposium Mathematical Aspects of Finite Elements in Partial Differential Equations, Madison, Wisconsin (Academic Press, New York, 1974). [12] P. Lesaint, Continuous and discontinuous finite element methods for solving the transport equation, The mathematics of finite elements and applic. II, Proc. 2nd Brunel Univ. Conf. Inst. Math. Appl., Uxbridge, 1975 (Academic Press, London, 1976) 151-161. [13] J.T. Oden and J.N. Reddy, An introduction to the mathematical theory of finite elements (Wiley, New York, 1976). [14] H.H Rachford and M.F. Wheeler, An H-’ Galerkin procedure for the two-point boundary value problem, in: C. de Boor (ed.), Mathematical aspects of finite element methods in partial differential equations (Academic Press, New York, 1975) 353-382. [15] W.H. Reed and T.R. Hill, Triangular mesh methods for the neutron transport equation, LA UR 73479 (Los Alamos Lab., 1973). [16] L.C. Wellford and J.T. Oden, Discontinuous finite element Galerkin approximations of shock waves in nonlinear elastic solids, TICOM Report 74-6 (The Texas Institute for Computational Mechanics, Univ. Texas at Austin, 1974). [17] L.C. Wellford and J.T. Oden, Discontinuous finite element approximations for the analysis of shock waves in nonlinear elastic materials, J. Comp. Phys. 19 (1975). [18] L.C. Wellford and J.T. Oden, A theory of discontinuous finite element Galerkin approximations of shock waves in nonlinear elastic solids -Part 1. Variational theory, Comp. Meths. Appl. Mech. Eng. 8 (1976) 1-16. [19] O.C. Zienkiewicz, The finite element method in engineering science (McGraw-Hill, London, 1971). [20] R. Courant and D. Hilbert, Methods of mathematical physics (Interscience, Wiley, New York, 1962). [21] P. Arminjon, Numerical solution of a boundary layer problem with continuous and discontinuous finite elements. To appear.