2004-Fez conference on Differential Equations and Mechanics Electronic Journal of Differential Equations, Conference 11, 2004, pp. 143–155. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu/ or http://ejde.math.unt.edu/ ftp ejde.math.txstate.edu (login: ftp)

ANALYTICAL NODAL METHODS FOR DIFFUSION EQUATIONS

NAJIB GUESSOUS, FOUZIA HADFAT

Abstract. We analyze the direct analytical nodal methods of index l (ANM) and show that the corresponding mathematical methods are equivalent to the physical methods when the components of the matrices are calculated by generalized Radau reduced integration. Some numerical examples applied to the diffusion equation, show that in case l = 0 the polynomial nodal methods (PNM) [7] obtained better results that those of ANM.

1. Introduction

Broadly speaking, nodal methods are fast and accurate methods which combine attractive features of the finite element method (f.m.e.) [3] as well as of the finite difference method (f.d.m.) [8]. With the f.e.m., they share the fact that the unknown function is approximated over a given coarse mesh by a piecewise continuous function, usually a polynomial one. As with the f.d.m., the resulting systems of algebraic equations are quite sparse and well structured, essentially because the basic parameters are cell and edge moments of the unknown function, which connect at most two adjacent cells. These methods were developed in the late 1970s for numerically solving partial differential equations (PDEs) in general, especially in relation with static and dynamic nuclear reactor calculation. Classical references are Finnemann and al. [4], Langenbuch et al. [11, 12], Shober et al [14], Frohlich [5], Dorning [2] and Wagner and Kobke [15].

In §2, we present the transverse integrated version of the analytical nodal methods of index l (ANM) [6]. In §3 we analyze the direct procedure of ANM introduced in [6] and we show that the corresponding mathematical ANM are equivalent to the physical ones when components of the matrices are calculated by generalized Radau reduced integration. To solve the arising system of discretisation which are sparse and symmetric, we use in §4, a Gauss-seidel iterative method. In§5, Some numerical examples applied to the diffusion equation one permitted to notice in case l = 0, that the polynomial nodal methods (PNM) [7] succeeds to the best results that those of the ANM.

2000 Mathematics Subject Classification. 74S30.

Key words and phrases. Partial differential equations of elliptic type; numerical analysis;

diffusion equations; finite element methods.

\begin{matrix}c \\ ◯ \end{matrix}2004 Texas State University - San Marcos.

Published October 15, 2004.

2. Basic analytical nodal formalism

Let us consider the elliptic PDE

−∇ · p∇u + qu = f ∀ r ∈ Ω ⊂ ℝ^{n}  (2.1)

where p = p(r) > 0 while q = q(r) ≥ 0. Moreover u is subject to boundary conditions on Γ = ∂Ω, that we shall take of the Dirichlet type for the sake of simplicity, namely

u = 0 on Γ  (2.2)

In typical neutronics applications, u would be the neutron flux, while p and q would be the diffusion coefficient and the removal cross section, respectively, usually assumed to be piecewise constant. To make things easier to visualize, we shall take n = 2 and assume that Ω is of the “union of rectangles” type. Ω can be completely covered by I vertical and J horizontal slices. The basic analytical nodal formalism consists first in “transverse-integrating” (2.1) over (n − 1) dimensions. This is done successively and leads to n 1D equations in the different directions (here x and y). Assume that Ω is the union of rectangular cells C_{ij} of size h × k (or more generally h_{i} × k_{j}) belonging to the intersection of the i^{th} vertical and j^{th} horizontal slices mentioned above. Now write (2.1) over on such node in dimensionless variables (x and y again) as

−\frac{4}{h^{2}} (pu_{x})_{x} − \frac{4}{k^{2}} (pu_{y})_{y} + qu = f ∀(x, y) ∈ Ĉ = [−1, 1]^{2}  (2.3)

Let us now introduce some notation which will be quite useful in the following for the i^{th}-order vertical moment of a function w of x and y, which is thus still a function of x, namely:

m_{y}^{i}(w; x) = \frac{1}{N_{i}} ∫_{−1}^{1} P_{i}(y)w(x, y)dy,

and similarly:

m_{x}^{i}(w; y) = \frac{1}{N_{i}} ∫_{−1}^{1} P_{i}(x)w(x, y)dx

where P_{i} is the normalized Legendre polynomial of degree i over [−1, 1] while

N_{i} = \frac{2}{2i + 1}

is a convenient normalization factor. In particular, we shall define the following edge moments of w:

m_{L}^{i}(w) = m_{y}^{i}(w; x = −1),
m_{R}^{i}(w) = m_{y}^{i}(w; x = +1),
m_{B}^{i}(w) = m_{x}^{i}(w; y = −1),
m_{T}^{i} (w) = m_{x}^{i}(w; y = +1),

where L, R, B and T stand for Left, Right, Bottom and Top, respectively. Cell moments of w can also be defined as

m_{C}^{ij}(w) = \frac{1}{N_{i}N_{j}} ∫_{−1}^{1} ∫_{−1}^{1} P_{i}(x)P_{j}(y)w(x, y)dx dy  (2.4)

As we said before, P_{i} will be the normalized Legendre polynomial of degree i on [−1, 1] with the well known properties:

P_{i}(+1) = 1 ,
P_{i}(−1) = (−1)^{i},
P_{i}(−x) = (−1)^{i}P_{i}(x),
∫_{−1}^{1} P_{i}(x)P_{j}(x)dx = δ_{ij}N_{i},
P_{ij}(x, y) = P_{i}(x)P_{j}(y) .

In the following, we shall frequently refer to some basic spaces of polynomials in 2D, namely

P_{i} = {x^{a}y^{b}/ 0 ≤ a + b ≤ i i ∈ ℕ},
Q_{i,j} = {x^{a}y^{b}/ 0 ≤ a ≤ i 0 ≤ b ≤ j i, j ∈ ℕ} ,
Q_{i} ≡ Q_{i,i} .

Multiplying (2.3) by P_{i}(y) and dividing by N_{i} (i = 0, …, l), we obtain after transverse integration in the y direction

−\frac{4}{h^{2}} \frac{d}{dx}(p\frac{d}{dx}m_{y}^{i}(u, x)) + qm_{y}^{i}(u, x) = m_{y}^{i}(f, x) + \frac{4}{k^{2}} l_{y}^{i}(u, x) i = 0, …, l (2.5)

where l_{y}^{i}(u, x) is a transverse leakage depending on u (and x) at the top and bottom of the cell C_{ij}:

l_{y}^{i}(u, x) = m_{y}^{i}(\frac{d}{dy} (p\frac{d}{dy} u(x, y))) i = 0, …, l

A similar equation can be obtained by transverse integrating in the x direction, namely

−\frac{4}{k^{2}} \frac{d}{dy} (p\frac{d}{dy} m_{x}^{i}(u, x)) + qm_{x}^{i}(u, y) = m_{x}^{i}(f, y) + \frac{4}{h^{2}} l_{x}^{i}(u, y) i = 0, …, l (2.6)

with

l_{x}^{i}(u, y) = m_{x}^{i}(\frac{d}{dx}(p\frac{d}{dx}u(x, y))) i = 0, …, l

and clearly compatibility conditions must be satisfied, requiring that m_{x}^{i} and m_{y}^{j} commute

m_{x}^{i}(m_{y}^{j}(w; x)) = m_{y}^{j}(m_{x}^{i}(w; y)) = m_{C}^{ij}(w) i, j = 0, …, l

The numerical solution of (2.5) for the J horizontal slices and of (2.6) for the I vertical ones is then intended as follows. To simplify the exposition, we shall first consider that each of them reduces to the following canonical form

−\frac{d}{dx}(p\frac{d}{dx})u + qu = f ∀x ∈ [a, b]  (2.7)

on the corresponding slice, where both p and q are piecewise constant and that due to (2.2), u is zero at the extremities of the slice

u(a) = u(b) = 0

Equation (2.7) is first replaced by a “weak” form obtained by multiplying it by some test function v(x) and integrating it over [a, b]. As with standard finite element methods, the original problem is equivalent to the following problem:

Find u ∈ H_{0}^{1}[a, b] such that a(u, v) = f(v)for all v ∈ H_{0}^{1}[a, b]  (2.8)

where

a(u, v) = ∫_{a}^{b} (p\frac{du}{dx} \frac{dv}{dx} + quv)dx f(v) = ∫_{a}^{b} fv dx  (2.9)

Some discretization must be introduced, namely a finite dimensional subspace S_{h} of H_{0}^{1}[a, b] is defined and (2.8) is replaced by: Find u_{h} ∈ S_{h} ⊂ H_{0}^{1}[a, b] such that

a(u_{h}, v_{h}) = f(v_{h}) ∀v_{h} ∈ S_{h} ⊂ H_{0}^{1}[a, b]  (2.10)

In our description of analytical nodal methods, the basis aspects of standard f.e.m.’s will all be preserved, except that the definition of D_{l} and S_{l} will be slightly different:

D_{l} = {m_{L,}m_{R,}m_{C}^{i}, i = 0, …, l}, l ∈ ℕ dim D_{l} = l + 3 .  (2.11)

If the particular element I_{i} = [x_{i−1}, x_{i}] of size h_{i} is mapped onto the reference interval Ĉ ≡ [−1, 1], equation (2.7) becomes in reduced variables [6]

−u_{xx} + λ_{i}^{2}u = f ∀x ∈ [−1, 1]  (2.12)

where λ_{i}^{2} = \frac{h_{i}^{2}q_{i}}{4p_{i}}, while f stands now for f = \frac{h_{i}^{2}f_{i}}{4p_{i}}. Over I, we thus have (with λ = λ_{i})

S_{l} = {exp(λx), exp(−λx), 1, …, x^{l}} l ∈ ℕ dim S_{l} = l + 3  (2.13)

Theorem 2.1 (m1). D_{l} is S_{l} unisolvent

Proof. With card D_{l} = dim S_{l}, it is sufficient to exhibit the corresponding basis functions let us call

Q_{i}(x) = P_{i}(x) i = 0, …, l  (2.14)

Q_{l+1}(x) = \begin{cases}cos h_{l}^{∗}(λx) & l is odd \\ sin h_{l}^{∗}(λx) & l is even \end{cases}

Q_{l+2}(x) = \begin{cases}sin h_{l}^{∗}(λx) & l is odd \\ cos h_{l}^{∗}(λx) & l is even, \end{cases}

where the star (*) indicates that the corresponding functions are normalized to +1 when x=1. Equivalently

Q_{l+i}(x) = \frac{exp_{l}(λx) + (−1)^{l+i+1} exp_{l}(−λx)}{exp_{l}(λ) + (−1)^{l+i+1} exp_{l}(−λ)} i = 1, 2  (2.15)

where

exp_{l}(λx) = exp(λx) − ∑_{i=0}^{l} \frac{1}{N_{i}} ∫_{−1}^{1} P_{i}(x) exp(λx)dx  (2.16)

the basis functions are

u_{L}(x) = −\frac{1}{2}(−1)^{l+1}(Q_{l+2}(x) − Q_{l+1}(x))  (2.17)
u_{R}(x) = +\frac{1}{2}(Q_{l+2}(x) + Q_{l+1}(x))  (2.18)
u_{C}^{i}(x) = Q_{i}(x) − Q_{l+m(i)}(x) i = 0, …, l ,  (2.19)

where m(i) = 1 or 2 is such that i and l + m(i) have the same party. □

3. Direct analytical nodal methods

In correspondence with polynomial nodal schemes of index l [7], the analytical analog schemes has been introduced by Hennart [6]. Let us define

D_{l} ≡ {m_{L}^{i}, m_{R}^{i}, m_{B}^{i}, m_{T}^{i} , i = 0, …, l; m_{C}^{ij}, i, j = 0, …, l}  (3.1)
S_{l} = {Q_{l,l}⊕ < exp_{l}(+λx)P_{i}(y), exp_{l}(−λx)P_{i}(y), P_{i}(x)exp_{l}(+λy),  (3.2)
P_{i}(x) exp_{l}(−λy), i = 0, …, l}  (3.3)
card D_{l} = dim S_{l}  (3.4)

Theorem 3.1 (m1). D_{l} is S_{l} unisolvent.

Proof. The basis functions are

u_{L}^{i}(x, y) = −\frac{1}{2}(−1)^{l+1}(Q_{l+2,i}(x, y) − Q_{l+1,i}(x, y)) i = 0, …, l  (3.5)
u_{R}^{i}(x, y) = +\frac{1}{2}(Q_{l+2,i}(x, y) + Q_{l+1,i}(x, y)) i = 0, …, l  (3.6)
u_{B}^{i}(y, x) = −\frac{1}{2}(−1)^{l+1}(Q_{i,l+2}(x, y) − Q_{i,l+1}(x, y)) i = 0, …, l  (3.7)
u_{T}^{i} (y, x) = +\frac{1}{2}(Q_{i,l+2}(x, y) + Q_{i,l+1}(x, y)) i = 0, …, l  (3.8)
u_{C}^{ij}(x, y) = Q_{ij}(x, y) − Q_{l+m(i),j}(x, y) − Q_{i,l+m(j)}(x, y) i, j = 0, …, l (3.9)

where

Q_{l+m(i),j}(x, y) = Q_{l+m(i)}(x)P_{j}(y)Q_{i,l+m(j)}(x, y) = P_{i}(x)Q_{l+m(j)}(y) i, j = 0, …, l

and where m(i) (resp m(j)) = 1 or 2 is such that i (resp j) and l + m(i) (resp l + m(j)) have the same party. □

We can write the basis functions as follow

u_{L}^{i}(x, y) = u_{L}(x)P_{i}(y) i = 0, …, l  (3.10)
u_{R}^{i}(x, y) = u_{R}(x)P_{i}(y) i = 0, …, l  (3.11)
u_{B}^{i}(x, y) = u_{B}(y)P_{i}(x) i = 0, …, l  (3.12)
u_{T}^{i} (x, y) = u_{T} (y)P_{i}(x) i = 0, …, l  (3.13)
\begin{matrix}u_{C}^{ij}(x, y) = P_{i}(x)P_{j}(y) − Q_{l+m(i)}(x)P_{j}(y) − P_{i}(x)Q_{l+m(j)}(y) \\ i, j = 0, …, l \end{matrix}  (3.14)

In [14], analytical nodal methods were developed which relied on more physical arguments. They will be called in the following physical analytical nodal methods (§3.1) in contrast with the mathematical analytical nodal methods (§3.2).

3.1. Physical analytical nodal method. In most papers dealing with nodal methods, the final nodal equations are derived in agreement with the physics of the problem. Using the same basis functions as above, but where the equations are derived through more physical arguments, namely:

- Cell balance equations: the first moments of the equation P_{ij,} i, j = 0, …, l over a cell are zero

∫_{K∈T_{h}} P_{i}(x)P_{j}(y)[Lu_{h} − f]dr = 0 i, j = 0, …, l  (3.15)

-Current continuity conditions: with the basis functions defined above, the first (l +1) moments of the function u_{h} are automatically continuous from one cell to its neighbors, the current continuity conditions simply ensure that the corresponding moments of p∇u are continuous through cell interfaces

∫_{K_{1}∩K_{2}} P_{i}[p∇u_{h}|_{K_{1}} − p∇u_{h}|_{K_{2}}].1_{n}ds = 0 i = 0, …, l  (3.16)

where 1_{n} is some arbitrary unit normal to Γ_{12} = K_{1} ∩ K_{2}

3.2. Mathematical analytical nodal method. A mathematical way to approximate the solution of the problem at hand consists in considering the weak form of (2.1) in other words, u ∈ H_{0}^{1}(Ω) is looked for such that

a(u, v) = f(v) ∀v ∈ H_{0}^{1}(Ω)

where a(., .) and f(.) are the usual bilinear and linear forms

a(u, v) = ∫_{Ω}(p∇u∇v + quv)dx, f(v) = ∫_{Ω} fv dx

In standard, i.e. conforming situations, a finite-dimenional S_{h} subspace of H_{0}^{1}(Ω) is considered and a discretized form of the above problem is to find u_{h} ∈ S_{h} such that

a(u_{h}, v_{h}) = f(v_{h}) ∀v_{h} ∈ S_{h}  (3.17)

if we choose to approximate H_{0}^{1}(Ω) by the finite-dimensional space S_{h} (or S_{l} element by element defined here above, S_{h} ⊄ H_{0}^{1}(Ω) because u_{h} is not continuous through the faces of K (only moments of order up to l are). We are thus in a nonconforming situation, where the discrete equations (3.17) can still be used provided a(., .) is replaced by a_{h}(., .) where

a_{h}(u, v) = ∑_{K∈T_{h}} ∫_{K}(p∇u∇v + quv)dx  (3.18)

since a(u, v) would have no meaning for u, v ∈ S_{h}. The convergence of this approximate solution depend on two components: the first one has to do with the approximation properties available in S_{h} globally or rather in S_{l} element by element. These approximation properties depend on the interpolation capablilities of S_{l}, and these in term rely on the highest index i such that P_{i} ⊂ S_{l}. The second component in the L^{2} error is a consistency error due to the nonconformity of S_{h}. More precisely, if P_{l+1} ⊂ S_{l} ∀l ∈ ℕ and that moreover moments of order up to l are common between two neighboring cells so that a Patch Test of order l is passed, the non-conforming schemes exibit convergence of O(h^{l+2}) in norme L^{2}.

General properties of the corresponding nonconforming schemes.

Basis functions associated to boundary values. Let us consider for instance the basis functions u_{L}^{j}(x, y) associated to the left (x = −1) boundary values m_{L}^{j}, j = 0, …, l. By (2.14), u_{L}^{j}(x, y) = 0 as soon as x is one of the (l + 2) generalized right Radau points x_{i}^{RR}, i = 1, …, l + 2 including x = +1 [6], or y is one of the j (assuming j > 0) Gauss y_{m}^{G} m = 1, …, j [7].

Basis functions associated to cell values. Let us consider basis functions associated to cell values. By (3.14) we have

\begin{aligned}u_{C}^{ij}(x, y) & = P_{ij}(x, y) − Q_{l+m(i)}(x)P_{j}(y) + P_{ij}(x, y) \\ & − P_{i}(x)Q_{l+m(j)}(y) − P_{ij}(x, y) \end{aligned}  (3.19)

combining (2.19) and (3.19), we obtain

u_{C}^{ij}(x, y) = u_{C}^{i}(x)P_{j}(y) + P_{i}(x)u_{C}^{j}(y) − P_{ij}(x, y) .  (3.20)

for example we have

u_{C}^{ij}(x, +1) = u_{C}^{i}(x) − P_{i}(x)  (3.21)

combining (2.19) and (3.21), we obtain

u_{C}^{ij}(x, +1) = −Q_{l+m(i)}(x)  (3.22)

denote T_{0} as follow

T_{0}(x) = −((−1)^{i}u_{L}(x) + u_{R}(x))  (3.23)

combining (3.5) and (3.6) we obtain

T_{0}(x) = −\frac{1}{2}[−(−1)^{l+i+1}(Q_{l+2}(x) − Q_{l+1}(x)) + (Q_{l+2}(x) + Q_{l+1}(x))] (3.24)

then

T_{0}(x) = \begin{cases}−Q_{l+1}(x) & if l + i + 1 is even \\ −Q_{l+2}(x) & if l + i + 1 is odd \end{cases}

in other words

T_{0}(x) = −Q_{l+m(i)}(x)  (3.25)

where m(i) = 1 or 2 is such that i (resp j) and l + m(i) (resp l + m(j)) have the same party. combining (3.22), (3.23) and (3.25) we have

u_{C}^{ij}(x, +1) = −((−1)^{i}u_{L}(x) + u_{R}(x))  (3.26)

Similarly, we have

u_{C}^{ij}(x, −1) = (−1)^{i}u_{L}(x) + u_{R}(x)  (3.27)
u_{C}^{ij}(+1, y) = −((−1)^{j}u_{B}(y) + u_{T} (y))  (3.28)
u_{C}^{ij}(−1, y) = (−1)^{j}u_{B}(y) + u_{T} (y)  (3.29)

(3.30)

Relationships between the boundary and cell basis functions. Denote T_{1} and T_{2} as follow

T_{1}(x, y) = (−1)^{i}u_{L}^{j}(x, y) + u_{R}^{j}(x, y)  (3.31)
T_{2}(x, y) = (−1)^{j}u_{B}^{i}(x, y) + u_{T}^{i} (x, y)  (3.32)

by (2.14) and (3.11)

T_{1}(x, y) = ((−1)^{i}u_{L}(x) + u_{R}(x))P_{j}(y)  (3.33)

combining (3.22), (3.33) we have

T_{1}(x, y) = P_{j}(y)Q_{l+m(i)}(x)  (3.34)

similarly

T_{2}(x, y) = P_{i}(x)Q_{l+m(j)}(y)  (3.35)

Using (3.31)-(3.35) and the definition of the boundary and cell basis functions, we obtain the following relationships

\begin{matrix}P_{ij}(x, y) = & (−1)^{i}u_{L}^{j}(x, y) + u_{R}^{j}(x, y) + (−1)^{j}u_{B}^{i}(x, y) \\ +u_{T}^{i} (x, y) + u_{C}^{ij}(x, y) i, j = 0, …, l \end{matrix}

(3.36)

Elementary mass and stiffness matrices. Let us introduce the elementary mass and stiffness matrices M^{e} and K^{e} of order N as

M^{e} = (m_{ij}^{e}), K^{e} = (k_{ij}^{e})

where

m_{ij}^{e} = ∫_{[−1,1]^{2}} u_{i}(x, y)u_{j}(x, y)dx dy  (3.37)
k_{ij}^{e} = ∫_{[−1,1]^{2}} ∇u_{i}(x, y)∇u_{j}(x, y)dxdy  (3.38)

assume now that the vector ^{t}u = [u_{1}, …, u_{N}] of basis functions is partitioned as

tu = [^{t}u_{H} ^{t}u_{C} ^{t}u_{V}]

where

tu_{H} = [u_{L}^{0}, …, u_{L}^{l}, u_{R}^{0}, …, u_{R}^{l}]  (3.39)
tu_{C} = [u_{C}^{00}, u_{C}^{10}, …, u_{C}^{ll}],  (3.40)
tu_{V} = [u_{B}^{0}, …, u_{B}^{l}, u_{T}^{0} , …, u_{T}^{l}]  (3.41)

with M^{e} partitioned accordingly. We have for instance

M^{e} = [\begin{matrix}M_{HH}^{e} & M_{HC}^{e} & M_{HV}^{e} \\ M_{CH}^{e} & M_{CC}^{e} & M_{CV}^{e} \\ M_{V H}^{e} & M_{V C}^{e} & M_{V V}^{e} \end{matrix} ]  (3.42)

if the matrix elements are evaluated exactly, it is easy to realize that

M_{HV}^{e} = M_{V H}^{e} = K_{HV}^{e} = K_{V H}^{e} = 0

so that the coupling between the H (for horizontal) and V (for vertical) components is only via the cell parameters.

The mass matrix can be greatly simplified if generalized reduced integration is used, each time a boundary basis function u_{E} where E is L, R, B or T appears, the corresponding Radau rule should have sampling abscissae or ordinates x^{RR}, x^{RL}, y^{RB} or y^{RT} respectively. As a result the elementary mass matrix M^{e} only has nonzero entries in M_{CC}^{e}, we have then thanks to (3.36)

(u_{C}^{ij}, u_{C}^{km}) = \frac{4δ_{ik}δ_{jm}}{(2i + 1)(2j + 1)}  (3.43)

we obtain a similar result for K_{e}, indeed : with u_{i} = u_{E}^{i}(E = L, R, B, T) and u_{j} = u_{G}^{j}(G = L, R, B, T, C),∇u_{E}^{i} is proportional to P_{i}(i = 1, …, l) or to (Q_{k+1} ± Q_{k+2})(x or y ). In the first case, ∇u_{E} is orthogonal to ∇u_{G}^{j} ; in the second case we use the generalized Radau reduced integration.

Theorem 3.2. Assuming that p and q are cellwise constant, the physical analytical nodal method is equivalent to the mathematical analytical nodal method with generalized Radau reduced integration.

Proof. a. The mathematical analytical nodal method consists of finding u_{h} ∈ S_{h} ⊄ H_{0}^{1}(Ω)

a_{h}(u_{h}, v_{h}) = f(v_{h}) ∀v_{h} ∈ S_{h} ⊄ H_{0}^{1}(Ω)  (3.44)

Let us take the cell basis functions u_{C}^{ij}, i, j = 0, …, l associated to some element K. Equation (3.44) becomes:

∫_{K}(p∇u_{h}^{T} ∇u_{C}^{ij} + qu_{h}u_{C}^{ij} − fu_{C}^{ij})dr = 0  (3.45)

or equivalently

∫_{∂K} u_{C}^{ij}(p∇u_{h}.1_{n})ds + ∫_{K} u_{C}^{ij}(Lu − f)dr = 0  (3.46)

where 1_{n} is some arbitrary unit normal to ∂K. If K = [−1, 1]^{2} then

\begin{matrix}∫_{∂K} u_{C}^{ij}(p∇u_{h}.1_{n})ds \\ = p ∫_{−1}^{1} u_{C}^{ij}(1, y)\frac{du_{h}}{dx} |_{x=1}dy − p ∫_{−1}^{1} u_{C}^{ij}(−1, y)\frac{du_{h}}{dx} |_{x=−1}dy \\ p ∫_{−1}^{1} u_{C}^{ij}(x, 1)\frac{du_{h}}{dy} |_{y=1}dx − p ∫_{−1}^{1} u_{C}^{ij}(x, −1)\frac{du_{h}}{dy} |_{x=−1}dx \end{matrix}  (3.47)

By (3.28) we have

p ∫_{−1}^{1} u_{C}^{ij}(1, y)\frac{du_{h}}{dx} |_{x=1}dy = −p ∫_{−1}^{1} u_{B}(y)\frac{du_{h}}{dx} |_{x=1}dy − p(−1)^{j} ∫_{−1}^{1} u_{T} (y)\frac{du_{h}}{dx} |_{x=1}dy

(3.48)

then, thanks to the quadrature rule based on the (l + 2) zero y_{BR}^{i} (resp y_{T R}^{i}) (i = 1, …, l + 2) of u_{B}(y) (resp u_{T} (y)) in [−1, 1]

p ∫_{−1}^{1} u_{C}^{ij}(1, y)\frac{du_{h}}{dx} |_{x=1}dy = 0  (3.49)

thus, it is easy to verify that the boundary term disappear in (3.46) with generalized Radau reduced integration. By (3.36) we have

u_{C}^{ij}(x, y) = P_{ij}(x, y) − (−1)^{i}u_{L}^{j}(x, y) − u_{R}^{j}(x, y) − (−1)^{j}u_{B}^{i}(x, y) − u_{T}^{i} (x, y) (3.50)

where i, j = 0, …, l. Equation (3.45) under generalized Radau reduced integration finally becomes

∫_{K} P_{ij}(x, y)(∇(p∇u_{h}) + qu_{h} − f)dr = 0 i, j = 0, …, l  (3.51)

b. Let us now take v_{h} the boundary basis functions u_{E}^{i}, i = 0, …, l where E = L, R, B or T. For simplicity, let us consider

u_{E}^{i} = \begin{cases}u_{R}^{i} & on K_{1} \\ u_{L}^{i} & on K_{2} \end{cases}

where Γ_{12} = K_{1} ∩ K_{2}. Then (3.44) becomes

∫_{K_{1}} (p∇u_{h}^{T} ∇u_{R}^{i} + qu_{h}u_{R}^{i} − f u_{R}^{i})dr + ∫_{K_{2}} (p∇u_{h}^{T} ∇u_{L}^{i} + qu_{h}u_{L}^{i} − f u_{L}^{i})dr = 0

(3.52)

or equivalently, after generalized Radau reduced integration

∫_{∂K_{1}} u_{R}^{i}(p∇u_{h}.1_{n_{1}})ds + ∫_{∂K_{2}} u_{L}^{i}(p∇u_{h}.1_{n_{2}})ds = 0  (3.53)

because of generalized Radau reduced integration, (3.53) reduces to current continuity conditions on Γ_{12},

∫_{Γ_{12}} P_{i}(y)(p∇u_{h}|_{K_{1}} − p∇u_{h}|_{K_{2}}).1_{n}ds = 0 i = 0, …, l  (3.54)
□

4. Iterative procedure

The numerical results with (l = 0) are tested using the diffusion equation

−Δu + u = f_{i} ∀r ∈ Ω = [−1, 1]^{2}  (4.1)
u = 0 in Γ .  (4.2)

assume that Ω is the union of rectangular cells K of size h × k belonging to the intersection of the I vertical and J horizontal slices. We number the cells of left on the right and of bottom in top. By using the physical nodal methods polynomial [7] or direct analytical (§3) version (without numerical integration), we get the following system:

AU + BV = 0
B^{t}U + EV + CW = F C^{t}V + ZW = 0

Assume that U, V and W are the vector of horizontal, center and vertical moments. After simplification, we get the following system

HV = F  (4.3)

where

H = E − B^{t}A^{−1}B − CZ^{−1}C^{t}
U = −A^{−1}BV
W = −Z^{−1}C^{t}V

where A and Z are tridiagonal matrices, B and C are two diagonals, the components are either zero or satisfy

Z_{ii} = A_{ii}
Z_{i,i+I} = A_{i,i+1}
B_{ii} = C_{ii}
B_{i,i+1} = C_{i,i+I}
E = cI_{d}

where c is a constant, I_{d} is an identity matrix.

Here H is a matrices of order I × J. The diagonal blocks of H are matrices of order I × I and corresponds to the J horizontal slices (for the sake of simplicity we suppose I = J) and take the form: (see fig 4.1)

• The blocks H_{ii}, i = 1, …, I are full symmetric, positive-definite.

• The blocks H_{ij} = α_{ij}I_{d}, i, j = 1, …, I, i ≠ j, α_{ij} is constant. Thus, to part the inversion of the H_{ii}, i = 1, …, I blocks, the resulted linear system is solved by the blocks Gauss-Seidel iterative method where the blocks of index ij, i ≠ j are identified to points. In addition, we must declared only A and B.

H = (\begin{matrix}(\begin{matrix}. & . & . \\ . & . & . \\ . & . & . \end{matrix}) & α_{12} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) & α_{13} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) \\ α_{21} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) & (\begin{matrix}. & . & . \\ . & . & . \\ . & . & . \end{matrix}) & α_{23} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) \\ α_{31} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) & α_{32} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) & (\begin{matrix}. & . & . \\ . & . & . \\ . & . & . \end{matrix}) \end{matrix})

fig 4.1: General shape of the H matrix for I = 3

5. Numerical results for l = 0

Consider three problems P_{i} i = 1, 2, 3 in equation (4.1), the first problem P_{1} has an exact solution

u_{1}(x, y) = (1 − x^{2})(1 − y^{2})
f_{1}(x, y) = 5 − 3(x^{2} + y^{2}) + x^{2}y^{2} .

The second problem, P_{2}, has an exact solution

u_{2}(x, y) = (1 − x^{4})(1 − y^{4})
f_{2}(x, y) = 12x^{2}(1 − y^{4}) + 12y^{2}(1 − x^{4}) − (x^{4} + y^{4}) + x^{4}y^{4} .

The last problem, P_{3}, has an exact solution

u_{3}(x, y) = (1 − \frac{chy}{ch1} )(1 − \frac{chx}{ch1})

f_{3}(x, y) = 1 − \frac{chxchy}{ch^{2}1} .

The numerical results obtained by direct physical analytical schemes are displayed in Table 1. Those obtained by physical polynomial schemes are displayed in Table 2. Here ε_{i} denote the error in discrete norm in L^{2} associated to the problem P_{i} i = 1, 2, 3 ε_{i}(c) denote the error in continuous norm in L^{2},

ε_{i} = ( ∑_{K∈T_{h}}(m_{C}(u_{i}) − m_{C}(u_{h_{i}}))^{2})^{1/2},  (5.1)

ε_{i}(c) = ( ∫_{Ω}(u_{i} − u_{h_{i}})^{2})^{1/2}  (5.2)

From these values of the error, we deduct the evaluation of the optimal order of convergence while calculating for every grid the numerical exponent α of

ε ∝ h^{α}  (5.3)

h being the size of the mesh. As passing the mesh i of size h_{i} to mesh i + 1 of size h_{i+1} = \frac{1}{2}h_{i}, we have

\frac{ε_{i}}{ε_{i+1}} = 2^{α}

Table 1. Direct physical analytical nodal methods for problems P_{1} and P_{3}

mesh  ε_{1}  α_{1}  ε_{3}  α_{3}
2x2  0.01300  0.06375
4x4  0.00303 2.07045 0.01440 2.1034
8x8  0.00074 2.02865 0.00347 2.0377
16x16 0.00018 2.00806 0.00086 2.0101

Table 2. Direct physical polynomial nodal methods for problems P_{1}, P_{2} and P_{3}

mesh ε_{1}(c) α_{1}(c) ε_{1} ε_{2}(c) α_{2}(c)  ε_{2}  α_{2}  ε_{3}  α_{3}
2 x2 0.4787  0 1,8839  0,0468  0.00632
4 x 4 0.1265 1,9435 0 0.4971 1,9466 0,0165 1.7768 0.00214 1.72
8 x 8 0.0321 1,9848 0 0,1272 1,9767 0,0045 1,9366 0.00057 1.94
16x16 0.0081 1,9961 0 0,0320 1.9933 0,0012 1,9837 0.00015 1.95

6. Conclusions

(1) - We analyzed the direct analytical nodal methods of index l, we showed that the corresponding mathematical methods are equivalent to the physical ones when components of the matrix are calculated by generalized Radau reduced integration (exact for polynomials of degree (2l + 2)).

(2) - We compare in the table2, the discrete (ε_{i}) and continuous (ε_{i}(c)) i = 1, 2, 3 error:

• ε_{1} = 0
• ε_{2}(c) ≪ ε_{2}

in the first case we have the similar result for mathematical method because of generalized Radau reduced integration is exact for polynomial of degree 2, in second case we have a superconvergence.

(3) - We notice that results of the polynomial method are better than the analytical ones (note that these last are O(h^{l+3−δ_{l0}}) order in norm L^{2}). Indeed we have :

• ε_{1} = 0 in table2 (polynomial method), in table1 (analytical methods), ε_{1} is small but far from being zero.

• ε_{3}(table 1) ≈ 6ε_{3}(table 2).

References

[1] Ciarlet P. G.; The finite element method for elliptic problems, North-Holland, Amsterdam, 1978.

[2] Dorning J. J.; Modern coarse-mesh methds. A development of the 70’s, Computational methds in Nuclear Engineering, vol. 1, pp. 3.1-3.31 ,American Nuclear Society, Williamsburg, Virginia (1979).

[3] Fedon magnand, Hennart J. P., and Lautard J. J.; On the relationship between some nodal schemes and the finite element method in static diffusion calculations, Advances in Reactor Computations, Vol. 2, Salt Lake City, UT, 1983, pp. 987-1000.

[4] Finnemann H., Bennewitz F., and Wagner M.R.; Interface current techniques for multidimensional reactor calculations, Atomkernenergie 30, 123- 128 (1977).

[5] Frohlich R.; Summary discussion and state of the art review for coarse-mesh computational methods, Atomkernenergie 13n 117-126 (1983).

[6] Hennart J. P.; On the numerical analysis of analytical nodal methods, Numerical methods for partial differential equations, 4, 233-254 (1988).

[7] Hennart J. P.; A general family of Nodal Schemes, SIAM J.Sci.Stat.Comput. 7.264 (1986).

[8] Hennart J.P.; A Unified view of finite differences, finite elements, and nodal schemes, National University of Mexico, IIMAS, Mexico City, 01000 Mexico 1986.

[9] Hennart J. P., On the numerical analysis of analytical nodal methods, Numerical methods for partial differential equations, 4, 233-254(1988).

[10] Hennart J.P.; Nodal method for the numerical solution of partial differential Equations, Numerical Analysis, Guanajuato, p. 175-190, 51984. J. P. Hennart, ed., Lecture Notes in Mathematics 1230 Springer-Verlag, Berlin, 1987.

[11] Langenbuch,Maurer W., and Werner W.; coarse-mesh flux-expansion method for the analysis of space-time effects in large light water reactor cores, Nucl. Sci.Engng., 63 (1977), pp. 437-456.

[12] Langenbuch S., W. Maurer, and Werner W.; High order shemes for neutron kinetics calculations, based on local polynomial approximation, Nucl.Sci.Engng, 64 (1977), pp. 508-516.

[13] Meade D., Del Valle E. and Hennart J.P.; Unconventional finite element methods for neutron group diffusion equations, IIMAS-UNAM (1986).

[14] Shober R. A., Simms R. C., and Henry A. F.; Tho Nodal Methods for Solving Time Dependent Group Diffusion Equations, Nucl. Sci. Engng.. 64. 582 (1977).

[15] Wagner M. R. and Kobke K.; Progress in nodal reactor analysis, Atomkernenergie 13,117-126 (1983).

Département de mathématiques et Informatique, Faculté des sciences Dhar-Mehraz, Université S.M. Ben Abdellah, B.P. 1796 Fès-Atlas, Fès, Morocco

E-mail address, N. Guessous: ngessous@yahoo.fr E-mail address, F. Hadfat: fhadfat@yahoo.fr