\documentclass[reqno]{amsart} \AtBeginDocument{{\noindent\small 2004-Fez conference on Differential Equations and Mechanics \newline {\em Electronic Journal of Differential Equations}, Conference 11, 2004, pp. 143--155.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu/ or http://ejde.math.unt.edu/ \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2004 Texas State University - San Marcos.} \vspace{9mm}} \setcounter{page}{143} \begin{document} \title[\hfilneg EJDE/Conf/11 \hfil analytical nodal methods for] {Analytical nodal methods for diffusion equations} \author[N. Guessous, F. Hadfat\hfil EJDE/Conf/11 \hfilneg] {Najib Guessous, Fouzia Hadfat} % in alphabetical order \address{D\'epartement de math\'ematiques et Informatique\\ Facult\'e des sciences Dhar-Mehraz\\ Universit\'e S.M. Ben Abdellah\\ B.P. 1796 F\`es-Atlas, F\`es, Morocco} \email[N. Guessous]{ngessous@yahoo.fr} \email[F. Hadfat]{fhadfat@yahoo.fr} \date{} \subjclass[2000]{74S30} \keywords{Partial differential equations of elliptic type; numerical analysis; \hfill\break\indent diffusion equations; finite element methods} \thanks{Published October 15, 2004.} \begin{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) \cite{h2} obtained better results that those of ANM. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}[theorem]{Remark} \section{Introduction} Broadly speaking, nodal methods are fast and accurate methods which combine attractive features of the finite element method (f.m.e.) \cite{f1} as well as of the finite difference method (f.d.m.) \cite{h3}. 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. \cite{f2}, Langenbuch et al. \cite{l1,l2}, Shober et al \cite{s1}, Frohlich \cite{f3}, Dorning \cite{d1} and Wagner and Kobke \cite{w1}. In \S 2, we present the transverse integrated version of the analytical nodal methods of index $l$ (ANM) \cite{h1}. In \S 3 we analyze the direct procedure of ANM introduced in \cite{h1} 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 \S 4, a Gauss-seidel iterative method. In\S 5, Some numerical examples applied to the diffusion equation one permitted to notice in case $l=0$, that the polynomial nodal methods (PNM) \cite{h2} succeeds to the best results that those of the ANM. \section{Basic analytical nodal formalism} Let us consider the elliptic PDE \begin{equation} -\nabla \cdot p\nabla u+qu = f \quad\forall \;r\in \Omega \subset \mathbb{R}^n\, \label{e1a} \end{equation} where $p = p (r) >0$ while $q =q (r) \ge 0$. Moreover $u$ is subject to boundary conditions on $\Gamma =\partial \Omega $, that we shall take of the Dirichlet type for the sake of simplicity, namely \begin{equation} u = 0 \quad \mbox{ on }\Gamma \label{e1b} \end{equation} 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 $\Omega $ is of the ``union of rectangles" type. $\Omega $ can be completely covered by $I$ vertical and $J$ horizontal slices. The basic analytical nodal formalism consists first in ``transverse-integrating" \eqref{e1a} over $(n-1)$ dimensions. This is done successively and leads to $n$ 1D equations in the different directions (here $x$ and $y$). Assume that $\Omega $ is the union of rectangular cells $C_{ij}$ of size $h \times k$ (or more generally $h_{i} \times k_{j})$ belonging to the intersection of the $i^{th}$ vertical and $j^{th}$ horizontal slices mentioned above. Now write \eqref{e1a} over on such node in dimensionless variables ($x$ and $y$ again) as \begin{equation} -\frac{4}{h^{2}}(pu_{x})_{x}-\frac{4}{k^{2}}(pu_{y})_{y}+qu=f \quad\forall (x,y)\in \hat{C}=[-1,1]^{2}\, \label{e2} \end{equation} 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 }\int_{-1}^1 P_{i} (y) w(x,y) dy, \] and similarly: \[ m_x^i(w;y) = \frac{1}{N_i}\int_{-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$: \begin{gather*} %2.7-2.10 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), \end{gather*} where $L$, $R$, $B$ and $T$ stand for Left, Right, Bottom and Top, respectively. Cell moments of $w$ can also be defined as \begin{equation} m_C^{ij}(w)=\frac{1}{N_i N_j} \int_{-1}^1\int_{-1}^1 P_{i}(x)P_{j} (y)w(x,y)dx\,dy \end{equation} As we said before, $P_{i}$ will be the normalized Legendre polynomial of degree i on $[-1,1]$ with the well known properties: \begin{gather*} P_{i}(+1) = 1\,,\\ P_{i}(-1) = (-1)^{i},\\ P_{i}(-x) =(-1)^{i}P_{i}(x) ,\\ \int_{-1}^1 P_{i} (x) P_{j} (x)dx = \delta _{ij} N_{i},\\ P_{ij}(x,y)=P_{i}(x) P_{j}(y)\,. \end{gather*} In the following, we shall frequently refer to some basic spaces of polynomials in 2D, namely \begin{gather*} \mathcal{P}_{i} = \{ x ^{a }y ^{b } / \quad 0\leq a + b \leq i \quad i \in \mathbb{N} \} ,\\ Q_{i,j} =\{ x^{a }y ^{b } / \quad 0 \le a \le i \quad 0 \le b \le j\quad i ,j \in \mathbb{N} \}\,,\\ Q_{i} \equiv Q_{i,i} \,. \end{gather*} Multiplying \eqref{e2} by $P_{i}(y)$ and dividing by $N_{i}$ ($i=0, \dots, l$), we obtain after transverse integration in the y direction \begin{equation} -\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) \quad i = 0, \dots , l \label{e3} \end{equation} 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}$: \begin{equation*} l_y^i (u, x) = m_y^i (\frac{d}{dy}(p\frac{d}{dy}u (x,y)))\quad i = 0, \dots , l \end{equation*} A similar equation can be obtained by transverse integrating in the $x$ direction, namely \begin{equation} -\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) \quad i = 0, \dots ,l \label{e4} \end{equation} with \begin{equation*} l_x^i (u, y) = m_x^i (\frac{d}{dx}(p\frac{d}{dx}u (x,y))) \quad i = 0, \dots , l \end{equation*} and clearly compatibility conditions must be satisfied, requiring that $m_x^i $ and $m_y^j $ commute \begin{equation*} m_x^i (m_y^j (w;x)) = m_y^j (m_x^i (w;y))=m_C^{ij} (w) \quad i, j=0, \dots , l \end{equation*} The numerical solution of \eqref{e3} for the $J$ horizontal slices and of \eqref{e4} 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 \begin{equation} -\frac{d}{dx}(p\frac{d}{dx})u + q u = f \quad \forall x\in [a,b] \label{e6} \end{equation} on the corresponding slice, where both $p$ and $q$ are piecewise constant and that due to \eqref{e1b}, $u$ is zero at the extremities of the slice \begin{equation*} u(a) = u(b) = 0\end{equation*} Equation \eqref{e6} 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: \begin{equation} \mbox{Find $u\in H_0^1 [a,b]$ such that $a(u,v)=f(v)$for all $v\in H_0^1 [a,b]$ }\label{e7} \end{equation} where \begin{equation} a(u,v) = \int_a^b {(p\frac{du}{dx}\frac{dv}{dx}} + quv ) dx \quad f(v) = \int_a^b fv\,dx\, \end{equation} Some discretization must be introduced, namely a finite dimensional subspace $S_{h}$ of $H_0^1 [a,b]$ is defined and \eqref{e7} is replaced by:\\ Find $u_{h}\in S_{h}\subset H_0^1[a,b]$ such that \begin{equation} a(u_{h},v_{h})=f(v_{h})\quad \forall v_{h}\in S_{h}\subset H_0^1 [a,b] \end{equation} 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: \begin{equation} D_{l }= \{ m_{L,} m_{R,} m_C^i , i = 0, \dots ,l\}, \quad l \in \mathbb{N}\quad \dim D_{l }=l + 3\,. \end{equation} If the particular element $I_i=[x_{i-1},x_i$] of size $h_i$ is mapped onto the reference interval $\hat{C}\equiv [-1,1]$, equation \eqref{e6} becomes in reduced variables \cite{h1} \begin{equation} -u _{xx}+\lambda ^{2}_{i} u = f \quad\forall x \in [-1,1]\, \label{e10} \end{equation} where $\lambda ^{2}_{i} = \frac{h^{2}_{i}q_{i}}{4p_{i}}$, while $f$ stands now for $f =\frac{ h^{2}_{i} f_{i}}{4p_{i}}$.\\ Over $I$, we thus have (with $\lambda =\lambda _{i})$ \begin{equation} S_{l}=\{ \exp(\lambda x), \exp (-\lambda x), 1, \dots , x^{l} \} \quad l \in \mathbb{N}\quad \dim S_{l }=l + 3\label{e11} \end{equation} \begin{theorem}[{m1}] $D_{l}$ is $S_{l}$ unisolvent \end{theorem} \begin{proof} With card $D_{l}=\dim S_{l}$, it is sufficient to exhibit the corresponding basis functions let us call \begin {equation} Q_{i}(x)=P_{i}(x) \quad i=0, \dots , l \label{e1.9} \end {equation} $$ Q_{l+1}(x)=\begin{cases} \cos h_l^*(\lambda x) & l \mbox{ is odd} \\ \sin h_l^*(\lambda x) & l \mbox{ is even} \end{cases} $$ $$ Q_{l+2}(x)= \begin{cases} \sin h_l^*(\lambda x) & l \mbox{ is odd} \\ \cos h_l^*(\lambda x) & l \mbox{ is even} , \end{cases} $$ where the star (*) indicates that the corresponding functions are normalized to +1 when x=1. Equivalently \begin{equation} Q_{l+i}(x) = \frac{\exp_{l}(\lambda x)+(-1)^{l+i+1} \exp_{l}(-\lambda x)} {\exp_{l}(\lambda ) +(-1)^{l+i+1}\exp_{l}(-\lambda )}\quad i=1,2 \label{e1.4} \end{equation} where \begin {equation} \exp_{l}(\lambda x) = \exp(\lambda x) - \sum_{i=0}^l \frac{1}{N_{i}}\int_{-1}^1 P_{i}(x) \exp(\lambda x)dx \label{e1.8} \end{equation} the basis functions are \begin {gather} u_{L}(x) = - \frac{1}{2}(-1)^{l+1} (Q_{l+2}(x) - Q_{l+1}(x))\label{e1.1} \\ u_{R}(x) = +\frac{1}{2} (Q_{l+2}(x) + Q_{l+1}(x)) \label{e1.2} \\ u_C^i (x) = Q_{i}(x) - Q_{l+m(i)}(x) \quad i = 0 , \dots , l\,, \label{e1.3} \end {gather} where $m(i) = 1$ or $2$ is such that $i$ and $l + m(i)$ have the same party. \end {proof} \section{ Direct analytical nodal methods} In correspondence with polynomial nodal schemes of index $l$ \cite{h2}, the analytical analog schemes has been introduced by Hennart \cite{h1}. Let us define \begin{gather} D_{l} \equiv \{m_L^i ,m_R^i ,m_B^i ,m_T^i ,i=0, \dots ,l; m_C^{ij} , \quad i, j= 0, \dots,l \} \label{e12}\\ S_{l}=\{ Q_{l, l} \oplus <\exp_{l}(+\lambda x)P_{i}(y), \exp_{l}(-\lambda x)P_{i}(y), P_{i}(x)exp_{l}(+\lambda y),\\ P_{i}(x) \exp_{l} (-\lambda y),\; i=0, \dots , l \} \label{e13}\\ \mathop{\rm card} D_{l}= \dim S_{l} \end{gather} \begin{theorem}[{m1}] $D_{l}$ is $S_{l}$ unisolvent. \end{theorem} \begin{proof} The basis functions are \begin {gather} u_L^i (x,y) =- \frac{1}{2} (-1)^{l+1} (Q_{l+2,i}(x,y) - Q_{l+1,i}(x,y))\quad i = 0, \dots ,l \label{e1.91}\\ u_R^i (x,y) =+ \frac{1}{2} (Q_{l+2,i}(x,y) + Q_{l+1,i}(x,y)) \quad i = 0, \dots ,l \label{e1.92}\\ u_B^i (y,x) =- \frac{1}{2} (-1)^{l+1} (Q_{i,l+2}(x,y) - Q_{i,l+1}(x,y)) \quad i = 0, \dots ,l \label{e1.93}\\ u_T^i (y,x) =+ \frac{1}{2} (Q_{i,l+2}(x,y) + Q_{i,l+1}(x,y)) \quad i = 0, \dots ,l \label{e1.94}\\ u_C^{ij}(x,y) = Q_{ij}(x,y) - Q_{l+m(i),j}(x,y) - Q_{i,l+m(j)}(x,y) \quad i, j = 0, \dots , l \label{e1.13b} \end{gather} where \begin{gather*} 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) \quad i,j=0, \dots , l \end{gather*} 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. \end{proof} We can write the basis functions as follow \begin{gather} u_L^i (x,y) =u_{L}(x) P_{i}(y)\quad i=0, \dots , l \label{e1.9b}\\ u_R^i (x,y) =u_{R}(x) P_{i}(y)\quad i=0, \dots , l\label{e1.10}\\ u_B^i (x,y) =u_{B}(y) P_{i}(x)\quad i=0, \dots , l\label{e1.11}\\ u_T^i (x,y) =u_{T}(y) P_{i}(x)\quad i=0, \dots , l\label{e1.12}\\ \begin{gathered} 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, \dots , l \end{gathered}\label{e1.13} \end {gather} In \cite{s1}, analytical nodal methods were developed which relied on more physical arguments. They will be called in the following physical analytical nodal methods (\S3.1) in contrast with the mathematical analytical nodal methods (\S3.2). \subsection{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, \dots , l$ over a cell are zero \begin{equation} \int_{{K}\in {T}_{h} } P_{i}(x)P_{j}(y)[Lu_{h}-f] dr = 0 \quad i,j=0, \dots , l \label{e28} \end{equation} -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\nabla $u are continuous through cell interfaces \begin{equation} \int_{{K}_{1} \cap {K}_{2} } P_{i} [p\nabla u_{h} | _{{K}_{1}} - p\nabla u_{h} | _{{K}_{2}} ]. 1_{n } ds = 0 \quad i=0, \dots , l\, \label{e29} \end{equation} where $1_{n }$ is some arbitrary unit normal to $\Gamma _{12}= K_{1} \cap K_{2}$ \subsection{Mathematical analytical nodal method} A mathematical way to approximate the solution of the problem at hand consists in considering the weak form of \eqref{e1a} in other words, $u\in H_0^1 (\Omega )$ is looked for such that \begin{equation*} a(u,v)=f(v) \quad \forall v\in H_0^1 (\Omega ) \end{equation*} where $a(.,.)$ and $f(.)$ are the usual bilinear and linear forms \begin{equation*} a(u,v)=\int_{\Omega}(p\nabla u \nabla v+quv)dx, \quad f(v)=\int_{\Omega}fv\,dx \end {equation*} In standard, i.e. conforming situations, a finite-dimenional $S_{h}$ subspace of $H_0^1 (\Omega )$ is considered and a discretized form of the above problem is to find $u_{h}\in S_{h}$ such that \begin{equation} a(u_{h},v_{h})=f(v_{h}) \quad \forall v_{h}\in S_{h} \label{e26} \end{equation} if we choose to approximate $H_0^1 (\Omega )$ by the finite-dimensional space $S_h$ (or $S_l$ element by element defined here above, $ S_{h}\not\subset H_0^1 (\Omega )$ 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 \eqref{e26} can still be used provided $a(.,.)$ is replaced by $a_h(.,.)$ where \begin{equation} a_{h}(u,v)=\sum_{K\in T_h } \int_{K}(p\nabla u \nabla v+quv)dx \end{equation} since $a(u,v)$ would have no meaning for $u, v\in 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 $\mathcal{P}_i\subset S_l$. The second component in the $L^2$ error is a consistency error due to the nonconformity of $S_h$. More precisely, if $\mathcal{P}_{l+1}\subset S_l \quad\forall l\in \mathbb{N}$ 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$. \subsection*{General properties of the corresponding nonconforming schemes} \subsubsection*{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, \dots , l$. By \eqref{e1.9}, $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, \dots , l+2$ including $x=+1$ \cite{h1}, or $y$ is one of the $j$ (assuming $j>0$) Gauss $y_m^G$ $m=1, \dots , j$ \cite{h2}. \subsubsection*{Basis functions associated to cell values} Let us consider basis functions associated to cell values. By \eqref{e1.13} we have \begin{equation} \begin{aligned} u_C^{ij}(x,y)&=P_{ij}(x,y)-Q_{l+m(i)}(x)P_{j}(y)+P_{ij}(x,y)\\ &\quad -P_{i} (x)Q_{l+m(j)}(y)-P_{ij}(x,y) \end{aligned}\label{e16} \end{equation} combining \eqref{e1.3} and \eqref{e16}, we obtain \begin{equation} u_C^{ij}(x,y) = u_C^i (x)P_{j}(y)+P_{i}(x)u_C^j (y)-P_{ij } (x,y) \,.\label{e17} \end{equation} for example we have \begin{equation} u_C^{ij}(x,+1) = u_C^i (x)-P_{i}(x) \label{e18}\end {equation} combining \eqref{e1.3} and \eqref{e18}, we obtain \begin{equation} u_C^{ij}(x,+1) = -Q_{l+m(i)}(x) \label{e49} \end {equation} denote $T_0$ as follow \begin {equation} T_0(x) =-((-1)^i u_L(x) + u_R(x)) \label{e1.98} \end{equation} combining \eqref{e1.91} and \eqref{e1.92} we obtain \begin {equation} 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))] \end{equation} then $$ T_0(x) = \begin{cases} -Q_{l+1}(x) & \mbox{if }l+i+1 \mbox{ is even} \\ -Q_{l+2}(x) & \mbox{if }l+i+1 \mbox{ is odd} \end{cases} $$ in other words \begin {equation} T_0(x) = -Q_{l+m(i)}(x) \label{e1.97} \end{equation} 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 \eqref{e49}, \eqref {e1.98} and \eqref{e1.97} we have \begin {equation} u_C^{ij}(x,+1) =-((-1)^i u_L(x) + u_R(x)) \label{e42} \end{equation} Similarly, we have \begin{gather} u_C^{ij}(x,-1) =(-1)^i u_L(x) + u_R(x) \label{e43} \\ u_C^{ij}(+1,y) =-((-1)^j u_B(y) + u_T(y)) \label{e40} \\ u_C^{ij}(-1,y) =(-1)^j u_B(y) + u_T(y) \label{e41} \\ \end{gather} \subsubsection*{Relationships between the boundary and cell basis functions} Denote $T_1$ and $T_2$ as follow \begin{gather} T_1(x,y)=(-1)^i u_L^j(x,y) + u_R^j(x,y)\label{e51}\\ T_2(x,y)=(-1)^j u_B^i(x,y) + u_T^i(x,y)\label{e52} \end{gather} by \eqref {e1.9} and \eqref{e1.10} \begin{equation} T_1(x,y)=((-1)^i u_L(x) + u_R(x))P_j(y) \label{e1.99} \end{equation} combining \eqref{e49}, \eqref{e1.99} we have \begin{equation} T_1(x,y)=P_j(y)Q_{l+m(i)}(x)\label{e50} \end{equation} similarly \begin{equation} T_2(x,y)=P_i(x)Q_{l+m(j)}(y)\label{e53} \end{equation} Using \eqref{e51}-\eqref{e53} and the definition of the boundary and cell basis functions, we obtain the following relationships \begin{equation} \begin{aligned} P_{ij}(x,y) =& (-1)^i u_L^j(x,y) + u_R^j(x,y) + (-1)^ju_B^i(x,y) \\ + u_T^i(x,y) + u_C^{ij}(x,y)\quad i,j=0, \dots , l \end{aligned}\label{e19} \end{equation} \subsubsection*{Elementary mass and stiffness matrices} Let us introduce the elementary mass and stiffness matrices $M^{e}$ and $K^{e}$ of order $N$ as \begin{equation*} M^{e}=(m_{ij}^e ), \quad K^{e}=(k_{ij}^e ) \end{equation*} where \begin{gather} m_{ij}^e = \int_{[-1,1]^{2}} u_{i}(x,y) u_{j}(x,y) dx \,dy \label{e20}\\ k_{ij}^e =\int_{[-1,1]^{2}} \nabla u_{i}(x,y) \nabla u_{j}(x,y) dx dy \label{e25} \end{gather} assume now that the vector $^{t}u=[u_{1}$, \dots , $u_{N}$] of basis functions is partitioned as \[ {}^{t}u=[{}^{t}u_{H} \;{}^{t}u_C \;{}^{t}u_{V}] \] where \begin{gather} ^{t}u_{H} = [u_{L}^{0} , \dots , u_{L}^{l},u_{R}^{0} ,\dots , u_{R}^{l} ] \label{e21} \\ ^tu_{C}=[u_C^{00}, u_C^{10} ,\dots ,u_C^{ll} ],\\ {}^tu_{V} = [u_B^0 , \dots , u_B^l ,u_T^0 ,\dots ,u_T^l ] \label{e22} \end{gather} with $ M^{e}$ partitioned accordingly. We have for instance \begin{equation} M^{e}=\begin{bmatrix} M_{HH}^e & M_{HC}^e & M_{HV}^e \\ M_{CH}^e & M_{CC}^e & M_{CV}^e \\ M_{VH}^e & M_{VC}^e & M_{VV}^e \\ \end{bmatrix} \label{e23} \end{equation} if the matrix elements are evaluated exactly, it is easy to realize that \begin{equation*} M_{HV}^e=M_{VH}^e=K_{HV}^e=K_{VH}^e=0 \end{equation*} 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 \eqref{e19} \begin{equation} (u_C^{ij} ,u_C^{km})= \frac{4\delta_{ik}\delta_{jm} }{(2i + 1)(2j+1)} \label{e24} \end{equation} 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)$,$\nabla u_{E}^{i} $ is proportional to $P_{i} (i=1, \dots ,l) $ or to $(Q_{k+1}\pm Q_{k+2})(x \;or \;y\;)$. In the first case, $\nabla u_{E}$ is orthogonal to $\nabla u_{G}^{j} $ ; in the second case we use the generalized Radau reduced integration. \begin{theorem} 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. \end{theorem} \begin{proof} a. The mathematical analytical nodal method consists of finding $u_{h}\in S_{h}\not\subset H_0^1 (\Omega )$ \begin{equation} a_{h}(u_{h},v_{h}) = f(v_{h})\quad \forall v_{h} \in S_{h}\not\subset H_0^1 (\Omega ) \label{e30} \end{equation} Let us take the cell basis functions $u_C^{ij}, i, j = 0, \dots ,l$ associated to some element $K$. Equation \eqref{e30} becomes: \begin{equation} \int_{K}(p\nabla u_{h}^{T} \nabla u_{C}^{ij} + qu_{h} u_{C}^{ij} -fu_{C}^{ij})dr=0 \label{e31} \end{equation} or equivalently \begin{equation} \int_{{\partial}{K}} u_C^{ij} (p\nabla u_h .1_{n})ds+\int_K {u_C^{ij} } (Lu - f )dr = 0 \label{e32} \end{equation} where $1_{n }$ is some arbitrary unit normal to $\partial{K}$. If $K=[-1,1]^2$ then \begin{equation} \label{e46} \begin{aligned} &\int_{{\partial}{K}} u_C^{ij} (p\nabla u_h .1_{n})ds\\ &= p \int_{-1}^1 u_C^{ij}(1,y)\frac{du_h}{dx}\big|_{x=1}dy -p \int_{-1}^1 u_C^{ij}(-1,y)\frac{du_h}{dx}\big|_{x=-1}dy\\ &\quad p \int_{-1}^1 u_C^{ij}(x,1)\frac{du_h}{dy}\big|_{y=1}dx -p \int_{-1}^1 u_C^{ij}(x,-1)\frac{du_h}{dy}\big|_{x=-1}dx \end{aligned} \end {equation} By \eqref{e40} we have \begin {equation} p \int_{-1}^1 u_C^{ij}(1,y)\frac{du_h}{dx}\big|_{x=1}dy= -p \int_{-1}^1 u_B(y)\frac{du_h}{dx}\big|_{x=1}dy-p(-1)^j \int_{-1}^1 u_T(y) \frac{du_h}{dx}\big|_{x=1}dy \label{e47} \end{equation} then, thanks to the quadrature rule based on the ($l+2$) zero $y_{BR}^{i}$ (resp $y_{TR}^{i}$) ($i=1,\dots, l+2$) of $u_{B}(y)$ (resp $u_{T}(y))$ in $[-1,1]$ \begin {equation} p \int_{-1}^1 u_C^{ij}(1,y)\frac{du_h}{dx}\big|_{x=1}dy= 0 \end{equation} thus, it is easy to verify that the boundary term disappear in \eqref{e32} with generalized Radau reduced integration. By \eqref{e19} we have \begin{equation} u_C^{ij}(x,y)=P_{ij}(x,y)-(-1)^i u_L^j(x,y) - u_R^j(x,y) - (-1)^ju_B^i(x,y)-u_T^i(x,y) \label{e34} \end{equation} where $i,j=0, \dots , l$. Equation \eqref{e31} under generalized Radau reduced integration finally becomes \begin{equation} \int_K {P_{ij}(x,y) (\nabla (p\nabla u_h) +\,qu_h } -f)dr=0 \quad i,j=0, \dots , l \label{e35} \end{equation} b. Let us now take $ v_{h }$ the boundary basis functions $u_E^i , i = 0, \dots , l$ where $E=L$, $R$, $B$ or $T$. For simplicity, let us consider $$ u_E^i = \begin{cases} u_R^i &\mbox{on } K_{1} \\ u_L^i &\mbox{on } K_{2} \end{cases} $$ where $\Gamma _{12}= K_{1} \cap K_{2}$. Then \eqref{e30} becomes \begin{equation} \int_{{K}_{1} } {{(p}\nabla {u}_{h}^{T} \nabla {u}_{R}^{i} +{qu}_{h} } {u}_{R}^{i} -{f}\,{u}_{R}^{i}){dr} +\int_{{K}_{2} } {{(p}\nabla {u}_{h}^{T} \nabla {u}_{L}^{i} \,+{qu}_{h} } {u}_{L}^{i} -{f}\,{u}_{L}^{i})dr=0 \label{e36} \end{equation} or equivalently, after generalized Radau reduced integration \begin{equation} \int_{{\partial}{K}_{1}}{u}_{R}^{i}(p\nabla u_h .1_{{n}_{1}})ds + \int_{{\partial}{K}_{2}}{u}_{L}^{i} (p\nabla u_h .1_{{n}_{2}})ds=0 \label{e37} \end{equation} because of generalized Radau reduced integration, \eqref{e37} reduces to current continuity conditions on $\Gamma _{12 }$, \begin{equation} \int_{\Gamma _{12} } P_{i} (y) ( p\nabla u_h |_{K_1 } - p\nabla u_h |_{K_2 } ).1_{n} ds = 0 \quad i = 0, \dots , l \label{e38} \end{equation} \end{proof} \section{Iterative procedure} The numerical results with $(l=0)$ are tested using the diffusion equation \begin{gather} -\Delta u+u = f_{i} \quad \forall r \in \Omega = [-1,1]^{2} \label{Pi}\\ u = 0 \quad \mbox {in } \Gamma \,. \end{gather} assume that $\Omega $ is the union of rectangular cells $K$ of size $h \times 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 \cite {h2} or direct analytical (\S 3) version (without numerical integration), we get the following system: \begin{gather*} AU + BV=0\\ B^{t}U+ EV + CW=F C^{t}V + ZW=0 \end{gather*} Assume that $U$, $V$ and $W$ are the vector of horizontal, center and vertical moments. After simplification, we get the following system \begin{equation} HV=F \end{equation} where \begin{gather*} H=E-B^{t}A^{-1}B-CZ^{-1}C^{t } \\ U=-A^{-1}BV \\ W=-Z^{-1}C^{t}V \end{gather*} where $A$ and $Z$ are tridiagonal matrices, $B$ and $C$ are two diagonals, the components are either zero or satisfy \begin{gather*} 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 } \end{gather*} where $c$ is a constant, $I_{d}$ is an identity matrix. Here $H$ is a matrices of order $I\times J$. The diagonal blocks of $H$ are matrices of order $I\times 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) \begin{itemize} \item The blocks $H_{ii},\quad i=1, \dots , I$ are full symmetric, positive-definite. \item The blocks $H_{ij}=\alpha _{ij}I_{d}$, $i,j=1, \dots , I$, $i\neq j$, $\alpha _{ij}$ is constant. \end{itemize} Thus, to part the inversion of the $H_{ii}$, $i=1, \dots , I$ blocks, the resulted linear system is solved by the blocks Gauss-Seidel iterative method where the blocks of index $ij$, $i\neq j$ are identified to points. In addition, we must declared only $A$ and $B$. \begin{gather*} H= \begin{pmatrix} \quad\begin{pmatrix} . & . & . \\ . & . & . \\ . & . & . \end{pmatrix} & \alpha _{12} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} & \alpha _{13}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \\ \alpha _{21}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} & \quad\begin{pmatrix} . & . & . \\ . & . & . \\ . & . & . \end{pmatrix} & \alpha _{23}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \\ \alpha_{31}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} & \alpha_{32}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} & \quad\begin{pmatrix} . & . & . \\ . & . & . \\ . & . & . \end{pmatrix} \end{pmatrix} \\[5pt] \mbox{fig 4.1: General shape of the $H$ matrix for $I=3$} \end{gather*} \section{Numerical results for $l=0$} Consider three problems $P_{i}$ $i=1,2,3$ in equation \eqref{Pi}, the first problem $P_{1}$ has an exact solution \begin{gather*} u_{1}(x,y)=(1-x^{2}) (1-y^{2}) \\ f_{1}(x,y)=5-3(x^{2}+y^{2})+x^{2}y^{2}\,. \end{gather*} The second problem, $P_{2}$, has an exact solution \begin {gather*} 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}\,. \end{gather*} The last problem, $P_{3}$, has an exact solution \begin{gather*} u_{3}(x,y)=(1-\frac{chy}{ch1})(1-\frac{chx}{ch1}) \\ f_{3}(x,y)=1-\frac{chxchy}{ch^{2}1}\,. \end{gather*} 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 $\varepsilon_i$ denote the error in discrete norm in $L^{2}$ associated to the problem $P_{i}$ $i=1,2,3$ $\varepsilon_i(c)$ denote the error in continuous norm in $L^{2}$, \begin{gather} \varepsilon_i =(\sum_{K\in T_h } (m_{C}(u_i)-m_{C}(u_{h_i}))^{2})^{1/2},\\ \varepsilon_i(c)=\big(\int_\Omega (u_i-u_{h_i}\big)^{2})^{1/2} \end{gather} From these values of the error, we deduct the evaluation of the optimal order of convergence while calculating for every grid the numerical exponent $\alpha$ of \begin {equation} \varepsilon \propto h^\alpha \end {equation} $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{\varepsilon_i}{\varepsilon_{i+1}}=2^\alpha $$ \begin{table}[ht] \caption{Direct physical analytical nodal methods for problems $P_{1}$ and $P_{3}$} \begin{center} \begin{tabular}{|l|l|l|l|l|} \hline mesh & $\quad\varepsilon _{1}$& $\quad\alpha _{1}$& $\quad\varepsilon _3 $& $\quad \alpha _3$ \\ \hline 2x2& 0.01300& & 0.06375& \\ \hline 4x4& 0.00303& 2.07045& 0.01440& 2.1034 \\ \hline 8x8& 0.00074& 2.02865& 0.00347& 2.0377 \\ \hline 16x16& 0.00018& 2.00806& 0.00086& 2.0101 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[ht] \caption{Direct physical polynomial nodal methods for problems $P_{1}$, $P_{2}$ and $P_{3}$} \begin{center} \begin{tabular}{|l|l|l|l|l|l|l|l|l|l|} \hline mesh& $\varepsilon _{1}(c)$& $\alpha _{1}(c)$ &$\varepsilon _{1}$&$\varepsilon _{2}(c)$& $\alpha_{2}(c)$ &$\quad\varepsilon_2$& $\quad\alpha_{2}$& $\quad\varepsilon_{3}$ &$\;\alpha_{3}$ \\ \hline 2 x2& 0.4787& & 0& 1,8839& & 0,0468& & 0.00632& \\ \hline 4 x 4& 0.1265& 1,9435& 0& 0.4971& 1,9466& 0,0165& 1.7768& 0.00214& 1.72 \\ \hline 8 x 8& 0.0321& 1,9848& 0& 0,1272& 1,9767& 0,0045& 1,9366& 0.00057& 1.94 \\ \hline 16x16& 0.0081& 1,9961& 0& 0,0320& 1.9933& 0,0012& 1,9837& 0.00015& 1.95 \\ \hline \end{tabular} \end{center} \end{table} \section {Conclusions} \noindent (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))$. \noindent (2) - We compare in the table2, the discrete $(\varepsilon _{i})$ and continuous $(\varepsilon _{i}(c)) \; i=1,2,3 $ error: \begin{itemize} \item $\varepsilon _{1}=0 $ \item $\varepsilon_{2}(c)\ll\varepsilon _{2}$ \end{itemize} 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. \noindent (3) - We notice that results of the polynomial method are better than the analytical ones (note that these last are $O(h^{l+{3-}\delta _{l{0}} })$ order in norm $L^2$). Indeed we have : \begin{itemize} \item $\varepsilon_{1}=0$ in $table2$ (polynomial method), in $table1$ (analytical methods), $\varepsilon _{1}$ is small but far from being zero. \item $\varepsilon _{3}(\mbox{table 1}) \approx 6\varepsilon_{3}(\mbox{table 2})$. \end{itemize} \begin{thebibliography}{00} \bibitem{c1} Ciarlet P. G.; \emph{The finite element method for elliptic problems}, North-Holland, Amsterdam, 1978. \bibitem{d1} Dorning J. J.; \emph{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). \bibitem{f1} Fedon magnand, Hennart J. P., and Lautard J. J.; \emph{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. \bibitem{f2} Finnemann H., Bennewitz F., and Wagner M.R.; \emph{Interface current techniques for multidimensional reactor calculations}, Atomkernenergie 30, 123- 128 (1977). \bibitem{f3} Frohlich R.; \emph{Summary discussion and state of the art review for coarse-mesh computational methods}, Atomkernenergie 13n 117-126 (1983). \bibitem{h1} Hennart J. P.; \emph{On the numerical analysis of analytical nodal methods}, Numerical methods for partial differential equations, 4, 233-254 (1988). \bibitem{h2} Hennart J. P.; \emph{A general family of Nodal Schemes}, SIAM J.Sci.Stat.Comput. 7.264 (1986). \bibitem{h3} Hennart J.P.; \emph{A Unified view of finite differences, finite elements, and nodal schemes}, National University of Mexico, IIMAS, Mexico City, 01000 Mexico 1986. \bibitem{h4} Hennart J. P., \emph{On the numerical analysis of analytical nodal methods}, Numerical methods for partial differential equations, 4, 233-254(1988). \bibitem{h5} Hennart J.P.; \emph{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. \bibitem{l1} Langenbuch,Maurer W., and Werner W.; \emph{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. \bibitem{l2} Langenbuch S., W. Maurer, and Werner W.; \emph{High order shemes for neutron kinetics calculations}, based on local polynomial approximation, Nucl.Sci.Engng, 64 (1977), pp. 508-516. \bibitem{m1} Meade D., Del Valle E. and Hennart J.P.; \emph{Unconventional finite element methods for neutron group diffusion equations}, IIMAS-UNAM (1986). \bibitem{s1} Shober R. A., Simms R. C., and Henry A. F.; \emph{Tho Nodal Methods for Solving Time Dependent Group Diffusion Equations}, Nucl. Sci. Engng.. 64. 582 (1977). \bibitem{w1} Wagner M. R. and Kobke K.; \emph{Progress in nodal reactor analysis}, Atomkernenergie 13,117-126 (1983). \end{thebibliography} \end{document}