- Email: [email protected]
Primal-dual variational problems by boundary and finite elements
Primal-dual variational problems by boundary and finite elements F. Sgallari Dipartimento di Matematica, Universit; di Bologna, Italy (Received August 1985; revised February 1985)
Via Vallescura
2, 40136 Bologna,
In this paper the primal-dual (or mixed) formulation is studied for selfadjoint elliptic problems coupled with a boundary integral equation. It is shown that, after introducing a suitable complementary variational principle, the problem is reduced to finding a stationarity point of a constrained functional. Some numerical examples are reported for a secondorder differential equation on unbounded domains. Key words: mathematical models, mixed variational elements, boundary elements
Increasing interest has recently been shown in so-called primal-dual (or mixed) variational formulations and the related mixed finite element methods.’ Moreover, in many engineering problems, it is convenient to combine boundary elements with finite element procedures.*-‘j Such a situation may arise when using finite elements in part of the domain and boundary elements in another part, and in ‘exterior’ problems. Classical ‘complementary’ variational principles’ provide an already established theoretical framework for the saddle-point or Lagrangian formulation of elliptic self-adjoint problems. In general, it is desirable to obtain a linking procedure between boundary and mixed finite elements, which will provide a symmetric (but in general non-definite) matrix and be capable of standard finite element assembly in conventional programs.8Y9 In this paper, the general case of self-adjoint elliptic problems is discussed. It is shown that, if a suitable complementary variational principle is introduced, the differential problem coupled with a boundary integral equation can be reduced to finding a stationarity point of a constrained functional. As far as numerical approximations are concerned, papers by BabuEka, r” Brezzi” and Nedelec12 have provided the basic convergence conditions. In this paper, some numerical examples are reported and discussed for a second-order differential equation related to a two-dimensional potential flow in an unbounded domain.
Variational principles for elliptic problems In the following, some basic definitions and results concerning elliptic self-adjoint problems are given. For a more detailed discussion, refer to Arthurs7 and Ciarlet.r3
246
Appl. Math. Modelling, 1985, Vol. 9, August
problems, finite
Let T denote a linear operator from one real Hilbert space (I’, (. , .)) to another space (IV, (. , .)). Thus one writes: TEL@,
W)
The adjoint operator
T$ = u T* EL(W,
(1)
V) is such that:
(u, T$) = (T*u, 6) + S(u, G) where S(u, 4) denotes boundary terms. An important class of boundary-value considered, described by the equations:
A@=f
(2) problems is now
in Q
with A = T*T an elliptic differential self-adjoint subject to the Dirichlet-Neumann conditions: o@ =
@$‘B
o*T@ + /!%#J = (s*uB
(3) operator,
on I?r
(4)
on I2
(5)
Here Q denotes a bounded open domain of II? with sufficiently regular boundary r = I’1 U I’* and u, u* are adjoint trace operators defined by: u: tr(V) -+ tr(W) (W o@)r = co*u, G)t-
(6)
Under the hypothesis that (@,A$) 2 k11@f, k E RR+, for any 4 # 0, the solution of equations (3)-(5) is unique and may be determined as the minimum of the functional: J(Q) = ;(TG, W - Cf,@ + it@,l%Vr,-(o*UB@)r,
(7)
with the constraint u$ = o@Bon Fr. This latter constraint is essential and needs to be explicitly imposed. 0307-904X/85/042246-07/$03.00 0 1985 Butterworth & Co. (Publishers)
Ltd
Primal-dual
If the ‘dual’ variable u = T@is of interest, the solution is found as the maximum of: I(u) = - 3@, u) + (u, u&)r, - ;(p-‘o*(u
- ug), a*@ -u&r,
W)
rr
on r
t=
x J/(E) dy(x) d-M
where F(E) represents the (known) contributions and u*T&. Denoting equation (14) as:
N[t4, ti, = (F, ti)
(14) of u&
Qll,
(15)
where the linear operator K is defined as:
domains
{[.I u*TG - Gu*T[.] )dy(x)
(16)
r
In many applications, it is common to encounter elliptic problems in unbounded domains, with suitable regularity conditions at infinity. Usually, the whole solution domain is divided into a bounded domain ai and an outer domain &, extending to infinity. The conditions at infinity are then brought back to the finite boundary r, dividing ai and S&. In the following, the discussion is restricted to this class of problem; however, the same approach can be used whenever it is necessary to deal with a differential problem coupled with an integral equation. In the inner domain al, the problem is:
u$ = undefined
integral
[qbu*TG - Gu*T$]
65) dy(E) +
K[.] = [.I -1
in ai
(13)
r
r
(9)
r
After substitution of (13) into (12) a boundary equation is obtained for:
=
It is well-known’ that the Lagrangian multiplier h can be identified with the ‘primal’ variable $ in Q (thus with 4 on I’a, if X is taken as regular as is needed here).
A$=T*T@=f
on
u*T# = u*T#, + u”T$i
+ (u*& - Q)> @jr,
Elliptic problems in unbounded
(12)
where $($) is a test function. Now, imposing the linking conditions:
h)
+ 104%&or,
F. Sgallari
w
d-i(x)d-r(t)
= - +(u, u) + (T*u, h) - (f, h) + (u, uhh,
problems:
(8)
with the constraint T*u = f in 52 (if /3E 0, u*u = u*ug on r,>. To remove the constraint T*u -f = 0, an appropriate space of Lagrange multipliersM and a Lagrangian d: : W x M + CRis considered in such a way that the unknown u is obtained as the first argument of the saddle-point (u, X) of the Lagrangian over the space IVx M. In the present case, it turns out that one may choose: ,c(u, A) = I(u) + (T*u -f,
variational
(10)
Constrained variational problems In the following, the problem of Figure 1 is considered, that is: A$ = T*Tq5 = f
in a1
u*T@ = 0
on
K[@l=F
on
r, rz
(17)
(18) (19)
For the sake of simplicity, essential homogeneous conditions have been chosen on rr ; otherwise, functional (9) has to be considered, where essential boundary conditions are to be explicitly imposed. The basic idea of the variational formulation is as follows: find the stationarity point
Obviously boundary conditions on r are needed. The solution in an2 is unknown, but it is possible to write a boundary integral equation given a fundamental solution G of the differential problem. The usual theory of fundamental solutions of partial differential equations, leading to the boundary integral formulation of the problem, must now be stated in terms of distributions. A complete exposition of distribution theory can be found in Guelfand and Chilov;r4 for boundary equations, see Alliney and Tralli.” Let @Edenote the solution in S& and assume: @E=@k+&
(11)
where @kis an assigned function which satisfies the first equation of (10); $i is infinitesimal at infinity together with its derivatives up to the (2m - 1)th order if A is of order 2m;16 f IS ’ assumed here with a support bounded in al. Therefore, the integral equation for & will involve only the values of the function on the boundary (contributions at infinity do actually vanish, due to the asymptotic behaviours of both G and @ii).Thus:
I r
@i(S)ti(t) WC-9=
Is
[h u*TG - Gu*T@I G(f)
rr
Figure
7
Domain ~2, of boundary
Appl.
value problem
Math. Modelling, 1985, Vol. 9, August
247
Primal-dual variational problems:
F. Sgallari
of the Lagrangian functional ,C(u, h) associated with the differential operator A, subject to the constraint given by the boundary integial equation. It should be noted that the problem in s2, has no boundary conditions explicitly prescribed on rZ. It follows that L(u, X) is simply: ‘c(u, h) = - &,
u) + (T*u -f,
X)
w9r,l
+ :[(o*u> @jr, + (u,
(20)
In the following previous results by Brezzi,“.17’18 Babushka” and Nedelec12 are used. First, details are required of the functional spaces related to the problem. Given an elliptic operator A of order 2m, the solution ti, h: of the variational problem for L(u, h) is to be found in the spaces: X = H(T*; L&) = {u E (L2(fll)r;
T*u EL’(&)}
(21)
tation of the boundary (further details in Ciarlet13 and Johnson and Nedelec19). Now, a set of basis functions u&x), hk(x) is chosen such that: VN(s2,) = Span {uk(x);
k = 1,2, . . , N} C X
(28)
WL(s21)=Span{hk(x);
k=
(29)
As is well-known, tions U&Z), hk(x) appropriate order with local support. described as:
1,2,...,L}CM
in numerical applications are real-valued piecewise but, here, not necessarily Thus, any function u E
the basis funcpolynomials (of continuous) X, h EM is
N
u(x) = c Y&(X) k=l
and
(30) M = L2(3 1)
(22)
Note that functions belonging to X do not need to be continuous, provided that a generalized Green’s theorem holds. As has already been pointed out, the boundary integral equation is to be considered in the weak sense; so equation (14) is considered which can be rewritten as: (23) where:
h(x) *
4
bkhk(X)
k=l
where the meaning of the parameters Yk, bk is strictly related to the finite element implementation. Some attention must be devoted to the approximation of the ‘trace’ spaces on the boundary. Recall that it is necessary to represent the contribution of the boundary terms of K1 and K2. One considers: T(r2) = tr(VN(al))
K,[-]
= [.I -
[.I
o*TG dy(x)
= Span{Uk(x)lrz;
k= 1,2,. . . ,J}
(31)
(24) S(r,)
I.2
= Span {@k(x); k = 1, 2,
. , I} C L2(r2)
(32)
Thus, one has:
and
J
K,[.] =
Ga*[.] dy(x)
(25)
s
u(x) N c YkUk(X) k=l
x E r2
rz
It should be noted that in the present ‘weak’ formulation q51r, does not correspond to hlr, (the ‘trace’ theorem fails); its meaning is simply that of an additional Lagrange multiplier. On the contrary, in the ‘strong’ formulation matching with the domain variable makes sense. As far as u is concerned, in equation (25) one can take the restriction u lr2. It is easy to see that both 4 and $ require only L2 summability on r2. Now, the problem (17)-(19) may be stated in abstract variational terms: find an element ti E X and an element fi EM such that:
L(ri,h)
VuEX
VhEM
(26)
ti, + W2b41, $, = UC $) V$ E L2(S12)
(27)
Under previous assumptions on the operator A and with the appropriate spaces X and M, the necessary and sufficient conditions (see Brezzill and Babufia” for (26) and Nedelec12 for (27)) for existence and uniqueness of the solution (ti, i) are satisfied.
Discrete approximation of the variational problem To find an approximate numerical solution, the discrete analoaue of variational problem (26) is considered. First, suppose a discretization over the-problem domain a1 has been established. For simplicity, it is assumed that the discretization does not introduce any error into the represen-
248
Appl.
Math.
Modelling,
1985,
Vol.
9, August
@(x) = i k=l
(33)
uk@k(X)
According to the previous remarks, functions G(x) themselves can be identified with @k(x). Using equations (30) and (33) one can rewrite functional (20) and the constraint (27) in compact form: ~~,b,a)“--fyT~Y+YTDb--BTf + &@%
+2-y]
(_T-M)a+_Ny=F
(34) (35)
Refer to the Appendix for detailed definitions of matrices and vectors. The constraint equation (35) can be solved2’ with respect to (I, to give: (36)
a=g+_4, where the following positions are introduced: g=(_T-MM)-‘F
(37)
p=-(T--_M)_‘_N
After substitution of (36) into the expression of the objective functional, one obtains: L(y, b) = - $y’cy + HYTST4
+ yTpb - bTf +yTsTg
+YTETsY +sTsyl or equivalently:
(38)
Primal-dual
+ &lyTzTg+gT$y]
-bTf
(39)
We have now an unconstrained saddle-point problem for L’b, b). Its solution is determined by the linear system: -Ey+a,+Db=-&S’g
(40)
DTb =f
problems:
F. Sgallari
Conditions on r2 are obtained by a boundary integral equation with fundamental solution C(x,, yo; x, u) = log(2R/r) + log(2R/r’), r = [(x -x0)’ + (y -yo)2]1’2, r’ = [(x - xo)2 + 0, + yo)2]1’2 and R is the radius of the circle r2. For a detailed exposition refer to Alliney.” Using the previous notation, one has:
qy,b)“-$yT~y+$yT[&~Tf+PTS)]y -+y’gb
variational
T* = - div
T = grad thus:
where the N x N matrix:
u=n
u = V@
(41)
u* = n
and :
resulting from the introduction of the boundary integral equation (considered as a constraint for the functional) turns out to be symmetric.
X=H(div;~;2)={uE(L2(~;2))2;divuEL2(~1)}
Example The procedure outlined in the previous sections, is well suited to the study of steady flows of ideal fluids past completely reflecting obstacles. It is assumed that these flows are irrotational and that there are no sources and sinks. Under the hypothesis stated, the determination of a flow in a region a can be reduced to determining a realvalued potential function $(x,y), such that u = V@. With reference to Figure 2, in the unbounded plane domain s2 the function Q satisfies Laplace’s equation: -V2$
inSZ=a1Ufi2
= 0
Functional
(20) can be rewritten
as:
~(~,A)=-illlu12dS2,-Il(divu)hdn, n,
a, -t$
(42)
where the constraint
5 rz
[u.n6+@uun]dy
(47)
is:
At the boundary rr (assuming that there is no absorption of the fluid through the walls) the flow is in the tangential direction; Thus, the boundary condition is:
a$
-0
on l?,
(43)
an-
In addition, when one determines the flow past an obstacle, one imposes the condition that the flow is homogeneous at infinity. This means that the potential @ has the following asymptotic behaviour: ~(x,Y)=~k(x,Y)+~i(x,Y)
+
rz rz
=
(44)
r-++m
(45)
IV&I
r++m
(46)
= W/r2)
with r = (x2 + _Y~)“~.
-Y
x
Figure 2
Definition
sketch for test problem
s
K9 9(t) WC;)
(48)
r2
Here, Gk(x, v) denotes an assigned function, which satisfies equation (42); physically, it represents the potential of the incident flow. The function &(x,v) is the perturbation potential due to the obstacle, and it is required that: $J&,Y) = 0(1/r)
W-G8 ; 6) d4.9d-r(x)d%9
The major difficulty, when using a dual, or primal-dual formulation, is to take appropriately into account the constraint T*u - f = 0, in this case div u -t-f = 0. There are essentially three ways of circumventing this difficulty. First, the constraint is approximated, in such a way that the discrete solution ii satisfies a relation of the form diva + f= 0 in !ZZ,where fis a typical finite element approximation of the function 5 By reference to elasticity problems, such methods, which are directly based on the dual formulation, are known as equilibrium methods.13 While the equilibrium methods are, by definition, based on a formulation where there is only one unknown (the gradient V@= u for the problem here, the stress tensor u for the elasticity problem, etc.), one may use the techniques of duality theory to get rid ot the constraint, a process which results in the addition of a second unknown, the Lagrangian multiplier. For a detailed exposition of such methods refer to Ciarlet13 and Mercier.21 For the discrete approximation of X, the mixed finite element method implementation of Johnson and Mercier” was used. Consider a triangulationyj over the set 52, i.e. the set Sz is subdivided into a finite number of triangles K, the index h representing the maximum diameter of the element K. For 2 and x the expansions:
Appl.
Math. Modelling,
1985,
Vol.
9, August
249
Primal-dual
variational
problems:
F. Sgallari
The numerical model presented above has been implemented on the VAX 1 l/780 computer of the Computing Centre at the Faculty of Engineering of Bologna University. To assess the numerical performance of the method, first consider a half-circular obstacle centred at the coordinate origin (Figure 4 gives the finite element mesh). For this particular case, the analytical solution is well known:
Figure 3 Triangular element for mixed implementation discussed in text. Degrees-of-freedom for dual variable are normal components of u at points on edges
6=
il Y&(X> Y>
x= i
(49)
bjXj(X,Y)
=-vx( 1+-&)
(56)
where --z, is the unperturbed incoming velocity and a is the radius of the obstacle. In Table 1 the results are presented for a = 1 and z1= 50. It should also be noted that with a coarse mesh there is excellent agreement with the exact solution. In the second example the obstacle in Figure 5 is depicted, with z, = 4 m/s. Figure 6 holds for any element a straight line segment of arbitrary length, which is tangential to the streamline in the centroid of the triangle.
j=1
were used, where the basis functions ui have local support on two adjacent triangles and $ have support on one triangle. Moreover, ui have a continuous normal component on the common edge of the two adjacent triangles and vanishing normal components on the other edges; thus ui are linear vector-valued functions on each triangle where they satisfy the ‘equilibrium’ equation divui .= g, for g piecewise constant. As degrees-of-freedom the normal components on six points were chosen (two for each triangle edge); at a single point ui. n = f 1 is imposed, with ui. n = 0 on the others. For hi consider constant functions on each triangle (corresponding to one degree-of-freedom). For the practical construction of the basis functions refer to Thomassetz3 and Johnson and Mercier22 where the implementation of the dual finite element is discussed for the more general case of elasticity problems. As far as error estimates are concerned one has: l/11- ii llL2 < c-h2
(51)
IIh - KllL2 < ch
(52)
40
30
20
IO
0 40 Figure 4
20
30
IO
00
20
IO
Y
30
Finite element
mesh for test example
Finite element
mesh for second obstacle problem
40
where c is a constant; for more details refer again to Ciarlet l3 Mercier,‘l Thomasset and Johnson and Mercier.22 The diicrete functional (47) and the constraint (48) become: 0 I(
.X(y,b,a)=-_yTEy+yTgb
Figure 5
+ + LvTsra + a’$y] (T-ma++
=F
(53) (54)
The constraint equation (54) can be solved with respect to a and, after substitution into (53) the objective functional LQ, b) is obtained as in (38). The solution is now determined by the linear system: --_Ey+_@+Db=-$sTg
-
gTb=Q
(55)
It is worth noting that the coefficient matrix is non-definite; however, one can build it up element-by-element using the standard assembling techniques and using a new iterative process9 that requires only minor modifications on Irons’ frontal solution program.8
250
Appl.
Math.
Modelling,
1985, Vol. 9, August
80 60
40 20 0 100
0.0
60
40
2.0
00
20
Figure 6 Straight length segments of arbitrary streamlines in centroid of each triangle
40
60
80
length tangent
10.0 to
Primal-dual
variational
problems:
F. Sgallari
Table 1
Element
x, y coordinates
Computed
1
1.491
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
1.236 1.191 1.368 1.168 0.877 0.703 0.688 0.313 0.114 2.237 1.978 1.900 2.047 1.762 1.415 1.110 1.020 0.484 0.199 2.983 2.720 2.609 2.726 2.356 1.953 1.517 1.351 0.655 0.285 3.483 3.713 3.534 3.159 2.789 2.708 2.009 1.522 0.826 0.456
-107.477 - 101.074 -99.155 -99.155 -85.122 - 72.702 -57.617 -49.846 -23.113 -9.877 -132.701 -122.715 -118.346 - 122.436 - 105.989 -88.816 -69.490 -61.211 -29.277 - 12.654 -162.120 -152.792 -146.814 - 151.875 -131.961 -111.360 -86.091 - 75.434 -36.468 -16.388 -185.199 -188.334 - 192.443 - 174.888 - 150.744 -149.319 - 106.059 -85.814 -41.946 - 26.235
39
40
0.101 0.159 0.325 0.559 0.905 0.868 1.015 1.308 1.444 1.229 0.144 0.246 0.537 0.853 1.346 1.376 1.633 1.969 2.164 1.964 0.188 0.332 0.749 1.146 1.788 1.884 2.251 2.630 2.884 2.699 0.188 0.419 1.086 1.396 2.038 2.517 3.104 3.100 3.354 3.669
Exact
References 1
2
3 4 5 6 I 8 9 10 11 12 13 14 15 16 17
Atluri, S. N., Gallagher, R. H. and Zienkiewicz, 0. C. (eds), ‘Hybrid and mixed finite element methods’, John Wiley, New York, 1983 Brebbia, C. A., Telles, J. and Wrobel, L. ‘Boundary element methods - theory and application in engineering’, SpringerVerlag, Berlin, 1983 Mei, C. C. and Cheu, H. S. Int. J. Num. Meth. Eng. 1976, 10, 1153 Yue, D. K. P., Cheu, H. S. and Mei, C. C. Znt. J. Num. Meth. Eng. 1978,12,245 Zienkiewicz efal. Int. J. Num. Meth. Eng. 1977, 11, 355 Allinev. S. Aunl. Math. Model& 1982, 6. 424 Arthur;, A. M: ‘Complementary &riational principles’, 2nd edn, Clarendon Press, Oxford, 1980 Irons, B. M. Int. J. Num. Meth. Eng. 1970, 2, 5 Alliney, S. ‘Iterative solution of large linear systems arising from saddle-point problems’ (submitted) Babuxka, I. Math. Comp. 1973,27, 221 Brezzi, F. R.A.I.R.O. Sk, Rouge 1974,3, 129 Nedelec, J. C. ComputerMeth. Appl. Mech. Eng. 1976,8,61 Ciarlet, P. G. ‘The finite element method for elliptic problems’, North-Holland, Amsterdam, 1978 Guelfand, L. M. and Chilov, C. E. ‘Les distributions’, Vol. 1, Dunod, Paris, 1972 Alliney, S. and Tralli, A. Appl. Math. Modelling 1984, 8, 75 Lions, J. L. and Magenes, E. ‘Problemes aux limites non homogenes et applications’, Vol. 1, Dunod, Paris, 1968 Brezzi, F., Johnson, C. and Mercier, B. Math. Comp. 1977, 31, 809
- 107.934 - 101.594 -98.617 -99.721 -85.180 - 72.666 -58.208 -50.174 -22.848 -9.436 -134.121 -123.811 -119.377 -123.185 -106.037 -88.939 -69.737 -61.376 -29.150 -12.533 - 165.869 -154.149 -148.163 -151.920 -131.276 - 110.939 -86.143 - 75.320 -36.523 -16.184 - 188.488 - 198.959 -189.657 -171.226 -151.148 -145.307 -107.819 -82.524 -44.790 -24.468
2.560 8.034 17.405 16.566 21.753 33.653 30.019 19.090 9.758 8.039 0.736 2.556 6.151 6.032 8.838 12.281 11.462 7.953 4.404 2.829 0.984 1.633 2.791 3.690 4.853 6.248 6.866 4.098 3.039 0.989 1.403 3.213 10.123 2.014 7.383 5.141 5.838 4.171 8.624 2.537
-26.900 -16.681 -20.611 -33.730 - 43.890 -47.120 -61.910 -63.971 - 71.357 -81.498 -38.505 -36.647 -38.594 -42.734 -47.485 -49.500 -55.089 -56.156 -59.054 -62.426 -40.928 -41.430 -43.272 -47.427 -49.403 -50.173 -53.154 -53.683 -54.539 -57.186 -36.388 -36.125 -52.087 -56.793 -47.588 -57.657 -43.402 -64.174 -59.124 -48.718
18
19 20 21 22 23
-27.821 -18.858 -21.788 -33.678 -44.276 -49.661 -61.526 -62.954 - 70.820 -82.214 -40.135 -37.804 -39.080 -42.843 -47.329 -49.642 -54.718 -55.863 -59.192 -62.570 -44.149 -43.141 -44.249 -46.002 -48.462 -49.756 -52.547 -53.327 -55.153 -56.634 -45.915 -46.509 -46.974 -47.179 -48.727 -49.734 -51.497 -52.561 -53.709 -53.545
3.027 8.155 16.678 16.566 22.150 32.771 30.679 18.857 9.482 7.024 1.281 3.080 6.716 7.215 9.805 12.814 11.919 8.300 4.332 2.576 0.752 1.604 3.600 4.083 5.503 6.780 6.288 4.646 2.469 1.416 0.442 0.799 2.053 3.098 3.991 3.646 3.336 3.315 1.946 0.895
Brezzi, F. and Raviart, P. A. ‘Mixed finite element methods for 4th order elliptic equations’, in ‘Topics in numerical analysis’ (ed. Miller, J. J. H.), Academic Press, London, 1976 Johnson, C. and Nedelec, J. C. Math. Comp. 1980,3.5, 1063 Albney, S. Appl. Math. Modelling 1982,6, 291 Mercier, B. ‘Topics in finite element solution of elliptic problems’, Springer-Verlag, Berlin, 1979 Johnson, C. and Mercier, B. Numer. Math. 1978,30, 103 Thomasset, F. ‘Implementation of finite element methods for Navier-Stokes equations’, Springer-Verlag, Berlin, 1981
Appendix The following are the definitions of matrices and vectors of the section covering the discrete approximation of the variational problem. One has:
E = [EiiI NxN Eii
= Eji
=
s
ui(x>q(x) d%(x)
a, D = Piil Dii =
Appl.
NxL
Xi(X) T*Ui(X)
Math.
(AlI
.
dCLl(X)
Modelling,
1985,
(43
Vol. 9, August
251
Primal-dual variational problems: F. Sgallari
As far as the boundary is concerned, the matrices --S, T, M, & with dimensions Ix J, Ix I, Ix I, I x J, respectively, are used, with elements: Sii =
=
s
$i (X)
s
@i
CJ*Ui (X)
dy(X)
Moreover, the following vectors were introduced:
(X)uj (XI dy(x)
6431
rz
Tii = Tji =
fi s r*
@i(X) $j (X> d+&)
Mii zzMji = ss rz rz
(A41
Appl.
Math.
I
f(x)
AjW d%(x)
(A7)
a2, F= [F1, . . ..FI]’
@i(x) o*TG(x,
[I @j(E) d?(x) dy(O
G48) WI
252
=
Modelling,
1985,
Vol.
9, August
Via Vallescura
2, 40136 Bologna,
In this paper the primal-dual (or mixed) formulation is studied for selfadjoint elliptic problems coupled with a boundary integral equation. It is shown that, after introducing a suitable complementary variational principle, the problem is reduced to finding a stationarity point of a constrained functional. Some numerical examples are reported for a secondorder differential equation on unbounded domains. Key words: mathematical models, mixed variational elements, boundary elements
Increasing interest has recently been shown in so-called primal-dual (or mixed) variational formulations and the related mixed finite element methods.’ Moreover, in many engineering problems, it is convenient to combine boundary elements with finite element procedures.*-‘j Such a situation may arise when using finite elements in part of the domain and boundary elements in another part, and in ‘exterior’ problems. Classical ‘complementary’ variational principles’ provide an already established theoretical framework for the saddle-point or Lagrangian formulation of elliptic self-adjoint problems. In general, it is desirable to obtain a linking procedure between boundary and mixed finite elements, which will provide a symmetric (but in general non-definite) matrix and be capable of standard finite element assembly in conventional programs.8Y9 In this paper, the general case of self-adjoint elliptic problems is discussed. It is shown that, if a suitable complementary variational principle is introduced, the differential problem coupled with a boundary integral equation can be reduced to finding a stationarity point of a constrained functional. As far as numerical approximations are concerned, papers by BabuEka, r” Brezzi” and Nedelec12 have provided the basic convergence conditions. In this paper, some numerical examples are reported and discussed for a second-order differential equation related to a two-dimensional potential flow in an unbounded domain.
Variational principles for elliptic problems In the following, some basic definitions and results concerning elliptic self-adjoint problems are given. For a more detailed discussion, refer to Arthurs7 and Ciarlet.r3
246
Appl. Math. Modelling, 1985, Vol. 9, August
problems, finite
Let T denote a linear operator from one real Hilbert space (I’, (. , .)) to another space (IV, (. , .)). Thus one writes: TEL@,
W)
The adjoint operator
T$ = u T* EL(W,
(1)
V) is such that:
(u, T$) = (T*u, 6) + S(u, G) where S(u, 4) denotes boundary terms. An important class of boundary-value considered, described by the equations:
A@=f
(2) problems is now
in Q
with A = T*T an elliptic differential self-adjoint subject to the Dirichlet-Neumann conditions: o@ =
@$‘B
o*T@ + /!%#J = (s*uB
(3) operator,
on I?r
(4)
on I2
(5)
Here Q denotes a bounded open domain of II? with sufficiently regular boundary r = I’1 U I’* and u, u* are adjoint trace operators defined by: u: tr(V) -+ tr(W) (W o@)r = co*u, G)t-
(6)
Under the hypothesis that (@,A$) 2 k11@f, k E RR+, for any 4 # 0, the solution of equations (3)-(5) is unique and may be determined as the minimum of the functional: J(Q) = ;(TG, W - Cf,@ + it@,l%Vr,-(o*UB@)r,
(7)
with the constraint u$ = o@Bon Fr. This latter constraint is essential and needs to be explicitly imposed. 0307-904X/85/042246-07/$03.00 0 1985 Butterworth & Co. (Publishers)
Ltd
Primal-dual
If the ‘dual’ variable u = T@is of interest, the solution is found as the maximum of: I(u) = - 3@, u) + (u, u&)r, - ;(p-‘o*(u
- ug), a*@ -u&r,
W)
rr
on r
t=
x J/(E) dy(x) d-M
where F(E) represents the (known) contributions and u*T&. Denoting equation (14) as:
N[t4, ti, = (F, ti)
(14) of u&
Qll,
(15)
where the linear operator K is defined as:
domains
{[.I u*TG - Gu*T[.] )dy(x)
(16)
r
In many applications, it is common to encounter elliptic problems in unbounded domains, with suitable regularity conditions at infinity. Usually, the whole solution domain is divided into a bounded domain ai and an outer domain &, extending to infinity. The conditions at infinity are then brought back to the finite boundary r, dividing ai and S&. In the following, the discussion is restricted to this class of problem; however, the same approach can be used whenever it is necessary to deal with a differential problem coupled with an integral equation. In the inner domain al, the problem is:
u$ = undefined
integral
[qbu*TG - Gu*T$]
65) dy(E) +
K[.] = [.I -1
in ai
(13)
r
r
(9)
r
After substitution of (13) into (12) a boundary equation is obtained for:
=
It is well-known’ that the Lagrangian multiplier h can be identified with the ‘primal’ variable $ in Q (thus with 4 on I’a, if X is taken as regular as is needed here).
A$=T*T@=f
on
u*T# = u*T#, + u”T$i
+ (u*& - Q)> @jr,
Elliptic problems in unbounded
(12)
where $($) is a test function. Now, imposing the linking conditions:
h)
+ 104%&or,
F. Sgallari
w
d-i(x)d-r(t)
= - +(u, u) + (T*u, h) - (f, h) + (u, uhh,
problems:
(8)
with the constraint T*u = f in 52 (if /3E 0, u*u = u*ug on r,>. To remove the constraint T*u -f = 0, an appropriate space of Lagrange multipliersM and a Lagrangian d: : W x M + CRis considered in such a way that the unknown u is obtained as the first argument of the saddle-point (u, X) of the Lagrangian over the space IVx M. In the present case, it turns out that one may choose: ,c(u, A) = I(u) + (T*u -f,
variational
(10)
Constrained variational problems In the following, the problem of Figure 1 is considered, that is: A$ = T*Tq5 = f
in a1
u*T@ = 0
on
K[@l=F
on
r, rz
(17)
(18) (19)
For the sake of simplicity, essential homogeneous conditions have been chosen on rr ; otherwise, functional (9) has to be considered, where essential boundary conditions are to be explicitly imposed. The basic idea of the variational formulation is as follows: find the stationarity point
Obviously boundary conditions on r are needed. The solution in an2 is unknown, but it is possible to write a boundary integral equation given a fundamental solution G of the differential problem. The usual theory of fundamental solutions of partial differential equations, leading to the boundary integral formulation of the problem, must now be stated in terms of distributions. A complete exposition of distribution theory can be found in Guelfand and Chilov;r4 for boundary equations, see Alliney and Tralli.” Let @Edenote the solution in S& and assume: @E=@k+&
(11)
where @kis an assigned function which satisfies the first equation of (10); $i is infinitesimal at infinity together with its derivatives up to the (2m - 1)th order if A is of order 2m;16 f IS ’ assumed here with a support bounded in al. Therefore, the integral equation for & will involve only the values of the function on the boundary (contributions at infinity do actually vanish, due to the asymptotic behaviours of both G and @ii).Thus:
I r
@i(S)ti(t) WC-9=
Is
[h u*TG - Gu*T@I G(f)
rr
Figure
7
Domain ~2, of boundary
Appl.
value problem
Math. Modelling, 1985, Vol. 9, August
247
Primal-dual variational problems:
F. Sgallari
of the Lagrangian functional ,C(u, h) associated with the differential operator A, subject to the constraint given by the boundary integial equation. It should be noted that the problem in s2, has no boundary conditions explicitly prescribed on rZ. It follows that L(u, X) is simply: ‘c(u, h) = - &,
u) + (T*u -f,
X)
w9r,l
+ :[(o*u> @jr, + (u,
(20)
In the following previous results by Brezzi,“.17’18 Babushka” and Nedelec12 are used. First, details are required of the functional spaces related to the problem. Given an elliptic operator A of order 2m, the solution ti, h: of the variational problem for L(u, h) is to be found in the spaces: X = H(T*; L&) = {u E (L2(fll)r;
T*u EL’(&)}
(21)
tation of the boundary (further details in Ciarlet13 and Johnson and Nedelec19). Now, a set of basis functions u&x), hk(x) is chosen such that: VN(s2,) = Span {uk(x);
k = 1,2, . . , N} C X
(28)
WL(s21)=Span{hk(x);
k=
(29)
As is well-known, tions U&Z), hk(x) appropriate order with local support. described as:
1,2,...,L}CM
in numerical applications are real-valued piecewise but, here, not necessarily Thus, any function u E
the basis funcpolynomials (of continuous) X, h EM is
N
u(x) = c Y&(X) k=l
and
(30) M = L2(3 1)
(22)
Note that functions belonging to X do not need to be continuous, provided that a generalized Green’s theorem holds. As has already been pointed out, the boundary integral equation is to be considered in the weak sense; so equation (14) is considered which can be rewritten as: (23) where:
h(x) *
4
bkhk(X)
k=l
where the meaning of the parameters Yk, bk is strictly related to the finite element implementation. Some attention must be devoted to the approximation of the ‘trace’ spaces on the boundary. Recall that it is necessary to represent the contribution of the boundary terms of K1 and K2. One considers: T(r2) = tr(VN(al))
K,[-]
= [.I -
[.I
o*TG dy(x)
= Span{Uk(x)lrz;
k= 1,2,. . . ,J}
(31)
(24) S(r,)
I.2
= Span {@k(x); k = 1, 2,
. , I} C L2(r2)
(32)
Thus, one has:
and
J
K,[.] =
Ga*[.] dy(x)
(25)
s
u(x) N c YkUk(X) k=l
x E r2
rz
It should be noted that in the present ‘weak’ formulation q51r, does not correspond to hlr, (the ‘trace’ theorem fails); its meaning is simply that of an additional Lagrange multiplier. On the contrary, in the ‘strong’ formulation matching with the domain variable makes sense. As far as u is concerned, in equation (25) one can take the restriction u lr2. It is easy to see that both 4 and $ require only L2 summability on r2. Now, the problem (17)-(19) may be stated in abstract variational terms: find an element ti E X and an element fi EM such that:
L(ri,h)
VuEX
VhEM
(26)
ti, + W2b41, $, = UC $) V$ E L2(S12)
(27)
Under previous assumptions on the operator A and with the appropriate spaces X and M, the necessary and sufficient conditions (see Brezzill and Babufia” for (26) and Nedelec12 for (27)) for existence and uniqueness of the solution (ti, i) are satisfied.
Discrete approximation of the variational problem To find an approximate numerical solution, the discrete analoaue of variational problem (26) is considered. First, suppose a discretization over the-problem domain a1 has been established. For simplicity, it is assumed that the discretization does not introduce any error into the represen-
248
Appl.
Math.
Modelling,
1985,
Vol.
9, August
@(x) = i k=l
(33)
uk@k(X)
According to the previous remarks, functions G(x) themselves can be identified with @k(x). Using equations (30) and (33) one can rewrite functional (20) and the constraint (27) in compact form: ~~,b,a)“--fyT~Y+YTDb--BTf + &@%
+2-y]
(_T-M)a+_Ny=F
(34) (35)
Refer to the Appendix for detailed definitions of matrices and vectors. The constraint equation (35) can be solved2’ with respect to (I, to give: (36)
a=g+_4, where the following positions are introduced: g=(_T-MM)-‘F
(37)
p=-(T--_M)_‘_N
After substitution of (36) into the expression of the objective functional, one obtains: L(y, b) = - $y’cy + HYTST4
+ yTpb - bTf +yTsTg
+YTETsY +sTsyl or equivalently:
(38)
Primal-dual
+ &lyTzTg+gT$y]
-bTf
(39)
We have now an unconstrained saddle-point problem for L’b, b). Its solution is determined by the linear system: -Ey+a,+Db=-&S’g
(40)
DTb =f
problems:
F. Sgallari
Conditions on r2 are obtained by a boundary integral equation with fundamental solution C(x,, yo; x, u) = log(2R/r) + log(2R/r’), r = [(x -x0)’ + (y -yo)2]1’2, r’ = [(x - xo)2 + 0, + yo)2]1’2 and R is the radius of the circle r2. For a detailed exposition refer to Alliney.” Using the previous notation, one has:
qy,b)“-$yT~y+$yT[&~Tf+PTS)]y -+y’gb
variational
T* = - div
T = grad thus:
where the N x N matrix:
u=n
u = V@
(41)
u* = n
and :
resulting from the introduction of the boundary integral equation (considered as a constraint for the functional) turns out to be symmetric.
X=H(div;~;2)={uE(L2(~;2))2;divuEL2(~1)}
Example The procedure outlined in the previous sections, is well suited to the study of steady flows of ideal fluids past completely reflecting obstacles. It is assumed that these flows are irrotational and that there are no sources and sinks. Under the hypothesis stated, the determination of a flow in a region a can be reduced to determining a realvalued potential function $(x,y), such that u = V@. With reference to Figure 2, in the unbounded plane domain s2 the function Q satisfies Laplace’s equation: -V2$
inSZ=a1Ufi2
= 0
Functional
(20) can be rewritten
as:
~(~,A)=-illlu12dS2,-Il(divu)hdn, n,
a, -t$
(42)
where the constraint
5 rz
[u.n6+@uun]dy
(47)
is:
At the boundary rr (assuming that there is no absorption of the fluid through the walls) the flow is in the tangential direction; Thus, the boundary condition is:
a$
-0
on l?,
(43)
an-
In addition, when one determines the flow past an obstacle, one imposes the condition that the flow is homogeneous at infinity. This means that the potential @ has the following asymptotic behaviour: ~(x,Y)=~k(x,Y)+~i(x,Y)
+
rz rz
=
(44)
r-++m
(45)
IV&I
r++m
(46)
= W/r2)
with r = (x2 + _Y~)“~.
-Y
x
Figure 2
Definition
sketch for test problem
s
K9 9(t) WC;)
(48)
r2
Here, Gk(x, v) denotes an assigned function, which satisfies equation (42); physically, it represents the potential of the incident flow. The function &(x,v) is the perturbation potential due to the obstacle, and it is required that: $J&,Y) = 0(1/r)
W-G8 ; 6) d4.9d-r(x)d%9
The major difficulty, when using a dual, or primal-dual formulation, is to take appropriately into account the constraint T*u - f = 0, in this case div u -t-f = 0. There are essentially three ways of circumventing this difficulty. First, the constraint is approximated, in such a way that the discrete solution ii satisfies a relation of the form diva + f= 0 in !ZZ,where fis a typical finite element approximation of the function 5 By reference to elasticity problems, such methods, which are directly based on the dual formulation, are known as equilibrium methods.13 While the equilibrium methods are, by definition, based on a formulation where there is only one unknown (the gradient V@= u for the problem here, the stress tensor u for the elasticity problem, etc.), one may use the techniques of duality theory to get rid ot the constraint, a process which results in the addition of a second unknown, the Lagrangian multiplier. For a detailed exposition of such methods refer to Ciarlet13 and Mercier.21 For the discrete approximation of X, the mixed finite element method implementation of Johnson and Mercier” was used. Consider a triangulationyj over the set 52, i.e. the set Sz is subdivided into a finite number of triangles K, the index h representing the maximum diameter of the element K. For 2 and x the expansions:
Appl.
Math. Modelling,
1985,
Vol.
9, August
249
Primal-dual
variational
problems:
F. Sgallari
The numerical model presented above has been implemented on the VAX 1 l/780 computer of the Computing Centre at the Faculty of Engineering of Bologna University. To assess the numerical performance of the method, first consider a half-circular obstacle centred at the coordinate origin (Figure 4 gives the finite element mesh). For this particular case, the analytical solution is well known:
Figure 3 Triangular element for mixed implementation discussed in text. Degrees-of-freedom for dual variable are normal components of u at points on edges
6=
il Y&(X> Y>
x= i
(49)
bjXj(X,Y)
=-vx( 1+-&)
(56)
where --z, is the unperturbed incoming velocity and a is the radius of the obstacle. In Table 1 the results are presented for a = 1 and z1= 50. It should also be noted that with a coarse mesh there is excellent agreement with the exact solution. In the second example the obstacle in Figure 5 is depicted, with z, = 4 m/s. Figure 6 holds for any element a straight line segment of arbitrary length, which is tangential to the streamline in the centroid of the triangle.
j=1
were used, where the basis functions ui have local support on two adjacent triangles and $ have support on one triangle. Moreover, ui have a continuous normal component on the common edge of the two adjacent triangles and vanishing normal components on the other edges; thus ui are linear vector-valued functions on each triangle where they satisfy the ‘equilibrium’ equation divui .= g, for g piecewise constant. As degrees-of-freedom the normal components on six points were chosen (two for each triangle edge); at a single point ui. n = f 1 is imposed, with ui. n = 0 on the others. For hi consider constant functions on each triangle (corresponding to one degree-of-freedom). For the practical construction of the basis functions refer to Thomassetz3 and Johnson and Mercier22 where the implementation of the dual finite element is discussed for the more general case of elasticity problems. As far as error estimates are concerned one has: l/11- ii llL2 < c-h2
(51)
IIh - KllL2 < ch
(52)
40
30
20
IO
0 40 Figure 4
20
30
IO
00
20
IO
Y
30
Finite element
mesh for test example
Finite element
mesh for second obstacle problem
40
where c is a constant; for more details refer again to Ciarlet l3 Mercier,‘l Thomasset and Johnson and Mercier.22 The diicrete functional (47) and the constraint (48) become: 0 I(
.X(y,b,a)=-_yTEy+yTgb
Figure 5
+ + LvTsra + a’$y] (T-ma++
=F
(53) (54)
The constraint equation (54) can be solved with respect to a and, after substitution into (53) the objective functional LQ, b) is obtained as in (38). The solution is now determined by the linear system: --_Ey+_@+Db=-$sTg
-
gTb=Q
(55)
It is worth noting that the coefficient matrix is non-definite; however, one can build it up element-by-element using the standard assembling techniques and using a new iterative process9 that requires only minor modifications on Irons’ frontal solution program.8
250
Appl.
Math.
Modelling,
1985, Vol. 9, August
80 60
40 20 0 100
0.0
60
40
2.0
00
20
Figure 6 Straight length segments of arbitrary streamlines in centroid of each triangle
40
60
80
length tangent
10.0 to
Primal-dual
variational
problems:
F. Sgallari
Table 1
Element
x, y coordinates
Computed
1
1.491
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
1.236 1.191 1.368 1.168 0.877 0.703 0.688 0.313 0.114 2.237 1.978 1.900 2.047 1.762 1.415 1.110 1.020 0.484 0.199 2.983 2.720 2.609 2.726 2.356 1.953 1.517 1.351 0.655 0.285 3.483 3.713 3.534 3.159 2.789 2.708 2.009 1.522 0.826 0.456
-107.477 - 101.074 -99.155 -99.155 -85.122 - 72.702 -57.617 -49.846 -23.113 -9.877 -132.701 -122.715 -118.346 - 122.436 - 105.989 -88.816 -69.490 -61.211 -29.277 - 12.654 -162.120 -152.792 -146.814 - 151.875 -131.961 -111.360 -86.091 - 75.434 -36.468 -16.388 -185.199 -188.334 - 192.443 - 174.888 - 150.744 -149.319 - 106.059 -85.814 -41.946 - 26.235
39
40
0.101 0.159 0.325 0.559 0.905 0.868 1.015 1.308 1.444 1.229 0.144 0.246 0.537 0.853 1.346 1.376 1.633 1.969 2.164 1.964 0.188 0.332 0.749 1.146 1.788 1.884 2.251 2.630 2.884 2.699 0.188 0.419 1.086 1.396 2.038 2.517 3.104 3.100 3.354 3.669
Exact
References 1
2
3 4 5 6 I 8 9 10 11 12 13 14 15 16 17
Atluri, S. N., Gallagher, R. H. and Zienkiewicz, 0. C. (eds), ‘Hybrid and mixed finite element methods’, John Wiley, New York, 1983 Brebbia, C. A., Telles, J. and Wrobel, L. ‘Boundary element methods - theory and application in engineering’, SpringerVerlag, Berlin, 1983 Mei, C. C. and Cheu, H. S. Int. J. Num. Meth. Eng. 1976, 10, 1153 Yue, D. K. P., Cheu, H. S. and Mei, C. C. Znt. J. Num. Meth. Eng. 1978,12,245 Zienkiewicz efal. Int. J. Num. Meth. Eng. 1977, 11, 355 Allinev. S. Aunl. Math. Model& 1982, 6. 424 Arthur;, A. M: ‘Complementary &riational principles’, 2nd edn, Clarendon Press, Oxford, 1980 Irons, B. M. Int. J. Num. Meth. Eng. 1970, 2, 5 Alliney, S. ‘Iterative solution of large linear systems arising from saddle-point problems’ (submitted) Babuxka, I. Math. Comp. 1973,27, 221 Brezzi, F. R.A.I.R.O. Sk, Rouge 1974,3, 129 Nedelec, J. C. ComputerMeth. Appl. Mech. Eng. 1976,8,61 Ciarlet, P. G. ‘The finite element method for elliptic problems’, North-Holland, Amsterdam, 1978 Guelfand, L. M. and Chilov, C. E. ‘Les distributions’, Vol. 1, Dunod, Paris, 1972 Alliney, S. and Tralli, A. Appl. Math. Modelling 1984, 8, 75 Lions, J. L. and Magenes, E. ‘Problemes aux limites non homogenes et applications’, Vol. 1, Dunod, Paris, 1968 Brezzi, F., Johnson, C. and Mercier, B. Math. Comp. 1977, 31, 809
- 107.934 - 101.594 -98.617 -99.721 -85.180 - 72.666 -58.208 -50.174 -22.848 -9.436 -134.121 -123.811 -119.377 -123.185 -106.037 -88.939 -69.737 -61.376 -29.150 -12.533 - 165.869 -154.149 -148.163 -151.920 -131.276 - 110.939 -86.143 - 75.320 -36.523 -16.184 - 188.488 - 198.959 -189.657 -171.226 -151.148 -145.307 -107.819 -82.524 -44.790 -24.468
2.560 8.034 17.405 16.566 21.753 33.653 30.019 19.090 9.758 8.039 0.736 2.556 6.151 6.032 8.838 12.281 11.462 7.953 4.404 2.829 0.984 1.633 2.791 3.690 4.853 6.248 6.866 4.098 3.039 0.989 1.403 3.213 10.123 2.014 7.383 5.141 5.838 4.171 8.624 2.537
-26.900 -16.681 -20.611 -33.730 - 43.890 -47.120 -61.910 -63.971 - 71.357 -81.498 -38.505 -36.647 -38.594 -42.734 -47.485 -49.500 -55.089 -56.156 -59.054 -62.426 -40.928 -41.430 -43.272 -47.427 -49.403 -50.173 -53.154 -53.683 -54.539 -57.186 -36.388 -36.125 -52.087 -56.793 -47.588 -57.657 -43.402 -64.174 -59.124 -48.718
18
19 20 21 22 23
-27.821 -18.858 -21.788 -33.678 -44.276 -49.661 -61.526 -62.954 - 70.820 -82.214 -40.135 -37.804 -39.080 -42.843 -47.329 -49.642 -54.718 -55.863 -59.192 -62.570 -44.149 -43.141 -44.249 -46.002 -48.462 -49.756 -52.547 -53.327 -55.153 -56.634 -45.915 -46.509 -46.974 -47.179 -48.727 -49.734 -51.497 -52.561 -53.709 -53.545
3.027 8.155 16.678 16.566 22.150 32.771 30.679 18.857 9.482 7.024 1.281 3.080 6.716 7.215 9.805 12.814 11.919 8.300 4.332 2.576 0.752 1.604 3.600 4.083 5.503 6.780 6.288 4.646 2.469 1.416 0.442 0.799 2.053 3.098 3.991 3.646 3.336 3.315 1.946 0.895
Brezzi, F. and Raviart, P. A. ‘Mixed finite element methods for 4th order elliptic equations’, in ‘Topics in numerical analysis’ (ed. Miller, J. J. H.), Academic Press, London, 1976 Johnson, C. and Nedelec, J. C. Math. Comp. 1980,3.5, 1063 Albney, S. Appl. Math. Modelling 1982,6, 291 Mercier, B. ‘Topics in finite element solution of elliptic problems’, Springer-Verlag, Berlin, 1979 Johnson, C. and Mercier, B. Numer. Math. 1978,30, 103 Thomasset, F. ‘Implementation of finite element methods for Navier-Stokes equations’, Springer-Verlag, Berlin, 1981
Appendix The following are the definitions of matrices and vectors of the section covering the discrete approximation of the variational problem. One has:
E = [EiiI NxN Eii
= Eji
=
s
ui(x>q(x) d%(x)
a, D = Piil Dii =
Appl.
NxL
Xi(X) T*Ui(X)
Math.
(AlI
.
dCLl(X)
Modelling,
1985,
(43
Vol. 9, August
251
Primal-dual variational problems: F. Sgallari
As far as the boundary is concerned, the matrices --S, T, M, & with dimensions Ix J, Ix I, Ix I, I x J, respectively, are used, with elements: Sii =
=
s
$i (X)
s
@i
CJ*Ui (X)
dy(X)
Moreover, the following vectors were introduced:
(X)uj (XI dy(x)
6431
rz
Tii = Tji =
fi s r*
@i(X) $j (X> d+&)
Mii zzMji = ss rz rz
(A41
Appl.
Math.
I
f(x)
AjW d%(x)
(A7)
a2, F= [F1, . . ..FI]’
@i(x) o*TG(x,
[I @j(E) d?(x) dy(O
G48) WI
252
=
Modelling,
1985,
Vol.
9, August