Tracking Mobile Objects with Several Kinect using HMMs and Component Labelling Amandine Dubois 1,2 and François Charpillet 2,1 Abstract-This paper proposes a markerless system whose purpose is to detect falls of elderly people at home. To track human movements we use Microsoft Kinect camera which allows to acquire at the same time a RGB image and a depth image. One contribution of our work is to propose a method for fusioning the information provided by several Kinects based on an occupancy grid. The observed space is tesselated into cells forming a 3D occupancy grid. We calculate a probability of occupation for each cell of the grid independently of its nearby cells. From this probability we distinguish whether the cells are occupied by a static object (wall) or by a mobile object (chair, human being) or whether the cells are empty. This categorization is realized in real-time using a simple three states HMM. The use of HMMs allows to deal with an aliasing problem since mobile objects result in the same observation as static objects. This paper also deals with the issue regarding the presence of several people in the field of view of camera by using the component labelling method. The approach is evaluated in simulation and in a real environment showing an efficient realtime discrimination between cells occupied by different mobile objects and cells occupied by static objects.

I. INTRODUCTION

During the last fifty years, the many progresses of medicine as well as the improvement of the quality of life have resulted in a longer life expectancy. The increase of the number of elderly people is a matter of public health because although elderly people can age in good health, old age also causes embrittlements in particular on the physical plan which can result in a loss of autonomy. If we look at elderly people (people of more than 65 years old), we realize that one of the main preoccupations at this age is fall. Indeed by looking at the figures of the INPES [START_REF] Bourdessol | Prévention des chutes chez les personnes âgées à domicile[END_REF] we notice that one elderly people out of three falls in the year. Our work is related to fall prevention and detection.

Many systems exist to detect falls. One of the categories consists in systems with sensors, that the person wears on her. These sensors are either accelerometer, gyroscopes or goniometers. These various sensors can be integrated in devices detecting the fall automatically, as shown in article of Bourke et al. [START_REF] Bourke | Evaluation of a threshold-based tri-axial accelerometer fall detection algorithm[END_REF]. There exists also systems made of an alarm button, in this case it's the person who must press herself on a button to alert after the fall. Another category of fall detection devices are the systems using 2D or 3D cameras. For example Jansen et al. [START_REF] Jansen | 3D human pose recognition for home monitoring of elderly[END_REF] used a 3D camera in order to classify the pose of a person among an a priori set of characteristic poses: "standing", "sitting" or "lying down". We chose a RGB depth camera, more precisely the Kinect camera. One of the disadvantages of the Kinect camera is its rather reduced field of view. The sensor can maintain tracking approximately 0.7-6 m. However we want to detect falls in parts of an apartment often larger than what the Kinect camera is able to see. Thus to cover a room entirely we should integrate several cameras in the same room.

This article is focused on the problem of human movements tracking with several Kinect. For this, three problems need to be solved. The first one is the fusion of the information provided by the different cameras. The second one is the discrimination between mobile objects and background. The last problem consists in gathering the mobile elements belonging to the same object. This stage makes it possible to identify several people as distinct in the same scene. In earlier work [START_REF] Dubois | Using hmms for discriminating mobile from static objects in a 3D occupancy grid[END_REF], we presented a method to detect pixels belonging to a mobile object adressing to the two first problems of human movements tracking. This paper entends previous work by adressing the third problem. We identify pixels belonging to the same object.

For making the fusion of the information provided by several 3D cameras, our method is to use a 3D occupancy grid. From the depth image, we create a spatial representation called an occupancy grid. In this representation, the space is divided into cells of a few centimeters with a probabilistic occupation state.

The second problem for tracking persons is to succeed at discriminating mobile objects from the background. Our approach is based on an extension of occupancy grids [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF] using hidden Markov models such that each voxel of the grid is determined by a three state model (the voxel belongs to the background, the voxel belongs to a mobile object, the voxel is not occupied). Our work is related to others. Among them, let us quote Yapo et al. [START_REF] Yapo | A probabilistic representation of LiDAR range data forefficient 3D object detection[END_REF] who proposed a method to detect 3D objects using LIDAR sensors. Their approach is also based on the concept of occupancy grids. From a probabilistic representation they determine if the voxels are free, occupied or hidden.

The third problem of identification of the mobile objects is carried out thanks to the component labelling method which consist in gathering the voxels belonging to the same mobile object.

This paper is organized as follows. Section II is dedicated to the data fusion of several Kinect cameras. Then, section III describes the method for categorizing the voxels (cells) occupation state. Section IV explains the method to identify the different mobile objects. Finally in section V, we present experimental results obtained with our method.

II. KINECT FUSION

We use an occupancy grid as a common coordinate system so as to fusion information from several Kinect cameras.

A. Occupancy grids

Occupancy grids such as defined in the article of Alberto Elfes [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF] consist in dividing into cells a 2D or in our case a 3D space. The grid provides a representation of the environment. For each cell C i (or voxels) we estimate its state, which can be either occupied or empty, from a probability of occupation P (C i |r). The probability of occupation is the probability of cell C i to be in state occupied given sensor reading r. For the sake of simplicity, each cell is estimated independently of its nearby cells.

B. Calibration of Kinect

The fusion information of several Kinect consists in determining the position and orientation of the different Kinect compared to the common coordinate system (the grid). For that we calibrate the Kinect cameras. The aim of camera calibration is to determine the transformation matrix (represented in Figure 1) between the RGB camera of each Kinect (K RgbRgb ), between RGB camera and the grid (K RgbGrid ), between RGB and depth camera of the same Kinect (K RgbDepth ) and between grid and depth camera (K GridDepth ). So with all these transformation matrix, the cells of the grid can be projected to the different camera coordinate systems (Kinect) and we will be able to define an observation function for each Kinect presented as in section III-B.

The transformation matrix K RgbRgb and K RgbGrid can be determined with the method of epipolar geometry [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] or chessboard calibration.

The transformation matrix K RgbDepth is known by manufacturer of Kinect (provided by openNI).

The transformation matrix K GridDepth is deduced from previous transformations with the following equation :

K GridDepth = K -1 RgbGrid × K RgbDepth

III. CLASSIFICATION OF CELLS

A. HMM models

We want to distinguish mobile objects (chair, human beings) from static objects (walls). We define a cell containing a mobile object as being an occupied cell that has previously been empty. Whereas cells containing a static object are cells that are occupied and that have never been empty. In the classical occupancy grid method the state "empty" or "occupied" is calculated from the probability of occupation P (r|C i ). We can consider the occupancy grid model as a two states HMM with no transitions between states as shown in Figure 2.

O E We modify the method of occupancy grid. For each cell, we use a three states HMM allowing to represent its dynamic and to determine its state. We remind that each cell is estimated independently of its nearby cells. The three states are:

1 1
• the state "O" meaning that the cell is occupied and has always been occupied; • the state "M" meaning that the cell is occupied but has already been empty at least once; • the state "E" meaning that the cell is not occupied. In other words we can write:

• for mobile objects:

C i = M • for statics objects: C i = O • for cells occupied by object: C i = M ∨ O
The representation of this HMM is shown in Figure 3. These probabilities respect the following assumptions :

O M E α 1 -α 1 -γ γ β 1 -β
• the detected cells in state O are regarded as being part of a wall. Consequently these cells have a weak probability of changing state. That's why the probability of transition of passing from O to another state is very weak. This probability is however not null because there is possibility that the model was mistaken since at the beginning it cannot make the difference between a wall and a mobile object which did not move yet (like a chair) ; • when a cell appeared at least once free, there is no possibility for this cell to return to O because it is not a cell belonging to a wall. If not, this cell would have never passed in a free state ; • the difference between a cell occupied in O and in M lies only in the passage or not in the state E that's why there is no transition between O and M.

B. Observation function

Each voxel C i is represented by its center of mass, defined by coordinates (x,y,z). We can obtain at which distance is located the voxel from the camera by projecting the voxel to the camera coordinate system using the camera transformation matrix K GridDepth . We denote as l this distance. The distance l of the voxel is compared to the depth, denoted as d, of the corresponding pixel provided by the Kinect camera (Figure 5). The observation r (see the section II-A) takes as value the error of distance (ε) between d and l calculated as ε = d -l. An observation function is built to evaluate the probability of occupation of the cell from the depth image Assuming that the information provided by the different camera is conditionally independent, we can multiply the different observation functions:

P (r|C i ) = f (ε). f (ε) is represented in Figure 4.
P (r 1 , ..., r N |C i ) = N j=1 f (ε j )
where N is the number of cameras and ε j the error of distance calculated with camera j.

C. Inference

To calculate the probability to be in one of the three states of HMM, we use to the Forward procedure [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. We don't use the procedure Backward because we want the result to be given online.

The three observation functions are given by :

P (r|C i = O) = f (ε) P (r|C i = M ) = f (ε) P (r|C i = E) = 1 -f (ε).
Cells are categorized by choosing the maximum a posteriori (MAP), that is to say the most likely state of the corresponding HMM.

We denote as po t , pm t , pe t the probability of a cell being respectively occupied, mobile or empty at time t. We fix the initial probability with po 0 = pe 0 = 0.5 and pm 0 = 0.0.

IV. IDENTIFICATION OF DIFFERENT MOBILE OBJECTS

In this part the aim is to gather the mobile cells belonging to the same object, so it will be possible to distinguish several persons in the same scene.

To gather the mobile cells we used the method «Component labelling» [START_REF] Suzuki | Linear-time connected-component labeling based on sequential local operations[END_REF]. This method consists in assigning a label (a number) to each cell detected as mobile (state M of HMM). In Figure 6a), the colored cells are cells in state M. Thus in Figure 6b), the algorithm assign a different number to each colored cells. Then the technique is to look for each cell p if one of its neighbors has a smaller number. In this case, one assigns the number of the smallest neighbor to the cell p. This operation is shown in figure 6c). The cells at step t+2 take the smallest number among their neighbors at step t+1. One carries out this operation until there is no more change in the assignment of the cells labels. The Figure 6d) is the last step of the algorithm because the cells labels can't change. Thus all the cells having a same number will be gathered as being a same mobile object. In Figure 6d) the algorithm has detected two different objects, one object with cell carrying number 1 and another object with cell having number 7.

One key question in this algorithm is what is most judicious distance to consider for defining cell neighbors? Should we limit neighborhood at juxtaposed cells? Or is it necessary it to consider a wider range?

One of the problems if we consider only juxtaposed neighbors of the cell (neighborhood of size 1) is that it is possible that these neighbors aren't detected as mobile cell. As a consequence the leg or the arm can be cut of the body and the leg or the arm will be identified as an object separated from the rest of the body. The second solution is to take farther cells in neighborhood, i.e. for cell C i if we defined a neighborhood of size 2, the algorithm takes the neighbor cells positioned at C i-2 , C i-1 , C i+1 and C i+2 . The problem with this method is the risk of integrating cells which don't belong to the same object.

For illustrating these two methods we can give the following example. Figure 6 shows the case where we considered only juxtaposed neighbors cells. Two object in 6d) were detected. On the other hand if we had considered the neighbor of size 3 in Figure 6d), the two blocks of mobile cells would have gathered in only one object. In section V-B.2 we discuss this problem.

V. RESULTS

A. Simulation

This section describes our method for evaluating the sensitivity and the specificity of the system in a simulated environment. The sensitivity is the capacity to detect mobile objects when they are present and the specificity is the capacity of the system to detect the absence of mobile objects when there is no mobile object.

Figure 7 shows the output of the simulator and the result of the system classification. In order to perform the evaluation, the output of the system should be compared to a reference image pixel by pixel. Since it is impossible to evaluate the system in real conditions due to the fact that we need to index real images, we propose to limit the quantitative evaluation to a simulated environment. We have recorded a simulated human activity in a virtual scene and used these images as a reference. We simulate a Kinect by generating depth and RGB images from the virtual scene. In addition, a reference image that index each pixel in the scene as static or mobile object is also generated. RGB and depth images are supplied to our system to perform the classification. Finally we compare the output of the algorithm to the reference image. Results show a sensitivity of 87.14% and a specificity of 98.52% for a total of 430 frames (73.8M pixels). In spite of good specificity we can notice that there are as many false-positives (pixel detected as mobile whereas the pixel is static) as of true-positives (pixel detected as mobile and pixel is mobile). The problem comes from the fact that in the reference images only 2% of the pixels corresponded to the moving person whereas the other 98% were static objects or background. However visually the mobile points are always near to mobile object. The model has a some inertia and so if a pixel is mobile it takes a certain time before to return a static pixel. Table I shows the number of pixels for static and mobile objects obtained from the simulation and detected by our system.

B. Behavior in a smart room

We have tested our algorithm in an experimental apartment. Results presented here are qualitative. In Figure 8(a) and 8(b) we see the RGB and depth images. The image sent back by the Kinect is illustrated by Figure 8(c). Black spaces correspond to non reconstructed zones.

1) Results with one camera and one mobile object: We have tested the algorithm with one camera and one person walking in front of the camera. As illustrated by Figure 9(a), walls, furnitures and the ground are correctly detected as static objects represented by green color and the person as a mobile object represented by blue color. We can see that the feet of the person in figure are detected as a static object, it's due to the size of the voxels (6 cm) and the uncertainty of the observation. The feet of the person are integrated in voxels representing the ground. We can notice that there is very limited noise on the background where a few badly positioned blue cubes remain. Moreover the tracking of mobile objects is fast enough to distinguish visually the members (leg, arm) of the person as it can be seen on Figure 9(a). A space without color is present above in the left of this figure. This is due to the size of the grid which is limited here to the perception range of the Kinect. The obstacles in a room don't disturb the discrimination between mobile and static objects as shown on Figure 9(b).

When we move a furniture, this furniture, previously detected as a static object, is recognized as a mobile object. Figure 10 shows a chair becoming a mobile object. This result is allowed by the transition (O → E) which models the fact that a furniture can be moved resulting in new empty space.

After a certain amount of time, it could be interesting to consider a furniture that has been moved as a new static object. This can be realized simply by adding a link (M → O) with a small probability γ 2 to the transition matrix as illustrated by Figure 11.

2) Results with one camera and several mobile objects:

To test component labelling we realized a situation where several persons were in the field of view of the Kinect camera. The results of this test are showned in Figure 12(a). In this figure we haven't represented the occupancy grid, nor the scene but only the points representative of the center of mass of each cell of the object. Each group of points detected as one object is represented by a different color.

We said in section IV that to consider neighbor of size 1 is not a judicious choice because the parts of the body can be detected as separate objects. Thus to avoid this defect we increase the neighborhood to 3. But as before said the risk for taking a depth too large is to gather distinct object as shown in Figure 12(b). In this case two persons are in the room but they are too close and are considered as a single object (represented by a same color). When there is more space between the two persons as Figure 12(a), the algorithm correctly detected the two persons as separated mobile objects (as we can see because each object is in a different color).

3) Results with two cameras: To finish the experiment was then realized with two cameras placed as illustrated in Figure 13(a). We can see that the fusion of several cameras allows to discover more space while increasing the noise around static objects as illustrated in Figure 13. The noise is

C. Behavior in realistic conditions

We have tested the algorithm, out of laboratory conditions, in a real apartment. One of the differences is on the level of lighting. The experimental apartment is located in a larger room with wall painted in black and not having windows. Thus lighting comes primarily from the artificial light. We wanted to test the algorithm in a more natural scene. We can see in Figure 14(c) that there are less noise compared to the experimental apartment. But we have noticed that when there is too much sun light on a white surface, the Kinect badly reconstructs the zone which is represented by black color in the lower right corner of Figure 14(a) which is the room where the test carryied out. We notice visually that the results are correct, the algorithm detects correctly the person.

VI. CONCLUSION

We presented a markerless system using Kinect cameras in the aim of tracking elderly people at home. First we proposed a system to merge several cameras by using a 3D occupancy grid. Secondly, compared to previous work on occupancy grid we proposed a method to allow the tracking of mobile objects. This method is based on a three states HMM: cell is empty, cell has always been occupied (static objects) and cell is occupied but has already been empty (mobile objects). This three state HMM is a simple yet elegant solution for solving a state aliasing problem (the observation for a static object is the same as the observation for a mobile object). Since each cell is updated independently one of the other, the process can be easily parallelized and implemented in a GPU allowing real-time (30 FPS) tracking with 2 cameras on a 1M cell grid. Thirdly, to solve the problem of distinguishing several persons in the field of view of the Kinect camera, we have implemented component labelling. This algorithm gather the cells belonging to the same object. Results in simulation allowed us to measure the quality of classification performed by the system in terms of sensitivity and specificity. Results on real images concerning the detection of cells occupied by mobile objects are visually satisfying. Concerning the detection of different persons the results are correct in most cases but the algorithm doesn't adress all the problems. For example, when two persons are too close the algorithm can't distinguish them. In the continuity of this work we will have to improve this part of the algorithm.

The aim of this project, as said in introduction, is to detect falls of elderly people. In this article, we presented the first

Fig. 1 .

 1 Fig. 1. Kinect grid calibration.

Fig. 2 .

 2 Fig. 2. Representation of the occupancy grid model (O: occupied; E: empty). Notice that P (O) = 1 -P (E).

Fig. 3 .

 3 Fig. 3. The same HMM is used to model the evolution of each cell. Probabilities of transition are α = 0.01, β = 0.1 and γ = 0.4.

Fig. 4 .

 4 Fig. 4. Representation of occupation probability.

Fig. 6

 6 Fig. 6. a) mobile cells are colored, b) assignment of numbers to mobile cells, c) cells with number higher than it neighbor are modified, d) final result.

  (a) Simulated apartment. (b) Distinction between static and mobile objects.

Fig. 7 .

 7 Fig. 7. Results in simulation.

Fig. 8 .

 8 Fig. 8. Image of Kinect camera.

Fig. 9 .

 9 Fig. 9. Green color: static objects. Blue color: mobile object.

  (a) The chair is considered as a static object. (b) The chair has been moved and is considered as a mobile object.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Chair becoming a mobile object.

  (a) Two persons are in the room and two objects are detected. (b) Two persons are in the room but only one object is detected.

Fig. 12 .

 12 Fig. 12. Results with component labelling.

  (a) Position of the two cameras in the apartment. (b) View of one of the cameras. (c) View of the other cameras. (d) Fusion of the two cameras.

Fig. 13 .

 13 Fig. 13. Test with two cameras.

  (a) View of camera Kinect. (b) All the objects are detected as static (view without texture).(c) A person is detected as a mobile object.(d) A person is detected as a mobile object (view without texture).

Fig. 14 .

 14 Fig. 14. The use of the algorithm in a real apartment.

TABLE I NUMBER

 I OF PIXELS IN EACH CATEGORY.

part of this project. It would be necessary in continuation of this work to learn characteristics of a person (as her size...) so as to be able to recognize her, track her and detect her activity (sitting, standing...). The purpose is to learn the habits of a person for thus detecting when an unusual behavior occurs (for example lying on the ground).

VII. ACKNOWLEDGEMENT

This work has been partly funded by Region Lorraine. The authors thanks Abdallah Dib for his contribution to this work and Cedric Rose for his valuable comments.


Ludivine Pradeleix
email: ludivine.pradeleix@irstea.fr

Bellon Véronique

P Maurel

O Roux

S Philippon

Bouarfa

Véronique Bellon Maurel

Philippe Roux

Olivier Philippon

Sami Bouarfa

Life Cycle Assessment at the regional scale: innovative insights based on the systems approach used for uncertainty characterization

Keywords: variability, LCI, irrigation, agrarian system diagnosis

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Rapidly increasing population growth and food requirements call for an increase in agricultural production, out of 40% is provided by irrigated areas. The International Water Management Institute has called for "more crop per drop" [START_REF] Molden | Water for Food, Water for Life: A Comprehensive Assessment of Water Management in Agriculture[END_REF]. And drip irrigation is endorsed as it is less water intensive par area than surface irrigation. But eventually this technique led to manifold indirect environmental impacts such as cropland extension, global input intensification and groundwater overexploitation due to individually managed tube wells . It is therefore essential to develop tools in order to assess the impact of various agricultural planning scenarios in irrigated areas.

Life Cycle Assessment could be a candidate tool. As public water management and decision making is carried out, Life Cycle Assessment (LCA) should be applied at regional scale, which is still a methodological challenge [START_REF] Guinée | Life Cycle Assessment: Past, Present, and Future[END_REF]. In spite of being product oriented, several authors underline LCA relevancy at regional scale, namely for a farming region, defined as a set of farms in a given geographical area [START_REF] Aubin | Méthodes D'évaluation. Panorama Critique. in duALIne -Durabilité De L'alimentation Face à De Nouveaux Enjeux. Questions à La Recherche[END_REF]Payraudeau and Van der Werf, 2005). Accounting for the variability of farming systems and management practices is one of the main challenges in agricultural LCA [START_REF] Nemecek | The Advantages of Generic LCA Tools for Agriculture1: Examples SALCAcrop and SALCAfarm[END_REF]. Unlike the case described by Payraudeau et al., (2005), the farm population in an irrigated area is too large to be surveyed one by one. Reviewing 70 LCA studies conducted in tropical and semi-arid locations, [START_REF] Basset-Mens | Is LCA-based Eco-labelling Reasonable1? The Issue of Tropical Food Products[END_REF], highlighted the failure to account for farming system diversity and the lack of specific data and data collection methods. Building up farm typologies is a way of dealing with the great variability of flows related to agricultural practices [START_REF] Dalgaard | Modelling Representative and Coherent Danish Farm Types Based on Farm Accountancy Data for Use in Environmental Assessments[END_REF]. Some studies already performed LCA at regional scale. At this scale, statistics e.g. the Farm Accounts Data Network (FADN) were used for building farm typologies [START_REF] Dalgaard | Modelling Representative and Coherent Danish Farm Types Based on Farm Accountancy Data for Use in Environmental Assessments[END_REF][START_REF] Mishima | Adaptation of Life Cycle Assessment (LCA) to Agricultural Production on a Regional Scale in Japan[END_REF], mainly because they are in line with LCA for being product oriented [START_REF] Meeusen | Agricultural Data for Life Cycle Assessments[END_REF]. Nonetheless, this approach cannot be widely applied: on the one hand, very few countries in the world offer agricultural statistics device and on the other hand these statistics are based on economic inputs and outputs (I/O) and thus present several drawbacks for LCA purposes. Agriculture is poorly described by I/O because it is mainly dependent on many self produced resources [START_REF] Haas | Life Cycle Assessment Framework in Agriculture on the Farm Level[END_REF] and also because data are expressed in monetary units; a large uncertainty [START_REF] Huijbregts | Application of Uncertainty and Variability in LCA. Part I: A General Framework for the Analysis of Uncertainty and Variability in Life Cycle Assessment[END_REF]) is then linked with major environmental impacting flows such as fertilizers and chemicals. Moreover, the European Commission, (2010) mentioned that when performing a LCA, data related to the foreground system, i.e. field scale agricultural production, should be inventory specific.

Facing a double challenge of data scarcity in the context of Southern countries and high variability of crop management practices in irrigated areas, we propose a new methodological framework to model activity data and build a regional model of agriculture, with the objective of reducing the uncertainty. In this paper, we propose to adapt the methodology of agrarian system diagnosis for modelling farming systems to create as accurate as possible Life Cycle Inventories. This methodological framework enables us to characterise the uncertainty linked to the "real world variability" and to imprecision also called epistemic uncertainty [START_REF] Huijbregts | Application of Uncertainty and Variability in LCA. Part I: A General Framework for the Analysis of Uncertainty and Variability in Life Cycle Assessment[END_REF]. Preliminary results of the agrarian system diagnosis conducted in Tunisia are presented. Finally relative pros and cons of this methodological framework for conducting a regional LCA are discussed.

Material and Methods

2.1 Adapting the methodology of Agrarian System Diagnosis for building a Life Cycle Inventory at regional scale Below is a brief description of the several steps of the Agrarian System Diagnosis (ASD) hereafter designated "diagnosis" and of adaptations made for LCA purposes (right column). This framework for building LC inventory should allow us to reduce the global uncertainty in the LCA; this demonstration is part of the results.

Step We propose to mobilize the Agrarian System Diagnosis (ASD), hereafter designated "diagnosis", to model the agricultural region for LCA purposes. ASD was initially designed for targeting farm diversity in development projects. For being systems oriented, ASD aims at understanding the diversity and complexity of regional agricultural production modes at different scales and then model them into a farming systems typology. Farming systems all together are interconnected and compose the agrarian system at the regional scale [START_REF] Cochet | L'agriculture Comparée[END_REF]. Each farming system is modelled as functionally representative of a set of comparable production units. These units carry out a given combination of cropping systems (crop rotation and associated cultivation practices) and livestock systems and rely on comparable resources and socio-economics constraints [START_REF] Cochet | Fonctionnement Et Performances Économiques Des Systèmes De Production Agricole1: Une Démarche à l 'échelle Régionale[END_REF][START_REF] Moreau | Reconciling Technical, Economic and Environmental Efficiency of Farming Systems in Vulnerable Areas[END_REF].

The modelling process is progressive but not linear, and iterative with several feedback procedures. Table 1 is a brief description of the several steps of the ASD and of some adaptations made for LCA purposes. Starting from a global standpoint by analysing landscape heterogeneity on maps (# 0) and in the field (#1a), several hypothesis about spatial distribution of cropping systems are formulated. Then, assumptions are checked during field surveys (#2) and cropping systems are modelled; other hypothesis on their combination are made into a pre-typology of farming systems. For each pre type, a set of representative farms is sampled (#3a) and in depth interviewed (#3b). Finally, an archetype is designed, whose agricultural practices and economical values are modelled for a "normal year", i.e. exceptional events are not modelled (#1b). The archetype is modelled for being for the most probable case according to the farm structure, its objective, opportunities and constraints. This approach is system-oriented: it uses triangulation for ensuring data reliability, cross checking structural, functional and historical information about farming systems. In the same vein, disciplinary viewpoints and scales of analysis enrich data consistency. Finally, technical and economic thresholds are calculated for each farming system for outliers identification. A restitution to surveyed people and local expert allows us checking data completeness and validating their reliability.

Characteristics of the irrigated plain of Kairouan, Tunisia

Located in central Tunisia, in semi arid to arid climate, area under study is an alluvial plain of 30 000 ha and comprises around 2 000 farms. Agriculture has much evolved with drip irrigation introduction, from sheep herding and rain fed cereals and low density olive groves to irrigated vegetables, fruit orchards and high density olive groves. Groundwater provides irrigation water and is overexploited. Nonetheless, economic profitability of irrigated crops led people to drill unauthorised boreholes. A pilot area of 6 000 ha out of 30 000 ha was selected for being a hotspot in terms of water exploitation and intensification of agricultural management practices, i.e. several crop cycles per year and numerous intercropping. Data were collected by two students during a three month stay.

Results

Hereafter, we demonstrate how the new framework based on ASD for LCI can support the characterization of uncertainty sources in a regional LCA and public decision making for land planning options. Uncertainty sources are manifold when aiming at modelling the Life Cycle Inventory of an agricultural region. They are usually separated into variability of the "real world" and uncertainty [START_REF] Huijbregts | Application of Uncertainty and Variability in LCA. Part I: A General Framework for the Analysis of Uncertainty and Variability in Life Cycle Assessment[END_REF].

Methodological output: the Agrarian System Diagnosis as a methodology to characterise uncertainty in agricultural Life Cycle Inventories

Figure 1 describes the sources of uncertainty in LCI of agricultural systems and the solutions proposed to characterise them via the ASD framework

In the upper part of the figure, uncertainty sources found at the regional scale are listed, in line with uncertainty classification proposed by [START_REF] Huijbregts | Application of Uncertainty and Variability in LCA. Part I: A General Framework for the Analysis of Uncertainty and Variability in Life Cycle Assessment[END_REF]. In the lower part the way the Agrarian System Diagnosis contributes to characterise each uncertainty source is explained. In the very bottom part we describe how regional LCA outcomes can support public land planning decision making.The Variability between Sources and Objects (VBSO) stands for the "differences in inputs and emissions of comparable processes in a product system"; parameter uncertainty is caused by inaccurate, unrepresentative, incomplete data e.g. chemicals specifications; uncertainty due to choices originates from choices made regarding allocation, Functional Unit, and the LCIA stage that is out of our scope; model uncertainty also occurs during the LCIA stage. Five out of the six categories of uncertainty sources are addressed by the ASD. Only model uncertainty is totally beyond our study scope.

As explained above, uncertainty related to variability of farming systems and management practices is addressed by building a functional typology of farming systems and cropping / livestock systems based on practices modelling. Practices are contextualised and different from standard technical guidelines.

The typology allows accounting for VBSO1 and spatial variability over the studied region. Then, a sample set of farms is in-depth interviewed for each identified type. Intra type variability (VBSO234should4 be less than inter-type variability (VBSO13. If not, a new type should be designed by splitting the type with high variability. Exceptional values due to temporal variability (e.g. climate hazards) for instance and looking inconsistent with regard to the functioning identified are discarded. The archetype built at farm scale is drawn after parameter uncertainty has been reduced and allocations have been made clear; This is done respectively by running several procedures of data consistency checking (cf. part 2.2) and by surveying the whole farming system and interconnections between farms.

Preliminary Results

Nine archetypes of farming systems and sub systems (cropping and livestock systems) were modelled in the pilot area, illustrating a high variability. A large range of cultivated species, crop sequences and associated crop management practices were identified. This is usual in irrigated areas because rainfall is no more a limiting factor for growing crops. Olive groves, a crop shared by most of farming systems, is part of very diverse cropping systems. Indeed, density, irrigation management and intercropping vary among olive-grove based cropping systems. Vegetables intercropping is widespread, mostly during unproductive stage of perennial crops. The smaller the cultivated area is, the more crops are intercropped. The level of intensification regarding crop inputs and water quantities vary from a factor one to five. Land productivity is high in case of several crop cycles per year or overlapping crop cycles. These results illustrate the degree of complexity of cropping systems in irrigated areas, especially in case of intercropping, and the need for having a typologybased inventory.

Indeed, the methodology proposed for archetype modelling enable us to design the typology, according to uncertainty source and also according to its magnitude in case of "Variability Between Sources and Objects".

Discussion

The objectives of a farming system typology for LCA purposes were: "to lower data variability, thereby allowing a better selection of representative farms for detailed research; better determine the marginal effects of a studied change." [START_REF] Lindeijer | Conclusions of the Working Group on Farm Typologies[END_REF].

The methodological framework we propose, namely the ASD-based LCA is a powerful method for LCAoriented data collection that accounts for farming systems and practices variability and provide LCA specific activity data "for a marginal supplementary effort." ASD is of higher interest since farming systems and sub systems modelling are based on in depth analysis of agricultural practices that are to be turned into LCI data. Consistency of each type is ensured by calculations of technical and economical thresholds. The consistency of activity data is also enhanced by collecting data at different scales, crossing field observations and farm and literature surveys.

The farm archetype that is modelled is neither a virtual farm average [START_REF] Basset-Mens | Uncertainty of Global Warming Potential for Milk Production on a New Zealand Farm and Implications for Decision Making[END_REF][START_REF] Dalgaard | Modelling Representative and Coherent Danish Farm Types Based on Farm Accountancy Data for Use in Environmental Assessments[END_REF] nor a single farm chosen by experts for being representative [START_REF] Haas | Estimation of Environmental Impact of Conversion to Organic Agriculture in Hamburg Using the Life-Cycle-Assessment Method[END_REF]). On the contrary, regarding temporal variability, the farm archetype is built for a "normal year" whereas statistics or farm accounting data would refer to a single year. Thus, yields were "modelled" to represent a normal year, instead of being averaged. For example, due to climate hazards yields of pepper ranged from 1.6 to 3.2t/ha; the "modelled" yield was 3t/ha and the average one 2.5t/ha, i.e. there is a 25% difference which would change significantly the LCA outcomes, especially if impacts are expressed per mass unit.

Parameter uncertainty can be large if used data are not specific. In ASD-based LCA we are able to calibrate our own data collection to LCA requirements and include critical flows expected to heavily impact conversely to statistics and FADN data [START_REF] Dalgaard | Modelling Representative and Coherent Danish Farm Types Based on Farm Accountancy Data for Use in Environmental Assessments[END_REF]. Indeed, fertilisers elements, agrochemicals properties and details about intercropping are crucially lacking into these generic databases. Extension services acknowledge that some values, such as cultivated areas in intercropping systems, can be registered twice in local statistics (personal comm.); this may be highly misleading since this practice is widespread on our field study.

ASD allows us to identify innovative cropping systems and current tendencies even if in minority, conversely to statistics already outdated when published. The holistic standpoint provides important insights regarding allocation rules. Indeed, statistics prevent the LCA practitioner from designing allocation rules among the numerous products of the farm, especially mixed ones which are highly represented in our study area. Efole [START_REF] Efole Ewoukem | Environmental Impacts of Farms Integrating Aquaculture and Agriculture in Cameroon[END_REF] highlighted that mixed up farming systems tend to make the most of their limited resources and allocate biomass among several productions, thus making allocation rules more complex. Other limitations of the statistic approach for building farming systems typology is that farm functioning cannot be described [START_REF] Dalgaard | Farm Types -How Can They Be Used to Structure, Model, and Generalise Farm Data?[END_REF], by-products and near-to-zero values are overlooked [START_REF] Lindeijer | Conclusions of the Working Group on Farm Typologies[END_REF], and data could be too much aggregated [START_REF] Cochet | Fonctionnement Et Performances Économiques Des Systèmes De Production Agricole1: Une Démarche à l 'échelle Régionale[END_REF].

Conclusion

ASD demonstrated its ability to help us designing and characterising typical farming systems and their crop and livestock components, in qualitative and quantitative terms. This method of typology is particularly relevant when data are scarce and key criteria for classifying the population of farms cannot be taken from statistics nor from expert knowledge [START_REF] Dalgaard | Modelling Representative and Coherent Danish Farm Types Based on Farm Accountancy Data for Use in Environmental Assessments[END_REF][START_REF] Haas | Estimation of Environmental Impact of Conversion to Organic Agriculture in Hamburg Using the Life-Cycle-Assessment Method[END_REF]. The agrarian system approach decreases uncertainty linked to the inherent variability of "real agriculture" (i.e. farming systems and management practices). Unlike statistics that process data and deduce mathematical correlations between variables, this approach make causalities clear within the frame of each farming system functioning.

Moreover, by revealing material and energy interconnection flows between farms or within a farm, it allows for clearer burden allocation rules among the different product systems. Double counting could also be avoided through the ASD holistic standpoint. Indeed, farming systems that are diversified are likely to support interconnection flows and thus would particularly benefit from this approach.

In addition, this work supports the identification of each farming system's room for manoeuvre to mitigate their environmental impact, within agro-ecological and technical values that define their range of existence.

Based on this methodology and the typology, our next objective is to complete the characterisation of farming systems archetypes, conduct LCA and lastly assess alternative land planning scenarios for agriculture based on the diagnosis outputs.

Figure 1 .

 1 Figure 1. Sources of uncertainty in LCI of agricultural systems and solutions proposed to characterise them within the Agrarian System Diagnosis (ASD) framework

Table 1 .

 1 Proposed methodology for adapting the Agrarian System Diagnosis (ASD) to Life Cycle Assessment.

		Major Steps of Agrarian Systems Diagnosis	Original	Agrarian	Systems	Agrarian Systems Diagnosis
	n°		Diagnosis			LCA oriented
	0	Choice of a pilot zone	Selection criteria: Most of knowl-	Selection criteria: Worst case
		Study of available information sources: Maps, pre-	edge			with regard to potential environ-
		vious studies (soil, slope, climate, water re-				mental impacts
		sources…)
	1a	Landscape analysis / identification of agro ecologi-	Identify cultivation dynamics, spa-	Identify vulnerable areas with
		cal units (soil, slope…) and pre-types of cropping	tial distribution		regard to major impacts
		and livestock systems
	1b	Historical analysis / interviews	For capturing past differentiation	Identify innovative systems
		Climate hazards frequency	processes.			Foresee potential evolutions
	2	Surveys of cropping and livestock systems	Focus on spatial distribution, crop	Investigate co-product destina-
		(diversity, varieties, soil fertility management, ani-	sequences, crop-animal interactions	tion, material flows between
		mal feeding calendar…)				farms
			Farms illustrating differentiation	Focus on potential environmental
			processes			impacts drivers (contrasted
						yields, fertilizers and agrochemi-
						cals).
						Select farms with most records.
	3b	In depth interviews of Farming Systems: techno-	Focus on Practices and Economy	Focus on input/output quantifica-
		economic characterization: cross-checking of quali-	Focus on farm strategy, opportuni-	tion (e.g. fertilisers & Agro-
		tative and quantitative information, iterative proc-	ties and bottlenecks linked to capi-	chemicals), including internal
		ess, systems triangulation procedure	tal, labour force…		flows
	4	Extrapolation to the whole area	Based on local knowledge about the	Only necessary for "snapshot"
			representation of each type in the	LCA, not for agricultural plan-
			whole area.			ning scenarios.

3a

Sampling design of farms to be surveyed for each farming system pre-type (steps 1 &2). Sampling criteria: maximise diversity, farms chosen according to criteria explaining diversity.

8th International Conference on LCA in the Agri-Food Sector, Rennes, France, 2-4 October 2012


Jean-Antoine Désidéri
email: jean-antoine.desideri@inria.fr

Manuel Bompard
email: manuel.bompard@onera.fr

Jacques Peter
email: jacques.peter@onera.fr

Hermitian Interpolation Subject to Uncertainties

. Various types of meta-models exist (interpolation polynomials, neural networks, Kriging models, etc.). Such meta-models are constructed by pre-calculation of a database of functional values by the costly high-fidelity model. In adjoint-based numerical methods, derivatives of the functional are also available at the same cost, although usually with poorer accuracy. Thus, a question arises: should the derivative information, available but known to belessaccurate,beused to construct the meta-model or be ignored? As the first step to investigate this issue, we consider the case of the Hermitian interpolation of a function of a single variable, when the function values are known exactly, and the derivatives only approximately, assuming a uniform upper bound ! on this approximation is known. The classical notion of best approximation is revisited in this context, and a criterion is introduced to define the best set of interpolation points. This set is identified by either analytical or numerical means. If n + 1i st h en u m b e ro fi n t e r p o l a t i o np o i n t s ,i ti s advantageous to account for the derivative information when ! ≤ ! 0 ,w h e r e! 0 de- creases with n,andthisisinfavorofpiecewise,low-degreeHermitianinterpolants. In all our numerical tests, we have found that the distribution of Chebyshev points is always close to optimal, and provides bounded approximants with close-to-least sensitivity to the uncertainties.

Introduction: the Classical Notion of Best Approximation

In this section, we review certain classical notions on polynomial interpolation, in particular the concept of "best approximation" or "Chebyshev economization". The literature contains numerous elementary and advanced textso nt h i sf u n d a m e n t a l issue, and we refer to [START_REF] Conte | Elementary numerical analysis: An algorithmic approach[END_REF][START_REF] Pearson | Handbook of applied mathematics: Selected results and methods[END_REF][START_REF] Quarteroni | Numerical mathematics[END_REF].

Let n be an integer and x 0 , x 1 ,...,x n be n + 1d i s t i n c tp o i n t so ft h en o r m a l i z e d interval [-1,1]. Let "(x) be the following polynomial of degree n + 1:

"(x)= n # i=0 (x -x i ) (1)
and consider the following n + 1polynomialsofdegreen:

L i (x)= n # j=0 j =i x -x j x i -x j (i = 0, 1,...,n) . (2)
Clearly ∀i, j ∈{0, 1,...,n} :

L i (x j )=$ i, j , (3)
where $ stands for Krönecker's symbol. Application of L'Hospital's rule yields the following compact formula:

L i (x)= "(x) " ′ (x i )(x -x i ) . (4)
Let f : [-1, 1] → R be a smooth function of the real variable x.Thepolynomial

P n (x)= n % i=0 f (x i ) L i (x) (5)
is of degree at most equal to n,a n di tc l e a r l ys a t i s fi e st h ef o l l o w i n gi n t e r p o l a t i o n conditions: ∀i ∈{0, 1,...,n} :

P n (x i )= f (x i ) . (6)
One such interpolant is unique among all polynomials of degree ≤ n.Thus,P n (x) is the Lagrange interpolation polynomial of f at the points {x i } 0≤i≤n .

It is well known that if f ∈ C n+1 ([-1, 1]),foranygivenx ∈ [a, b],theinterpolation error is given by

e n (x)= f (x) -P n (x)= f (n+1) (& ) (n + 1)! "(x) (7)
for some & ∈ [-1, 1].

Proof. Let x ∈ [-1, 1] be given. If x = x i for some i,theresultistri viallysatisfied. Otherwise, "(x) = 0; then let ' = e n (x) "(x) so that f (x)=P n (x)+'"(x)

and define the function

( (t)= f (t) -P n (t) -'"(t) t ∈ [-1, 1].
The function ( (t) is of class C n+1 ([-1, 1]) and it admits a nonempty set of roots in the interval [-1, 1] that includes X = {x 0 , x 1 ,...,x n , x}.T h en + 2e le m e n tso fX are distinct and can be arranged as the elements of a strictly increasing sequence {x 0 i } 0≤i≤n+1 whose precise definition depends on the position of x w.r.t. the interpolation points {x i } 0≤i≤n .ByapplicationofRolle'stheoremto( (t)=( (0) (t) over the subinterval [x 0 i , x 0 i+1 ], i = 0, 1,...,n,itfollowsthat( ′ (t) admits a root x 1 i in the open interval ]x 0 i , x 0 i+1 [,a n dt h i s ,f o re a c hi.I nt h i sw a yw ei d e n t i f yas t r i c t l y -i n c r e a s i n g sequence of n + 1r o o t so f( ′ (t), {x 1 i } 0≤i≤n .T h e nR o l l e ' st h e o r e mc a nb ea p p l i e d in a similar way, this time to ( ′ (t),andsoontothesuccessi v ederi v ati v esof( (t). We conclude that in general ( (k) (t) admits at least n + 2 -k distinct roots in [- 1,1], {x k i } 0≤i≤n+1-k ,0≤ k ≤ n + 1. In particular, for k = n + 1, ( (n+1) (t) admits at least one root, x n+1 0 ,hereafterdenoted& for simplicity. But since

P (n+1) n (& )=0and " (n+1) (& )=(n + 1)!, one gets ' = f (n+1) (& ) (n + 1)!
and the conclusion follows.

Hence, if

K = 1 (n + 1)! max x∈[-1,1] f (n+1) (x) (8) we have ∀x ∈ [-1, 1] : |e n (x)| ≤ K |"(x)| . (9)
Therefore, a natural way to optimize the choice of interpolation points apriori,that is, independently of f ,istosolvethefollowingclassicalmin-maxproblem:

min {x i } 0≤i≤n x i ∈[-1,1], ∀i max x∈[-1,1] |"(x)| . (10)
In view of this, the problem is to find among all polynomials "(x) whose highestdegree term is precisely x n+1 ,a n dt h a ta d m i tn

+ 1d i s t i n c tr o o t si nt h ei n t e r v a l [-1, 1],anelement,uniqueornot,thatminimizesthesup-normover[-1, 1].
The solution of this problem is given by the n + 1rootsoftheChebyshevpolynomial of degree n + 1. Before recalling the proof of this, let us establish some useful auxiliary results. Let k be an arbitrary integer and T k (x) denote the Chebyshev poly-

nomial of degree k.Recallthatforx ∈ [-1, 1] T k (x)=cos k cos -1 x , k ∈ N, (11) so that, for k ≥ 1andx ∈ [-1, 1], T k+1 (x)+T k-1 (x)=cos k + 1) + cos k -1) = 2cos(k) ) cos ) = 2xT k (x), (12)
where one has let

x = cos ) , 0 ≤ ) ≤ ". (13)
Thus, if the leading term in T k (x) is say a k x k ,thefollowingrecursionapplies:

a k+1 = 2a k , k ≥ 1, (14)
and, since a 0 = a 1 = 1, it follows that

a k = 2 k-1 , k ≥ 1. ( 15
)
Therefore, an admissible candidate solution for the min-maxp r o b l e m ,( 1 0 ) ,i st h e polynomial

" ⋆ (x)= 1 2 n T n+1 (x). (16)
It remains to establish that " ⋆ (x) is the best choice among all admissible polynomi- als, and its roots the best possible interpolation points. Toa r r i v ea tt h i s ,w ec l a i m the following lemma:

Lemma 1. Fo r a l l a d m i s s i b l e p o l y n o m i a l q (x) one has

" ⋆ ≤ q , ( 17
)
where is the sup-norm over [-1, 1].

Proof. Assume otherwise that an admissible polynomial q(x) of a strictly smaller sup-norm over [-1, 1] exists: q < " ⋆ .

(18)

Let r(x)=" ⋆ (x)q(x).Sincetheadmissiblepolynomials" ⋆ (x) and q(x) have the same leading term, x n+1 ,thepolynomialr(x) is of degree at most n.Letusexamine the sign of this polynomial at the n + 2points

* i = cos i" n + 1 , i = 0, 1,...,n + 1, (19)
at which " ⋆ (x) as well as T n+1 (x) achieve a local extremum. At such a point,

|" ⋆ (* i )| = 1 2 n = " ⋆ > q = max x∈[-1,1] |q(x)| ≥ |q(* i )| (20)
and r(* i ) is nonzero and has the sign of the strictly dominant term " ⋆ (* i )=

1 2 n T n+1 (* i )= (-1) i 2 n .
T h e r e f o r e ,r(x) admits at least n + 1s i g na l t e r n a t i o n sa n da s many distinct roots. But this is in contradiction with the degree of this polynomial. The contradiction is removed by rejecting the assumption made on q .

Consequently, in (10), the min-max is achieved by the roots of T n+1 (x):

x ⋆ i = cos (2i + 1)" 2(n + 1) , i = 0, 1,...,n, (21)
and the value of the min-max is 1 2 n .

Best Hermitian Approximation

Now assume that the points {x i } 0≤i≤n are used as a support to interpolate the function values {y i = f (x i )} 0≤i≤n ,b u ta l s ot h ed e r i v a t i v e s{y ′ i = f ′ (x i )} 0≤i≤n ,t h a ti sa set of 2(n + 1) data. Thus, we anticipate that the polynomial of least degreec o mplying with these interpolation conditions, say H 2n+1 (x),isofdegreeatmostequal to 2n + 1. One such polynomial is necessarily of the form

H 2n+1 (x)=P n (x)+"(x) • Q(x), ( 22
)
where the quotient Q(x) should be adjusted to comply with the interpolation conditions on the derivatives. These conditions are

y ′ i = H ′ 2n+1 (x i )=P ′ n (x i )+" ′ i • Q(x i ), i = 0, 1,...,n, (23)
where

" ′ i = " ′ (x i )= # j=0 j =i (x i -x j ) = 0, (24)
and since "(x i )=0. Thus Q(x) is solely constrained by the following n + 1interpolation conditions:

Q i = Q(x i )= y ′ i -P ′ n (x i ) " ′ i , i = 0, 1,...,n. (25)
Therefore, the solution of least degree is obtained when Q(x) is the Lagrange interpolation polynomial associated with the above function values:

Q(x)= n % i=0 Q i L i (x). (26)
The corresponding solution is thus unique, and we will refer to it as the global Hermitian interpolant.

The interpolation error associated with the above global Hermitian interpolant H 2n+1 (x) is given by the following result valid when

f ∈ C 2n+2 ([-1, 1]): ∀x ∈ [-1, 1], ∃* ∈ [-1, 1] : f (x)=H 2n+1 (x)+ f (2n+2) (*) (2n + 2)! " 2 (x). ( 27
)
Proof. Let x ∈ [-1, 1] be given. If x = x i for some i,theresultistri viallysatisfied. Otherwise, "(x) = 0; then let i of + ′ (t) exists. But the interpolation points, at which the derivative also is fitted, are themselves n + 1otherdistinctrootsof+ ′ (t). Thus we get a total of at least 2n + 2r o o t sf o r+ ′ (t),a n db yi n d u c t i o n ,2 n + 1f o r + ′′ (t),andsoon,andfinallyone,say*,for+ (2n+2) (t).Now ,sinceH

µ = f (x) -H 2n+1 (x) " 2 (x) so that f (x)=H 2n+1 (x)+µ" 2 (x) . Let +(t)= f (t) -H 2n+1 (t) -µ" 2 (t), t ∈ [-1, 1]. The function +(t) is of class C 2n+2 ([-1, 1]),
(2n+2) 2n+1 (*)=0 because the interpolant is of degree 2n + 1a tm o s t ,a n ds i n c e(" 2 (t)) (2n+2) (t)= (2n + 2)!, one gets 0 = f (2n+2) (*) -0 -µ (2n + 2)!
which yields the final result.

As a consequence of (27), the formulation of the best approximation problem for the global Hermitian interpolant is as follows:

min {x i } 0≤i≤n x i ∈[-1,1], ∀i max x∈[-1,1] " 2 (x). (28)
But, if we define the following functions of R n+1 → R:

     P(x 0 , x 1 ,...,x n )= max x∈[-1,1] " 2 (x), p(x 0 , x 1 ,...,x n )= max x∈[-1,1] |"(x)|, (29)
it is obvious that

∀x 0 , x 1 ,...,x n : P(x 0 , x 1 ,...,x n )=p 2 (x 0 , x 1 ,...,x n ). (30)
Hence the functions P and p achieve their minimums for the same sequence of interpolation points, and

min {x i } 0≤i≤n x i ∈[-1,1], ∀i P(x 0 , x 1 ,...,x n )=    min {x i } 0≤i≤n x i ∈[-1,1], ∀i p(x 0 , x 1 ,...,x n )    2 = 1 4 n . ( 31
)
Therefore the best Hermitian interpolation is achieved for the same set of interpolation points as the best Lagrangian interpolation, that is, the roots, {x ⋆ i } 0≤i≤n of (21), of the Chebyshev polynomial T n+1 (x).

Best Inexact Hermitian Approximation

In PDE-constrained global optimization [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF], it is often useful to model the functional criterion to be optimized by a function f (x) of the design variable x ∈ R n ,w h o s e values are computed through the discrete numerical integration of a PDE, and the derivative f ′ (x),orgradientvector ,bymeansofanadjointequation.Adatabase of function values and derivatives is compiled by Design of Experiment, and the surrogate model, or meta-model is constructed from it. This meta-model is then used in some way in the numerical optimization algorithm with the objective of improving computational efficiency (see, for example, [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]). The success of such a strategy depends on the accuracy of the meta-model f (x) to represent the dependency on x of the actual functional criterion. If all the data were exact, and properly used, the accuracy would undoubtedly improve by the addition of the derivative information. However, in practice, since the PDE is solved discretely, thederi v ati vesarealmost inevitably computed with inferior accuracy. Therefore it isnotclearthataccounting for the derivatives is definitely advantageous if the corresponding accuracy of the data is poor. Should special precautions be taken to guarantee it?

In order to initiate a preliminary analysis of this problem, we examine the simple one-dimensional situation of a function f (x),w h e nx is scalar (x ∈ R), and consider the case of a Hermitian interpolation meta-model basedo ni n e x a c ti n f o r m ation. Specifically, we assume that the function values {y i } 0≤i≤n are known, whereas only approximations { y ′ i } 0≤i≤n of the derivatives {y ′ i } 0≤i≤n are available, and we let:

$ y ′ i = y ′ i -y ′ i := ! i , i = 0, 1,...,n. (32)
Hence the computed interpolant is H 2n+1 (x) instead of H 2n+1 (x),andinviewofthe definitions ( 22)-( 25), we have

             $ H 2n+1 (x)=H 2n+1 (x) -H 2n+1 (x)="(x) $ Q(x), $ Q(x)= n % i=0 $ Q i L i (x), $ Q i = $ y ′ i " ′ i = ! i " ′ i . (33)
Now, suppose an upper bound ! on the errors ! i 's is known:

|! i | ≤ !, 0 ≤ i ≤ n. ( 34
)
The following questions arise:

1. What is the corresponding upper bound on max x∈[-1,1] |$ H 2n+1 (x)|? 2.
Can we choose the sequence of interpolation points {x i } 0≤i≤n to minimize this upper bound? 3. Is the known sequence of the Chebyshev points a good approximation of the optimum sequence for this new problem?

This article attempts to bring some answers to these questions. Presently, we try to identify how the interpolation points should be defined to minimize or reduce the effect on the meta-model accuracy of uncertainties on the derivatives only. It follows from (33) that

$ H 2n+1 (x)="(x) n % i=0 ! i " ′ i L i (x) (35)
which by virtue of (4) simplifies as follows:

$ H 2n+1 (x)=" 2 (x) n % i=0 ! i " ′ i 2 (x -x i ) . (36)
Thus if (34) holds, we have

|$ H 2n+1 (x)| ≤ !,(x) (37)
where

, (x)=" 2 (x) n % i=0 1 " ′ i 2 |x -x i | . ( 38
)
These considerations have led us to analyze the min-max problem applied to the new criterion , (x) in place of " 2 (x).Insummary,thesolutionofthemin-maxprob- lem associated with the criterion , (x) minimizes the effect of uncertainties in the derivatives on the identification of the global Hermitian interpolant. In the subsequent sections, this solution is identified formally, or numerically, and compared with the Chebyshev distribution of points, which is optimal w.r.t. the interpolation error. Lastly, the corresponding interpolants are comparedbynumericalexperiment.

Formal or Numerical Treatment of the Min-Max-, Problem

We wish to compare three particular relevant distributions of interpolation points in terms of performance w.r.t. the criterion , (x).Thesethreedistributionsaresymmetrical w.r.t. 0, and recall that the total number of interpolation points is n + 1. Thus, we let

n + 1 = 2m + - (39)
and when n is odd (-= 0; n = 2m -1 ≥ 1),

{x i } 0≤i≤n = {±& 1 , ±& 2 ,...,±& m }, (40)
where

0 < & 1 < & 2 <...<& m ( 41
)
and m ≥ 1. Otherwise, when n is even (-= 1; n = 2m ≥ 0), we adjoin to the list & 0 = 0(once).W econsiderspecifically:

1. The uniform distribution:

n = 2m: & u 0 = 0associatedwiththeinterpolationpointx 0 = & 0 = 0, and & u k = k m ,1≤ k ≤ m,associatedwith2interpolationpoints±& u k . n = 2m -1: & u k = 2k-1 n ,1≤ k ≤ m. 2. The Chebyshev distribution: & ⋆ k = x ⋆ m-k = cos 2(m -k)+1 n + 1 " 2 , 1 ≤ k ≤ m. ( 42
)
3. The optimal distribution:

& = arg min & max x∈[0,1] , (x; & ), (43)
where & =(& 1 , & 2 ,...,& m ) denotes the vector of adjustable parameters defining, along with & 0 = 0ifn is even, the distribution of interpolation points, and the dependence of the criterion , on & is here indicated explicitly for clarity. (Note that due to symmetry, the interval for x has been reduced to [0, 1] without incidence on the result.)

To these three distributions are associated the corresponding values of the maximum of , (x; & ) over x ∈ [0, 1];thesemaximumsaredenoted, u , , ⋆ and , ,respectively.

As a result of these definitions, the polynomial "(x) is expressed as follows:

"(x)=x - m # k=1 x 2 -& 2 k , n + 1 = 2m + -; -= 0or1, (44)
and for x > 0, the criterion , (x) becomes

, (x)=" 2 (x) n % i=0 1 " ′ i 2 |x -x i | = " 2 (x) - " ′ 0 2 1 x + m % k=1 1 " ′ k 2 1 x + & k + 1 |x -& k | .
(45) Then, given x,let j be the index for which

& j-1 ≤ x < & j (46) so that x -& k ≥ 0 ⇐⇒ k ≤ j -1 . ( 47
)
As a result,

, (x)=" 2 (x) - " ′ 0 2 1 x + j-1 % k=1 1 " ′ k 2 2x x 2 -& 2 k + m % k= j 1 " ′ k 2 2& k & 2 k -x 2 . ( 48
)
Calculation of the derivatives " ′ k :

First, for -= 0, "(x) is an even polynomial, and " ′ 0 = 0. Otherwise, for -= 1,

" ′ 0 = lim x→0 "(x) x = m # k=1 -& 2 k , -= 1; n = 2m. ( 49
)
Regardless -,fork ≥ 1

"(x)=x -(x -& k )(x + & k ) # i=1 i =k x 2 -& 2 i ( 50
)
so that:

" ′ k = lim x→& k "(x) x -& k = 2& -+1 k # i=1 i =k & 2 k -& 2 i . ( 51
)

Expression of $

0 Applicable Whenever n + 1 is Even (-= 0; n + 1 = 2m) Suppose that 0 < x < & 1 .Then j = 1and(48)reducesto , (x)=" 2 (x) m % k=1 1 " ′ k 2 2& k & 2 k -x 2 . ( 52
)
But,

" 2 (x)= m # i=1 x 2 -& 2 i 2 (53) so that , (x)= m % k=1 2& k " ′ k 2 × m # i=1 i =k & 2 i -x 2 2 × & 2 k -x 2 . ( 54
)
All the terms in this sum are composed of three factors: a positive constant and two positive factors that are monotone decreasing as x varies from 0 to & 1 .Hence, (x) is monotone decreasing and

$ 0 = max x∈[0,& 1 ] , (x)=, (0) = m % k=1 2& 3 k " ′ k 2 m # i=1 i =k & 4 i = 2 m # i=1 & 4 i m % k=1 1 & k " ′ k 2 , -= 0; n + 1 = 2m. ( 55
)

Expression of

$ 1 Applicable Whenever & m < 1 Suppose that & m < x < 1. Then j = m + 1and(48)reducesto , (x)=" 2 (x) - " ′ 0 2 1 x + m % k=1 1 " ′ k 2 2x x 2 -& 2 k = - " ′ 0 2 x m # i=1 x 2 -& 2 i 2 + m % k=1 2 " ′ k 2 x 2-+1 x 2 -& 2 k m # i=1 i =k x 2 -& 2 i 2 . ( 56
)
All the terms in , (x) are products of positive factors that are monotone-increasing with x.Consequently ,themaximumisachievedatx = 1. But

, (1)= - " ′ 0 2 m # i=1 1 -& 2 i 2 + m % k=1 2 " ′ k 2 1 -& 2 k m # i=1 i =k 1 -& 2 i 2 . ( 57
)
This gives

$ 1 = max x∈[& m ,1] , (x)=, (1)= m # i=1 1 -& 2 i 2 - " ′ 0 2 + 2 m % k=1 1 " ′ k 2 1 -& 2 k .
(58)

Special Cases

For n = 0, 1, 2, the formal treatment is given in [START_REF] Bompard | Best Hermitian interpolation in presence of uncertainties[END_REF].

For n = 3, the four interpolation points form a symmetrical set {±& , ±*}.Hence the optimization is made over two parameters & and *.Thus,thefunction

, (x) (59)
is defined discretely. The iso-value contours of this function are indicated in Fig. 1, which permits after refinement to identify the optimum & ≈ 0.351 and * ≈ 0.926.

For n = 4, the five interpolation points form a symmetrical set 0, {±& , ±*}. Hence the optimization is again made over two parameters & and *.Thus,thefunction z(*, *) is again defined by (59) and evaluated discretely. The iso-value contours of this function are indicated in Fig. 2, which permits after refinement to identify the optimum & ≈ 0.571 and * ≈ 0.948.

General Results (n > 4)

The min-max-, problem has been solved by either analytical or numerical means for values of n in the range from 0 to 40. The results are collected in Table 1 inwhich the first column indicates the number of interpolation points n + 1, the second gives the definition of the Chebyshev points & ⋆ (n ≤ 4), the third provides the definition of the optimal distribution & ,andthefourthacomparisonofperformancebygiving, when available, the values of 2. , ⋆ = max x , (x, & ⋆ ),t h eu p p e rb o u n do n, (x) corresponding to the approximately optimal distribution & = & ⋆ of interpolation points (the Chebyshev dis- tribution); 3. , u = max x , (x, & u ),t h eu p p e rb o u n do n, (x) corresponding to the uniform distribution & = & u of interpolation points.

The analytical results are related to the cases for which n ≤ 4, and have been outlined in a previous subsection.

For n + 1 ≥ 10, the distribution & (not given here) has been identified by a numerical minimization realized by a particle-swarm (PSO) algorithm. The table indicates the corresponding values of , .

From these results, one observes that the upper bound , achieved when the distribution of interpolation points is optimized, is not only bounded, but it even diminishes with increasing n.T h eC h e b y s h e vd i s t r i b u t i o nh a sa na l m o s te q u i v a l e n t performance. Inversely, the uniform distribution yields a value of the upper bound , u that is unbounded with n.Inconclusion,usingtheChebyshevdistribution,which is known explicitly, is highly recommended in practice.

5G e n e r a l i z e dH e r m i t i a nI n t e r p o l a t i o n

In this section, we generalize the notions introduced in the first three to the situation where one wishes to construct a(the) low(est)-degree polynomial interpolant of the values, as well as the derivatives up to order, say p (p ∈ N), of a given smooth

; , = max x , (x, & ), , ⋆ = max x , (x, & ⋆ ) and , u = max x , (x, & u ), where & u i = -1 + 2 n-1 (i -1), i = 1,...,
y (k) i = f (k) (x i ), k = 0, 1,...,p; i = 0, 1,...,n. ( 60
)
The interpolation polynomial is denoted H n,p (x) and it is now constrained to the following (p + 1)(n + 1) interpolation conditions: ∀k ∈{0, 1,...,p}, ∀i ∈{0, 1,...,n} :

H (k) n,p (x i )=y (k) i . ( 61
)
We associate such kind of interpolation with the expression "generalized Hermitian interpolation".

Existence and Uniqueness

We first establish existence and uniqueness by the following:

Theorem 1. There exists a unique polynomial H n,p (x) of degree at most equal to (p + 1)(n + 1) -1 satisfying the generalized interpolation conditions (61).

Proof. By recurrence on p.F o rp = 0, the generalized Hermitian interpolation reduces to the classical Lagrange interpolation, whose solution is indeed unique among polynomials of degree at most equal to (p + 1)(n + 1) -1 = n:

H n,0 (x)=P n (x). ( 62
)
For p ≥ 1, assume H n,p-1 (x) exists and is unique among polynomials of degree at most equal to p(n + 1) -1. This polynomial, by assumption, satisfies the following interpolation conditions: ∀k ∈{0, 1,...,p -1}, ∀i ∈{0, 1,...,n} :

H (k) n,p-1 (x i )=y (k) i . (63)
Hence, by seeking H n,p (x) in the form

H n,p (x)=H n,p-1 (x)+R(x), (64)
one finds that R(x) should be of degree at most equal to (p + 1)(n + 1) -1a n d satisfy ∀k ∈{0, 1,...,p -1}, ∀i ∈{0, 1,...,n} :

R (k) (x i )=0( 6 5 )
and ∀i ∈{0, 1,...,n} :

R (p) (x i )=y (p) i -H (p) n,p-1 (x i ). (66)
Now, (65) is equivalent to saying that R(x) is of the form

R(x)= n # i=0 (x -x i ) p • Q(x)="(x) p Q(x) (67)
for some quotient Q(x).Then,thederivativeoforderp of R(x) at x = x i is calculated by the Leibniz formula applied to the product u(x)v(x) where

u(x)=(x -x i ) p , v(x)= n # j=0 j =i (x -x j ) p • Q(x). (68)
This gives

R (p) (x i )= p % k=0 p k u (k) (x i ) v (p-k) (x i ). ( 69
) But, u (k) (x i )=0forallk except k = p yielding R (p) (x i )=p! v(x i )=p! n # j=0 j =i (x i -x j ) p Q(x i )=p! " ′ (x i ) p Q(x i ) . (70)
Thus, all the interpolation conditions are satisfied iff the polynomial Q(x) fits the following interpolation conditions: ∀i ∈{0, 1,...,n} :

Q(x i )=Q i = R (p) (x i ) p!" ′ (x i ) p = y (p) i -H (p) n,p-1 (x i ) p!" ′ (x i ) p . (71)
Therefore, solutions exist, and the lowest-degree solutionisuniquelyobtainedwhen Q(x) is the Lagrange interpolation polynomial associated with the above function values. This polynomial is of degree at most equal to n.Hence,R(x) and H n,p (x) are of degree at most equal to p(n + 1)+n =(p + 1)(n + 1) -1.

Interpolation Error and Best Approximation

We have the following:

Theorem 2 (Interpolation error associated with the generalized Hermitian interpolant). Assuming that f

∈ C (p+1)(n+1) ([-1, 1]),wehave ∀x ∈ [-1, 1], ∃& ∈ [-1, 1] : f (x)=H n,p (x)+"(x) p+1 f ([(p+1)(n+1)]) (& ) [(p + 1)(n + 1)]! . (72) Proof. Let x ∈ [-1, 1] be fixed. If x = x i for some i ∈{0, 1,...,n}, f (x)=H n,p (x)
and "(x)=0, and the statement is trivial. Hence, assume now otherwise that x = x i for any i ∈{0, 1,...,n}.Then,definetheconstant

. = f (x) -H n,p (x) "(x) p+1 (73) so that f (x)=H n,p (x)+."(x) p+1 . (74)
Now using the symbol t for the independent variable, one considers the function

) (t)= f (t) -H n,p (t) -."(t) p+1 . (75)
By virtue of the interpolation conditions satisfied by the polynomial H n,p (t), ∀k ∈{0, 1,...,p}, ∀i ∈{0, 1,...,n} :

) (k) (x i )=0( 7 6 )
but, additionally, by the choice of the constant .,wealsohave

) (x)=0. (77)
This makes n + 2distinctzeroesfor) (x) : x 0 , x 1 ,...,x n and x.Thus,byapplication of Rolle's theorem in each of the n + 1consecutiveintervalsthatthesen + 2points once arranged in an increasing order define, a zero of ) ′ (t) exists, yielding n + 1 distinct zeroes for ) ′ (t),t ow h i c h( 7 6 )a d d sn + 1d i s t i n c ta n dd i f f e r e n to n e s ,f o r at o t a lo f2 (n + 1)=2n + 2z e r o e s .S t r i c t l yb e t w e e nt h e s e ,o n efi n d s2 (n + 1) -1 zeroes of ) ′′ (t) to which (76) adds n + 1d i s t i n c ta n dd i f f e r e n to n e s ,f o rat o t a lo f 3(n + 1) -1 = 3n + 2zeroes.Thus,foreverynewderivative,wefindonelesszero in every subinterval, but n + 1m o r eb yv i r t u eo f( 7 6 ) ,f o rat o t a lo fn more, and this as long as (76) applies. Hence we get that ) (p) (t) admits at least (p + 1)n + 2 distinct zeroes. For derivatives of higher order, the numberofzeroesisonelessfor every new one; hence, (p + 1)n + 1f o r) (p+1) (t),a n ds oo n .W efi n a l l yg e tt h a t ) ([p+(p+1)n+1]) (t)=) ([(p+1)(n+1)]) (t) admits at least one zero & ,thatis

0 = f ([(p+1)(n+1)]) (& ) -.[(p + 1)(n + 1)]!( 7 8 )
because H ([(p+1)(n+1)]) (& )=0sin ceth ede g reeo fH n,p (t) is at most equal to (p + 1)(n + 1) -1, and the conclusion follows.

As a consequence of this result, it is clear that the best generalized Hermitian approximation is achieved by the Chebyshev distribution of interpolation points again.

Best Inexact Generalized Hermitian Interpolation

Now suppose that all the data on f and its successive derivatives are exact, except for the derivatives of the highest order, {y (p) i } that are subject to uncertainties {! i } i=0,1,...,n .Then,theuncertaintiesonthevalues{Q i } i=0,1,...,n of the quotient Q(x) are the following:

$ Q i = ! i p!" ′ (x i ) p ;( 7 9 )
on the quotient itself the following:

$ Q(x)= n % i=0 ! i p!" ′ (x i ) p L i (x)="(x) n % i=0 ! i p!" ′ (x i ) p+1 (x -x i )
;( 8 0 ) and, finally, the uncertainty on the generalized Hermitian interpolant H n,p (x) the following:

$ H n,p (x)="(x) p+1 n % i=0 ! i p!" ′ (x i ) p+1 (x -x i ) . (81)
In conclusion, for situations in which the uncertainties {! i } i=0,1,...,n are bounded by the same number !,t h ec r i t e r i o nt h a to n es h o u l dc o n s i d e rt oc o n d u c tt h em i nmax optimization of the interpolation points {x i } i=0,1,...,n is now the following one to replace the former , (x):

For given n and p,thisexpressionisunboundedinx.Thus,(theboundon)theerror is inevitably degraded in the order of magnitude due to presence of uncertainties.

However, the actual dilemma of interest is somewhat different. The first interpolant differs from f (x) by the sole interpolation error, µ n,p-1 |"(x)| p . The second interpolant is associated with the higher-order interpolation error, µ n,p |"(x)| p+1 ,butissubjecttotheuncertaintyterm!, (p) (x),whichisdominant,as we have just seen. Thus, the decision of whether to include derivatives or not should be guided by the ratio of the uncertainty term, !, (p) (x),tothelo werinterpolation error, µ n,p-1 |"(x)| p .Thisratioisequalto

! µ n,p-1 |"(x)| n % i=0 1 p! |" ′ (x i )| p+1 |x -x i | (89)
and it admits the bound

!B n,p µ n,p-1 , (90)
where the bound

B n,p = max x∈[-1,1] |"(x)| n % i=0 1 p! |" ′ (x i )| p+1 |x -x i | (91)
exists since, in the above, the function over which the max applies is piecewise polynomial for fixed n and p. Hermitian interpolation is definitely preferable whenever the expression in (90) is less than 1. This criterion permits us to identify trends as !, n and p vary, but is not very practical in general since the factors ! and µ n,p-1 are problem-dependent and out of control. The variation with n of the bound B n,p has been plotted in Fig. 3 for p = 1, 2 and 3. Visibly, the bound B n,p can be large unless p and n are small. Therefore, unsurprisingly, unless n and p,a sw e l la st h eu n c e r t a i n t yl e v e l!,a r e small enough, the criterion in (90) is larger than 1, and the interpolant of the sole exactly known values is likely to be the more accurate one.

To appreciate this in a practical case, we have considered thecaseoftheinterpolation of the function

f (x)= f ' (x)= 1 1 + ' x 2 (92)
over the interval [-1, 1] for p = 0( L a g r a n g ei n t e r p o l a t i o n )a n dp = 1( H e r m i t i a n interpolation). This smooth function is bounded by 1, and itsm a x i m u md e r i v a t i v e increases with ' .For' = 64/27, this maximum is equal to 1. For ' = 256/27, this maximum is equal to 2.

In the first experiment (Figs. 4 and5), ' = 64/27 and n = 5. The Lagrange interpolant is fairly inaccurate, mostly near the endpoints. Thus the error distribution indicates that the approximate Hermitian interpolant is preferable even for a fairly high level of uncertainty on the derivatives (20 % is acceptable).

In the second experiment (Figs. 6 and7), the interpolated function is the same, but the number n is doubled (n = 10). Consequently, the Lagrange interpolant of the sole exact function values is very accurate. The approximateHermitianinterpolant can only surpass it if the level of uncertainty on the derivatives is small (less than 5%).

Lastly, with the same number of interpolation points (n = 10), we have considered the case of a function with larger derivatives (' = 256/27). As a result (see Figs. 8 and9), the accuracy of the Lagrange interpolation hasb e e ns e v e r e l yd egraded. Then again, the approximate Hermitian interpolation is found superior for higher levels of uncertainty in the derivatives (the switch is between 20 and 50 %).

6C o n c l u s i o n s

Recalling that the Chebyshev distribution of interpolationpointsisoptimalw.r.t.the minimization of the (known bound on the) interpolation error, we have proposed an alternate criterion to be subject to the min-max optimization. The new criterion to be minimized aims at reducing the sensitivity of the Hermitian interpolant of function values and derivatives, to uncertainties assumed to be present in the derivatives only. We have found by analytical developments and numericale x p e r i m e n t st h a t the Chebyshev distribution is close to be optimum w.r.t. thisnewcriterionalso,thus giving the stability of the corresponding approximation a somewhat larger sense.

We have also considered the generalized Hermitian interpolation problem in which the derivatives up to some order p (p > 1) are fitted. For this problem we have derived the existence and uniqueness result, as well as the expression of the interpolation error, and also the definition that one could use for the criterion to be subject to the min-max optimization to reduce the sensitivity of the interpolant to uncertainties in the derivatives of the highest order. We conjectured from the derived expression that the corresponding optimal distribution of interpolation points converges to the Chebyshev distribution as p → /.

Lastly, we have made actual interpolation experiments in cases of a function bounded by 1, whose derivative is either bounded by 1 or 2. These experiments have confirmed that the approximate Hermitian interpolant was superior to the Lagrange interpolant of the sole exact function values, when the uncertainty on the derivatives is below a certain critical value which decreases when n is increased.

In perspective, we intend to examine a much more complicated case in which the meta-model depends nonlinearly on the adjustable parameters, by means of semiformal or numerical analysis tools.

Fig. 1

 1 Fig. 1 Contour plot of function z(& , *) of (59) for n + 1 = 4: the plot is symmetrical w.r.t. the bisecting line * = & , and one location of the minimum is & = 0.351, * = 0.926

1 .Fig. 2

 12 Fig. 2 Contour plot of function z(& , *) of (59) for n + 1 = 5: the plot is symmetrical w.r.t. the bisecting line * = & , and one location of the minimum is & = 0.571, * = 0.948

Fig. 3

 3 Fig. 3 Coefficient B n,p as a function of n for p = 1, 2 and 3

Fig. 4 Fig. 5 Fig. 6 Fig. 7

 4567 Fig. 4 Case ' = 64/27 (max x | f ' (x)| = max x f ′ ' (x) = 1); function f ' (x) and various interpolapolynomials (n = 5)

Table 1

 1 Variation of the criterion max x , (x, & ), related to Hermitian interpolation with uncertain derivatives, for different choices of the set & = {& i }, i = 1,...,n,ofinterpolationpointsin[-1, 1], and different degrees, 2n + 1

  n.

	Number interpol. pts.: n + 1 Degree of interpolant: 2n + 1		Chebyshev points: & ⋆	& = arg min & max x ,	Performance: ,, ⋆ , u
	1 1				00	11 1
	2 3	± 1 √ 2	. = ±0.7071	±0.7549	0.3774	0.5	0.5
	3 5	±	√ 3 2	0 = ±0.8660 .	0 ±0.8677	0.3258	0.3755	0.3333
	4 7	± 1 2 -1 √ 8 ± 1 2 + 1 √ 8	. = ±0.3827 . = ±0.9239	±0.351 ±0.926	0.282	0.439	0.299
	5 9	± 5-8 √ 5 ± 5+ √ 5 8	0 . = ±0.5878 . = ±0.9511	0 ±0.571 ±0.948	0.249	0.652	0.262
	10							0.164	0.179
	19								39.
	11							0.154	0.167
	21								111.
	20							0.103	0.112
	39							3.9 × 10 6
	21							0.100	0.108
	41							1.3 × 10 7
	function f (x) over [-1, 1].Theinterpolationpointsareagaindenoted{x i } i=0,1,...,n ,
	and we use the notation

  It is the following:given the values {y i , y ′ i ,...,y

	(p-1) i (p) of the higher derivative {y i },whichofthefollo wingtwointerpolantsis(guaran-},0≤ i ≤ n,andcorrespondingly ,approximations
	teed to be) more accurate:
	1. the Hermitian interpolant of the sole exact values: {y i , y ′ i ,...,y (p-1) i or	},0≤ i ≤ n,
	2. the Hermitian interpolant of the entire data set?

, (p)

or, equivalently,

We note that this expression is a homogeneous function of "(x) of degree 0. We conjecture that the variations of the above criterion, as p → /,a r ed o m inated by those of the factor |"(x)|.Hence,inthislimit,theoptimaldistrib utionof interpolation points should approach the Chebyshev distribution.

Overall Bound on the Approximation Error

The quantity !, (p) (x) is an absolute bound on the error committed in the com- putation of the generalized Hermitian interpolant based on function and derivative values in presence of uncertainties on the derivatives of theh i g h e s to r d e r ,p,o n l y , when these are uniformly bounded by !:

where H n,p (x) represents the actually computed approximation.

On the other hand, the interpolation error is the difference between the actual function value, f (x),andthe"true"interpolant,H n,p (x),thatcouldbecomputedif all function and derivative information was known. The interpolation error satisfies

(85) where one has let

Consequently, we have

Now, examining the precise expression for , (p) (x),thatis(82),weseethatthe ratio of the second term to the first on the right of the above inequality is equal to