Finite volume method

• Prof. Robert Eymard, Université Paris-Est, France

• Prof. Herbin, Université Aix Marseille 1

• Prof. Thierry Gallouët, Centre de Mathématiques et Informatique, Université de Provence, France

The Finite Volume Method (FVM) is a discretization method for the approximation of a single or a system of partial differential equations expressing the conservation, or balance, of one or more quantities. These partial differential equations (PDEs) are often called conservation laws; they may be of different nature, e.g. elliptic, parabolic or hyperbolic, and they are used as models in a wide number of fields, including physics, biophysics, chemistry, image processing, finance, dynamic reliability. They describe the relations between partial derivatives of unknown fields such as temperature, concentration, pressure, molar fraction, density of electrons or probability density function, with respect to variables within the domain (space, time,...) under consideration.

As in the finite element method, a mesh is constructed, which consists in a partition of the domain where the space variable lives. The elements of the mesh are called control volumes. The integration of the PDE over each control volume results in a balance equation. The set of balance equations is then discretized with respect to a set of discrete unknowns. The main issue is the discretization of the fluxes at the boundaries of each control volume: in order for the FVM to be efficient, the numerical fluxes are generally

• conservative, i.e. the flux entering a control volume from its neighbour must be the opposite of the one entering the neighbour from the control volume,

• consistent, i.e. the numerical flux of a regular function interpolation tends to the continuous flux as the mesh size vanishes.

It is sometimes possible to discretize the fluxes at the boundaries of the control volume by the finite difference method (FDM). In this case, the method has often been referred to as a finite difference method or conservative finite difference method (see Samarskii 2001). The specificity of the FVM with respect to the FDM is that the discretization is performed on the local balance equations, rather than the PDE: the fluxes on boundaries of the control volumes are discretized, rather than the continuous differential operator.

The resulting system of discrete equations depends on a discrete (finite) set of unknowns, and may be either linear or non linear, depending on the original problem itself; this system is then solved exactly or approximately, using for example direct or iterative solvers in the case of linear equations and fixed point or Newton type methods in the case of nonlinear equations.

Let us emphasize that there is no systematic way to derive the discrete system from the continuous one; in fact, the physical or engineering knowledge of the problem is in many cases the key ingredient in order to obtain a scheme which satisfies the expected properties of accuracy, robustness, and computational cost. Let us mention that in some cases, the engineering problem is initially set at a discrete level: the model consists of a set of discrete balance equations and therefore written as a finite volume scheme, and the continuous model of the problem is not explicitly given.

Many of these issues have often been addressed in computational fluid mechanics. A review of the schemes for the discretization of the compressible and incompressible Navier-Stokes and Euler equations may be found in e.g. Wesseling 2001, Feistauer et al. 2003. The FVM has also been very early introduced in the framework of oil reservoir engineering (see Peaceman 1991).

Fundamental principles

Consider the following PDE under conservative form\[\tag{1} {\partial_t A}(x,t) + \nabla\cdot F(x,t) = S(x,t),\]

where the space variable \(x\) belongs to the domain \(\Omega \subset {\mathbb R}^d\) (\(d\) is the space dimension, greater or equal to 1), and the time variable \(t\) belongs to some time interval \([0,T]\ ,\) with \(T>0\ .\) The scalar function \(A\ ,\) defined in \(\Omega\times [0,T]\ ,\) expresses the density of some quantity, and \(\partial_t A\) denotes its time derivative. The function \(F\ ,\) defined in \(\Omega\times [0,T]\) and valued in \(\mathbb R^d\ ,\) expresses the flux of this quantity, and \(\nabla \cdot F\) denotes its space divergence. The function \(S\ ,\) defined in \(\Omega\times [0,T]\) and valued in \(\mathbb R\ ,\) denotes some source term. Some initial condition \(A(x,0) = A_{\rm ini}(x)\) for \(x\in\Omega\) is imposed, where the function \(A_{\rm ini}\) is defined in \(\Omega\) and valued in \(\mathbb R\ ,\) as well as some boundary conditions, which depend on the considered equation. Note that the functions \(A\) and \(F\) are not necessarily regular, so that the derivatives involved in (1) may be weak derivatives.

These functions \(A\ ,\) \(F\ ,\) \(S\) are assumed to be related to a set of unknown fields \((u_j)_{j=1,\ldots,N}\ ,\) where \(u_j\) is an unknown function defined from \(\Omega\times [0,T]\) to \(\mathbb R\ .\) Such relationships may for instance be written as\[\tag{2} \begin{matrix} A(x,t) = {\mathcal A}(u_1(x,t),\ldots,u_N(x,t),x,t), \\ S(x,t) = {\mathcal S}(u_1(x,t),\ldots,u_N(x,t),x,t), \\ F(x,t) = {\mathcal F}(u_1(x,t),\ldots,u_N(x,t),\nabla u_1(x,t),\ldots,\nabla u_N(x,t),x,t), \end{matrix}\]

where the functions \({\mathcal A}\ ,\) \({\mathcal S}\) and \({\mathcal F}\) are given. Classical examples with \(N=1\) (we then denote by \(u\) the unique component \(u_1\)), are:

• the heat equation which corresponds to \(A(x,t) = u(x,t)\ ,\) and \(u\) is the temperature; \(F(x,t) = -\Lambda(x)\nabla u(x,t)\) (where for any \(x \in \mathbb R^d\ ,\) \(\Lambda(x)\) is a linear operator from \(\mathbb R^d\) to \(\mathbb R^d\)), \(S(x,t) = {\mathcal S}(x,t)\ ;\)
• the Poisson equation, which is the steady version of the heat equations \(A(x,t) = 0\ ,\) \(F(x,t) = -\Lambda(x)\nabla u(x)\ ,\) \(S(x,t) = {\mathcal S}(x)\ ;\)
• the linear convection equation \(A(x,t) = u(x,t)\ ,\) \(F(x,t) = u(x,t)V(x,t)\) where \(V\) is a function defined from \(\overline\Omega\times[0,T]\) to \({\mathbb R}^d\ ,\) \(u\) is, say, a concentration, and \(V\) its transport velocity.

The domain \(\Omega\) is partitioned into a mesh, denoted by \({\mathcal M}\ .\) In general, the mesh \(\mathcal M\) is a finite family of nonoverlapping subsets of \(\Omega\ ,\) which satisfy a number of regularity properties; some of these properties are needed for the definition of the scheme, while others are required to establish the precision and the convergence of the method, depending on the type of term that needs to be discretized. The elements of \({\mathcal M}\ ,\) denoted by \(K\ ,\) \(L\ ,\) are called the control volumes; the measure of a control volume \(K\) (its length if \(d=1\ ,\) area if \(d=2\ ,\) volume if \(d=3\)) is denoted by \(|K|\ .\) The boundary \(\partial K\) of each control volume \(K\) is partitioned into the finite set \({\mathcal E}_K\ ,\) called the set of the faces (\(d=3\)) or edges (\(d=2\)) or extremities (\(d=1\)) of \(K\ .\) An element \(\sigma\) of \({\mathcal E}_K\) is either located on the boundary of \(\Omega\ ,\) or belongs to \({\mathcal E}_K\cap {\mathcal E}_L\) where \(K\) and \(L\) are two adjacent control volumes.

Figure 1: Two control volumes with a common face

The \((d-1)\)-dimensional measure of a face \(\sigma\in{\mathcal E}_K\) (its length if \(d=2\ ,\) area if \(d=3\)) is denoted by \(|\sigma|\ .\) A strictly increasing finite sequence \(t^{(0)} = 0 <t^{(1)} <\ldots <t^{(N)} = T\) is defined for the time discretization, and one denotes \(\delta t^{(n)} = t^{(n)} - t^{(n-1)}\ ,\) for \(n=1,\ldots,N\ .\) The classical finite volume approximation of (1) relies on the approximation of the balance equations on the control volumes between time \(t^{(n-1)}\) and \(t^{(n)}\) (in fact, the balance equations usually precede the continuous equation (1) in the derivation of the model; this is the reason why, in several engineering applications, the FVM is applied directly to the discrete balance equations without any knowledge of the continuous equations). These balance equations may be obtained by integrating (1) on \(K\times [t^{(n-1)},t^{(n)}]\) and using the divergence formula\[ \int_K ( A (x,t^{(n)})- A (x,t^{(n-1)})) {\mathrm d}x + \sum_{\sigma\in {\mathcal E}_K}\int_{t^{(n-1)}}^{t^{(n)}} \int_{\sigma} F(x,t) \cdot \boldsymbol n_{K,\sigma}(x) {\mathrm d}s(x) {\mathrm d}t = \int_{t^{(n-1)}}^{t^{(n)}} \int_K S(x,t) {\mathrm d}x {\mathrm d}t, \]

where \({\mathrm d}s(x)\) denotes the integration with respect to the \((d-1)\)-dimensional measure of the boundary and \(\boldsymbol n_{K,\sigma}(x)\) the unit vector normal to \(\sigma\) at point \(x\ ,\) outward to \(K\ .\) The FV scheme is a discretization of the above balance equation, and may be written as\[ |K|(A_K^{(n)} - A_K^{(n-1)}) + \delta t^{(n)} \sum_{\sigma\in {\mathcal E}_K} \vert \sigma \vert \,F_{K,\sigma}^{(n)} = \delta t^{(n)} |K| \, S_K^{(n)}, \]

where\[ A_K^{(n)} \, (\hbox{ resp. } A_K^{(0)} \hbox{) is an approximation of } \frac 1 {|K|}\int_K A(x,t^{(n)}) {\rm d}x\hbox{ (resp. } \frac 1 {|K|}\int_K A_{\rm ini}(x){\rm d}x\hbox{)}, \]

\( F_{K,\sigma}^{(n)} \hbox{ is an approximation of } \phi^{(n)}_{K,\sigma} = \frac 1 {\delta t^{(n)} \vert \sigma \vert }\int_{t^{(n-1)}}^{t^{(n)}}\int_\sigma F(x,t)\cdot \boldsymbol n_{K,\sigma}(x) {\mathrm d}s(x){\mathrm d}t, \)

and

\( S_K^{(n)} \hbox{ is an approximation of } \frac 1 {\delta t^{(n)}|K|} \int_{t^{(n-1)}}^{t^{(n)}}\int_K S(x,t){\rm d}x{\rm d}t. \)

For a face \(\sigma\) common to two control volumes \(K\) and \(L\ ,\) since \(\boldsymbol n_{K,\sigma}(x) +\boldsymbol n_{L,\sigma}(x) =0\) for \(x\in\sigma\ ,\) the normal fluxes in the original problem satisfy the conservation property \(\phi^{(n)}_{K,\sigma} + \phi^{(n)}_{L,\sigma} = 0\ .\) A discrete version of this property of conservation is then required\[ F_{K,\sigma}^{(n)} + F_{L,\sigma}^{(n)}= 0. \] This relation is one of the pillars of the classical FVM; moreover, it is crucial in the mathematical analysis for proving convergence properties.

Discrete expressions in terms of discrete unknowns must then be provided for the terms \(A_K^{(n)}\ ,\) \(F_{K,\sigma}^{(n)}\) and \(S_K^{(n)}\ .\) These discrete unknowns are usually expected to be approximations of the unknown functions \((u_j)_{j=1,\ldots,N}\) at different locations (grid blocks, faces, vertices,...) and times: hence \(u_{j,K}^{(n)}\) (resp. \(u_{j,\sigma}^{(n)}\)) denotes the approximation of function \(u_j\) in control volume \(K\) (resp. at face \(\sigma\)) at time \(t^{(n)}\ .\) For example, considering the case given in (2), it is generally possible to simply use the function \({\mathcal A}\) to define \(A_K^{(n)}\ ,\) using an expression such as

\( A_K^{(n)} = \frac 1 {|K|} \int_K{\mathcal A}(u_{1,K}^{(n)},\ldots,u_{N,K}^{(n)},x,t^{(n)}){\rm d}x. \)

It is generally less simple to derive from (2) an expression for \(F_{K,\sigma}^{(n)}\) with respect to the discrete variables, since there is no systematic method to use the function \({\mathcal F}\) which ensures accuracy and stability on general grids; the discretization of diffusion and convection fluxes are discussed below.

The arguments for the discrete expression of \(S_K^{(n)}\) are generally \(u_{j,K}^{(m)}\ ,\)..., \(m=n\) or \(n-1\ .\)

Discretization of diffusion fluxes

We consider here a diffusive flux \(F(x,t)\) of the form \(F(x,t) = -\Lambda(x,t) \nabla u(x,t)\ ,\) where \(\Lambda(x,t)\) is a linear operator from \({\mathbb R}^d\) to \({\mathbb R}^d\ ,\) and \(u\) is one of the continuous unknowns \((u_j)_{j=1,\ldots,N}\) of the problem. Such an expression of the flux is encountered when using e.g. Fourier's law (heat conduction), Darcy's law (flow in porous media) or Fick's law (chemical diffusion).

In some cases, a simple, conservative and consistent discrete expression for \(F_{K,\sigma}^{(n)}\) may be obtained. This is the case when \(\Lambda(x,t)\) is the identity operator and when the grid satisfies the following "orthogonal property": in each control volume \(K\in{\mathcal M}\ ,\) there exists a particular point \(x_K\ ,\) called the "center" (or centroid) of the control volume, such that, for any adjacent control volume \(L\) to \(K\) with common face \(\sigma\ ,\) the straight line \((x_K,x_L)\) is orthogonal to \(\sigma\ .\)

Figure 2: Two control volumes with a common face, satisfying an orthogonality property

This orthogonality property is satisfied by several types of meshes, such as triangles, rectangles, Voronoï boxes. Then a simple finite difference expression for the approximation of \(\phi_{K,\sigma}^{(n)} =-\frac 1 {\delta t^{(n)} \vert \sigma \vert }\int_{t^{(n-1)}}^{t^{(n)}}\int_\sigma \nabla u(x,t)\cdot n_{K,\sigma}(x) {\rm ds}(x){\rm d}t\) is given by

\(
F_{K,\sigma}^{(n)} = \frac {u_K^{(m)} - u_L^{(m)}} {d(x_K,x_L)},
\)

where \(d(x_K,x_L)\) denotes the Euclidean distance between the points \(x_K\) and \(x_L\ ,\) and provides an implicit (resp. explicit) approximation if \(m=n\) (resp. \(m=n-1\)). This is the famous "two-point flux approximation" scheme, which has been widely used. It is simple, and also much appreciated because it respects a discrete version of the maximum principle; in particular, the approximate solution does not oscillate and stays within its physical bounds. The convergence analysis of the FVM for diffusion equations in this particular case is thoroughly performed in Eymard et al. 2000. The error analysis relies on the conservativity and consistency of the numerical flux.

In the general case of unstructured meshes and general (possibly anisotropic and heterogeneous) linear operators \(\Lambda(x,t)\ ,\) a two point flux formula is not consistent in general, and yields an approximation which might not converge to the appropriate solution.

There are several ways to derive expressions for \(F_{K,\sigma}^{(n)}\) in the general case (intensive research is going on for the analysis of these schemes and their implementation, see the proceedings of the conferences on Finite Volume for Complex Applications):

• The flux may for instance be approximated by a consistent formula using several points and respecting the conservativity property, such as in the MPFA (multiple point flux approximation) scheme.
• Expressions for \(F_{K,\sigma}^{(n)}\) may also be obtained thanks to an identification procedure, using a variational expression and discrete gradients. Such a technique is used in mixed finite volume, hybrid finite volumes (and the SUSHI scheme) and the mimetic method. These three methods can in fact be viewed as elements of one unified family of schemes.
• Another way to express \(F_{K,\sigma}^{(n)}\) in the anisotropic case is to partition the domain \(\Omega\) in more than one mesh (typically two in 2D and three in 3D), and use the differences of the unknowns in the different meshes to reconstruct an expression for \(F_{K,\sigma}^{(n)}\ .\) This is the principle of the discrete duality FVM.
• The expression of \(F_{K,\sigma}^{(n)}\) may also be inspired by the finite element framework, as in the example of Finite Element-Finite Volume methods (see Feistauer et al. 2003).
• \(\ldots\)

Unfortunately, all the linear schemes devised in the general case of unstructured meshes and general linear operators do not seem to satisfy a discrete version of the maximum principle; for severely distorted meshes or for highly anisotropic operators, the risk of oscillations and violation of the physical bounds is high. Ongoing research using in particular nonlinear schemes for linear problems is under way to avoid these drawbacks (see the proceedings of Finite Volume for Complex Applications V).

Approximation of convection terms

Nonlinear scalar convection equations

Consider the partial differential equation\[ {\partial_t u} (x,t) + \nabla\cdot {\mathcal F}(u(x,t),x,t) = 0. \]

This is the conservation equation (1) with \(A(x,t) = u(x,t) \ ,\) \(F(x,t) = {\mathcal F}(u(x,t),x,t)\) and \(S(x,t) =0\ .\) The linear transport equation is a particular case, with \({\mathcal F}(u,x,t) = u V(x,t)\ ,\) and \(V\) is defined from \(\Omega\times[0,T]\) to \({\mathbb R}^d\ .\) The discrete unknowns are the values \(u_K^{(n)}\ ,\) which are expected to be approximations of \(u\) in the control volumes \(K\in {\mathcal M}\) at time \(t^{(n)}\ .\) Since the convection flux is of order 0 (no derivative involved) its approximation does not require any assumption on the mesh, contrary to the diffusion flux. Nevertheless, the fluxes must be carefully approximated in order to ensure the stability (and convergence) of the scheme. Consider the flux \(\phi^{(n)}_{K,\sigma}\) through an interface \(\sigma\) at time \(t^{(n)}\ ,\) seen as a function of the solution \(u\ :\) \(\phi^{(n)}_{K,\sigma}(u) = \frac 1 {\delta t^{(n)}\vert \sigma \vert}\int_{t^{(n-1)}}^{t^{(n)}}\int_\sigma {\mathcal F}(u(x,t),x,t)\cdot\boldsymbol n_{K,\sigma}(x) {\rm ds}(x){\rm d}t\ .\)

A numerical monotone flux  associated to \(\phi^{(n)}_{K,\sigma}(u)\) is a function \(G_{K,\sigma}^{(n)}\) from \({\mathbb R}^2\) to \(\mathbb R\) such that:
*\(G_{K,\sigma}^{(n)}(a,b) = - G_{L,\sigma}^{(n)}(b,a)\ ,\) for an interface \(\sigma\in {\mathcal E}_K\cap {\mathcal E}_L\ ,\)
*\(G_{K,\sigma}^{(n)}\) is  non-decreasing with respect to  its first variable  and non-increasing with respect to
its second variable,  
*\(G_{K,\sigma}^{(n)}\) is locally Lipschitz-continuous with respect to its variables,  
*\(G_{K,\sigma}^{(n)}(a,a) =  \phi_{K,\sigma}^{(n)}(a),\) for any \(a\) respecting the bounds of the solution.  

Then FV schemes are obtained by approximating \( \phi^{(n)}_{K,\sigma}(u)\) by

• \(F_{K,\sigma}^{(n)} = G_{K,\sigma}^{(n)}(u_K^{(m)}, u_L^{(m)})\ ,\) for an interface \(\sigma\in {\mathcal E}_K\cap {\mathcal E}_L\ ,\)
• \(F_{K,\sigma}^{(n)} = G_{K,\sigma}^{(n)}(u_K^{(m)}, u_\sigma^{(n)})\ ,\) for an interface \(\sigma\subset\partial\Omega\ ,\) where \(u_\sigma^{(n)}\) is a suitable approximation of the boundary condition on \(\sigma\times]t^{(n-1)},t^{(n)}[\ ,\)

with \(m = n-1\) (explicit scheme) or \(m = n\) (implicit scheme).

Provided that, for an explicit scheme, the so-called CFL condition be satisfied (this condition prescribes a bound for the time step with respect to the size of the mesh), these schemes satisfy the following properties (see Eymard et al. 2000, Godlewski et al. 1996):

• they are \(L^\infty\)-stable, in the sense that the bounds of the initial and boundary data are preserved,
• they are Total Variation Diminishing (TVD), in the case of structured meshes.

Hyperbolic systems

We now consider a system of \(N\) one-dimensional conservation laws \({\partial_t A_j} (x,t) + {\partial_x F_j}(x,t) = 0\ ,\) where, for \(j=1,\ldots,N,\) \(A_j(x,t) = u_j(x,t)\) and the scalar flux \(F_j(x,t)\) is given by \(F_j(x,t) = {\mathcal F}_j(u_1(x,t), u_2(x,t),\ldots,u_N(x,t),x,t)\ .\) The continuous problem itself is not easy to study from a mathematical point of view, and open questions include existence and uniqueness of a physically relevant solution in the general system case. The continuous fluxes \((\phi_j)^{(n)}_{K,\sigma}\) between times \(t^{(n-1)}\) and \(t^{(n)}\) through an interface \(\sigma\) (which is a real point in the one-dimensional framework) read\[ (\phi_j)^{(n)}_{K,\sigma} = \frac {n_{K,\sigma}} {\delta t^{(n)} }\int_{t^{(n-1)}}^{t^{(n)}} {\mathcal F}_j(u_1(\sigma,t),\ldots,u_N(\sigma,t),\sigma,t){\rm d}t,\ j=1,\ldots,N, \]

with \(n_{K,\sigma} = 1\) or \(n_{K,\sigma} = -1\ ,\) depending on the position of \(\sigma\) with respect to \(K\ .\) A time-explicit FVM approximation reads\[ |K|(u_{j,K}^{(n)} - u_{j,K}^{(n-1)}) + \delta t^{(n)} \sum_{\sigma\in {\mathcal E}_K} (F_j)_{K,\sigma}^{(n)} = 0, \, j = 1, \ldots,N, \] where \((F_j)_{K,\sigma}^{(n)}\) is an approximation of \( (\phi_j)^{(n)}_{K,\sigma}.\) Such an approximation may be performed by the very famous Godunov scheme (presented here in its explicit version, which is the most often used). The idea is, for given \(K,\sigma,n\ ,\) to approximate the flux \((\phi_j)^{(n)}_{K,\sigma}\) by \((F_j)_{K,\sigma}^{(n)} = f_j(w^R_1,\ldots,w^R_N),\) where \(f_j\) is the function defined from \(\mathbb R^N\) to \(\mathbb R\) by\[ f_j(v_1,\ldots,v_N) = \frac {n_{K,\sigma}} {\delta t^{(n)} }\int_{t^{(n-1)}}^{t^{(n)}} {\mathcal F}_j(v_1,\ldots,v_N,\sigma,t) {\rm d}t,\ j=1,\ldots, N, \, \forall (v_1,\ldots,v_N) \in \mathbb R^N, \] and \((w^R_1,\ldots,w^R_N)\) is solution at \(\xi=0\) for all \(\tau>0\) of the so-called one-dimensional Riemann problem, defined by\[\begin{matrix} {\partial_\tau} w_j (\xi,\tau) + \partial_\xi f_j(w_1(\xi,\tau), w_2(\xi,\tau),\ldots,w_N(\xi,\tau)) = 0,& \\ w_j(\xi,0)= u_{j,K}^{(n-1)}, \hbox{ if } \xi < 0,& \\ w_j(\xi,0)= u_{j,L}^{(n-1)}, \hbox{ if } \xi > 0,&\ j=1,\ldots,N. \end{matrix}\]

The Godunov scheme is a very efficient scheme when the above solution of the Riemann problem is easily constructed. When the solution of the Riemann problem becomes too complex or expensive, alternative schemes based on approximate Riemann solvers can be used (see Eymard et al. 2000, Feistauer et al. 2003, Godlewski et al. 1996, Wesseling 2001). See also these references for some extensions to multi-dimensional problems.

References

R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In P.G. Ciarlet and J.L. Lions, editors, Techniques of Scientific Computing, Part III, Handbook of Numerical Analysis, VII, pages 713--1020. North-Holland, Amsterdam, 2000.

R. Eymard and J.-M. Hérard, editors. Finite volumes for complex applications V. ISTE, London, 2008.

M. Feistauer, J. Felcman, and I. Straskraba. Mathematical and computational methods for compressible flow. Numerical Mathematics and Scientific Computation. The Clarendon Press Oxford University Press, Oxford, 2003.

E. Godlewski and P.-A. Raviart. Numerical approximation of hyperbolic systems of conservation laws, volume 118 of Applied Mathematical Sciences. Springer-Verlag, New York, 1996.

D.W. Peaceman. Fundamentals of Numerical Reservoir Simulation. Elsevier Science Inc., New York, NY, USA, 1991.

A.A. Samarskii. The theory of difference schemes, volume 240 of Monographs and Textbooks in Pure and Applied Mathematics. Marcel Dekker Inc., New York, 2001.

P. Wesseling. Principles of computational fluid dynamics, volume 29 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2001.

Sponsored by: Prof. Barbara Zubik-Kowal, Department of Mathematics, Boise State University, Idaho, USA

Reviewed by: Anonymous

Accepted on: 2010-06-23 19:44:18 GMT