Convergence of difference schemes of the reference-operator method to generalized solutions of the axisymmetric Poisson equation

u.S.S.R.comput.Matks.Matk.Pkys.,V01.30,N0.5,pp.140-147,1990 Printed in Great Britain

0041~5553/w $10.00+0.00 &X991 Pergamon Press plc

CONVERGENCE OF DIFFERENCE SCHEMES OF THE REFERENCE-OPERATOR METHOD TO GENERALIZED SOLUTIONS OF THE AXISYMMETRIC POISSON EQUATION'k A.A. DENISOV, A.V. KOLDOBA and YU.A. POVESHCHENKO The convergence of difference schemes of the reference-operatormethod for the axisymmetric Poisson equation to generalized solutions and the accuracy of such schemes are investigated. Sufficient conditions for convergence are formulated in terms of the energy norm. This paper is devoted to an investigationof difference schemes for the axisymmetric Poisson equation, based on the reference-operatormethod 11-31 with irregular meshes, in the context of generalized solutions. The convergence-ofdifference schemes for elliptic equations, including the axisymmetric Poisson equation, has been studied in numerous publications; a summary of the results may be found in [41. Convergence has been studied in weighted Sobolev spaces, which appear naturally when Cartesian coordinates in three-dimensionalspace are replaced by cylindrical coordinates r, z. However, only rectangular meshes have been considered. In [5] conditions were derived which guarantee convergence to first order of difference schemes of the reference-operatormethod for the Poisson equation in plane geometry, for minimal restrictions on the mesh. It was shown in [61 that these conditions are sufficient for the convergence to generalized solutions inbvZ*< Our object is to extend the results obtained in [6] to the case of the axisymmetric Poisson equation. It will be shown that the sufficient conditions for convergence obtained in [S, 63 remain valid. We shall be concerned with the Neumann,and Dirichlet problems in a bounded three-dimensionalregion 0 with smooth boundary. Unlike [4], we shall consider Sobolev spaces wk,p with the norms defined in the standard manner. As in [Sl, the method of proof is based on minimizing the error with respect to all mesh functions of a certain type in a suitable natural norm. We shall use theorems concerning embedding and equivalence of quotient norms on W*+'*~/9)* (4 is the wk+1.* (see [7, Chapter 3)). space of polynomials of degree at most k) and seminorms on 1. The family of difference schemes Let us consider a three-dimensionalregion 0 with boundary 80 and in it the Poisson equation, written as a system divs=f(z),

q=-gradu

(1.1)

with certain boundary conditions; here fix) is a given function, and s~(z,,zp,x~) is a point in three-dimensionalspace. We shall assume that the problem possesses axial s~metry, i.e. 0 is obtained by rotating a two-dimensionalregionb about an axis of symmetry, the functions u, f are independent of the azimuthal coordinate& depending only on the variables (r, 2). In addition, we shall assume that the input data (the region, right-hand side f(r) and boundary conditions) are such that problem (1.1) has a solution u(z)~W*.'(0), Let us introduce a family of irregular meshes in the regions in the (r, z)-plane. Each mesh consists of cells(Q) and edges(o).We single out from the set (a)(see Fig. 1) those elements that lie on the boundary of the computational region and denote them by(&). In addition, we introduce a set of fictitious boundary cells(B.),whichwill be used in imposing the boundary conditions. The mesh in the three-dimensionalregion 0 is then obtained by rotating the twodimensional mesh through unit angle about the axis of symmetry. The bodies and surfaces thus swept out by the plane cells and edges will be called cells and faces (Fig. 2). We shall consider the mesh to be connected, in the sense that one can enter any cell from any other by passing through faces. The fineness with which the mesh subdivides the region will be characterized by a parameter k; one possibility is h=(mesa/.'V)'", where N is the number of cells in the mesh. We now consider the linear span in [W1*'(01)13 of the following elements of this space: 1) the constant functions q = const; 2) the azimuthally symmetric functions, i.e. vectors q(r,z,h) obtained from vectors q(r, 2, 0) by rotation through an angle h about the axis symmetry. An operator projecting functions of the manifold thus obtained onto the mesh may be

141

FIG. 1

FIG. 3

FIG. 2

defined as follows. Let o(h) be the oriented surface swept out in f%%Xs)

-space by rota-

ting the edge a through the angle h (see Fig. 3, where this surface is indicated by hatching), and define j;(k)=Jqds. NM The surface may be oriented in an arbitrary manner; to fix our ideas, let us assume that the boundary faces have the same orientation as the boundary ofo.The integral is into L’ over an arbitrary smooth surface meaningful by virtue of the embedding of WiJ(0) in 0‘ (we are thus assuming that the boundary of 0‘ isC').The projection of a vector-valued function q(X) onto the mesh is now defined as the mesh function q defined on (a)by the formula q=(d^q/dX)+o. The projection of a function q = const will obviously have the form q"= (Q, s").

The vectors sa will be called the areas of the faces. For axisymmetric func-

tions,

(I”= j qds. O(1) As already noted, the axisymmetric Poisson equation will be considered as a system (1.1) in the three-dimensionalregion C? with a function u(X) independent of the azimuthal variable.From this point of view the axis of symmetry is not Part of theboundary, and therefore edges lying on it are not boundary faces and do not require adjacent fictitious cells. Since the fluxes through such edges vanish, all terms containing these fluxes as factors will also vanish. We may therefore assume that there are no edges on the axis, i.e. we shall exclude them from the start from the set(u).The space of mesh functions defined on (cf will be denoted by 3. Under this discretization of vector fields, it is natural to approximate the divergence operator and the first of equations (1.1) as follows:

P

(I

where the summation ranges over all faces of the cellP and o=ki according as the relevant side of the face is exterior or interior relative to% Since we are interested only in the error of the difference scheme due to the irregularity of the mesh, we shall assume that the integral in (1.2) is evaluated accurately. According to the method of reference operators, the difference analogue of the operator grad is constructed in such a way as to satisfy the difference analogue of the integral identity

where 0, is a region measuring one unit in the direction of the azimuthal coordinate, and 00, is the corresponding part of the boundary%Y. It is then natural to define the mesh analogue of the function u(X) on (a)U(Q.) and approximate the integral

s

(P, rl) dV

(1.3)

by a bilinear form (1.4)

142

Here g,=g..,(8)are the elements of a certain symmetric matrix, whose actual form is not important for the sequel; P, Q are difference analogues of the operators f tie, 0

Sqda. (i

The integral on the left is approximated by a sum

we have the identity (1.5) The first term on the right of (1.5) approximates

while the second, therefore, approximates

- .f (P, gradu)dV - .f (p,n)N. 6‘

6%

On the other hand, the integral (1.3) is approximated by the expression (1.41, whence, since P is arbitrary, we have z g~(a)Q"+g~(a')Q"=a(n)u(a)+o(B')U(O'). ‘

(1.6)

The operator on the left of (l.S), applied to the mesh function Q, will henceforth be denoted by G and called the metric operator. Since U- il, the right-hand side of (1.6) is the expression -AU=U(a)-U(n'), where B ’ is the cell for which the selected side of the face u is interior, whereas for B it is exterior. Thus, the operator method yields the relation GQ = -AU. To each face a we assign an integer m(u) such that the part of the face in a sector with angle 211/mcan be enclosed in a sphere of radius O(h). Define the diagonal operator m:?if-+i% with elementsm(U). The matrix gab is assumed to satisfy the following estimate uniformly with respect to meshes, cells and faces: (1.7)

zb IgtiIGTl(hm), while the operators G satisfy the following estimate, again uniformly: G>Bm-l h

(1.8)

’

We also impose restrictions on the geometrical dimensions of the mesh elements: b,hP4s(a)/m9

b,h’;

(1.9)

the number m in the first of these inequalities is that corresponding to any face of the cell.a. The metric operator G yields a natural definition of the scalar product in%:

and of a norm [X{p=. We will now decomposeI into the direct sum of subspaces& anda.The elements of& are the functions G-'AT, where T is an arbitrary function on (sl)U(sZ.), while a is the set of all functions orthogonal to those of &in the sense of the above scalar product. Let us investigatethe structure of the elements ofa. LetXE.!$,YES. Then

143

Hence, since T is arbitrary, it follows that the elements of g are functions such that

Isu.Y"=O on

(8),

Y"=O on

(8s).

2. Estimate of the rate of convergence of difference scheme of the reference-operator method for the Neumann problem Let us consider the second boundary-value problem for equations (1.1). We shall assume that the problem has a solution ~(2. Let u~W'*'(a); then ~~@'h*(o)) component by component. Since we are considering wLV2(o)) as embedded in L' over an arbitrary smooth manifold in0,the function q,(z) = (q, n), regarded as a function on the face o (n is the normal to the face at a~),is an element of L‘(a), i.e. the mesh function introduced in Section 1 is defined. Let us compare the mesh function q = q"thus obtained with the mesh function Q which is the solution of the difference problem E o.Q"=F on D

(a), CQ=--AU on (a), Q"=Q. on (aa),

(2.1)

where G is the metric operator, F=

JfdV,

Q.= -

$ds.

P

The system of equations (2.1) has a solution provided that

(0)

c*w

The following theorem was proved in [6]: Theorem 1. Under the above conditions, we have

IkQII= id Ik~ll.

(2.‘)

The conditions for convergence of difference schemes of the reference-operatormethod for the axisymmetric Poisson equation, as applied to generalized solutions, are analogous to those obtained in [61 in two dimensions. In order to formulate these conditions, we consider the mesh function s defined above on (0). We will require this function to satisfy the following conditions. Condition 1. s is an element of the subspace&,i.e. Gs,= Axi (component by component). Obviously, the vector-valued function z:iis defined on (Q)U(a.)apart from a constant vector. In each cell we specify some point, say the centre of gravity (in the fictitious cells this will be the midpoint of the corresponding face), whose coordinates are denoted by Z‘(n). Condition 2. There exists xi(n) which satisfies Condition 1 and whose distance from f‘(9). in terms of the norm

ll41.== SUP lz(Q)I, ,Q,“C”., is of the same order of magnitude as h, i.e. inIl]x-Zlj,=0(h). We shall show that under these conditions, if there is a solution in IV*'(a),the difference scheme converges to it; we shall also determine the rate of convergence. Theorem 2. If Conditions 1 and 2 are satisfied and a solution exists in W'~'(fY), then the difference scheme (2.1) converges to a solution at a rate O(k), in the sense that the following estimate is valid:

II~-QII~C4l4. Proof.

2.B.

(2.3)

According to [8, pp.129-1901, a function ~(z)EIVJ(~) may be continued to a

144

FIG. 4

FIG. 5

function of kV~(IR').We denote this function again by u(x), noting the inequality Ilull~,t.k< BIMln, 2, 8. Since w*'([R')can be embedded in C(w'),we can form the mesh function u(~)=u(x 69). Put p=-6YAu in (2.2). Then II~-QI~z~~~~+G-lA.~~z~~~m~GpfAu~'. D Let us estimate the quantity IGp-+ Au\ for a single face. To that end, consider a sector with angle 2nfm; let S be a sphere of radius p, whose interior contains the parts cut out by the sector from the face in question and from the cells adjacent to it, as well as the centres z(Q) of these cells (see Fig. 4). By Condition 2, we can choose p
By the Hahn-Banach Theorem, any such funcwhere M,=O(f) and it is the same for all p. tional p can be continued to the entire space IV* (S,) without increasing its norm. Using the embedding of W8*2 (S,f into C(S,1, we obtain the following chain of inequalities:

Condition 1 means that the functionalcp vanishes on polynomials of degree 1. Therefore, by the equivalence of the quotient norm of Wc2/9, and the seminorm of W’*a we obtain ~~lczMJ(.~,cy+~*)Ivl*,2.S‘ Returning now to r-coordinates, we obtain l~pl~2Ms(~,cl+Ma)~'"luIz,2,.s, where S is the preimage of S, in (s,, x2,.za)-space. The spheres S thus constructed cover a certain neighbourhood of the half-plane h=O. In order to cover the entire regionO,we express ImI'as

where qa is the functional constructed in the same way as rp,exceptthat q=(Jq/dh)r_,..,,,,. Obviously, in problems with axial symmetry rp,-_cp. We thus have

145

here the summation ranges over spheres covering the whole of 0. Since there exists a number n such that any point in (.z,,x~,z,)-space is covered by at most n spheres S) it follows that

so that

finally,

where C=BBM,(1M,C+M,)(cn/S)'~. II~-QII~C~~IUII~,~. O.

3. Estimate of the rate of convergence for the Dirichlet problem In this section we shall examine the convergence of difference schemes of the reference-operator method for the Dirichlet problem, namely:

I2@“=F

on (9),GQ--AU

(3.1)

on (a),U=U. on (8.),

where G is the metric operator. It is easy to show that under our assumptions system (3.1) has a unique solution. We shall assume that u, = u(x+), where x, is some point of the corresponding boundary face. As in Section 2, we express8 as a direct sum &@g, but the subspaces .@,9 are defined differently:.@ is the set of all functions x=G-'AT, where T(8) is an arbitrary function such that,T(Q.)=O, and 9! is the set of all functions Y such that '&o.p =0 As in Section 2, we obtain

'&a.(q-Q)"=O, i.e. q -Q=B,

on(Q).

and so (see 161):

(3.2) IkQII= inf IIP-PII, P-1=& -where p=-G-'AU, U being an arbitrary function such that O(Q.)=U.. For the sequel we shall need the following mesh norm. Let y be a function defined on (Q)U(Q.) with y(Q.)=O. We define a mesh analogue of the norm in ah'(C) by the formula (IY~I,*=IIG-'AY~~~. AS in Section 2, the function u(x)E~~'(O‘)can be continued to the whole of (sl,%Z~) -space and we define a mesh function u(Q)=u(z(Q)) on (Q)U(Q.). Conditions 1 and 2 are modified as follows. Condition I’. There exists a function x(Q) such that Gs,=Asr (component by component). Condition 2 ‘. The function s(a) satisfies the conditions x

(8.)

=x.,

supI&2)-fl=Q(h). (0,

Theorem 3. If Conditions 1' and 2' are satisfied and the problem has a solution in then the difference scheme (3.1) converges to such a solution at a rate O(h); specifically, one has the estimates W’s * (0)

lI~-Qli~W4l~, 2.or

Il~-~ll~~w41*.

2.0.

(3.3)

Proof. The first of inequalities (3.3) is proved in exactly the same way as in Theorem 2. We will therefore only prove the second. Put <=-G-'AU. The error y=U(Q)-u(Q) satisfies the system of equations -~~.(G-'AY)"=CU~(Q-~)~=~~~(~-~)~. Li D S

(3.4)

In addition, it is clear that y(8.)=0. Multiplying (3.4) by y and summing over(Q),we obtain (IG-'Ayll'=
146

4. The difference scheme. Sample investigation To illustrate the approach presented above, we will construct a metric operator that satisfies the convergence conditions 1 and 2 and demonstrate the validity of estimates (1.7) and (1.8) for it. Consider a mesh all of whose cells are triangles in the (r, a) plane. The vertices of each cell are numbered by an indexIp.To each vertex we assign a volume V, and consider the metric operator based on a matrix gkz of the form

g,, -

v&J,‘, et%. zl P

(4.1)

The summation index 9,runs over all vertices of a cell Q; the indices k, 7, in the expression (S*', S1Yr take values U,Z corresponding to the indices of the faces that meet in the vertex cp.The vectors {s.',s,'} form a basis dual to{s',s'}.Obviously, the vector S,'is determined not only by the face a itself but also by the vertex which it forms (there are four such vertices for each interior face and two for each boundary face), i.e. by the pair {s",sT) forming the basis associated with the vertexcp.Thus (s;,81').is the Gram matrix constructed from the vectors {S.',s,'} and inverse to the matrix (s',s').. The index cp indicates that k, 1 run over the values a, 7. Each term of the form Vq(sl’r s,‘)@ "fills" the rows and columns of the matrix of the operator G indexed a,~. With this choice of gk2, we have

I:V VA‘.

Gs-

(4.2)

The summation here extends over all vertices formed by the face a,%' is the vector paired with 8" at the vertex cp. Consider a face forming the vertices cp=A,B,A’,B’ (Fig. 5). It is easy to see that s*’ =

r.d--rc SAW2(r*+r.)

’

s*’ =

where sABC is the area of the triangle ABC in the (P, 2) plane, PA, rB are the distances from the indicated vertices to the s axis, and r,, rB, rc are the radius-vectorsof the vertices of the triangle ABC. Set V~=~.&MK& V~=~E&S& in (4.1). Then

--.

3

Obviously, the end of the vector -- 2rr4-r,

’ =s3(r,fr,) ‘* +

r4+2rB

3(rA+rs)

‘a

lies on the side AB, R=(r_,+r&r,)/3 is the radius-vector of the centre of mass of the triangle ABC, and the vector p-R connects these two points. Carrying out similar constructions for the triangle A’B’C’ and noting that p depends only on the position of the vertices A,

B,

we see

that the vector (4.2) connects the centres of mass of the

triangles

ABC and

A'B'C' . Thus, for z+(Q) in Condition 1 we must take the coordinates of the centres of mass of the triangles. For fictitious cells, zi(Q) will be the coordinates of the endpoint of the vector p. Condition 2 is obviously satisfied. We will now verify the validity of estimates (1.7) and (1.8) for the operator G based on matrices gkt of the form (4.1). Consider an arbitrary cell&J.Its vertices will be denoted again by A,B,C, and the sides opposite these vertices by a,b,c, respectively. With this notation, the elements of the row of gkz corresponding to the side a will be (the other elements are evaluated in similar fashion)

147

Recall that the number m (see Section 2), by definition, is chosen in such a way that all the linear dimensions of the part of the three-dimensionalcell enclosed in a sector with angle 2n/m are of the order of h (see (1.9)). In view of this, we have

If we now require the angles of the cells to satisfy the condition m>B>G, the expression in braces will be bounded by 4sinAa@. Condition (1.7) will then hold with ~=8xn,(b,*sin*B)-*. It may also be shown (estimating the minimimum eigenvalue of the matrix gkt) that condition (1.8) hold; dith ~=4a,x/(Sba*). To summarize: the difference scheme of the reference-operatormethod based on a triangular mesh with metric operator (4.1) converges to first order in the mesh norm ofWL2. REFERENCES 1.

2. 3. 4. 5.

6.

7. 8.

SAMARSKII A.A., TISHKIN V.F., FAVORSKII A.P. and SHASHKOV M.YU., Difference analogues of the fundamental differential operators of first order. Preprint No. 8, Inst. Prikl. Mat., Akad. Nauk SSSR, Moscow, 1981. SAMARSKII A.A., TISHKIN V.F., FAVORSK~I A.P. and SHAS~OV M.YU., Operator schemes. Preprint No. 9, Inst. Prikl. Mat., Akad. Nauk SSSR, Moscow, 1981. SAMARSKII A.A., TISHKIN V.F., FAVORSKII A.P. and SHASHKOV M.YU., Operator difference schemes, Differents. Uravn., 17, 1317-1321, 1981. SAMARSKII A.A., LAZAROV R.D. and MAKAROV V.L., Difference Schemes for Differential Equations with Generalized Solutions., Nauka, Moscow, 1987. KOLDOBA A.V., KUZNETSOV O.A., POVESHCHENKO YU.A. and POPOV YU.P., On the convergence of difference schemes of the reference-operatormethod for Poisson's equation. Preprint No. 85, Inst. Prikl. Mat., Akad. Nauk SSSR, Moscow, 1987. DENISOV A.A., KOLDOBA A-V. and POVESHCHENKO YU.A., On the convergence of difference schemes of the reference-operatormethod for Poisson's equation in the case of generalized solutions., Zh. Vychisl. Mat. mat. Fiz., 29, 3, 371-381, 1989. CIARLET P., The Finite Element Method for Elliptic Problems., North-Holland, Amsterdam, 1977. FADDEYEV D.K., VULIKH B.Z., URAL'TSEV N.N. et al., Selected Chapters of Analysis and Higher Algebra. Izd. Leningrad. Gos. Univ., Leningrad, 1981.

Translated by D.L.