Daniel A Cervantes Cabrera

Pedro González-Casanova

L Héctor

Rafael Reséndiz

Vector field approximation using radial basis functions

Keywords: Vector field approximation, radial basis functions, Runge phenomenon

In this article, we investigate the performance of RBF-PDE methods to approximate solenoidal fields. It is well know that global RBF collocations methods present a trade-off principle, which means that smoothness implies high convergence order plus ill-condition. On the other hand, local methods have recently appeared in the literature to solve this problem. In this paper, we perform a numerical investigation between RBF global and local methods in order to investigate the possible advantage of local methods for the approximation of vector fields. More precisely, we compare the local Hermite interpolation technique using inverse multiquadrics versus the non symmetric collocation method of Kansa.

Introduction

Radial basis functions methods have proved to be highly effective for the solution of problems both in the field of approximation theory and the solution of PDEs problems. Among other elements, this is due to their capacity of handling complex geometries in higher dimensions for scattered nodes distributions, as well as the possibility of having spectral convergence. Due to the importance of approximating vector fields in subjects such as electromagnetism, meteorology, scientific visualization among others, we explore the possibility of these methods to solve such problems at great scale.

Despite their success, it is well known that one of the major limitations of radial basis functions (RBF) methods to solve PDEs problems is that they present high instabilities when the number of data is large or when the mesh parameter is large. This problem, known in the literature as the uncertainty principle of Schaback [START_REF] Schaback | Error estimates and condition numbers for radial basis functions interpolation[END_REF][START_REF] Schaback | Multivaraite interpolation and approximation by translates of a basis funciton[END_REF], means that when we have exponential convergence, the corresponding condition number of the Gram matrices increases in an exponential way. Another limitation, which has recently been incorporated to this discussion, is the so called Runge phenomenon. This important limitation is implied by the fact that Gaussian RBFs interpolants in 1D, converge to a polynomial as the shape parameter tends to zero (see [START_REF] Fornberg | Some observations regarding interpolants in the limit of flat radial basis functions[END_REF][START_REF] Fornberg | The Runge Phenomenon and Spatially Variable Shape Parameters in RBF Interpolation[END_REF][START_REF] Fornberg | Stable Computations With Gaussian Radial Basis Functions[END_REF] for more details).

We want to mention the following methods linked to this problem for global RBFs: Fornberg et al. [START_REF] Fornberg | A Stable Algorithm for Flat Radial Basis Functions on a Sphere[END_REF][START_REF] Fornberg | Stable Computations With Gaussian Radial Basis Functions[END_REF] propose an expansion of the RBFs in terms of spherical harmonics, a method called QR-RBF, allow the approximation, in a stable way, of the solution for small values of the shape parameter and large number of nodes. This, however, only applies to periodical domains, such as disks or spheres, where the Runge phenomenon does not appear [START_REF] Fornberg | The Runge Phenomenon and Spatially Variable Shape Parameters in RBF Interpolation[END_REF]. On the other hand, for complex domains on the plane or in space, to our knowledge, this algorithm has only been developed for non trivial node distributions [START_REF] Fornberg | Stable Computations With Gaussian Radial Basis Functions[END_REF]. However, we stress that these techniques are currently under research.

Platte and Driscoll [START_REF] Platte | Eigenvalue stability of radial basis function discretizations for time-dependent problems[END_REF] obtain node distributions making possible to calculate stable solutions for temporal problems in one dimension, or using tensor product techniques in two dimensions. They also remark that these distributions can not be obtained for problems defined on general domains and suggest to use minimum least square techniques of an over determined linear system built by using twice the number of collocation points relative to the PDEs centers.

Elsewhere, Kee, Mao and Liu [START_REF] Kee | A least-square radial point collocation method for adaptive analysis in linear elasticity[END_REF][START_REF] Liu | A stabilized least-squares radial point collocation method (LS-RPCM) for adaptive analysis[END_REF][START_REF] Mao | Least-square-based radial basis collocation method for solving inverse problems of Laplace equation from noisy data[END_REF][START_REF] Zhang | Least-squares collocation meshless method[END_REF] use a similar approach for elasticity problems as well as inverse problems for Laplace equation by means of a local approximating scheme and node adaptive techniques. In this case, the over determination of the local systems is obtained by solving the differential equations and the boundary conditions over a set of additional nodes on the boundary.

Stevens, Power et al. [START_REF] Stevens | A local hermitian RBF meshless numerical method for the solution of multi-zone problems[END_REF] use a local approximating scheme, based on Hermite symmetric interpolation over local domains applied to the partial differential operator of the problem, resulting in a sparse linear system of equations for the values of the solution in each center of the local domain.

In Le Guyader et al. [START_REF] Guyader | Spline approximation of gradient fields: applications to wind velocity fields[END_REF], the authors propose an alternative approach: they study a spline-based approximation of vector fields in the conservative case. In the modeling, they introduce a minimization problem on an Hilbert space: they introduce a regularized least-square problem defined on a space of potentials (real-valued functions) to fit the vector field data-set. For any > 0, they consider the functional J defined as follows:

J :

H m+1 (Ω, R) → R v → ρ(∇v) -u 2 N + |v| 2 m+1,Ω,R , where u = (u 1 , • • • , u N ) T ∈ (R n )
N is the vector field data-set and | • | m+1,Ω,R , the semi-norm on the Sobolev space H m+1 (Ω, R). The authors establish the existence and uniqueness of the solution, they also give convergence result in the introduced Sobolev space using norm equivalence and compactness arguments. Elsewhere, we can also mention several other approaches to study this problem: finite element approximation (see Dzhabrailov et al. [START_REF] Dzhabrailov | The finite element approximation of vector fields in curvilinear coordinates[END_REF]), PDE-based methods (see Amodei and Benbourhim [START_REF] Amodei | A Vector Spline Approximation[END_REF]), spline and RBF approximations (see Awanou and Lai [START_REF] Awanou | Trivariate spline approximations of 3D Navier-Stokes equations[END_REF], Dodu and Rabut [START_REF] Dodu | Vectorial interpolation using radial-basis-like functions[END_REF], Benbourhim and Bouhamidi [START_REF] Benbourhim | Approximation of vector fields by thin plate splines with tension[END_REF]- [START_REF] Benbourhim | Pseudo-polyharmonic vectorial approximation for div-curl and elastic semi-norms[END_REF], Ettl-Lowitsch et al. [START_REF] Ettl-Lowitzsch | Shape reconstruction from gradient data[END_REF][START_REF] Lowitzsch | Approximation and Interpolation employing divergence-free radial basis functions with applications[END_REF][START_REF] Lowitzsch | Error estimate for matrix-valued radial basis functions interpolation[END_REF]).

Although we have emphasized the actual limitations of global RBF approximation for general domains, we stress that they are currently under intensive research and that several techniques, like domain decomposition methods [START_REF] González-Casanova | Node Adaptive Domain Decomposition Method by Radial Basis Functions[END_REF] among others, are under way to overcome these limitations [START_REF] Boyd | Six strategies for defeating the Runge Phenomenon in Gaussian radial basis functions on a finite interval[END_REF]. In this work, we propose local methods to get a solution of the considered vector field approximation problem. Let us note that the resulting linear system of equations is sparse.

It is important to stress that the main advantage of Steven's method is that by solving a set of small local systems of equations, a global sparse matrix is build. This makes it possible, unlike Kansa's method whose condition number can grow in an exponential way, to solve large problems. Moreover, the local systems do not increases their size as the number of nodes increases, which makes it possible to solve them in parallel. In section 4 we further investigate the running time of Stevens' algorithm for our porblem.

We recall that compactly supported RBFs methods, see Wendlad [START_REF] Wendland | Multiscale analysis in Sobolev spaces on bounded domains[END_REF], which have been intensively studied, are a highly important alternative to these problems and in particular to vector field interpolation, but here we focus on RBFs spectral type convergence methods.

Vector field approximation: modelling

In this section, given an initial vector field u 0 , which can be a prescribed function or a given set of values, we describe the mathematical formulation to approximate a vector field u as closed as possible to u 0 , such that u • n = 0 in Γ N and ∇ • u = 0.

We formulate the problem by means of a variational approach, which led to the solution of the corresponding Euler-Lagrange equations [START_REF] Flores | Algorithms for Vector Field Generation in Mass Consistent Models[END_REF][START_REF] Ratto | Mass-consisten models for wind fields over complex terrain: The estate of arte[END_REF]. In what follows we introduce the energy functional and the corresponding Euler Lagrange equations.

Let Ω ⊂ R d (d = 2 or 3) be an open, simply connected and bounded set with boundary ∂Ω = Γ N ∪ Γ D , where Γ N = ∅ is the part of the boundary associated to the surface terrain (topography), Γ D = ∅ is the rest of the boundary and n is an exterior normal vector to Ω ( see Figure 1). So, according to the proposed method by Sasaki [29]: given an initial field u 0 , the exact field of divergence free is the result of minimizing the following functional

L(u, λ) = 1 2 Ω {S(u -u 0 ) • (u -u 0 ) + λ[∇u]}dV
where λ is a Lagrange multiplier and S is a diagonal matrix related to the scales of the components of the velocity, that is

S =   α 2 1 0 0 0 α 2 2 0 0 0 α 2 3  
where the α i are weight parameters called Gaussian precision moduli, related to the scales of the respective components of the velocity field. According to Ratto et al. [START_REF] Ratto | Mass-consisten models for wind fields over complex terrain: The estate of arte[END_REF] the Euler-Lagrange equations of the Lagrangian are

u = u 0 + S -1 ∇λ (1) λu • n = 0 (2)
Finally, in [START_REF] Girault | Finite Element Methods for the Navier-Stokes Equations: Theory and Algorithms[END_REF][START_REF] Flores | Algorithms for Vector Field Generation in Mass Consistent Models[END_REF], the authors show that the Lagrange multipliers satisfy the next elliptic PDE:

-∇ • (S -1 ∇λ) = ∇ • u 0 inΩ, (3)
B 1 λ = g 1 onΓ D (4) B 2 λ = g 2 onΓ N (5)
where B 1 and B 2 are boundary conditions and g 1 and g 2 are given functions. Once λ is calculated, the approximated field u can be recovered from (1).

RBF solution

In this section, we use the classical Kansa's collocation method, with inverse multiquadric kernel φ(r) = (r 2 + c 2 ) -1/2 and we apply it on the reconstruction of a 2 dimensional solenoidal vector field u.

Non-Symmetric Collocation Method

Here, we get a radial basis approximation, so that, we can obtain the adjusted vector field u. To do that, we use Kansa's approach [START_REF] Kansa | Multiquadrics -a scattered data approximation scheme with applications to computational fluid dynamics II[END_REF] to solve numerically the PDE system (3)-( 5) for the Lagrange multiplier λ and then we obtain u from [START_REF] Amodei | A Vector Spline Approximation[END_REF].

In order to do that, let us introduce the following sets of points Ξ = {ξ i } N i=1 and χ = {x i } M i=1 in R d named centers and collocations points respectively (in this case we will take Ξ = χ). Then, a radial approximation λ of the Lagrange multiplier λ, could be defined as follow:

λ(x) = N j=1 β j φ(r j ) = N j=1 β j φ( x -ξ j ) (6)
where

x = x 2 1 + • • • + x 2 d , φ : R d → R
is the radial basis function, and the set {β} N j=1 are the coefficients that we must determinate. To simplify the notations, we define the next differential operator

L := -∇ • (S -1 ∇) (7)
Now, after substituting ( 6) in ( 3)-( 5) we get

N j=1 β j Lφ( x -ξ j )| x=xi = ∇ • u 0 (x i ) x i ∈ Ω, i = 1 . . . , N I (8)
N j=1

β j B 1 φ( x i -ξ j ) = g 1 (x i ) x i ∈ Γ D , i = N I + 1, . . . , N I + N D (9) N j=1 β j B 2 φ( x -ξ j )| x-xi = g 2 (x i ) x i ∈ Γ N , i = N I + N D + 1, . . . , N (10)
where N I ,N D and N N are the number of points in Ω, Γ D and Γ N respectively (N = N I + N D + N N ). Equations ( 8)-( 10), imply the algebraic system

Aβ = b (11)
where the Gram matrix

A =    A L A B1 A B2   
is defined as follow: for j = 1 . . . , N ,

( A L ) i,j = Lφ( x -ξ j )| x=xi f ori = 1, . . . , N I (12)
( A B1 ) i,j = B 1 φ( x i -ξ j ), f ori = N I + 1, . . . , N I + N D ( A B2 ) i,j = B 2 φ( x -ξ )| x=xi f or i = N I + N D + 1, . . . , N (13)
To illustrate this approach, we now give an example which was studied in [START_REF] Flores | Algorithms for Vector Field Generation in Mass Consistent Models[END_REF].

Example 1. We consider the two dimensional solenoidal vector field u(x, z) = (x, -z) defined in Ω = (1, 2) × (0, 1). Assuming that u 0 (x, z) = (x, 0) as an initial horizontal wind field and α 1 = 1, α 3 = 0.001, we consider the following boundary conditions in ( 4) and ( 5)

B 1 λ = λ = g 1 = 0 onΓ D B 2 λ = -S -1 ∇(λ) • n = g 2 = u 0 • n onΓ N
In Table 1, we show the numerical results obtained with inverse multiquadric as radial basis function in [START_REF] Chen | Scientific Computing with Radial Basis Functions[END_REF], with different numbers of nodes and shape parameters.

Note that the condition number grows as the number of node increases or the shape parameter decreases in agreement with the uncertainty principle of Schaback [START_REF] Schaback | Error estimates and condition numbers for radial basis functions interpolation[END_REF][START_REF] Schaback | Multivaraite interpolation and approximation by translates of a basis funciton[END_REF]. In order to reduce the influence of ill-conditioning in the inversion of the Gram matrix, the truncated SVD decomposition can be applied to obtain a stable solution by means of a change of basis ( [START_REF] Schaback | Multivariate interpolation by polynomials and radial basis functions[END_REF]). The results using this decomposition are displayed in Table 1. Note that the L2 error, for the SVD case, is good enough and that the values of the errors are reasonable when N grows or the shape parameter decreases. It is important to note that although the truncated SVD method improves our results, this algorithm can not be applied in general to large scale problems, because the technique that eliminates the singular values, simultaneously remove their respective singular eigenvectors, degrading the basis of the space. Moreover, it should be noted that the computation of this decomposition is of order O(2N 3 ) for square matrices, and due to this fact, some alternatives to these techniques, like domain decomposition methods or techniques based on local approximations have been developed.

In the following tables we use the notation: N number of nodes, κ(A) number of condition the matrix A, ∇ • û divergence of the approximated field û and the L2 error given by: u-û 2 û 2 .

The Runge Phenomena

In the previous section, we noted that once the condition number is controlled, it is possible to improve the relative error when the shape parameter tends to zero, so we can expect further reducing (bounded to the machine epsilon) of this parameter, then it might be possible to obtain lower relative errors. But as shown in Table 2, this trend is preserved only for values greater than or equal to c = 10 -4 (getting a relative error of order O(10 -7 )), since for lower values, the relative error increases. This behavior which has been reported in [29,[START_REF] Fornberg | The Runge Phenomenon and Spatially Variable Shape Parameters in RBF Interpolation[END_REF], is due to the Runge phenomena which is well known in polynomial interpolation for equidistant partitions. The existence of the Runge phenomena for RBF approximations, was first proved by Platte and Driscoll [29], the authors prove that Gaussian kernels interpolants, for equispaced nodes in one dimension, converge to a polynomial when the shape parameter tends to zero or the number of nodes increases.

c N κ(A) ∇ • û L2

Local Hermitian Interpolation Method

Due to the limitations observed for global collocation methods, i.e. the uncertainty principle of Schaback and the Runge phenomena, in this section we propose to apply the local algorithm developed by Stevens, Power et al. [START_REF] Stevens | A local hermitian RBF meshless numerical method for the solution of multi-zone problems[END_REF] to the reconstruction of solenoidal fields. We aim to show that unlike Kansa's technique, the condition number of the global matrix corresponding to the Local Hermitian Interpolation (LHI) method does not increases in an exponential way as the number of nodes and or the shape parameter increases. In what follows we formulate the LHI method for the vector field problem and discuss whether or not it is capable of attaining exponential order of convergence. That is, we want to study if it is competitive with respect to the global collocation method of Kansa. We shall do this by means of different numerical examples.

In this approach, the analytic solution u of a PDE system

Lu(x) = f (x) ∈ Ω (14) Bu(x) = g(x) ∈ ∂Ω (15)
is approximated in a set of nodes Ω sc = {x 1 , . . . , x Nsc } ⊂ Ω (called solution centers), by means of N sc local sub-systems defined for each x i as follows:

u(x i ) = h i x i ∈ Ω sc (16
)

Lu(x) = f (x) x ∈ Ω pdec (17) Bu(x) = g(x) x ∈ ∂Ω f c ( 18
)
where h i are the unknown parameters, Ω pdec = {x 1 , . . . , x N pdec } ⊂ Ω is a set of interior nodes related to the differential operator L and ∂Ω f c = {x 1 , . . . , x N f c } ⊂ ∂Ω the boundary nodes. This procedure generates a set of local linear systems given by

A (k) β (k) = d (k) k = 1, . . . , N sc (19)
which are obtained by substituting the Stevens radial Anzatz given by

û(k) (x) = Nsc j=1 β (k) j φ j (r) + Nsc+N f c j=Nsc+1 β (k) j B ξ φ j (r)+ Nsc+N f c +N pdec j=Nsc+N f c +1 β (k) j L ξ φ j (r) + p m k (20)
in ( 16)-( 17), for the local domains displayed (in circles) in Figure 2. In [START_REF] Kee | A least-square radial point collocation method for adaptive analysis in linear elasticity[END_REF], k is the local system index, N sc the number of solution centers, N pdec the number of PDE centers in the local system and N f c the number of boundary centers in the local system. Moreover,

L ξ φ j (r) := Lφ( x -ξ )| ξ=ξj , B ξ φ j (r) := Bφ( x -ξ )| ξ=ξj , φ(r)
the inverse multiquadric and p m k a polynomial in R d and degree m, which is defined by the null space of ( 14) and ( 15), i.e. for d = 2 and m = 2, the polynomial is given by p(x, y) = a 0 + a 1 x + a 2 y + a 3 xy + a 4 x 2 + a 5 y 2 .

By using [START_REF] Kansa | Multiquadrics -a scattered data approximation scheme with applications to computational fluid dynamics II[END_REF], we obtain the following Gram matrix and the right hand vector Thus, each approximation û(k) (x) has the form

A (k) =       Φ ij L ξ [Φ ij ] B ξ [Φ ij ] P ij L x [Φ ij ] L x L ξ [Φ ij ] L x B ξ [Φ ij ] L x [P ij ] Φ ij L ξ [Φ ij ] B ξ [Φ ij ] P ij B x [Φ ij ] B x L ξ [Φ ij ] B x B ξ [Φ ij ] B x [P ij ]       and d (k) =     h i f i 0 g i   
û(k) (x) = H(x)β (k) (h 1 , . . . , h Ns ) (21)
where

H (k) (x) = φ( x -ξ ) L ξ φ( x -ξ ) φ( x -ξ ) B ξ φ(x -ξ)
Finally, the local systems [START_REF] Kansa | Multiquadrics -a scattered data approximation scheme with applications to computational fluid dynamics II[END_REF] are coupled in a global sparse linear system whose unknowns are the h i values. This is performed by applying the differential operator L to each local interpolator [START_REF] Liu | An introduction to meshfree methods and their programming[END_REF] and evaluating it at the corresponding solution center, i.e.

L[û (k) (x)]| x=xi = f (x i ) (22)
which is the discretization of [START_REF] Flores | Algorithms for Vector Field Generation in Mass Consistent Models[END_REF]. Once the values of the unknowns are obtained, the approximation of u in each solution center (see [START_REF] Guyader | Spline approximation of gradient fields: applications to wind velocity fields[END_REF]) can be determined.

In our problem, i.e the approximation of vector fields, as it has been mentioned, it is necessary to calculate an approximation of the Lagrange multiplier λ by solving the corresponding Euler-Lagrange equations and by using equation [START_REF] Amodei | A Vector Spline Approximation[END_REF].

It should be noted, however, that in this case it is necessary to approximate the gradient of λ, so the Stevens algorithm can not be applied directly. This is because in their case they approximate the solution, whereas in our case we need to estimate the components of the gradient of the Lagrange multiplier.

To achieve this result, we need to change the original Stevens Anzatz (20), by incorporating the the partial derivatives in x and y in the following way

λ(k) x (x) = Nsc j=1 β (k) j ∂φ j (r) ∂x + Nsc+N B 1 j=Nsc+1 β (k) j B ξ 1 φ(r) + Nsc+N B 1 +N B 2 j=Nsc+N B 1 +1 β (k) j B ξ 2 φ(r)+ (23)
+

Nsc+N f c +N pdec j=Nsc+N f c +1 β (k) j L ξ φ j (r) + p m k,x (x) λ(k) y (x) = Nsc j=1 β (k) j ∂φ j (r) ∂y + Nsc+N B 1 j=Nsc+1 β (k) j B ξ 1 φ(r) + Nsc+N B 2 j=Nsc+1 β (k) j B ξ 2 φ(r) (24)
+

Nsc+N f c +N pdec j=Nsc+N f c +1 β (k) j L ξ φ j (r) + p m k,y (x)
Where N B1 are the number of boundary nodes in Γ D and N B1 in Γ N respectively and N f c = N B1 + N B2 . The corresponding local systems are thus given by

         ∂ ∂x λ(x)| x=xi = h i x i ∈ Ω sc Lλ(x) = f (x) x ∈ Ω pdec B 1 λ(x) = g 1 (x) x ∈ Γ D B 2 λ(x) = g 2 (x) x ∈ Γ N (25)          ∂ ∂y λ(x)| x=xi = h i x i ∈ Ω sc Lλ(x) = f (x) x ∈ Ω pdec B 1 λ(x) = g 1 (x) x ∈ Γ D B 2 λ(x) = g 2 (x) x ∈ Γ N ( 26
)
with L defined in [START_REF] Dodu | Vectorial interpolation using radial-basis-like functions[END_REF] and B 1 and B 2 the boundary conditions, so that we can calculate ∇λ = (λ x (x), λ y (x)). Note that p m k,x and p m k,y for our case are defined by p(x, y) This is, in our opinion, a significant modification, which yields an effective methodology to compute numerical solutions of the problem, according to the results that we will show bellow. We conjecture that this type of modifications may be useful to approximate, not only the gradient, but also other type of differential operators like the rotational, when solving partial differential equations.

In order to illustrate the effectiveness of the LHI method, we display in tables (3) and ( 4) the CPU time of the Kansa's collocation method using SVD and the LHI method of Stevens, Power et al. ( [START_REF] Stevens | A local hermitian RBF meshless numerical method for the solution of multi-zone problems[END_REF]) for example number 1.

We can appreciate from these tables that the CPU time (in a Intel core 2, 2.13GHz) of the LHI technique decreases at least one order of magnitude with respect to the SVD-Kansa's method. Note also that for the LHI method, a quad-tree algorith of numerical complexity O(n log(n)) plus the invesion of the local matices corresponding to each solution center, must be taken into acount. We denote A x and A y the global Gram matrices resulting from applying the differential operator L to each Anzatz ( 23) and ( 24), and evaluating them at the solution centers in agreement to equation [START_REF] Liu | A stabilized least-squares radial point collocation method (LS-RPCM) for adaptive analysis[END_REF]. Simultaneously, A Loc x and A Loc y are the matrices which corresponds to the local systems arising form substituting ( 23) and [START_REF] Lowitzsch | Error estimate for matrix-valued radial basis functions interpolation[END_REF] in equation ( 25) and ( 26) respectively.

In what follows, we shall consider first the case where the boundary conditions are of Dirichlet type and then the case where Neumann boundary conditions are taken into account. We aim to investigate whether the LHI scheme is capable to attain exponential orders of convergence, with low condition numbers for the global matrices, and if it can improve the performance of the global collocation method of Kansa. We also analyze the limitations of this method.

For this purpose, in the following numerical experiments, we display, in the first two columns of the following tables, the condition number and the L2 error of the global matrices. Our intention is to examine whether exponential convergence can be reached and a low condition number can be attained simultaneously. On the other hand we display, in columns four and five, the condition numbers of the local systems to analyze whether this values influence the size of the condition number of the global matrix and the L2 error.

Example 2. Approximation with Dirichlet boundary conditions.

In this example we solve the same problem as in Example 1, but in this case we take Dirichlet boundary conditions i.e.

B 1 λ(x) = λ(x) = g 1 (x) B 2 λ(x) = λ(x) = g 2 (x) with g 1 (1, z) = g1(2, z) = α 2 3 2 (1-z 2 ), g 1 (x, 1) = 0 and g 2 (x, 0) = α 2 3
2 . The results are showed in the Tables [START_REF] Boyd | Six strategies for defeating the Runge Phenomenon in Gaussian radial basis functions on a finite interval[END_REF], [START_REF] Chen | Scientific Computing with Radial Basis Functions[END_REF] and Figure 3.

B 1 λ = λ = g 1 = 0 onΓ D B 2 λ = -S -1 ∇(λ) • n = g 2 = u 0 • n onΓ N
The results are showed in the Tables [START_REF] Dodu | Vectorial interpolation using radial-basis-like functions[END_REF], [START_REF] Wendland | Multiscale analysis in Sobolev spaces on bounded domains[END_REF] In tables (5, 6, 7 and 8) only the condition number of A y is showed while A x is omitted, the reason is that both values have a very similar numerical behavior.

In Tables [START_REF] Boyd | Six strategies for defeating the Runge Phenomenon in Gaussian radial basis functions on a finite interval[END_REF] and ( 6) we display, for Dirichlet boundary conditions, the numerical results for two different values of the shape parameter, namely for c = 0.1 and c = 0.01. It can be observed, that the condition number of the global matrix A y , increases with the number of nodes contained in the circle. This happen within an interval and then the condition number slows down within a region around the value 0.6. In this region it can be observed that both the condition number and the error, decreases. In fact, the error is of the order of 10 -6 , which is of the same order obtained with the global method of Kansa. These results are displayed in Figure 3. Also, it can be appreciated that the numerical approximation of the divergence is reduced within this region and its value remains constant having an order of 10 -3 . It can also be observed that it does not exists a direct relation between the maximum (Max local κ(A Loc y )) and minimum (Min local κ(A Loc y )) value of the condition number of the local matrices A Loc y . In tables ( 5) and (6) we display the numerical results when the boundary conditions are of Dirichlet-Neumann type. In this case it can be observed that the behavior is similar to the Dirichlet case, the main difference is that the error is of the order of 10 -5 instead of 10 -6 . These results are display in Figure 4.

We stress that in both cases, Dirichlet and Dirichlet-Neumann boundary conditions, there is a region for the value for the radius, i. e. for the number of local nodes, where the approximation is nearly as good as for the global collocation method, i.e of spectral type.

Conclusions

In this work we have studied numerically the approximation of solenoidal vector fields. Our main objective was to compare global radial basis functions techniques, of spectral type, with local methods for inverse multiquadric kernels. Our numerical results indicate that unlike global methods, local techniques and in particular the LHI method may solve problems of great scale since the global systems of equations has a global sparse matrix with a corresponding low condition number.

The Anzatz of the LHI is flexible and it is possible to incorporate differential operators that depend to the studied problem. However, in our case, this condition was not strong enough and it was necessary to modify the Anzatz itself in order to approximate the gradient of the Lagrange multipliers. This modification of the Anzatz is a significant innovation and confirm that this technique produces good results for problems like the one considered here.

We point that in the case of Dirichlet boundary conditions, local methods achieve an error equivalent to the error obtained by the global approximation collocation techniques for this problem. The numerical experiments for vector fields were done for Cartesian meshes as well as for random data (Halton nodes) verifying that LHI scheme produced similar error in both cases. This LHI local method is highly adequate for parallel computing because the computation of the inverse of each local matrix, corresponding to each solution center, can be performed independently for multiple processors.

On the other hand, we observed that for Dirichlet-Neumann boundary conditions the LHI method produces errors which are less accurate than global collocation schemes. This agree both, with the numerical experiments performed by Stevens,Power et al. ([27]) as well as the the results for least square methods reported by Liu ([20]). In both cases the authors found that Neumann boundary conditions produce instabilities with local methods. This behavior is a subject of current research at the scientific community.

The numerical results obtained in this work, showed that for the LHI method there is a "tread-off" principle similar to the one that exists for the compact support RBF method developed by Wendlad ( [START_REF] Wendland | Multiscale analysis in Sobolev spaces on bounded domains[END_REF]). More precisely, when the support of the local basis functions increases, the error of the approximated solution decreases and the condition number increases. On the other hand, if the support is small, the condition number is low but the error increases. Moreover, as it is displayed in Figures 3,4, there is a region where the sizes of the supports of the local functions, presents small -in fact spectral-errors, and this is true for different values of the shape parameter c. Simultaneously, in this same region it can be observed, see Figure 3, that the condition number of the global matrices reaches minimum values.

Although the maximum and minimum condition number of the local basis functions are high, our results indicate that the increase of the error and the condition number of the global matrix do not depends on these values. The study of this behavior is currently a topic of our current research.

Figure 1 :

 1 Figure 1: The studied domain Ω with boundary conditions (Neumann and Dirichlet).

  h so lu t io n c e n t e r PD E c e n t e r * So l u t i o n c e n t e r B o u n d a r y c e n t e r

Figure 2 :

 2 Figure 2: Centers and local subdomains for the LHI method.

Figure 3 :

 3 Figure 3: Top: Radius v.s. condition number and Radius v.s. L2 Error for c = 0.1 (table 5), Down: Radius v.s. condition number and Radius v.s. L2 Error for c = 0.01 (table6).

Table 1 :

 1 Inverse Gram matrix: SVD.

	c	N	κ(A)	∇ • û	L2 Error
	1	20	2.5e+10 2.3e-06	6.6e-01
	1	60	1.1e+15 6.6e-06	7.0e-01
	1	200 3.0e+23 2.0e-06	1.8e-03
	0.1	20	3.1e+19 2.4e-06	4.9e-01
	0.1	60	2.8e+23 1.7e-06	2.2e-05
	0.1	200 4.1e+22 1.6e-06	3.7e-06
	0.01	20	2.3e+19 1.7e-06	3.0e-04
	0.01	60	2.7e+20 1.8e-06	2.8e-06
	0.01 200 9.2e+20 1.1e-06	1.7e-06

Table 2 :

 2 Runge phenomena for Kansa scheme.

	Error

  = a 1 y + a 2 xy + a 3 xy + a 3

	2 k,x = a 1 x + a 2	x 2 2	+ a 3 xy + a 4	x 2 y 2	+ a 5	x 3 3	+ a 6 xy 2	(27)
	p(x, y) 2 k,y y 2 2	+ a 5	xy 2 3	+ a 6	y 3 3	(28)

Table 3 :

 3 Inverse Gram matrix: SVD.

	c	N	κ(A)	L2 Error	CPU	κ(A)	L2 Error	CPU
	0.1	60	2.8e+23	2.2e-05	0.00181	2.8e+13	2.4e-06	0.000210
	0.1	200	4.1e+22	3.7e-06	0.032	2.1e+14	9.99e-03	0.00556
	0.1	400	1.4e+24	1.4e-06	0.1873	1.26e+14 6.85e-03	0.0171
	0.01	60	2.7e+20	2.8e-06	0.00175	8.1e+13	4.4e-05	0.000204
	0.01 200	9.2e+20	1.7e-06	0.0273	1.3e+14	2.4e-04	0.00164
	0.01 400	4.8e+25	5.7e-06	0.1821	2.41+14	2.5e-03	0.00663

Table 4 :

 4 Inverse Global matrix: LHI.

table 5

 5

	Radius Global κ(Ay)	∇ • û	L2 Error Max local κ(A Loc y ) Min local κ(A Loc y )
	.2	3.8e+06	7.8e+00	1.4e-01	3.4e+27	1.2e+12
	.3	1.3e+12	4.1e-02	4.9e-01	2.9e+28	1.1e+17
	.4	7.4e+15	1.5e-01	4.0e+00	1.7e+29	1.0e+18
	.5	1.1e+17	3.9e-03	4.4e-01	1.1e+30	1.0e+27
	.6	1.9e+12	3.1e-03	2.1e-06	2.2e+30	1.6e+27
	.7	5.9e+12	4.1e-03	4.2e-06	5.6e+30	1.6e+27
	.8	2.7e+11	4.0e-03	3.7e-06	1.0e+30	2.6e+27
	.9	1.3e+17	3.4e-03	3.5e-01	3.5e+30	4.7e+27
	Table 5: LHI approximation with Dirichlet boundary conditions, number of interior nodes =
	100 (50 solution centers and 50 PDE centers) and 20 boundary nodes, with shape parameter
	c =0.1.

), Down: Radius v.s. condition number and Radius v.s. L2 Error for c = 0.01 (table

6

).

Table 6 :

 6 LHI approximation with Dirichlet boundary conditions, number of interior nodes = 100 (50 solution centers and 50 PDE centers) and 20 boundary nodes, with shape parameter c =0.01. In this example we solve the same problem as in Example 1 but this time with the same Dirichlet-Neumann boundary conditions i.e.

	Example 3. Approximation with Dirichlet-Neumann boundary condi-
	tions.

  and Figure 4.

	Radius Global κ(Ay)	∇û	L2 Error Max local κ(A Loc y ) Min local κ(A Loc y )
	.2	4.6e+10	4.3e+00	3.7e+00	1.2e+22	3.8e+18
	.3	1.0e+16	1.9e-01	5.4e+00	3.7e+28	3.3e+19
	.4	3.6e+20	3.1e-01	5.6e-01	1.6e+29	4.9e+19
	.5	3.1e+16	6.2e-01	5.4e-01	4.1e+29	8.0e+19
	.6	6.8e+14	8.2e-03	1.0e-05	4.7e+30	8.1e+26
	.7	3.0e+14	6.3e-03	7.7e-03	2.8e+30	1.3e+27
	.8	1.6e+15	4.4e-03	1.2e-01	5.1e+29	3.7e+27
	.9	2.9e+20	4.4e-01	3.1e-01	7.6e+29	4.4e+27
	Table 7: LHI approximation with Dirichlet-Neumann boundary conditions, number of interior
	nodes = 200 (100 solution centers and 100 PDE centers) and 40 boundary nodes, with shape
	parameter c =0.1.

Table 8 :

 8 LHI approximation with Dirichlet-Neumann boundary conditions; number of interior nodes = 200 (100 solution centers and 100 PDE centers) and 40 boundary nodes, with shape parameter c =0.01.

	Radius Global κ(Ay)	∇û	L2 Error Max local κ(A Loc y ) Min local κ(A Loc y )
	.2	1.1e+14	1.5e+02	6.7e+01	1.7e+28	3.9e+16
	.3	7.1e+21	2.8e-01	5.7e-01	2.0e+27	8.7e+19
	.4	2.9e+21	5.9e-01	5.3e-01	1.1e+27	3.9e+20
	.5	6.8e+17	1.4e-01	6.4e-01	5.5e+27	1.1e+21
	.6	2.4e+14	3.6e-03	4.5e-03	1.8e+28	1.3e+25
	.7	2.3e+15	1.7e-02	5.2e-03	3.3e+28	2.4e+25
	.8	3.8e+12	3.4e-02	1.1e-01	1.0e+28	5.7e+25
	.9	2.7e+22	4.3e-01	4.1e-01	2.7e+28	5.3e+25

Figure 4: Top: Radius v.s. condition number and Radius v.s. L2 Error for c = 0.1 (table 7), Down: Radius v.s. condition number and Radius v.s. L2 Error for c = 0.01 (table

8

).


Ben Faker

Belgacem

Naïma Débit

Henda El Fekih

Souad Khiari

A Finite Element Method for the Boundary Data Recovery in an Oxygen-Balance Dispersion Model

Keywords: Dispersive transport, Data completion, Ill-posedness, Mixed finite elements, Differential algebraic equations

The inverse problem under investigation consists of the boundary data completion in a deoxygenation-reaeration model in stream-waters. The unidimensional transport model we deal with is based on the one introduced by Streeter and Phelps, augmented by Taylor dispersion terms. The missing boundary condition is the load or/and the flux of the biochemical oxygen demand indicator at the outfall point. The counterpart is the availability of two boundary conditions on the dissolved oxygen tracer at the same point. The major consequences of these non-standard boundary conditions is that dispersive transport equations on both oxygen tracers are strongly coupled and the resulting system becomes ill-posed. The main purpose is a finite element space-discretization of the variational problem put under a non-symmetric mixed form. Combining analytical calculations, numerical computations and theoretical justifications, we try to elucidate the characteristics related to the ill-posedness of this data completion dynamical problem and understand its mathematical structure.

Introduction and setting of the inverse problem

The analysis of the oxygen household conditions is central in the water quality modeling. Aquatic life and ecosystem are strongly dependent on the oxygen content of stream-waters. The oxygen balance is governed by the deoxygenation-reaeration problem which is composed of two dispersive transport scalar equations. One involves the biochemical oxygen demand (BOD) concentration and the other is the dissolved oxygen (DO) density. After the mixing length is over practitioners may have good reasons for considering one dimensional models. Our investigation is actually based on oneD problem. Throughout, we denote I = (0, L) the curvilinear representation of the streamwater and T > 0 is a final instant. We use the symbol x for the curvilinear abscissa and t stands for the time variable. The linear dispersive transport system describing the oxygen balance reads therfore as

∂ t b -(db ′ + vb) ′ + rb = f in I × (0, T ) ∂ t c -(dc ′ + vc) ′ + r * c + rb = g in I × (0, T ) db ′ (L, t) = dc ′ (L, t) = 0 in (0, T ) b(x, 0) = b 0 (x), c(x, 0) = c 0 (x) in I. (1)
The longitudinal dispersion coefficient d, the average flow velocity v and the parameters of reaction rates (r, r * ) are all in L ∞ (I). Moreover d, r and r * are positive and are assumed to be piecewise continuous on I. The dispersion parameter is also bounded away from zero. There exists a constant d > 0 such that d(x) > d. We assume also that r(x) > r > 0. The space-time dependent function b denotes the biochemical oxygen demand concentration and c is the dissolved oxygen.

One dimensional models are currently used when one is interested in the downstream dispersion and propagation of an accidental pollution (see [START_REF] Okubo | Diffusion and Ecological Problems: Mathematical Models[END_REF][START_REF] Brown | The Enhanced Stream Water Quality Models QUAL2E and QUAL2E-UNCAS: Documentation and User Manual[END_REF][START_REF] Jolánkai | Basic river water quality models : computer aided learning (CAL) programme on water quality modelling (WQMCAL)[END_REF]). The term rb appearing in the transport equation on c is the depletion of oxygen. This oxygen fraction is necessary to the micro-organisms to oxidize the excessive biodegradable matter in the water. Let us point out that if the dispersion terms are canceled in both transport equations, we obtain the famous Streeter-Phelps model (see [START_REF] Streeter | A study of the pollution and natural purification of the Ohio river[END_REF]).

Zero Neumann boundary conditions at point x = L say that there is no polluting flux downstream.

When every thing is known, suitable boundary conditions are given at the outfall point (x = 0) on both concentrations b and c. They may have the following form b(0, t) = η(t), -dc ′ (0, t) = 0 in (0, T ).

(

) 2
No oxygen supply occurs at the outfall point x = 0 and the polluting load at that same point is known and fixed to η(•). Most often they are of Neumann type, -db ′ (0, t) = γ(t), -dc ′ (0, t) = 0 in (0, T ).

(

) 3
This tells that a polluting flux γ is taking place at (x = 0). That flux may have various causes.

It may for instance be originated from an accidental spill upstream of the river portion under investigation or may be the result of drainage of fertilizers from farm fields. Notice that we are interested in the health state of the river portion downstream. From the modeling point of view, it sounds reasonable, although not necessary from the mathematical view, to consider that the advection velocity is such that v < 0 at least at the vicinity of x = 0. The river flow then inside I at x = 0.

The direct system (1) and (2) (or (3)) is weakly coupled. It is in fact triangular and can be studied using well known tools from the theory of parabolic differential equations (see [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF][START_REF] Lions | Problèmes aux Limites Non Homogènes et Applications[END_REF]). The specific fact here is that, in real-life situations, measurements on the polluting load η or flux γ are too hard to obtain. In the contrary, recording the values of c at the border is easy and may be realized instantaneously. Conditions (2) (or (3)) are therefore replaced by both Dirichlet and Neumann

conditions on c, c(0, t) = α(t) -dc ′ (0, t) = 0 in (0, T ). ( 4
)
The effect of the modification of the boundary conditions on the problem is dramatic. The nature of the problem is entirely altered. With (3), it was well posed, with (4) it becomes ill-posed. It was weakly coupled, it becomes strongly coupled through the boundary conditions. It was a direct problem it becomes an inverse problem. One has for instance to reconstruct the polluting flux γ from a supplementary Dirichlet condition on c at x = 0.

The study realized for this inverse problem in [START_REF] Ben Belgacem | Uniqueness for an ill-posed reaction-dispersion model: Application to organic pollution in stream-waters[END_REF] provides a brief illustration of the instability of the (hypothetic) solution with respect to the data. We propose here to conduct a deeper investigation, in Section 2 to accurately quantify the ill-posedness degree as defined in [START_REF] Wahba | Ill posed problems: Numerical and statistical methods for mildly, moderately and severely ill posed problems with noisy data[END_REF]. Understanding the mathematical structure of this dynamical system becomes troublesome. An identifiability result has been proved in [START_REF] Ben Belgacem | Uniqueness for an ill-posed reaction-dispersion model: Application to organic pollution in stream-waters[END_REF]. Calling for the open map theorem in a suitable framework, one concludes that a by-product of the instability and the uniqueness is that the existence of the solution can not succeed for arbitrary data. The purpose here is to realize a space approximation of the inverse problem ( 1) and ( 4) by a finite element method. The semi-discrete problem turns out to be a differential algebraic equation. We check out how it is affected by the ill-posedness. In particular, the existence issue is tackled frontally owing to the Weierstrass-Kronecker canonical decomposition. This is the subject of Sections 3 and 4. In Section 5, we suggest a full discretization by advancing in time using the Euler scheme. We discuss some issues about the properties of the full discrete problem. To end, we present in Section 6 some numerical experiences to assess the computing devices we use in the simulation of the inverse problem we deal with.

Ill-posedness degree and Uniqueness

A wide literature exists dealing with the direct problem (1) and (3) in both theoretical and computational registers. To our opinion, putting aside the non-linearities not provided here, the most important feature that has been studied has to do with the influence of the Péclet number on the numerical simulations (see, eg, [START_REF] Rubin | Environmental Fluid Mechanic[END_REF]). The specific issue we are confronted to in this work is related to the redundancy of the boundary condition on the concentration c at x = 0 to compensate the lack of conditions on the density b. With this respect, the dispersion plays a preponderant role in the complications arisen in the mathematical handling of the inverse problem. We will therefore exclusively focus in this effect. As shown in [START_REF] Ben Belgacem | Uniqueness for an ill-posed reaction-dispersion model: Application to organic pollution in stream-waters[END_REF] and this will be confirmed along this work, the study is far from being easy to achieve and we thus prefer to switch off any further causes of technical difficulties so as to stress on the central issue.

The ill-posedness of the problem has been briefly illustrated in [START_REF] Ben Belgacem | Uniqueness for an ill-posed reaction-dispersion model: Application to organic pollution in stream-waters[END_REF]. We here intend to accurately evaluate the ill-posedness degree, as defined in [START_REF] Wahba | Ill posed problems: Numerical and statistical methods for mildly, moderately and severely ill posed problems with noisy data[END_REF], at the cost of a more complete investigation.

We accomplish this study when the physical parameters are all constants. According to the remark made above, the main point is the incidence of the dispersion. To simplify, one may get rid of the advection. Indeed, as well known in the theory of Sturm-Liouville operators, there is a suitable transformation that allows to hide the advective term. In addition, the most important feature is concerned with the sensitivity of the reconstructed flux γ(•) in (3) with respect to the measures α(•) in (4). To alleviate the exposition, we fix (f, g) and (b 0 , c 0 ) to zero. This is not a restriction to generality. Henceforth, we assume that v = 0 in both transport scalar equations on b and c. Since we intend to use Fourier series we put L = π to alleviate the calculations. We therefore consider the pure dispersion-reaction problem

∂ t b γ -db ′′ γ + rb γ = 0 in I × (0, T ) ∂ t c γ -dc ′′ γ + r * c γ + rb γ = 0 in I × (0, T ) db ′ γ (π, t) = dc ′ γ (π, t) = 0 in (0, T ) -db ′ γ (0, t) = γ(t), -dc ′ γ (0, t) = 0 in (0, T ) b γ (x, 0) = 0, c γ (x, 0) = 0 in I. (5)
The polluting agent is the flux γ(•). The right one to recover is the one that yields the observations α(•). It is thus solution of the following reduced equation

(Sγ)(t) := c γ (0, t) = α(t) in (0, T ). ( 6
)
If γ(•) is given then the direct system (5) is triangular and has a unique solution (b γ , c γ ). The linear operator S is well defined and maps continuously L 2 (0, T ) into L 2 (0, T ). Using Fourier series in (5), S can be expressed as a convolution operator. Indeed, assembling the dispersion-reaction equation, the boundary and the initial conditions on b γ , then solving the well-posed resulting problem, one obtain that

b γ (x, t) = 1 π t 0 γ(s)e -r(t-s) ds + 2 π k≥1 t 0 γ(s)e -λ k (t-s) ds cos(kx).
The sequence (λ k ) k is given by λ k = dk 2 + r. In the other hand side, gathering the equation on c

and due to the Neumann boundary and initial conditions it may be stated that

c γ (x, t) = r π t 0 γ(s)e -r * (t-s) 1 -e -(r-r * )(t-s) r -r * ds + 2r π k∈N t 0 γ(s)e -µ k (t-s) 1 -e -(r-r * )(t-s) r -r * ds cos(kx).
Here, we have µ k = dk 2 + r * . Now, considering c γ (0, t), this formula yields the convolution form of

S, that is (Sγ)(t) = t 0 K(t -s)γ(s) ds ∀t ∈ (0, T ).
The convolution kernel K(•) is given by

K(s) = r π(r -r * ) (1 -e -(r-r * )s ) e -r * s + 2 k≥1 e -µ k s ∀s ∈ (0, T ).
It is easily seen that K is indefinitely smooth away from the origin. The qualitative properties of S are thus tightly dependent on the behavior of this kernel at the vicinity of zero. By the way, a useful information one can have right away is the ill posedness of problem (6) because of K(0) = 0.

In particular, the operator S is compact. We aim to know more about the ill-posedness degree of ( 6). This may be ensued from the asymptotics of the singular values of the operator S (see [START_REF] Gripenberg | Volterra Integral and Functional Equations[END_REF]). We are hence required to find out the positive real number a > 0 so that K(s) = O(s a ). This pertains to the limit toward zero of the infinite sum defining K(•). Comparison with integrals allows to state that

K(s) ∼ - r 2 √ πd √ s,
at the vicinity of zero. Applying Theorem 3.1 proved in [START_REF] Dostanic | Asymptotic behavior of singular values of certain integral operators[END_REF], concludes to the fact that the sequence of the singular values of the operator S decays towards zero like Ck -3/2 . This tells that the illposedness degree is 3/2 and ascertains of the mild ill-posedness of the reduced problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. As a result, the resolvent of the inverse problem ( 1) and ( 4) cannot be continuous and that problem suffers thus from instability.

Remark 2.1 One may wonder why the data completion problem we consider is mildly ill-posedness.

A fundamental part of the answer is the fact that the reconstruction of the missing pollution flux γ(•) is pursued at the same location where the observations on the oxygen concentration c are made.

It may occur of course that collecting measures on c is not possible at that same point x = 0 while they are available at a different location. We refer to Figure 1 for an illustration

Let us then have a look at case where the abundant measures on c are obtained at the other extreme point x = π. The same Fourier computations can be achieved that yields the convolution kernel The other important issue for the inverse problem has to do with the uniqueness or in other words it is related to the injectivity of S. The result has been established in [START_REF] Ben Belgacem | Uniqueness for an ill-posed reaction-dispersion model: Application to organic pollution in stream-waters[END_REF]Theorem 2.7] in the general context of space-varying parameters. The tools used there come from functional analysis of the saddle point problems and are long to summarize. Nevertheless, it may be briefly discussed in some particular cases. We intend to provide an illustration of that uniqueness result using a simple argument. We then consider the same problem as above. The parameters are all constant and we make the additional assumption that T = ∞. We are consequently enabled to use Laplace's transform as a bounded operator in L 2 (0, ∞). The reduced problem, where α(•) is put to zero, can then be written as

H(t) = r π(r -r * ) (1 -e -(r-r * )s ) e -r * s + 2 k≥1 (-1) k e -µ k s ∀s ∈ (0, T ). db'(0) = ? dc'(0) = 0, c(0) = α db'(L) = 0 dc'(L) = 0 dc'(L) = 0, c(L) = α db'(L) = 0 db'(0) = ? dc'(0) = 0
Sγ(p) = K(p)γ(p) = 0 ∀p ∈ (0, ∞).
The symbol γ(•) is for the Laplace transform of the function γ(•). Easy computation yields that

K(p) = 2r π k∈N 1 (p + µ k )(p + λ k ) ∀p ∈ (0, ∞).
It can be checked out that K(p) = 0, ∀p ∈ (0, ∞), we derive that γ(•) vanishes identically. On account of the injectivity of the Laplace transform we conclude to γ(•) = 0. The uniqueness is assessed by linearity.

Before closing this section let us say some words about the surjectivity of S which has a narrow link to the existence for problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. Let us first indicate that the adjoint operator S * is an integral operator. Following the same arguments as for S, it can be proved that S * is injective. The range of S is thus dense in L 2 (0, ∞) and cannot be closed. Otherwise, applying the open map theorem yields that the inverse of S would be bounded, but we known that it is not. Hence, surjectivity fails for S and the existence of a solution for (6) may not be ensured. The same conclusion may be extended to problem (1) and ( 4).

Finite element approximation

A space finite element approximation of the problem ( 1) and ( 4) is the purpose here. As we will see later on this sheds more clarity about the existence issue which was addressed through an indirect argument for the continuous problem. We still keep the assumption, v = 0. Moreover, handling non-homogeneous Dirichlet conditions does not arise any particular obstacle. We assume thus that α(•) is also zero. Adding the advection in the inverse model or/and taking into account non-vanishing boundary condition do not introduce more complexity from the theoretical point of view, even though they are undeniably important in the physical and computational grounds.

Let I i = (x i , x i+1 ) 1≤i≤n-1 be a subdivision of the interval I with x 1 = 0 and x n = π. The grid-size h is equal to max 1≤i≤n-1 |I i |. We refer to [START_REF] Brenner | Mathematical Theory of Finite Element Methods[END_REF][START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] for the basics of the finite element method. Next, we define the discrete spaces V h and H h

V h = ϕ h ∈ C (I); ∀i, (ϕ h ) |I i ∈ P 1 , H h = ψ h ∈ C (I); ∀i, (ψ h ) |I i ∈ P 1 , ψ h (0) = 0 ,
where P 1 stands for the space of affine functions. The space approximation relies on a mixed variational form of the problem. The framework selected here reproduces the one used in [START_REF] Ben Belgacem | Uniqueness for an ill-posed reaction-dispersion model: Application to organic pollution in stream-waters[END_REF], adapted to the discrete level. That framework is the only one liable to provide the identifiability result for the exact model ( 1) and ( 4). It sounds natural to select it for our finite element approximation.

The semi-discrete problem is hence constructed following Ritz-Galerkin's procedure. It reads as:

find (b h , c h ) in V h × H h fulfilling ∀ψ h ∈ H h , (∂ t b h (t), ψ h ) L 2 (I) + k (ψ h , b h (t)) = (f (t), ψ h ) L 2 (I) , ∀ϕ h ∈ V h , (∂ t c h (t), ϕ h ) L 2 (I) + k * (c h (t), ϕ h ) + a r (b h (t), ϕ h ) = (g(t), ϕ h ) L 2 (I) , b h (•, 0) = b 0,h , c h (•, 0) = c 0,h . (7)
Using (•, •) L 2 (I) marks the inner product in the space L 2 (I). Different bilinear forms k(•, •) and

k * (•, •) and a r (•, •) are expressed by, ∀(ψ h , ϕ h ) ∈ H h × V h , k(ψ h , ϕ h ) = I (db ′ h ψ ′ h + rb h ψ h ) dx, ∀(ψ h , ϕ h ) ∈ H h × V h , k * (ψ h , ϕ h ) = I (dc ′ h ψ ′ h + r * c h ψ h ) dx, ∀(ψ h , ϕ h ) ∈ V h × V h , a r (χ h , ϕ h ) = I rχ h ϕ h dx.
To where c h vanishes. Problem ( 7) can be expressed under a matrix form,

M 0 0 M T ∂ t b(t) c(t) + K 0 A r K T * b(t) c(t) = f (t) g(t) , in (0,
M ∂ t Y (t) + KY (t) = G(t), in (0, T ) Y (0) = Y 0 . (8)
The data G contains the contribution of the input (f , g) and Y 0 contains the initial states (b 0 , c 0 ).

Both matrices M and K are square. Problem ( 8) is an implicit differential equation; it is not an ordinary differential equation because M is not invertible. The incidence of this singularity is that existence for (8) may fail for arbitrary data.

Understanding the particular structure of equation ( 8) is tightly connected to the matrix M . Let us therefore check out the important statement

The kernel N (M ) is a subspace of dimension one. The matrix M is therefore singular.

Proof : Let Y be given in N (M ). We consider (b h , c h ) ∈ V h ×H h as the finite element functions whose degrees of freedom form the coefficients of Y . According to the construction of N (M ), the following double orthogonality holds

I b h ψ h dx = 0, ∀ψ h ∈ H h ; I c h ϕ h dx = 0, ∀ϕ h ∈ V h . (9)
We deduce directly that c h = 0. The function b h considered alone has n degrees of freedom. It is subjected to (n -1) constraints that are linearly independent. As a consequence, the set of these functions determines a vector space of dimension one. The value b h (x 1 ) can be considered as the only degree of freedom. The dimension of the kernel N (M ) is thus one. This results in the non-invertibility of M . As a consequence of the non-invertibility of M , problem [START_REF] Griepentrog | Differential-Algebraic Equations and Their Numerical Treatment[END_REF], apparently a differential equation, is actually hiding an algebraic and a differential (sub)equations or components. It is a differential algebraic equation. An exposition of these of equations may be found in some manuscripts (see [START_REF] Kunkel | Differential-Algebraic Equations. Analysis and Numerical Solution[END_REF][START_REF] Lamour | Differential-Algebraic Equations: A Projector Based Analysis[END_REF][START_REF] Riaza | Differential-Algebraic Systems: Analytical Aspects and Circuit Applications[END_REF][START_REF] Brenan | Numerical Solution of Initial-Value Probin Differential Algebraic Equations[END_REF]) where specific tools are elaborated to cope with them. Before proceeding with this methodology, it is useful to know more about the stiffness matrix K. It turns out to be invertible under a mild assumption on the grid-size. This condition on the mesh is liable to ensure the discrete positivity principle. Observe first that the kernel of the matrix K can be identified to the subspace

N (K) = ϕ h ∈ V h ; k(ψ h , ϕ h ) = 0, ∀ψ h ∈ Q h .
The subspace N (K * ) in linked to k * (•, •) in the same way. Under the condition on the grid-size

h 2 max (sup x∈I r(x), sup x∈I r * (x)) ≤ 6 inf x∈I d(x), (10)
the following statement holds true (see [START_REF] Vejchodský | Necessary and sufficient condition for the validity of the discrete maximum principle[END_REF])

ϕ h ∈ N (K) ∩ N (K * ), ϕ h (0) ≥ 0 implies that ϕ h (x) ≥ 0, ∀x ∈ I.
We have thence the invertibility result, Assume condition (10) on the grid-size is fulfilled. Then, the stiffness matrix K is invertible.

Proof : Let Y be in N (K) and is the representation of (b h , c h ). It is direct that b h ∈ N (K).

Moreover, it can be easily checked out that

a r (b h , ϕ h ) + k * (ϕ h , c h ) = a r (b h , ϕ h ) = 0, ∀ϕ h ∈ N (K * ).
Choosing ϕ h such that ϕ h (x 1 ) = b h (x 1 ), then, by the positivity principle both b h and ϕ h are simultaneously positive or negative. The fact that a r (b h , ϕ h ) = 0 results in b h = ϕ h = 0. Deducing that c h = 0 becomes straightforward from the injectivity of (K * ) T . We conclude thus to Y = 0, hence the proof.

Remark 4.2 Arguing analogously, one can also obtain that (λM + K) is invertible for any λ ≥ 0.

To reveal useful properties of problem (8), we call for the algebraic transformation introduced by Kronecker and Weierstrass (see [START_REF] Kronecker | Gesammelte werke, Reduktion der scharen bilinearer formen[END_REF][START_REF] Weierstrass | Gesammelte werke, Zur theorie der bilinearen und quadratischen formen[END_REF]). The chief idea is to put the problem under a canonical 'Weierstrass-Kronecker' form that enables to uncouple the differential equation and the algebraic equation and to find out their respective dimension. Then, each of the two subproblems can be studied in its own context.

Realizing this requires the regularity in the sense of 'Weierstrass-Kronecker'. The matrix pencil formed by the pair of matrices (M , K) is regular because the characteristic polynomial p(λ) = det(λM + K) does not vanish identically. Indeed, p(0) = 0, as K is invertible. Due to the singularity of M , the degree of the polynomial p is obviously lower than n. The following result holds (see [START_REF] Lamour | Differential-Algebraic Equations: A Projector Based Analysis[END_REF]Proposition 1.3] and [12, Theorem 2

.7])

There exist two linear bijective transformations Z = LY = (u, v) T and R = HG = (q, r) T which yield the 'Weierstrass-Kronecker' canonical form given by

I n-m 0 0 N ∂ t u(t) v(t) + W 0 0 I m u(t) v(t) = q(t) r(t) , in (0, T ). ( 11
)
The vectors u, q and v, r are of dimension (n-m) and m ≥ 1 respectively. The vector function u is the differential variable while v the algebraic variable. The sub-matrix W is square of dimension (nm) and the bloc N is square of dimension m with N nilpotent. Recall that the nilpotency index µ of N is the smallest integer such that N µ = 0. The symbol I m is for the unity matrix of dimension m. Equation (11) can also be displayed under a condensed form

N ∂ t Z(t) + W Z(t) = R(t), in (0, T ),
with the following relations on the matrices, HN = M L and HW = KL.

Z 0 = LY 0 = (u 0 , v 0 ) T . Then, we derive in one hand side the ordinary differential problem

∂ t u(t) + W u(t) = q(t) in (0, T ), u(0) = u 0 , (12)
and we obtain in the other hand side the algebraic equation

N ∂ t v(t) + v(t) = r(t) in (0, T ), v(0) = v 0 . ( 13
)
The ordinary differential problem ( 12) is well posed and has unique solution provided by

u(t) = e -W t u 0 + t 0 e -W (t-s) q(s) ds, in (0, T ).
Letting aside the initial condition for a while, the solution of the algebraic equation ( 13) may be computed by induction. There comes out the following expression

v(t) = 1≤j≤µ-1 (-1) j N j r (j) (t), in (0, T ).
This formula indicates that the full algebraic equation may not have a solution unless the initial state v 0 and the source data r are consistent. Existence is guaranteed by the fulfillment of the constraint

v 0 = 1≤j≤µ-1
(-1) j N j r (j) (0).

In the pratice, computing the Weirstrauss-Kronecker decomposition of equation ( 8) may be expensive and may not be pursued. Actually, solving (8) is highly sensitive to the nilpotency index µ of the matrix N . It is called the Kronecker index of the matrix pencil (M , K). This Kronecker index µ and the dimension m of N are both related to the spectral decomposition of the matrix K -1 M . We already know that 0 is eigenvalue. If ℓ is its algebraic multiplicity then dim N = m coincides with the dimension of the characteristic subspace N ((K -1 M ) ℓ ). The Kronecker index µ is lower or equal to m. Knowing the value of the index µ for the linear differential algebraic equations is so important. It is a strong indication of the complexity of integrating them. Indexes higher or equal to two are causes of instability and generate tough theoretical and practical difficulties. If that index is one, the principal effect is that the initial value cannot be fixed freely in the algebraic (sub)equation. Notice that the result [START_REF] Kronecker | Gesammelte werke, Reduktion der scharen bilinearer formen[END_REF] does not provide the Kronecker index nor the dimension of N . This is rather supplied in

The Kronecker index of the pencil (M , K) is one, µ = 1. As a consequence, we have N = 0.

Moreover, the dimension of N is equal to one, m = 1.

Proof : Following [8, Appendix A, Theorem 13], the index µ is one if we show that

N (M ) ∩ K -1 (R(M )) = {0}.
The symbol R(M ) is for the range of M . Let Y be in the intersection. On one hand side we have that M Y = 0 and on the other hand side there exists Z such that KY = M Z. Now, denote (b h , c h ) and (β h , χ h ) the finite element functions in V h × H h associated to Y and Z, respectively.

As established above, that Y lies in N (M ) yields the orthogonality [START_REF] Gripenberg | Volterra Integral and Functional Equations[END_REF]. As seen there, this gives in particular that c h = 0. Next, the formula KY = M Z may be converted into a functional form

∀ϕ h ∈ V h , k * (c h , ϕ h ) + a r (b h , ϕ h ) = (χ h , ϕ h ).
Owing to the fact that c h vanishes identically and choosing

ϕ h = b h ends to a r (b h , b h ) = (χ h , b h ).
Now, since χ h ∈ H h , calling for (9), we derive that a r (b h , b h ) = 0. This provides that b h = 0 and concludes to Y = 0. An immediate consequencer is that N = 0. Finally, the dimension m of the nilpotent matrix N in (8) coincides with the dimension of N (M ) which is one. The proof is complete.

The analysis concludes to the fact that the differential algebraic problem (8) obtained by the space discretization of the inverse problem ( 1) and (4) may not be solved for arbitrary data Y 0 and G. The solvability requires a particular compatibility between these data. The initial state Y 0 and the source data G(0) should be consistent. However, when it comes to the identifiability, the following uniqueness result holds

The semi-discrete problem (7) has at most one solution. Existence is secured only for consistent data (b 0,h , c 0,h ) and (f, g).

5 Euler Time Scheme.

Numerical solvers for differential algebraic equation, commonly called descriptor space-state systems, is a sensitive issue especially for high (Kroenecker) indexed equations. Time schemes carelessly used for these equations may converge poorly or not at all. A specialized literature can be found in (see [START_REF] Lamour | Differential-Algebraic Equations: A Projector Based Analysis[END_REF][START_REF] Kunkel | Differential-Algebraic Equations. Analysis and Numerical Solution[END_REF][START_REF] Brenan | Numerical Solution of Initial-Value Probin Differential Algebraic Equations[END_REF]). The situation is not that complicated in our case. The Kroenecker index of problem ( 7) is low as it equals one and the algebraic part of problem ( 8) is of dimension one. So to say, we profoundly doubt that it conceals any worthy information. Its main influence seems to maintain the ill-posedness for the semi-discrete problem. Thus, all the care should be paid for the approximation of the ordinary differential (part of the equation. The concern for practitioners is possibly the consistency of the initial value and the source data which is not an easy matter because of the implicit form of differential algebraic equation. Because of the one-index of our problem it can be handled numerically using modified standard integration methods.

When the full approximation of variational problem [START_REF] Dostanic | Asymptotic behavior of singular values of certain integral operators[END_REF] is pursued, a time marching scheme has to be combined with the finite element method. Let then τ be time step, T = p * τ and p is the current time level. We consider ((b h ) p , (c h ) p ), the space/time approximation of (b(pτ, •), c(pτ, •)). BDF methods , for Backward Differential Formula, are well suited to differential algebraic equations. We opt for the one step scheme which coincides with the better backward Euler first order time scheme.

The variational problem to solve is thus expressed as follows: find the sequence

((b h ) p , (c h ) p ) 0≤p≤p *
that satisfy the induction

∀ψ h ∈ H h , (b h ) p+1 -(b h ) p τ , ψ h + k(ψ h , (b h ) p+1 ) = (f p+1 , ψ h ), ∀ϕ h ∈ V h , (c h ) p+1 -(c h ) p τ , ϕ h + k * ((c h ) p+1 , ϕ h ) + a r ((b h ) p+1 , ϕ h ) = (g p+1 , ϕ h ), (b h ) 0 = b 0,h , (c h ) 0 = c 0,h . (14)
The well definiteness of the sequence ((b h ) p , (c h ) p ) 0≤p≤p * remains questionable. We need to know if the induction defines without unambiguity ((b h ) p , (c h ) p ) 0≤p≤p * . To clarify this point, we need to look how the time scheme is echoed on the differential algebraic equation ( 8),

(λM + K)Y p+1 = G p+1 + λM Y p , Y 0 = Y 0 . (15)
The parameter λ is the inverse of the time step, that is λ = τ -1 > 0. Following Remark 4.2, the matrix (λM + K) is invertible and the sequence ((b h ) p , (c h ) p ) 0≤p≤p * is therefore uniquely constructed. The resulting algebraic system has to be solved repeatedly and we choose to use a direct algorithm. Moreover, if one is concerned with the accuracy of the scheme, we recall that given that the Kroenecker index one, then the Euler time scheme is consistent with order one provided that the data are consistent.

Numerical Investigation

We investigate some features of the semi-discrete and full-discrete versions of the oxygen balance inverse problem (1) and (4). We first address the issue of determining the Kronecker index. Then, we assess through some judicious examples the reliability of the Euler scheme/mixed finite element method discussed here. Computations are realized in Matlab. Scripts for one dimensional finite elements are specifically implemented for our purpose.

Eigenvalues of the pencil (M , K)

As proved above the Kronecker index and the dimension of the matrix N in formula ( 11) are equal to one. We assess these results numerically. Recall that the spectral decomposition of the matrix (K -1 M ) provides necessary knowledge to learn about these indices.

We The generalized eigenvalue of M with respect to K are therefore evaluated by means of the Matlab procedure eig. These eigenvalues coincides with those of (K -1 M ). The first thing we learn from the computations is that (K -1 M ) is diagonalizable in C. Some eigenvalues are real, others are complex. The set of eigenvalues is symmetric with respect to the real axis which is expected given that (K -1 M ) is a real matrix. These eigenvalues are accumulated around zero.

This can be explained by the fact that K is the finite element representation of a second order operator while M is mass matrix that describes a zero-degree operator. The matrix (K -1 M ) may be viewed as the discretization of a compact operator. Now, we come to the key point of the Kronecker index. Results by Matlab confirms that 0 is an eigenvalues and it is simple. To be convinced we display in We refer for instance to [START_REF] Tyagi | Mathematical modelling of stream do-bod accounting for settleable bod and periodically varying bod source[END_REF] for very close values to these.

We start the computational runs with a noise-free DO profile α(•), the one represented in Fig- Noise is expected to ruin the numerical solution. That is what exactly happen; facts can be checked out in the plots of Figure 6. If the magnitude of the white noise grows higher then the computed solution seems to behave in completely erratic way and is far away from the exact 'real' polluting load as a time dependent function. A possible remedy to restore a reasonable load is to apply a-posteriori Gaussian filter to dampen the harmful effect of noise on the load function (see [START_REF] Blinchikoff | Filtering in the time and frequency domains[END_REF]). The a-priori unrecognizable and apparently inaccurate computed function η D (•) conceals an effective signal. Indeed, that hidden signal can then be reconstructed after operating the filter on the approximated profile. Figure 7 depicts the filtered signals. They reproduce the main features of the exact profile in particular the number of peaks and dips. These are interesting information to hydrological engineers, especially in facing accidental spills in the stream waters.

The last issue assessed here is the impact on computations of noise presence on the initial condition c 0 (•). An example is run in the case where c 0 (•) is contaminated by white Gaussian noise of variance σ = 0.5. Both noisy initial data and the computed BOD load are portrayed in Figure 8.

Letting aside a visible effect at the initial instant, the profile of η(•) seems to be nicely approximated. The observation may be put in relation with the inconsistency of the noisy initial data wit the boundary condtions at the outfall point. Using a Gaussian filter is liable to bring some correction to the signal and repair that deviation a the vicinity of t = 0.

Conclusion

The relevance of the mixed finite element method for the boundary data completion in the dispersive BOD-DO model is the main purpose we pursue in both mathematical and computational grounds.

The semi-discrete problem we obtain turns out to be a differential algebraic equation (DAE). Tools specifically developed for implicit DAE enable a careful analysis of the semi-discrete equations.

As result, we highlight the reason why existence fails for some (inconsistent) data. Prospective numerical experiments presented show the liability of such a finite element approximation when combined with a well adapted time marching scheme in solving the inverse problem we cope with.

Further research in many directions is under consideration. Extension to higher dimension, using less naive filters, testing and sorting different regularization strategies are all in expected in the forthcoming work. The analysis of the inverse problem where measurements on the DO density are not accessible at the outfall point is also under consideration.

Figure 1 :

 1 Figure 1: Data can be collected either at of the outfall point or downstream.

Figure 2 :

 2 Figure 2: K(•) (left) and H(•) (right) so as their first derivatives on the time interval (0, 1).

  put the mixed problem under a matrix form, we denote by b and c the vectors of unknowns for the concentration b h and the density c h , respectively. The degrees of freedom in b are the values of b h at all grid vertices and those of c made of the values of c h at all the vertices except at x 1 = 0

  T ).The symbol T is for the transposition. The matrix A r , of type (n, n), is linked to a r (•, •) and coincides with the mass matrix with respect to the weight r. It is symmetric positive definite. The notation M is for the rectangular mass matrix of type (n, n -1). It is extracted from the full mass matrix representing the L 2 -inner product in V h . K and K * stand for the stiffness matrices related to the mixed bilinear forms k(•, •) and k * (•, •), respectively. They are rectangular, of type (n, n -1).Moreover, the obvious coerciveness of both k(•, •) and k * (•, •) on the space H h × H h yields that the matrices K and K * are full rank matrices. Consequently their transposes K T and (K * ) T are injective. To put the discrete problem under a compact form we list all the degrees of freedom in a single vector called Y . This vector contains therefore (2n -1) degrees of freedom and is solution of the initial value problem

Remark 4 . 1

 41 The dimension of N (M ) depends in fact upon the boundary points where Cauchy's conditions are enforced on the density c. If they are also prescribed at point x = π, the dimension of N (M ) becomes two.

Figure 3 :

 3 Figure 3: Space-varying parameters d (black), r (green) and r * (red). Non-uniform grids.

6. 2 2 m 2 s - 1 .

 2221 The inverse problemWe use the dissolved oxygen balance model for the study of the state of a portion of a streamwater with a length L = 2500 meters. This portion is contaminated by some organic polluting load η(t) = b(0, t) at the outfall point x = 0. Ultimately we need to quantify this time-varying density η(•). The way to achieve this consists in collecting measurements on the dissolved oxygen at x = 0 during two days, that is T = 48 × 3600 seconds. Then we tackle computationally the inverse problem exposed above to reconstruct the BOD load η(t) using the observations of the DO profile α(t) = c(0, t) at the same extreme-point x = 0. Let us remark that in this prospecting work the 'observed' function α(•) is in reality synthesized. The BOD load η(•) being known, the function α(t) is obtained after numerically solving the well posed direct BOD-DO problem (1) and[START_REF] Blinchikoff | Filtering in the time and frequency domains[END_REF]. A contrario, we will use this 'observed function' for the re-constitution of η(t) by dealing with the inverse problem. In Figure4, both functions are depicted, one, η(•), in the left diagram and the other, α(•), in the right plot. Throughout the simulations, our choice for the physical parameters is as follows r = 0.3 × 10 -4 s -1 , r * = 0.6 × 10 -4 s -1 , v = -0.1 ms -1 , d(x) = 1 40 50 + 7.5 × 10 -4 x 1.0 + 1.5 × 10 -4 x Symbols m and s are for meter and second. The space dependent dispersion d = d(x) is varying increasingly in the interval [34.2975, 35.5837].

ure 4 .Figure 4 :Figure 5 :

 445 Figure 4: BOD load η(•) (left) and DO profile α(•) (right).

Figure 6 :

 6 Figure 6: Computed BOD loads for noisy data.

Figure 7 :Figure 8 :

 78 Figure 7: Computed BOD loads at the outfall point after applying a Gaussian filter.

Table 1 ,

 1 the modulus of the eigenvalue that is the most close to zero. The algebraic multiplicity is thus one and so is the geometric multiplicity. As a result the dimension of nilpotent matrix N in the canonical Weierstrass-Kronecker form is one; actually it is identically zero. The Kronecker index is hence one, µ = 1.

	n =	10	20	30	40	50	Grids
	(d, r, r * ) = (1, 1, 1) 8.3	2.0	0.9	0.5	0.3	Uniform
	(d, r, r * ) = (1, 1, 1) 1.1 0.72	0.14	0.045 0.018 Non-uniform
	Varying (d, r, r * )	0.75 0.048 0.0095 0.0030 0.0012 Non-uniform

Table 1 :

 1 The modulus of the smallest non-vanishing eigenvalue (×10 3 ).