Formal verification of a static analyzer: abstract interpretation in type theory

Xavier Leroy

Inria Paris-Rocquencourt xavier.leroy@inria.fr (Joint work with David Pichardie, Sandrine Blazy, Jacques-Henri Jourdan, and Vincent Laporte.)

Abstract

Static analysis is the automatic inference and checking of simple properties of all executions of a program. Initially developed in the context of compilers to support code optimization, static analysis is very successful today for the formal verification of safety properties of critical software, owing to its good scalability. As is the case for all tools involved in the production and verification of critical software (compilers, code generators, program provers, model checkers), confidence in the results of a static analysis tool requires evidence that the tool is sound and correctly over-approximates all possible executions of the program. Such evidence can take the form of a soundness proof mechanized using a proof assistant [START_REF] Leroy | Formal verification of a realistic compiler[END_REF][START_REF] Herms | A certified multi-prover verification condition generator[END_REF].

Abstract interpretation [START_REF] Cousot | Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] is an elegant, powerful mathematical framework to define and reason about static analyses. In particular, it is not limited to so-called "non-relational" analyses (inferring properties of a single value or variable) and works naturally for "relational" analyses (inferring relations between several variables, such as linear inequalities). The classic presentation of abstract interpretation involves Galois connections. It has the advantage that, once the meaning of abstract data is chosen via a Galois connection, the abstract operators used by the static analyzer can, in principle, be derived systematically from the concrete semantics, in a way that is not only sound by construction, but also relatively optimal.

However, the theory of Galois connections is resolutely set-theoretical, involving non-computable functions and equational reasoning over set comprehensions, making it very hard to express in type theory and to use in a proof assistant such as Coq. To overcome this difficulty, Pichardie et al [START_REF] Pichardie | Interprétation abstraite en logique intuitionniste: extraction d'analyseurs Java certifiés[END_REF][START_REF] Cachera | Extracting a data flow analyser in constructive logic[END_REF] developed and mechanized an alternative presentation of abstract interpretation, using only the "γ" (concretization) part of Galois connections, viewed as relations "concrete-datum ∈ abstract-datum". The calculational style is lost, and relative optimality is no longer guaranteed, but soundness proofs are easily conducted with a proof assistant.

In the context of the Verasco project, we are currently trying to scale Pichardie's approach all the way to the development and Coq verification of a realistic static analyzer based on abstract interpretation for the CompCert subset of the C language. Proper modular decomposition is crucial to build the appropriate abstractions as a hierarchy combining numerical and memory abstract domains. While the general interface of a non-relational domain is well known, giving such an interface for relational domains is more challenging, and so is formulating generic composition operators (such as reduced products) between such domains. Another enabling technique is the opportunistic use of validation a posteriori to obviate the need to prove complicated algorithms such as fixpoint iteration with widening and narrowing, or operations over polyhedra for relational domains of linear inequalities [START_REF] Fouilhé | Efficient generation of correctness certificates for the abstract domain of polyhedra[END_REF].

Acknowledgments This work is supported by the Verasco project (ANR-11-INSE-003) of Agence Nationale de la Recherche (ANR).


Michel Grangeat
email: michel.grangeat@ujf-grenoble.fr

What is the engine of teacher development? CPD programmes vs teacher experience.*

This paper tackles the question of factors which underpin teacher competences development; it is drawn both upon a set of research from the European project S-TEAM 1 , and upon data collected in science classrooms. We will figure out to what extent teachers' approaches and practices with respect to inquiry based science teaching [IBST] could be altered by teacher education and training programmes or by teaching experiences. These two ways for enhancing teacher competences do exist within international expectations towards teachers. Currently, we could be concerned with a trend in many European countries, and specifically in France, which leads to maximize the role of teachers experience and so to minimize the effects of formal CPD programmes.

A continuum of teaching approaches and practices

The first problem we have to handle with is the way of identifying differences among teacher practices and conceptions.

Most studies about teachers' practices use two broad and opposite categories: the teacher-centred and the learner-centred approaches. The former is seen as a way of teaching in which students are more or less passive and the latter as a way that facilitates student's learning processes. Nevertheless, numerous researchers had shown variations within these broad categories and they had conceived them as two poles of a unique continuum [START_REF] Hudson | Comparing different traditions of teaching and learning: what can we learn about teaching and learning[END_REF].

The model we are elaborating will be drawn upon continuums of IBST conceptions and practices. One pole of these continuums will emphasise the teacher and content-centred approaches. The other pole will focus on students and learning-oriented approaches. These poles are not exclusive and opposite. We hypothesize that teachers and teacher educators locate themselves between those according to different factors that are included within the teaching context (available resources, exchanges within their department, etc.). Most of time, this location is quite unconscious but we might assume that the more teachers are aware of it the better their competences improve.

A system which organise professional knowledge

The second problem we need to address is the way teachers and teacher educators, as other actors in other professional sectors, elaborate knowledge about their activities.

Studies in vocational fields showed that actors elaborate knowledge about their work as a synthesis of knowledge that arises both from their education or training and from their lived experiences. This synthesis occurs when actors can overcome the contradictions between what standards and instructions imply, what theoretical knowledge indicates, and what professional reality shows. This knowledge concerns not only current actions but also the work process itself: thus it is called as work process knowledge [WPK] [START_REF] Fisher | Work process knowledge: origins of the concept and current development[END_REF].

1 Science-Teacher Education Advanced Methods was funded by EU FP7 SiS (http://www.s-teamproject.eu/)

Others studies, from a similar framework, focus on the organization of the professional knowledge. Researchers show that this knowledge is not fragmented but linked in a system that allows actors to act rapidly, relevantly, and efficiently. This activity system [AS] is oriented by the activity's object, and it depends of the work organisation, of the professional community's culture, of the available resources and tools within the situation, and also of the repertoire of actions which are mastered by the actors [START_REF] Engeström | Activity theory as a framework for analysing and redesigning work[END_REF].

According to these outcomes, we had specified what teacher knowledge is and how to grasp it [START_REF] Grangeat | Factors influencing teachers' professional competence development[END_REF]. It results firstly that teachers' activities are oriented by knowledge which overlaps their current tasks and concerns the schooling process itself (i.e. school project, curriculum, etc). It results secondly that the set of teacher knowledge is organized as a system according to the main dimensions of the activities. Teacher knowledge may be identified according of four elements:

1. Goal: the teachers' purpose that is held individually or collectively 2. Clue: the piece of information that is picked out from the teaching situation, that is seen as relevant by teachers, and that will activate specific teaching practices 3. Repertoire of actions: the teaching practices that had been triggered by the clue and was orientated by the goal 4. Reference knowledge: the set of knowledge that enable one situation to be matched to another in order to define and to justify an action This knowledge is adapted to the actual situation of each teacher. Within this situation, many factors alter the activities and the set of knowledge that is necessary for monitoring them: the grade level of the class, the specificities of the students, the content to be taught, the way the teacher masters this content, etc. Nevertheless, within each specific situation the set of knowledge used by a teacher could be identified. Henceforth, we will call this cognitive unit Teacher Work Process Knowledge [TWPK]. Teachers gain to organize these clusters of knowledge within a system according to the crucial aspects of the situation, as it happens for expert teachers. Henceforth, we will call this set of knowledge a Teacher Activity System [TAS].

This TAS may develop through two main factors: the individual history and experience which provide opportunities to enlarge the repertoire of actions; the education programmes which lead to renew the set of reference knowledge. Research findings are coherent and lead to conceive this development of professional competences as an evolution of conceptions and practices between two poles: the first one is centred on the core of the activities (e.g. as a teacher I have only to expose specific contents as clearly as possible to a whole class); the last one is open to the variability of the situations (e.g. as a teacher I have to consider the school collective projects and to act with equity towards all the students). This development could be identified through different modes from novice to expert but each mode doesn't consist of a definitive position: facing a new content or an unknown situation, an experienced teacher could act like a novice. These modes are not exclusive since expertise consists in playing on the whole continuum, being more or less aware of this.

Elaborating a model for a better understanding of science teacher activities

We had set out the main dimensions of IBST methods by analysing the S-TEAM products and we found out six dimensions: origin of questioning, nature of the problem, students' responsibility in managing the inquiry, teacher's acknowledgment of students' diversity, role of argumentation, and explanation of teacher's goals. To specify the four modes of each dimension, we observed and Figure 1 : IBST 6six dimensions model diagramme interviewed 20 science teachers which are novice, experienced or expert ones. They teach maths, sciences or technology in lower secondary schools. The TWPK of these teachers were used to specify each mode on the six dimensions (see appendix 1).

The researchers asked each ST to carry out a lesson which he or she considered as an IBST lesson. The lesson was videotaped from the back of the classroom. Afterwards, each ST was interviewed about the video. The ST was asked to stop the video when there was an event which had involved him or her in choosing among different alternatives. Doing that, STs were led to make explicit the goals which underlay the observed action, and the other components of their TWPK. Interviews had been fully transcribed and a text analysis had been used for extracting the TWPK and TAS from the transcriptions. On each dimension, each ST is given a score.

In order to test the model's validity, three groups amongst the 20 teachers will be compared. The commitment science teachers group [CST] consists of 8 sciences teachers. They are involved within CPD programmes based on teacher collaboration. The new science teachers group [NST] consists of 6 beginner teachers. During their first year as teacher, they intended to a specific training programme based on cooperation and controversy. The third group gathers 6 experienced science teachers [EST] who could only apply on their experience. Groups' medians were compared.

Comparing science teachers approaches and practices of IBST methods

The results show that CSTs' group reaches the upper modes of the model for 4/6 dimensions, NSTs' group for 3/6 dimensions, EST group reaches only for 2/6 dimensions (see Figure 1). Nevertheless, through the 3 groups, 'questioning' and 'explanation' are the dimensions which seem the most teacher-centred.

Consequently, within this sample, teacher experience cannot be evaluated as sufficient in order to develop IBST methods. Conversely, the development of complex methods in classroom, such as IBST, is supported by collaboration among teachers and with external professionals, such as teacher educators or researchers.

This strengthens the research methodology and the validity of the IBST 6 dimensions model.

Content-centred approaches

Learner-focused approaches Dimension 1: Who initiates the questioning? Dimension 3: What is the students' level of responsibility within the inquiry process? Dimension 4: How is managed the students' diversity of knowledge, needs and wills?

4.1-Teachers cope with specific pupils' behaviour in order to involve them within the inquiry.

4.2-Teachers adapt the task in order to maintain specific students' involvement.

4.3-Each student's team receives teachers' supervision.

4.4-Students with specific needs rely on specific adaptation of the inquiry situation.

Dimension 5: Which role is given to argumentation?

5.1-Teachers facilitate interactions through students within each team.

5.2-Teachers communicate each team's propositions to the whole class 5.3-Students are encouraged to consider their mates' assumptions, results, and conclusions.

5.4-Students are asked to justify their conclusions with respect to knowledge or evidence.

Dimension 6: What is the level of explanation of teacher's goals?


Michel Grangeat
email: michel.grangeat@ujf-grenoble.fr

Teacher Educators' Conceptualisations and Practices towards Inquiry Based Science Teaching: a renew research framework

Keywords: teacher development, teacher collaboration, teaching approaches, science education

Currently European teachers have to develop inquiry based science teaching [IBST] in order to enhance pupils' learning and motivation towards scientific objects. Such teaching approaches are learning-focused, and appear more complex than content-centred methods. This paper focuses on 8 proficient science teachers who are involved both within professional networks and with part time in service teacher education about IBST. It aims to identify their conceptualisations and practices, and particularly about their way to cope with pupils' diversity.

The lack of teacher confidence for IBST methods

According to many European projects1 , IBST approach combines four main specifications:

• authentic and problem based learning activities where there may not be a correct answer • a certain amount of experimental procedures, experiments and "hands on" activities, including searching for information

• self-regulated learning sequences where student autonomy is emphasised • discursive argumentation and communication with peers ("talking science").

(Doris Jorde & Olsen Moberg, 2010, p. 3) A round of seminars within the 15 European countries organized within S-TEAM project shows that in any national curriculum, such specifications consist of rules which could be narrowly applied by teachers. In fact, they combine plural elements which are interwoven and consist only of directions for actions which offer a wild scope of autonomy for teachers and students. Overcoming such complexity causes both an enthusiastic involvement for few teachers and a lack of confidence for a large amount of others. In order to disseminate IBST we need to make teachers more confident in their own capacity to be open to changes in the direction of science teaching methods. This paper addresses a methodology for identifying and assessing IBST methods. It consists in elaborating a renewed research framework which considers the IBST genuine complexity. It aims to allow teachers and teacher educators to situate their activities with respect to an accurate model. The point is to enhance actors' confidence by providing them with a model which is accurate and could support teaching improvement.

A 6 dimensions model describing IBST methods

The 4 stages of professional development

Research about professional development shows that it consists of progressive and repeated reorganisations of the conceptualisations and practices, of the ways to think about activity and to act professionally [START_REF] Boreham | A Socio-cultural analysis of Organisational Learning[END_REF][START_REF] Engeström | Expansive learning at work: toward an activity theoretical reconceptualization[END_REF][START_REF] Grangeat | Factors influencing teachers' professional competence development[END_REF]. Four stages are often identified as a continuum which underlines professional development: novice, advanced beginner, proficient performer, and expert (H. L. [START_REF] Dreyfus | Mind over Machine: The Power of Human Intuition and Expertise in the Era of the Computer[END_REF].

Michel Grangeat ECER 2011 2

During the two first stages, actors endeavour to elaborate meaningful units which gather and synthetize the different elements resulting from previous experience, discussions with other actors and knowledge resulting from education or training programs. Such units are identified as work process knowledge [WPK] within the industrial field by [START_REF] Boreham | Work process knowledge in technological and organizational development[END_REF]. Such WPK organizes and orders 4 elements: actors' objects, clues picked up from the situation in order to diagnose the case, set of actions which results from this diagnostic, and professional knowledge used as reference to justify these actions [START_REF] Grangeat | Effects of Expert Teachers' Collaboration on their Conceptualisations and practices towards Inquiry Based Methods in Science Teaching[END_REF]. These set of 4 elements could be called 'teacher work process-knowledge' [TWPK]. These actors are centred on their own activity.

During the two following stages, actors are able to match situations which appear similar and necessitate the same way to a ct. Thus, they can beneficiate of previous professional knowledge in order to achieve more challenging tasks and to cope with variations occurring amongst these. Consequently they reinforce the set of their TWPK by gathering similarities and distinguishing specificities through new wide cognitive units we could call 'teacher activity system' [TAS] following [START_REF] Engeström | Activity theory as a framework for analyzing and redesigning work[END_REF]. These actors are centered both on pupils and colleagues activity; they connect their own conceptualizations and practices both with learners' ways to act in school and with their colleagues' teaching approaches.

Consequently, the following model must include the main dimensions of IBST. On each dimension it needs to identify 4 stages: the first two stages would be content and teacher centred and the last ones would be learning focused. We assume that the latter approach includes the former: the learners-centred approach is more sophisticated and complex than the teacher-centred [START_REF] Grangeat | Teachers' knowledge: a synthesis between personal goals, collective culture and conceptual knowledge[END_REF][START_REF] Hudson | Comparing different traditions of teaching and learning: what can we learn about teaching and learning[END_REF].

An IBST 6 dimensions model

A 6 dimensions model results from the IBST state of the art conducted by S-TEAM project.

Origin of questioning

The first dimension is about the origin of questioning which triggers the inquiry. On the first stage, teachers could elaborate the questioning on their own and manage to transform it into a learners questioning. On the last stage this questioning could be elaborate by pupils' group.

Nature of problem

The second dimension tackles the nature of the problem which supports and orientates the investigation. On one hand, the problem is closed and pupils have to follow a narrow protocol. On the other hand, the problem is open-ended and the learners need to elaborate their own hypothesis and protocol.

Promoting pupils' responsibility within the inquiry process

The third dimension addresses the responsibility of the leaners to carry out the inquiry process. On the first stages, teachers steer quite narrowly the investigation, for example by leading all the pupils' team towards a same experimentation. On the last, teachers manage to let a part of their responsibility to the learners, for example through providing them self-assessment documents.

Considering positively pupils' diversity

This dimension addresses the way teachers take into account of and cope with pupils' diversity of knowledge, needs, and motivation with respect of science learning. With respect of the current model, the first stages which are teacher-centred aim to cope with specific pupils' behaviour in order to involve them within the inquiry. The following stages learners-centred consist in adapting teaching strategy to the learners specificity.

Development of argumentation

The fifth dimension tackles the argumentation development. The first stages aim to facilitate discussion amongst pupils within the teams. The following stages aim to enable pupils to justify their conclusion with respect of experimental or database searching results.

Explanation of teachers' goals and strategies

This sixth and last dimension aims to make explicit the teacher's goals and specifically the learning outcomes which are expecting through IBST. The first stages consist in communicating to the students the

Michel Grangeat

ECER 2011 3 teacher's expectations for the current lesson. The last stages aim to make explicit knowledge and metaknowledge which result from the current session and will be useful within further situation and problem.

Understanding IBST through a 6 dimensions model

This model could allow us to describe and compare teachers' practices and approaches. Nevertheless, this paper's main purpose consists in overcoming IBST complexity through a model both referring to current research about science education. This 6 dimensions model endeavours to achieve these goals in order to able teachers to be more confident facing the changes resulting from IBST development within class-room.

Comparing teachers' activities through the IBST 6 dimensions model

Two groups of science teachers

In order to understand the effects of CPD 14 science teachers [ST] split in two groups are compared. They teach mathematics, biology, physics and chemistry, or technology.

The first group gathers 8 proficient science teachers [PST] who are involved within CPD program through two activities: they frequently meet together and with local authorities in order both to improve their teaching approaches and practices, and to design and to carry out CPD program about IBST. The second group gathers 4 ordinary science teachers [OST] who could only apply on their own experience since they are involved neither within professional networks, nor within CPD program about IBST.

Data collecting

The data are collected through videotaped IBST sessions and interviews with each ST about the video. Firstly the research team ask each ST to carry out a lesson which he or she considers as IBST session. The lesson lasts for about 55 minutes. Afterwards, each ST is interviewed about the 20 last minutes of the video. The ST is asked to stop the video when an event which had conducted him or her to choose amongst different alternatives occurs. Thus, the ST is asked to make explicit both the event which challenged her teaching strategy, and the goals which steer the observed action; mostly, teachers explain professional knowledge which had underpinned their choices. The interviewer could also stop the video and ask questions. Such a methodology makes reachable the set of TWPK of each ST. Each interview lasts about 60 minutes. All lessons and interviews are transcribed and analysed in order to find out the 4 elements of TWPK which define the conceptualizations and practices system of each ST. This TAS is built on the IBST 6 dimensions model.

Expected findings

Two results are expected. Firstly, proficient science teachers [PST] reach the upper levels of IBST 6 dimensions model. Secondly, PSTs who collaborate through professional networks and attend to specific CPD program about IBST reach upper levels than OSTs which can apply only on their own experience.

Results: teacher collaboration lead to learning focused teaching approaches

Proficient teachers involved within teacher collaboration through CPD

For all the dimensions of the model, the 8 PSTs reach at least level 2 of the continuum. All PSTs [8/8] reach levels 3 or 4 for half of the 6 dimensions and most of PSTs [5/8] for a least 4 dimensions. The lower mean [2.38/4] occurs for dimension 2 (nature of the problem). The higher mean [3.50/4] occurs for dimensions 4 (pupils diversity) and 5 (argumentation). Thus they are quite learner-centred.

Ordinary teachers applying only on their experience

All OSTs [4/4] are ranking at level 1 for at least 1 dimension of the model, mostly for dimension 1 (origin of questioning). Half OSTs [2/4] reach levels 3, which are learner centred, for half of the 6 dimensions. The other one [2/4] reach this level for 2 dimensions which are either dimension 4 (diversity), 5 (argumentation), or 6 (explanation). The lower means occurs for dimension 1 (origin of questioning) [1.50/4] and dimension 2 (nature of the problem) [1.75/4]. The higher means are achieved by dimension 4 (pupils diversity) [3.25/4] and dimension 3 (pupils responsibility) [3.00/4]. Thus they are globally teacher-centred.

Mastering learned-centred approaches and practices through IBST development

The mean of PSTs group is up than 2 for all the 6 dimensions contrarily of OSTs' which reach this threshold for only 3 dimensions. Thus results confirm what was expected: CPD based on teacher collaboration could contribute to improve the mastering of IBST complexity. Nevertheless OSTs are learner-centred for some Michel Grangeat ECER 2011 4 crucial IBST dimensions: enhancing pupils' responsibility, considering positively their diversity, and promoting argumentation within inquiry. Thus, IBST by itself could promote more autonomous learning. The graph which synthetized these results (see Figure 1) could make teachers, teacher educators or inspectors able to locate their own activity or an observed lesson. It could be an instrument used by actors to self-assess their IBST approaches and practices and to improve them, following specific directions and goals according to their own situation, needs and wills. The fact that teachers don't need to act at each upper level of each dimension could support their confidence in developing IBST. This was a goal of this paper.

Proficient science teachers' activity system

Within this paper length, it is impossible to describe the TAS of each ST. Thus this paper focuses of the activity system of PSTs' group with respect of the fourth dimension (diversity) of the IBST model. This could promote this aspect of IBST towards new and ordinary teachers.

Here is the main teaching work process knowledge [TWPK] for the 4 levels of this dimension.

3.2.1.a. Level 1: Coping with some pupils' behaviour in order to involve them within inquiry

All the PSTs act at least at this first level [8/8]. Their goals range from 2 to 6 (mean=4). These 4 goals are [within bracket the number of PSTs who aim this goal]:

 To support a pupil who encounters difficulties and to help him or her to enter within the inquiry process The TWPK which corresponds to the goal which is more shared is (see Table 1):

Goal

To support a pupil who encounters difficulties and to help him or her to enter within the inquiry process Clue Some pupils are bogged down Actions I restore their confidence I support them really Reference knowledge I think it is crucial to begin by restoring their confidence: during peer-work they will understand that they are allowed and capable to propose their own questions and solutions. After that, learning and teaching could run in easiest ways. The TWPK which corresponds to the more shared goal is (see Table 2):

Goal To maintain the activity of pupils who have completed the task before the others Clue Some pupils are very advanced and there is a gap with other teams. Actions I propose another activity which complements the first. Reference knowledge I ever endeavour to prepare complementary activities for the teams who will achieve the inquiry before the others. The TWPK which corresponds to the more shared goal is (see Table 3):

Goal

To create pupils' teams.

Clue

When I want to promote exchange and discussion amongst pupils Actions I let pupils create the teams with their preferred mates.

Reference knowledge

The class-room is noisier with this way to create teams but the inquiry runs well. When I create teams by myself, it doesn't work as good. One of the common TWPK is (see Table 4):

Goal Letting low achiever pupils express their opinions Clue When a low achiever pupil is trying to express his or her opinion Actions I give time to finish his or her explanation Reference knowledge This kid is boring when other pupils are speaking, and the others are boring when this kid is speaking… But during IBST it's good to make an effort for this kid could express opinions.

Conclusion

The IBST 6 dimensions model drawn from S-TEAM products and deliverables appears as a relevant mean to understand and assess science teachers' conceptualizations and practices when they carry out an inquirybased lesson. The study had compared 8 teachers involved within specific CPD program with 4 ordinary teachers. It shows that teachers involved within CPD programs which promote teacher collaboration through exchange and controversy amongst attenders could achieve the highest levels of the IBST 6 dimensions model, which correspond to learning-focused approaches and practices.

Nevertheless, the fact that the sixth dimension of the model appears always amongst the weakest ask a crucial question: if teachers remain reluctant to make explicit the goals they want that pupils achieve through inquiry, we could think that most of these pupils cannot benefit of the inquiry to master more scientific knowledge. Here is a crucial challenge for improving IBST.

Figure 1 :

 1 Figure 1: Locating teachers' strategies on IBST 6 dimension model

  [7].  To choose the right class-room with which developing IBST session [6].  To gain pupils' attention and to support them to enter within the inquiry process [5].  To adapt IBST lesson to the global class-room level [2].

Table 1 :

 1 Supporting low achievers Half of PSTs acts at this second level [4/8]. They share two goals:  To maintain the activity of pupils who have completed the task before the others [3].  To give time to all pupils to complete the inquiry [2].

		Questioning
		4,00
		3,00
	Explanation	2,00	Problem
		1,00	PST
		0,00
			OST
	Argumentation		Responsability
		Diversity

Table 2 :

 2 Proposing complementary activity for advanced teams 3.2.1.c. Level 3: Assessing the problem understanding within each pupils' team All the PSTs act at this third level [8/8]. Their goals range from 1 to 3 (mean=2). Most frequent goals are:  To create pupils' teams [4].  To help the teams which encounter difficulties [3].  To choose the team which need teacher' help [3].  To identify the product of each pupils' team [3].

Table 3 :

 3 Creating pupils' teams 3.2.1.d. Level 4: Adapting the situation to specific pupils' particularities.

	Half of PSTs acts at this fourth level [4/8]. They have 1 or 2 goals:
	 Letting low achiever pupils express their opinions [2].
	 Adapting the situation to disabled pupils [2].

Table 4 :

 4 Letting low achiever pupils express their opinions

Mind The Gap and S-TEAM, for instance

Michel Grangeat

ECER 2011 6


Ivan Breslavsky

Marco Amabili
email: marco.amabili@mcgill.ca

Mathias Legrand

Nonlinear vibrations of thin hyperelastic plates

Keywords:

Static deflection as well as free and forced nonlinear vibration of thin square plates made of hyperelastic materials are investigated. Two types of materials, namely rubber and soft biological tissues, are considered. The involved physical nonlinearities are described through the Neo-Hookean, Mooney-Rivlin, and Ogden hyperelastic laws; geometrical nonlinearities are modeled by the Novozhilov nonlinear shell theory. Dynamic local models are first built in the vicinity of a static configuration of interest that has been previously calculated. This gives rise to the approximation of the plate's behavior in the form of a system of ordinary differential equations with quadratic and cubic nonlinear terms in displacement. Numerical results are compared and validated in the static case via a commercial finite element software package: they are found to be accurate for deflections reaching 100 times the thickness of the plate. The frequency shift between low-and large-amplitude vibrations weakens with an increased initial deflection.

Introduction

Thin-walled structures made of hyperelastic materials, such as rubbers and biomaterials, are common in mechanical and biomechanical engineering applications. Such structures are frequently subjected to significant static and dynamic loadings potentially yielding large deflections and deformations. Associated previous works reported in the literature employ, in a vast majority, a simplifying assumption that considers as a priori known the shape of the deformed structure; for instance, a cylindrical shell remains cylindrical after deformation. A general review on these works is available in [START_REF] Ivan | Physically and geometrically non-linear vibrations of thin rectangular plates[END_REF]. The present investigation discards this limitation.

One possible way to account for complicated deformed shapes is to use the finite element (FE) approach. For instance, Einstein [START_REF] Daniel | Dynamic finite element implementation of nonlinear, anisotropic hyperelastic biological membranes[END_REF] developed a finite element strategy to analyze the dynamic behavior of hyperelastic and viscoelastic membranes with amplitudes of the order of their thickness. With the proposed method, the dynamic response of a hemisphere made of a Mooney-Rivlin material under impact pressure load is studied. A relatively similar problem concerning the dynamic inflation of a Mooney-Rivlin spherical membrane is explored in [START_REF] Erwan | Dynamic inflation of non-linear elastic and viscoelastic rubber-like membranes[END_REF] also through the FE method.

Another technique lies in the approximation of the deformed configuration as a truncated series of continuous functions satisfying the essential geometric boundary conditions. This approach is implemented for circular membranes in [START_REF] Paulo | Nonlinear vibrations of a radially stretched circular hyperelastic membrane[END_REF][START_REF] Soares | Nonlinear vibrations and instabilities of a stretched hyperelastic annular membrane[END_REF]. The sensitivity of the respective vibratory behavior to the pre-stretch configuration of the membranes is carried out in these scientific contributions. The incompressible hyperelastic material is either described by a Neo-Hookean model [START_REF] Soares | Nonlinear vibrations and instabilities of a stretched hyperelastic annular membrane[END_REF] or by Mooney-Rivlin, Yeoh, Ogden, and Arruda-Boyce models in [START_REF] Paulo | Nonlinear vibrations of a radially stretched circular hyperelastic membrane[END_REF]. In this latter work, the frequencyamplitude relationships for various hyperelastic models are found to be analogous. The authors come to the conclusion that the higher the stretch, the closer the nonlinear forced response to its linear counterpart.

In [START_REF] Ivan | Physically and geometrically non-linear vibrations of thin rectangular plates[END_REF], the static deflection and vibration around a deformed configuration of a rectangular Neo-Hookean plate are investigated accounting for both geometrical (the nonlinearity of straindisplacements relations) and physical nonlinearities (the material nonlinearity of stress-strain relations). Is proposed a method which systematically builds approximate local models (LMM) in the form of polynomial expansions of the non-polynomial strain energy densities: this greatly simplifies the description of the plate's dynamic behavior. The present paper extends this investigation to various hyperelastic laws by using a more sophisticated nonlinear plate theory. The two targeted types of material (rubber and biomaterial) are described through (i) Neo-Hookean, (ii) Mooney-Rivlin, and (iii) Ogden hyperelastic laws as well as (iv) a physically linear material law. The comparison of static force-deflection curves and frequency responses obtained with these models is carried out and the contribution of in-plane nonlinearities in the sought displacement is estimated. The numerical solutions are verified through direct comparison with results from a commercial FE software package.

Kinematics and energy equations

The geometrical nonlinearity in the kinematics of the investigated flexible plate is described with the help of the Novozhilov nonlinear plate theory [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF], which stands as a limit case of the Novozhilov nonlinear shell theory. It is governed by the following strain-displacement relationships:

" 1 D u x C 1 2 .w 2 x C OEu 2 x C v 2 x / zw xx " 2 D v y C 1 2 .w 2 y C OEu 2 y C v 2 y / zw yy " 12 D u y C v x C w x w y C OEu x u y C v x v y 2zw xy (1)
where " 1 , " 2 , and " 12 are the components of the Green-Lagrange strain tensor for thin plates. The Novozhilov nonlinear plate theory, commonly recognized as the finest classical plate theory, reduces to the well-known von Kármán nonlinear plate theory [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF] when the in-plane terms in square brackets in (1) are neglected.

Although it neglects shear deformations and rotary inertia, it is accurate for thin plates. Also, it is known that in-plane nonlinearities play a major role in large deflections, and they are accordingly retained in the present work to improve accuracy. Lagrange equations are utilized to derive the dynamics of the plate, that is: where L D T … is Lagrange's functional, T is the kinetic energy of the plate, … is the potential elastic deformation energy, and Q n are the generalized forces. The potential and kinetic energies are expressed as follows [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF]:

d dt @L @ P q n C @L @q n D Q n ; n D 1; : : : ; N (2)
… D • V W dV (3)
T D h 2 " S . P u 2 C P v 2 C P w 2 / dS ( 4
)
where W is the strain energy density (SED), V is the volume of the plate, S is the surface of the middle plane of the plate, is the mass-density of the plate material, h is the thickness of the plate, and u; v; w are the displacements along the axes of the rectangular coordinate system x; y; z, respectively. The dot stands for differentiation with respect to time.

The displacements are expanded into corresponding truncated series involving the appropriate generalized coordinates q n :

w.x; y; t/ D N W X i D1 q i .t/W i .x; y/ u.x; y; t/ D N U X i D1 q i CN W .t/U i .x; y/ v.x; y; t/ D N V X i D1 q iCN U CN W .t/V i .x; y/ (5)
In Eq. ( 5), quantities W i ; U i ; V i are the admissible functions that satisfy the homogeneous boundary conditions (i.e. the geometric constraints) of the problem. The linear modes of vibration, which form a complete set properly capturing the dynamics of a structure, are eligible admissible functions and are selected in the present work. The total number of degrees of freedom is

N D N W C N U C N V .

Hyperelastic relations

Strain energy density

Usually, nonlinear elasticity of rubbers and soft biomaterials is described by hyperelastic laws and in most cases, such materials are assumed to be incompressible [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Bower | Applied Mechanics of Solids[END_REF]. Three hyperelastic laws together with their associated SED W are considered: Neo-Hookean:

W D E 6 .I 1 3/ (6)
Mooney-Rivlin:

W D 1 2 .I 1 3/ C 2 2 .I 2 3/ (7)
Ogden:

W D N T X iD1 i ˛i . ˛i 1 C ˛i 2 C ˛i 3 3/ (8)
Also, the strain energy density for a linear material reads:

W D E 6 4." 2 1 C " 2 2 C " 1 " 2 / C " 2 12 (9)
The following notation is used: I 1 is the first invariant of the right Cauchy-Green deformation tensor C; E is Young's modulus of the plate's material; I 2 stands for the second invariant of the right Cauchy-Green deformation tensor; 1 ; 2 ; 3 are the principal stretches of the plate; i ; ˛i denote the material parameters. Since SED ( 9) is polynomial in strains, spatial and temporal components of the solution can be uncoupled and the Lagrange equations reduce to simple ordinary differential equations with quadratic and cubic nonlinearities [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF]:

R q n C 2 n n P q n C N X i D1 k ni q i C N X i;j D1 k nij q i q j C N X i;j;`D1 k nij `qi q j q `D Q n ; n D 1; : : : ; N ( 10
)
where n is the natural frequency of mode n and n is the corresponding damping ratio; k ni ; k nij ; k nij `are known coefficients that result from integration in space.

2

Cauchy-Green tensor invariants and principal stretches

In order to derive the expressions of the invariants of the right Cauchy-Green strain tensor C in terms of displacements, the Green-Lagrange strain tensor:

E D 1 2 0 @ 2" 1 " 12 0 " 12 2" 2 0 0 0 2" 3 1 A (11)
is used, where the expressions of " 1 , " 2 , and " 12 (but not " 3 ) are given in [START_REF] Ivan | Physically and geometrically non-linear vibrations of thin rectangular plates[END_REF]. The right Cauchy-Green deformation tensor C is then defined as [START_REF] Bower | Applied Mechanics of Solids[END_REF]:

C D 2E C I D 0 @ 2" 1 C 1 " 12 0 " 12 2" 2 C 1 0 0 0 2" 3 C 1 1 A (12)
and its three invariants are:

I 1 D tr C D 2." 1 C " 2 C " 3 / C 3 (13)
I 2 D 1 2 .tr C 2 tr.C 2 // D 4." 1 C " 2 C " 3 C " 1 " 2 C " 1 " 3 C " 2 " 3 / " 2 12 C 3 (14) J 2 D jCj D .2" 3 C 1/ .2" 1 C 1/.2" 2 C 1/ " 2 12 ( 15
)
The third invariant J is used to reflect the incompressibility condition through J D 1 [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF] and the principal stretches are the square roots of the eigenvalues of C [8]:

1;2 D r 1 C " 1 C " 2 ˙q." 1 " 2 / 2 C " 2 12 ( 16
) 3 D p 1 C 2" 3 (17)

Transverse normal strain

The transverse normal strain " 3 expressed in terms of u, v, and w shall be inserted into expressions [START_REF] Eusebius | AUTO 97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont)[END_REF] to (17) and the incompressibility condition ( 15) is enforced to retrieve the expression of " 3 :

" 3 D 1 2 .2" 1 C 1/.2" 2 C 1/ " 2 12 1 2 (18)
Expression ( 18) is then introduced in ( 13)-( 15), (17) as well as in SEDs ( 6)-( 8).

Strain energy density local expansion

Expressions ( 6)-( 8) together with " 3 from (18) are not polynomials in strains, which essentially complicates the investigation of the plate's behavior. The analysis is thus simplified by introducing a transformation of SEDs ( 6)-( 8) in order to derive approximate governing equations in the form of ordinary differential equations with nonlinearities of order not higher than three.

The corresponding local model is reliable only in the vicinity of a configuration of interest around which the SED is expanded into a series in the generalized coordinates truncated at order 4.

To reach highly deformed configurations, successive local models have to be constructed. Further details on the local models method (LMM) can be found in [START_REF] Ivan | Physically and geometrically non-linear vibrations of thin rectangular plates[END_REF].

Exact low-dimensional models with both material and geometric nonlinearities

A distinct additional approach is implemented to measure the accuracy of the LMM. It involves the numerical solution of Equa-tion (2) in its static version, that is:

• V @W @q n dV D Q n ; n D 1; : : : ; N (19)
with SEDs ( 6)-( 8) and expression (18) previously substituted into it. Since the latter is a smooth function in the generalized coordinates, integration and differentiation operators commute, that is:

@ @q n • V W dV D • V @W @q n dV (20)
System (19) of nonlinear algebraic equations can be solved numerically through the Newton-Raphson iterative technique; this approach is named the "exact solution" hereinafter.

Numerical example: static and dynamic bending of a rubber plate 6.1 Hyperelastic models parameter identification

Experimental data for 8% sulfur rubber obtained by Treloar [START_REF] Leslie | Stress-strain data for vulcanised rubber under various types of deformation[END_REF] for uniaxial and equibiaxial tensions, as well as for pure shear, are exploited. This common rubber is chosen because experimental data for uniaxial and multiaxial loads is available in the literature: a good approximation for the strain energy density is then possible for this material. Corresponding details on the approximation of the experimental data are provided in Appendix A. From the given experimental data, Young's modulus is E D 1 247 060:2 Pa, which is in almost perfect agreement with Ogden's one, E D 1 242 992:9 Pa [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF]. Other relevant parameters are listed in Table 1. ,2, and 3 display the strain-stress relationships for the three hyperelastic laws together with the experimental points. The nominal stress S 1 versus engineering strain " is shown-see Appendix A for details. It is worthy to note that the Neo-Hookean and Mooney-Rivlin laws are very close and both laws satisfactorily approximate the rubber behavior at strains not higher than 30% (that is 0:3 in Figures 1 to 3). However, the approximation of the rubber behavior offered by the Ogden law is much more accurate on the full range of experimental strains.

Model Parameters

Neo-Hookean

Problem description

A simply supported square rubber plate illustrated in Fig. 4 is considered. It is defined on the following domain:

V D ˚x 2 OE0; a; y 2 OE0; b; z 2 OE h=2; h=2 « (21)
with a D 0:1 m, b D 0:1 m, h D 0:0005 m. The material characteristics are listed in Table 1 and the rubber mass density is D 1100 kg m 3 . The plate is simply supported with immovable edges yielding the following boundary conditions [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF]:

wj @S D M j @S D uj @S D vj @S D 0 ( 22
)
where @S denotes the boundary of the plate's middle surface. The bending moment per unit length M [6] reads:

M D D @ 2 w @n 2 C @ 2 w @ 2 (23)
and quantities n and are the outer normal and tangent directions to @S , respectively, as shown in Fig. The problem is made non-dimensional by introducing a nondimensional time D 1 t , where 1 is the circular frequency of the first natural mode of the deformed (pressure loaded) plate. Also, a change of notation has been introduced; the two-subscript generalized coordinates (time functions) w 2nC1;2mC1 , u 2n;2mC1 , and v 2nC1;2m are divided by the plate thickness h and replaced by the single-subscript generalized coordinates q i with corresponding increasing numbering.

Static analysis

First, the problem of static plate bending under uniform pressure is explored. As illustrated in Figure 5, a quick convergence analysis shows that a 12 degree-of-freedom (DOF) model stands as a convincing compromise between prediction capabilities and computational cost for the force-deflection curve of the Neo-Hookean plate. The generalized coordinates plugged into the expansion of the displacements are given in Table 2. Responses are calculated with the LMM, starting from the configuration with w 1;1 D 20h captured by the less sophisticated model involving geometrical nonlinearity only. Accordingly, the curves start at a central deformation of 20h in Figure 5. The maximal discrepancy between the 12 DOF and 27 DOF models in the deflection range OE0; 100h is 2.3%. Similar results are found for other material models.

Size Participating eigenmodes 3 DOF w 1;1 , u 2;1 , v 1;2 12 DOF w i;j i; j D 1; 3 u i;j , v j;i i D 2; 4, j D 1; 3 27 DOF w i;j i; j D 1; 3; 5 u i;j , v j;i i D 2; 4; 6, j D 1; 3; 5 Force-deflection relationships for different materials with 12 DOF are depicted in Figure 6. The displacement sensitivity to physical nonlinearities is highlighted by plotting the solution for linear elastic material and geometric nonlinearities only. Neo-Hookean and Mooney-Rivlin results almost coincide, while the Ogden curve slightly deviates for large deflections. Also, it is worthy to note that the LMM results agree very well with the exact solution.

Material nonlinearity for deflections smaller than 20h can be neglected as shown in Figure 6. For deflections up to 20h the model with only geometrical nonlinearity is a very good approximation. As a consequence, the configuration computed with the linear elastic material is used as an initial guess for LMM numerical iterations for hyperelastic materials [START_REF] Ivan | Physically and geometrically non-linear vibrations of thin rectangular plates[END_REF]. An additional validation of the results is carried out. The static deflection of the hyperelastic plate is also explored with the commercial finite element solution ANSYS [START_REF] Saeed | Finite Element Analysis: Theory and Application with ANSYS[END_REF] for the Neo-Hookean and Ogden materials. The central deflection with respect to the applied external force is displayed in Figure 7: FEA and LMM predictions are in good agreement.

The contribution of the in-plane nonlinearities in expressions ( 1) is now estimated. As an illustrative example, the pressure-deflection curves for the Neo-Hookean law with and without in-plane nonlinearities are shown in Figure 8. The maximal difference is 2%. Other materials feature discrepancies of the same order. Despite of their relatively small influence, in-plane nonlinearities are incorporated in all numerical investigations except in Figure 8. ] and without [ ] in-plane nonlinearities

Nonlinear vibration analysis

At small strains, the effect of physical nonlinearities is negligible (Figure 6), and the vibrations of the initially flat plate can be explored by only retaining geometrical nonlinearities in the model [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF]. Here, we focus on the investigation of a more challenging configuration where both nonlinearities are participating, i.e. the vibrations around a pre-loaded state. The deflection with principal generalized coordinate w 1;1 D 80h is chosen as an initial deformed configuration. Comparison with the exact static solution shows that the local model around this deformed configuration is accurate for deflections up to 10h, which therefore stands as an upper limit in the dynamic analysis. From previous analysis (Figure 5), it follows that the 12-DOF model is sufficiently accurate and can be selected in order to reduce computation costs. The harmonic balance method [START_REF] Thomas | Practical Numerical Algorithms for Chaotic Systems[END_REF] is implemented to find the sought periodic solutions to system (10) through a Fourier expansion of the generalized coordinates in time:

q n D A n0 C N h X j D1
A nj cos.j /; n D 1; : : : ; N

where is the non-dimensional frequency, normalized with respect to 1 . Expression (27) considers cosine terms only since the external force and damping in [START_REF] Ivan | Two modes nonresonant interaction for rectangular plate with geometrical nonlinearity[END_REF] are ignored in the free vibration analysis. Coefficients A nj are determined from the system of nonlinear algebraic equations resulting from balancing the coefficients associated to the same harmonics in equations [START_REF] Ivan | Two modes nonresonant interaction for rectangular plate with geometrical nonlinearity[END_REF] where expression (27) has been previously substituted. A conver- the deformed plate are given in Table 3 for an initial deflection w 1;1 D 80h. Figure 9 displays the backbone curves of the plate's free vibrations, with frequencies close to the first eigenfrequency of the pre-loaded plate and normalized with respect to the natural frequency 1 of the corresponding deflected plate. Comparison with the exact static solution shows that the Neo-Hookean and Mooney-Rivlin local models are accurate for deflections up to 10h, while the Ogden model can be considered for deflections not higher than 8h only. Again, the Neo-Hookean and Mooney-Rivlin models yield identical results. This is true for both the eigenfrequency of the deformed plate and the backbone curve. The Ogden backbone curve exhibits a slightly weaker nonlinearity. For all materials, the nonlinear nature of the deformed plate, as opposed to its flat counterpart [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF], is very weak: the frequency of large-amplitude vibrations is almost identical to the natural frequency. In order to show the sensitivity of the backbone curves to the initial deflection, they are constructed around various initial configurations in Figure 10. The Neo-Hookean material is selected since other materials exhibit similar features. Two effects are observed: the first one consists of an increased range of vibra- tory amplitudes where the backbone curve displays a softening behavior. This effect is well-known and might be attributed to the contribution of the quadratic terms in expression [START_REF] Ivan | Two modes nonresonant interaction for rectangular plate with geometrical nonlinearity[END_REF] with respect to the initial deflection amplitude [START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF]. The second effect is the attenuation of the nonlinear nature of the model with an increased initial deflection. A similar behavior has been previously reported for bended plates [START_REF] Ivan | Physically and geometrically non-linear vibrations of thin rectangular plates[END_REF] and stretched membranes [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF] and is related to the large in-plane stretching associated to the initial deflection of the loaded plate. The forced vibrations are also explored with the AUTO package [START_REF] Eusebius | AUTO 97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont)[END_REF]. External forcing derives from a time-dependent periodic pressure having an harmonic component P d D 4:53 Pa and a constant mean-value corresponding to the deformed configuration w 1;1 D 80h, around which the local model is built. A modal damping ratio n D D 0:001 is adopted. The frequency response for the principal bending coordinate w 1;1 is shown in Figure 11. The vibratory response around a highly deformed configuration is almost identical to the linear response, in contrast to the response around a moderately deflected plate for which a significant softening behavior turning to hardening for vibration amplitude around 5h is found as shown in Figure 12.

Static and dynamic plate bending for biologicaltype material 7.1 Hyperelastic models parameters

The experimental data for uniaxial tension available in [START_REF] Gerhard | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF] are exploited. The present study is intended to model only one key feature of biological materials, i.e. a sharp increase in stiffness after a given strain threshold is reached. To this end, the experimental data corresponding to the tunica adventitia of a human aorta in [START_REF] Gerhard | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF] are selected. The procedure described in Appendix A is implemented with the difference that only the uniaxial test data are available. Corresponding material parameters are listed in Table 4. The associated stress-strain relationships are shown in Figures 13 to 15. The difference between the Neo-Hookean and Mooney-Rivlin laws is clearly distinguishable for this type of material. However, both laws are in agreement with the experimental points for strains less than 8% and both are unable to reproduce the increase in stiffness, as opposed to the Ogden law.

Static analysis

The problem detailed in subsection 6.2 is solved for the hyperelastic materials listed in Table 4. Attention is paid to the static bending of the plate (Figure 4) with the help of the 12 DOF model (Table 2). The exact pressure-deflection response is depicted in Figure 16. The Neo-Hookean and Mooney-Rivlin results are almost identical, while the Ogden model starts to depart from these curves for amplitudes about 80h, as displayed in Figure 16, where it becomes much stiffer. Strains at the point where the Ogden and the Neo-Hookean curves cross each other in Figure 16 are also analyzed through the strain intensity, given by [START_REF] Lazar | Fundamentals of Theory of Plasticity[END_REF]: At the crossing point in the center of the plate, the strain intensity is " INT D 42:5% (Figure 16) near the values of the crossing points in the stress-strain curves: 40% for uniaxial test and 39.6% for equibiaxial test (Figures 13 and14).

" INT D r 2 3 r ." 1 " 2 / 2 C ." 1 " 3 / 2 C ." 2 " 3 / 2 C 3 2 " 12 (28)

Vibration analysis

We investigate the free nonlinear vibrations around the deformed response with principal generalized coordinate w plate is listed in Table 5 for the three hyperelastic laws. The mass density is D 1380 kg m 3 , which is the density of Polyethylene terephthalate, frequently used for artificial arteries [START_REF] Jaroslav | Blood vessel replacement: 50 years of development and tissue engineering paradigms in vascular surgery[END_REF]. We notice that the Neo-Hookean and Mooney-Rivlin results are again in good agreement for the full investigated range of deflections. The Ogden plate's eigenfrequencies are also similar for central deflections up to 80h but then become higher.

Similarly to the rubber plate example, the behavior of a deformed plate made of biomaterial is weakly nonlinear for all materials and the nonlinearity is of softening type. However, certain features differ. Although static results for the Neo-Hookean and Mooney-Rivlin models are almost identical, the Mooney-Rivlin backbone curve shows a softer behavior. Still, the Ogden curve is the softest one and is presented for vibration amplitude up to 3h in Figure 17.

Conclusions

Vibration of plates made of rubber and biological materials are explored with a dedicated local models method. Local models allow for the analysis of static bending as well as physically and geometrically large-amplitude nonlinear vibrations of an initially distorted plate. Static results are validated through a systematic comparison with available exact solutions as well as against commercial finite element software results. It is found that the local models method provides accurate predictions for a wide range of deflections. However, the sharp increase in stiffness peculiar to biological materials limits the range of achievable vibration amplitudes.

In most cases, it is found that the Mooney-Rivlin and Neo-Hookean materials exhibit similar static and vibratory behaviors. Corresponding constitutive laws properly capture the behavior of the actual material at moderate strains (30% for rubber and 8% for biomaterial). The best approximation is provided by the Ogden model. The latter correctly reproduces the behavior at high strains, including the well-known sharp increase in stiffness. However, the Ogden model has significant drawbacks. Due to the complicated form of the strain energy density and its formulation in terms of principal strains (unlike the Neo-Hookean and Mooney-Rivlin models which allow for a formulation in terms of strain invariants), it is computationally much more expensive. It is also shown that the pre-loaded plate exhibits very weak dynamic nonlinearity, i.e. the frequencies of the oscillations around the deformed configuration are almost identical to the associated eigenfrequencies.

The sensitivity of the free vibration backbone curves to the initial deflection is also discussed. It is shown that the higher the initial static deflection, the higher the range of amplitudes at which the backbone curve displays a softening behavior. The nonlinear nature of the system of interest (the frequency shift between low-and large-amplitude vibrations) weakens with an increased initial deflection.

A. Hyperelastic material parameter identification

The following nominal stress expressions for the Neo-Hookean law are used in order to fit the experimental data (uniaxial, equibiaxial tensions and pure shear, respectively) [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Bower | Applied Mechanics of Solids[END_REF]: while for the Ogden law, they are [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF]: The procedure for the determination of the material parameters is as follows. First, Young's modulus is determined using the expressions for the Neo-Hookean law. For this purpose, the material is assumed to be Neo-Hookean and E is determined through a least squares technique by using only one point (corresponding to the smallest stress) for each stress-strain test. For the sake of consistency, the Young's modulus is assumed to be identical for all the hyperelastic laws. This yields 3. 1 C 2 / D E for the Mooney-Rivlin law and 3 2 P N T i D1 i ˛i D E for the Ogden law [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Bower | Applied Mechanics of Solids[END_REF]. Once E is known, the Neo-Hookean strain energy density is completely defined. The Mooney-Rivlin law contains one unknown parameter and the Ogden law involves 2N T 1 unknowns. These parameters are determined by the least squares fitting of expressions (32) to (37) to the experimental data.

S
S 1 D N T X iD1 i ." C 1/ ˛i 1 C ." C 1/ ˛i =2 1 (35) S 1 D N T X iD1 i ." C 1/ ˛i 1 C ." C 1/ 2˛i 1 (36)

  E D 1 247 060:2 Pa Mooney-Rivlin 1 D 416 185:5 Pa 2 D 498:8 Pa Ogden 1 D 161:2 Pa 2 D 1831 Pa 3 D 781 111 Pa ˛1 D 7:295 ˛2 D 2:729 ˛3 D 1:056

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Rubber uniaxial tension: strain-stress curves for various material laws and corresponding experimental points [ ]. Neo-Hookean [ ]; Mooney-Rivlin [ ]; Ogden [ ]

Figure 5 .

 5 Figure 5. Neo-Hookean plate pressure-deflection response. Deflection measured at the center of the plate and normalized with respect to the plate thickness h

Figure 6 .

 6 Figure 6. Pressure-displacement response. Local models as lines and exact results as symbols: Neo-Hookean [ , ]; Mooney-Rivlin [ , ]; Ogden [ , ]. Linear elastic material [ ]

Figure 7 .Figure 8 .

 78 Figure 7. Comparison between ANSYS and exact solutions. Exact Neo-Hookean [ ]; exact Ogden [ ]; ANSYS with Neo-Hookean model [ ]; ANSYS with Ogden model [ ]

Figure 10 .

 10 Figure 10. Plate's free vibration backbone curves in the vicinity of various initially deformed configurations; Neo-Hookean model

Figure 9 .

 9 Figure 9. Rubber plate's free vibration dimensionless backbone curves around deformed configuration w 1;1 D 80h for various material laws: Neo-Hookean [ ]; Mooney-Rivlin [ ]; Ogden [ ]

Figure 12 .

 12 Figure 12. Rubber plate's forced frequency response around deformed configuration w 1;1 D 5h; P d D 0:063 Pa, D 0:001. Stable [ ] and unstable [ ] solutions

Figure 13 .Figure 14 .

 1314 Figure 13. Biomaterial uniaxial tension: strain-stress curves for various material laws and corresponding experimental points OE . Neo-Hookean [ ]; Mooney-Rivlin [ ]; Ogden [ ]

Figure 15 .Figure 16 .

 1516 Figure 15. Biomaterial pure shear: strain-stress curves for various material laws. Neo-Hookean [ ]; Mooney-Rivlin [ ]; Ogden [ ]

Figure 17 .

 17 Figure 17. Biomaterial plate's free vibration backbone curves around deformed configuration w 1;1 D 70h for various material laws: Neo-Hookean [ ]; Mooney-Rivlin [ ]; Ogden [ ]

C 1 /

 1 ˛i 1 C ." C 1/ ˛i 1 (37)

Table 1 .

 1 Hyperelastic law parameters for rubberFigures 1

  Rectangular plate and coordinate system in this study; accordingly, due to straightforward considerations on the symmetry of the plate's geometry, on the external applied forces, and on the boundary conditions, it is legitimate to expand the displacements in series of sine[START_REF] Marco | Nonlinear Vibrations and Stability of Shells and Plates[END_REF][START_REF] Ivan | Two modes nonresonant interaction for rectangular plate with geometrical nonlinearity[END_REF]:

	z
					h
						b
						y
	a					n
	x
	n
	Figure 4. w.x; y; t/ D X n;m2N w 2nC1;2mC1 .t /sin	.2n C 1/x a	sin	.2m C 1/y b	(24)
	u.x; y; t/ D X n;m2N u 2n;2mC1 .t / sin	2nx a	sin	.2m C 1/y b	(25)
	v.x; y; t/ D X n;m2N v 2nC1;2m .t / sin	.2n C 1/x a	sin	2my b	(26)

[START_REF] Paulo | Nonlinear vibrations of a radially stretched circular hyperelastic membrane[END_REF]

. Static and dynamic deflections under a uniformly distributed pressure are targeted

Table 2 .

 2 Generalized coordinates utilized in the models

Table 3 .

 3 First eigenfrequency of the deformed rubber plate. Initial deformation w 1;1 D 80h

  1;1 D 70h. The comparison with the exact static solution shows that the Neo-Hookean and Mooney-Rivlin local models are accurate for deflections up to 10h whereas the Ogden model is limited to deflections not larger 3h only. The corresponding backbone curves are obtained through the harmonic balance method (27) and are depicted in Figure17for vibrations around the first eigenfrequency of the pre-loaded plate. The first eigenfrequency of the deformed

	Model Natural Frequency
	Neo-Hookean 183 rad/s
	Mooney-Rivlin 177.88 rad/s
	Ogden 181.96 rad/s

Table 5 .

 5 First eigenfrequency of the deformed biomaterial plate. Initial deformation w 1;1 D 70h

  For the Mooney-Rivlin law, these expressions become[START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Bower | Applied Mechanics of Solids[END_REF]:S 1 D 1 ." C 1/ ." C 1/ 2 C 2 1 ." C 1/ 3 (32) S 1 D 1 ."C1/ ."C1/ 5 C 2 ."C1/ 3 ."C1/ 3 (33) S 1 D . 1 C 2 / ." C 1/ ." C 1/ 3 (34)

	1 D	E 3	." C 1/ ." C 1/ 2	(29)
	S 1 D	E 3	." C 1/ ." C 1/ 5	(30)
	S 1 D	E 3	." C 1/ ." C 1/ 3	(31)

Acknowledgements

The authors acknowledge the financial support of NSERC Discovery Grant, Canada Research Chair, Canada Foundation for Innovation (LOF) programs of Canada and the PSR-SIIRI program of Québec.


Michel Broniatowski

Giorgio Celant
email: giorgio.celant@stat.unipd.it

Optimal extrapolation design for the Chebyshev regression

Keywords: Chebyshev system, optimal design, extrapolation design, Borel-Chebyshev Theorem

This paper introduces optiomal designs in the context of a regression model when the regression function is assumed to be generated by a Chebyshev system of functions. The criterion for optimality is the variance of a Gauss Markov estimator for an extrapolated value.

Introduction

This paper deals with a natural extension of the Hoel Levine optimal extrapolation design, as described in [START_REF] Hoel | [END_REF]. We recall that this classical design results as a consequence of the following fact.

A design is de…ned as a discrete probability measure on a set of measurements points x 0 ; ::; x g 1 which for notational convenience belong to the observable environmental set [ 1; 1] , denoting n i =n := (x i ) the frequency of replications of the experiment to be performed at point x i ; 0 i g 1, where the n i 's satisfy n 0 + :: + n g 1 = n: The points x i 's are the nodes of the design, and (x i ) is the so-called frequency of the design at node x i : Recall that the model writes

Y (x) = f (x) + (x)
for x in [ 1; 1] , the real valued function f is unknown but belongs to a speci…ed class of functions, and the random variable (x) is centered, with a …nite variance, in the present context. Observations are performed under the design, with the constraint n 1 + :: + n g 1 = n on the global budget of the design. Replications of the n i measurements Y j (x i ) , 1 j n i are independent. Independence also holds from node to node, which is to assume that all measurements errors due to the r.v.'s (x) are independent. The model is supposed to be homoscedasticity; hence the variance of (x) may not depend on x:

For a given c not in [ 1; 1] consider an estimate of f (c) with smallest variance among all unbiased estimators of f (c) which are linear functions of the observations Y j (x i ) , 1 j n i ;,0 i g 1; hence under a given design : An optimal design achieves the minimal variance among all such designs. This design is achieved by the Hoel Levine design when the function f is assumed to belong to the class of all polynomials de…ned on R with degree less or equal g 1, hence to the span of the class of monomials f1; x; ::; x g 1 g :

The main mathematical argument in order to obtain the Hoel Levine design lies in the solution of the following basic question: …nd a polynomial with equioscillations in g + 1 points in [ 1; 1] which assumes maximal absolute values all equal to 1 at those points. Up to a multiplicative constant such a polynomial results as the best polynomial approximation of the null function on [ 1; 1] by polynomials with degree g 1: Existence and uniqueness of this polynomial follows from the Borel-Chebyshev Theorem. We refer to [START_REF] Dzyadyk | [END_REF] for details and derivation of these results.

The aim is now to provide a larger context for similar questions, assuming that the function f may belong to some other functional class, still in a …nitely generated set of functions. De…nition 1 The system of functions (' 0 ; :::; ' g 1 ) in C (R) is a Chebyshev (or Haar) system on [ 1; 1] when 1) (' 0 ; :::; ' g 1 ) are linearly independent 2) Any equation

a 0 ' 0 (x) + ::: + a g 1 ' g 1 (x) = 0
with (a 0 ; ::; a g 1 ) 6 = (0; ::; 0) has at most g roots in [ 1; 1] :

Denote V := span f' 0 ; :::; ' g 1 g C ([ 1; 1])
the linear space generated by the Chebyshev system (' 0 ; :::; ' g 1 ) :

Haar Theorem (see [START_REF] Dzyadyk | [END_REF] ) states that the two following assertions are equivalent: a) f' 0 ; :::

; ' g 1 g is a Chebyshev system in C ([ 1; 1]) b) for any f in C ([ 1; 1]
) there exists a unique best uniform approximation in V:

In the sequel we assume that the system f' 0 ; :::; ' g 1 g is a Chebyshev system in C ([ 1; 1]) and in C ([ 1; c]) with c > 1: This implies that no non null linear combination of the ' i 's may have roots in (1; c] :

We also make use of the following result.The following properties are equivalent Proposition 2 1) f' 0 ; :::; ' g 1 g is a Chebyshev system;

2) for any set of g points (x 0 ; :::; x g 1 ) in [ 1; 1] such that x i 6 = x j , and for any (y 0 ; :::; y g 1 ) in R g , there exists a unique function g in V such that g (x k ) = y k ;

3) for any g points (x 0 ; :::; x g 1 ) in [ 1; 1] such that x i 6 = x j , the determinant

:= det G; G := 0 B B B B @ ' 0 (x 0 ) : ' 0 (x j ) : ' 0 (x g 1 ) : : : : : ' i (x 0 ) : ' i (x j ) : ' i (x g 1 ) : : :: : : ' g 1 (x 0 ) : ' g 1 (x j ) : ' g 1 (x g 1 ) 1 C C C C A :
does not equal 0:

Proof. Assume 3) holds. With the set of g points (x 0 ; :::; x g 1 ) in [ 1; 1] such that x i 6 = x j , = 0 i¤ the matrix G is not invertible, which is to say that the system of equations de…ned through 0 = P g 1 i=0 a i g i (x j ) ; j = 0; :::; g 1; admits a solution a := a 0 ; :::; a g 1 di¤erent from (0; :::; 0) in R g : De…ne g := P g 1 i=0 a i g i , an element in V which is not the function x ! 0: Since , P g 1 i=0 a i g i (x) = 0 for x in fx 0 ; :::; x g 1 g it follows that g has g distinct roots in [ 1; 1] : It follows that whenever = 0; f' 0 ; :::; ' g 1 g is not a Chebyshev system. It follows that 3) is equivalent to 1): Now 2) is equivalent to 3): Indeed when G is invertible then for any (y 0 ; :::; y g 1 ) in R g 1 the system P g 1 i=0 a i g i (x j ) = y j ; j = 0; :::; g 1, has a unique solution , which means that there is a unique g in V with g (x j ) = y j for all j:

We therefore introduce the basic de…nition

De…nition 3 A regression model Y (x) = f (x) + (x)
is a Chebyshev regression model i¤ f belongs to V := span f' 0 ; :::; ' g 1 g where (' 0 ; :::; ' g 1 ) is Chebyshev system (or Haar system) of functions in C ([ 1; 1]) :

The following result stands as a generalization of the Borel Chebyshev Theorem and improves on the Haar Theorem Theorem 4 .(Generalization of Chebyshev-Borel Theorem) Let f' 0 ; :::; ' g 1 g be a Chebyshev system on [ 1; 1] , and g is any function in C ([ 1; 1]) : Then there exists a unique function h in V := span f' 0 ; :::; ' g 1 g de…ned on [ 1; 1] , which achieves

sup x2[ 1;1] jg(x) h(x)j = inf f 2V sup x2[ 1;1] jg(x) f (x)j :
Furthermore h is the only function in V such that p := g h attains its unique maximal values in at least g + 1 points in [ 1; 1]; the sign of p on those points alternates.

Proof. See [Achieser, 1992].

Remark 5

The above function h plays a similar role as the function T g 1 (Chebyshev polynomial of the …rst kind) in the polynomial regression case; see [START_REF] Broniatowski | [END_REF].

The notation M d ([ 1; 1]) designates the class of all discrete probability measures with support in [ 1; 1] :

The aim of this paper is to present the contribution of Hoel [START_REF] Hoel | [END_REF] to the construction of optimal designs for the extrapolated value of the regression function as treated by Kiefer and Wolfowitz [Kiefer and Wolfowitz, 1965]. The model and the Gauss Markov estimator are de…ned in the next Section. An orthogonalization procedure allows to express the extrapolated value as a parameter in an adequate regression model. Finally the support of the optimal design will be obtained through geometrical arguments; the number of replications of the experiments on the nodes will then be deduced.

The model and Gauss Markov estimator

We consider a Chebyshev system on [ 1; 1] f' 0 ; :::; ' g 1 g :

For any x 2 [ 1; 1] we assume that we may observe a r.v. y (x), such that, denoting := ( 0 ; ::

; g 1 ) 0 f (x) := E (Y (x)) = g 1 X j=0 j ' j (x) = (X (x)) 0 :
(1)

We notice that the function

f : R ! R; x 7 ! f (x)
is continuous on R: Indeed since the system of the g equations in 8 < :

f (x 0 ) = P g 1 j=0 j ' j (x 0 ) :::::::::::::::::::::::::::::::: f (x g 1 ) = P g 1 j=0 j ' j (x g 1 ) has a unique solution whenever (f (x 0 ) ; :::; f (x g 1 )) 0 is known, for any (x 0 ; :::; x g 1 ) 0 2 [ 1; 1] g with 1

x 0 < ::: < x g 1 1; the function f can be extended on R; this extension is continuous since so are the ' i 's.

Recall that the measurements can be performed only on [ 1; 1], and not for jxj > 1:

Examples of Chebyshev systems

Here is a short list of classical chebyshev systems. We refer to the classical treaties of [Karlin and Studden, 1966] for a extensive study of those systems and their applications in analysis and in statistics.

a) f' 0 (x) = 1; ' 1 (x) = x 3 g is a Chebyshev system on whole R, b) ' 0 (x) = 1; ' 1 (x) = 1
x 3 is a Chebyshev system on (0; +1), c) f1; cos x; cos 2x; :::; cos nxgis a Chebyshev system on [0; ), d)f1; sin x; cos x; sin 2x; cos 2x:::; sin nx; cos nxg is a Chebyshev system on R 2 e) fsin x; sin 2x; :::; sin nxg is a Chebyshev system on [0; ),

f) f' 0 (x) = x 2 x; ' 1 (x) = x 2 + x; ' 2 (x) = x 2 + 1gis a Chebyshev system on R;
g)fx a 0 ; :::; x an ; where 0 = a 0 < ::: < a n g is a Chebyshev system on [0; +1), h) fe a 0 x ; :::; e anx ; where 0 = a 0 < ::: < a n gis a Chebyshev system on R, i) f1; sinh x; cosh x; :::; sinh nx; :::; cosh nxg is a Chebyshev system on R, j) (x + a 0 ) 1 ; :::; (x + a n ) 1 ; where 0 = a 0 < ::: < a n is a Chebyshev system on [0; +1), k) f1; log x; x; x log x; x 2 ; x 2 log x; :::; x n ; x n log xgis a Chebyshev system on (0; 1), .... . Finally note that being a Chebyshev system is a linear property; indeed if (' 0 ; :::; ' g 1 ) is a Chebyshev system then any other basis of span f' 0 ; :::; ' g 1 g is a Chebyshev system.

Description of the dataset coming from the experiment

Given the set of nodes 1 x 0 < ::: < x g 1 1, the experiment is described through the following measurements

8 > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > :
Y 1 (x 0 ) = P g 1 j=0 j ' j (x 0 ) + " 1 (x 0 ) :::::::::::::::::::::::::::::::::::: Y n 0 (x 0 ) = P g 1 j=0 j ' j (x 0 ) + " n 1 (x 0 ) :::::::::::::::::::::::::::::::::::: Y 1 (x i ) = P g 1 j=0 j ' j (x i ) + " 1 (x i ) :::::::::::::::::::::::::::::::::::: Y n i (x i ) = P g 1 j=0 j ' j (x i ) + " n i (x i ) :::::::::::::::::::::::::::::::::::: Y 1 (x g 1 ) =

P g 1 j=0 j ' j (x g 1 ) + " 1 (x g 1 ) :::::::::::::::::::::::::::::::::::: Y n g 1 (x g 1 ) = P g 1 j=0 j ' j (x g 1 ) + " n g 1 (x g 1 ) or through the more synthetic form, with

Y (x i ) := Y 1 (x i ) ; :::; Y n j (x i ) 0 , X (x i ) := (' 0 (x i ) ; :::; ' g 1 (x i )) 0 , " (x i ) := (" 1 (x i ) ; :::; " n i (x i )) 0 , Y (x i ) = 0 @ ' 0 (x i ) ::: ' g 1 (x i ) ::: ::: ::: ' 0 (x i ) ::: ' g 1 (x i ) 1 A 0 @ 0 : g 1 1 A +" (x i ) , i = 0; :::; g 1.

Denote,

X i := 0 @
' 0 (x i ) ::: ' g 1 (x i ) ::: ::: ::: ' 0 (x i ) ::: ' g 1 (x i ) 1 A :

The matrix X i has n i lines and g columns. All lines of X equal

(X (x i )) 0 . Denote H := Im X (x) := fX (x) 2 R g : x 2 [ 1; 1]g : (2)
The set H is called the regression range.

It may be at time convenient to attribute distinct indices to the same x j when repeated n j times.

The discrete measure de…ned through z }| { x 0 ; :::; x 0 (n 0 times); :::; z }| { x j ; :::; x j (n j times ); :::; z }| { x g 1 ; :::; x g 1 (n g 1 times) with n 0 + ::: + n g 1 = n will hence be written as

t 1 ; :::; t n (3)
with t 1 = t 2 = :: = t n 0 = x 0 ; :::; t n 0 +::+n g 2 +1 = t n 0 +::+n g 2 +2 = ::: = t n 0 +::+n g 2 +n g 1 = x g 1 ; hence t 1 ; :::; t n 0 indicates the same point x 0 repeated n 0 times, etc.

The system which describes the n observations writes therefore as

Y = C + "
where

Y := 0 @ Y 1 : Y n 1 A C := 0 B B B B @
' 0 (t 1 ) ::: ' g 1 (t 1 ) :: ::: ::: ' 0 (t i ) ::: ' g 1 (t i ) ::: :: :::

' 0 (t n ) ::: ' g 1 (t n ) 1 C C C C A , := 0 @ 0 : g 1 1 A , " := 0 @ " 1 : " n 1 A ; E (Y) = C , var (") = 2 I n ;
and I n is the Identity matrix of order n:

The Gauss Markov estimator of f (x) = E (y (x)) is the solution of the linear system

X 0 i X i = X 0 i Y (x i ) , i = 0; :::; g 1: It holds X 0 i X i = 0 @ ' 0 (x i )
::: ' 0 (x i ) :::

::: :::

' g 1 (x i ) ::: ' g 1 (x i ) 1 A 0 @
' 0 (x i ) ::: ' g 1 (x i ) ::: ::: :::

' 0 (x i ) ::: ' g 1 (x i ) 1 A = n i M i
where

M i := 0 B B B B B B @ (' 0 (x i )) 2 ::: ' 0 (x i ) ' k (x i )
::: ' 0 (x i ) ' g 1 (x i ) :::::: :::: :::: :::: :::::

' h (x i ) ' 0 (x i )
.:::

' h (x i ) ' k (x i )
.::: ' h (x i ) ' g 1 (x i ) :::::

.::: .::: .::: .:::

' g 2 (x i ) ' 0 (x i ) ::::: ' g 2 (x i ) ' k (x i ) ::::: ' g 2 (x i ) ' g 1 (x i ) ' g 1 (x i ) ' 0 (x i ) ::::: ' g 1 (x i ) ' k (x i ) ::::: (' g 1 (x i )) 2 1 C C C C C C A .
We have

M i = X (x i ) X 0 (x i ) : In X (x i ) X 0 (x i ) = X 0 i Y (x i
) , i = 0; :::; g 1 sum both sides with respect to i to obtain

g 1 X i=0 X i X 0 i = g 1 X i=0 X 0 i Y (x i ) : Therefore n g 1 X i=0 n i n M i ! = g 1 X i=0 X 0 i Y (x i ) . Denote i := (x) := n i n if x = x i 0 if x = 2
fx 0 ; :::; x g 1 g :

The matrix

M ( ) := g 1 X i=0 n i n M i = g 1 X i=0 i M i (4)
is the moment matrix of the measure . By de…nition supp ( ) = fx 0 ; :::; x g 1 g .

Since

M i = X (x i ) X 0 (x i )
we may write

M ( ) = g 1 X i=0 i M i = g 1 X i=0 i X (x i ) X 0 (x i ) = Z [ 1;1] X (x) X 0 (x) d (x) .
Speci…c study of this matrix is needed for the estimation of linear forms of the coe¢cients i 's. This area has been developed by Elfving see e.g. [Pukelsheim, 1993]), out of the scope of the present paper.

3 An expression of the extrapolated value through an orthogonalization procedure

We will consider an alternative way, developed by Kiefer and Wolfowitz [Kiefer and Wolfowitz, 1965] as follows. It has the main advantage that up to a coe¢cient g 1 which depends on the values of f on the x 0 j s, the estimate of f (c) is ' g 1 (c): It follows that only the coe¢cient g 1 has to be estimated, a clear advantage. Recall that c does not belong to [ 1; 1] :

It is more convenient, at this stage, to introduce the following notation. It will be assumed that n measurements of Y are performed, namely Y (t 1 ); ::; Y (t n )

where the t 0 i s belong to [ 1; 1] : The points of measurement t 1 ; ::; t n might be distinct or not, as de…ned in (3). Obviously when de…ning the optimal design with nodes x 0 ; ::; x g 1 , then n j values of the t 0 i s coincide on x j for 0 j g 1: In order to de…ne the estimator, and not the design, it is however more convenient to di¤erentiate between all the measurements Y (t i ); 1 i n: This allows to inherit from the classical geometric least square formalism.

We consider the basis of V de…ned as follows: Set for all j between 0 and g 2

h j (x) := ' j (x) ' j (c) ' g 1 (c) ' g 1 (x) (5)
and

h g 1 (x) := ' g 1 (x):
Clearly (h 0 ; ::; h g 1 ) generate V . Also (h 0 ; ::; h g 1 ) is a Chebyshev system on [ 1; c] : Denote ( 0 ; ::; g 1 ) the coordinates of f on (h 0 ; ::; h g 1 ) , namely

f (x) = g 1 X j=0 j h j (x):
We evaluate the coe¢cients j with respect to the 0 k s de…ned in (1). It holds j := j for j = 0; :::; g 2 and

g 1 := P g 1 j=0 j ' j (c) ' g 1 (c)
assuming ' g 1 (c) 6 = 0; and obviously we have

f (x) = g 1 X j=0 j h j (x) = g 1 X j=0 j ' j (x) :
In x = c we get

f (c) := g 1 X j=0 j h j (c) = g 1 X j=0 j ' j (c) :
By the de…nition of g 1 we have

g 1 := P g 1 j=0 j ' j (c) ' g 1 (c)
and therefore we have proved

Lemma 6 f (c) = g 1 ' g 1 (c) : (6)
4 The Gauss Markov estimator of the extrapolated value

It holds

f (c) = g 1 X i=1 i ' i (c)
where the i 's are de…ned through g equations of the form

f (x j ) = g 1 X i=0 i ' i (x j )
with 1 x j 1 for all 0 j g 1: Replace f (x j ) by its estimate

[ f (x j ) := 1 n j n j X i=1 Y i (x j ):
Under the present model, [ f (x j ) is an unbiased estimate of f (x j ): Determine b i though the system de…ned by

[ f (x j ) = g 1 X i=0 b i ' i (x j ):
The resulting b i 's are unbiased and so is

d f (c) = g 1 X i=0 b i ' i (c):
The natural optimality criterion associated to this procedure is the variance of the estimate d f (c) which depends on the location of the nodes and on the weights n j 's.

We now write the above Gauss Markov estimator of f (c) on the new basis (h 0 ; ::; h g 1 ) : Substituting the function f by its expansion on the basis (h 0 ; ::; h g 1 ) the model write as

8 > > > > < > > > > :
Y (t 1 ) = 0 h 0 (t 1 ) + ::: + g 2 h g 2 (t 1 ) + g 1 ' g 1 (t 1 ) + " 1 ::::::::::::::::::::::::::::::::::::::::::::::: Y (t i ) = 0 h 0 (t i ) + ::: + g 2 h g 2 (t i ) + g 1 ' g 1 (t i ) + " i :::::::::::::::::::::::::::::::::::::::::::::::::

Y (t n ) = 0 h 0 (t n ) + ::: + g 2 h g 2 (t n ) + g 1 ' g 1 (t n ) + " n : because of (5), Y (t) = T + "
where t := (t 1 ; ::; t n ) 0

T := 0 B B B B @
h 0 (t 1 ) : : h g 2 (t 1 ) ' g 1 (t 1 ) : : : : : h 0 (t i ) : : h g 2 (t i ) ' g 1 (t i )

: : :

: : h 0 (t n ) : : h g 2 (t n ) ' g 1 (t n ) 1 C C C C A ; := 0 B B @ 0 : g 2 g 1 1 C C A ; " := 0 @ " 1 : " n 1 A :
Recall that we intend to estimate g 1 : We make a further change of the basis of V: We introduce a vector G g 1 , which together with h 0 ; :::; h g 2 will produce a basis (h 0 ; :::; h g 2 ; G g 1 ) for which the vector G g 1 is orthogonal to any of the h j , 0 j g 2: The aim of this construction is to express f (c) as a linear combination of the components of G g 1 : Since G g 1 belongs to V = span (h 0 ; ::; h g 1 ) we write

G g 1 (t i ) := h g 1 (t i ) g 2 X j=0 j h j (t i )
for some vector := ( 0 ; ::; g 1 ) 0 :

We impose the following condition * 0

@ G g 1 (t 1 ) : G g 1 (t n ) 1 A ; 0 @ h j (t 1 ) : h j (t n ) 1 A +
= 0; for all j = 0; :::; g 2 ; where the above symbol <; > is the inner product in R n :The j 's in R are to be chosen now.

The linear system n X i=1 G (t i ) h j (t i ) = 0; for j = 0; :::; g 2 with g 1 equations has g 1 unknown variables j .

Once obtained the solution j ; j = 0; :::; g 2 ; and since

h g 1 (t) = G g 1 (t) + g 2 X j=0 j h j (t) ;
we may write f (t) for any t

f (t) = g 1 X j=0 j h j (t) = 0 h 0 (t) + ::: + g 2 h g 2 (t)
+ g 1 G g 1 (t) + g 1 0 h 0 (t) + :::

+ g 1 g 2 h g 2 (t) = ( 0 + g 1 0 ) h 0 (t) + ::: + ( g 2 + g 1 g 2 ) h g 2 (t) + g 1 G g 1 (t) = 0 h 0 (t) + ::: + g 2 h g 2 (t) + g 1 G g 1 (t) ;
where the 0 j s are de…ned by j := j + g 1 j for j = 0; :::; g 2 g 1 for j = g 1

:

The point is that g 1 appears as the coe¢cient of G g 1 , namely the last term in the regression of f (t) on the regressors (h 0 ; ::; h g 2 ; G g 1 ) : Furthermore G g 1 is orthogonal to the other regressors. The system which describes the data is now written by

Y (t) = e T e + "
where

e T := 0 B B B B @ h 0 (t 1 ) : : h g 2 (t 1 ) G g 1 (t 1 ) : : : : : h 0 (t i ) : : h g 2 (t i ) G g 1 (t i ) : : : : : h 0 (t n ) : : h g 2 (t n ) G g 1 (t n ) 1 C C C C A ; e := 0 B B @ 0 : g 2 g 1 1 C C A :
The minimum least square estimation of g 1 is obtained through the normal equations imposing

Y (t) e T b e 2 V ?
where b e hence designates the least square estimator of the vector of coe¢cients e ; and where V ? is the orthogonal linear space of V:

We have, denoting d g 1 the least square estimator of g 1 ; and noting that V = span fh 0 ; :::; h g 2 ; G g 1 g

* 0 B B @ Y (t 1 ) P g 2 j=0 j h j (t 1 ) d g 1 G g 1 (t 1 ) : : Y (t n ) P g 2 j=0 j h j (t n ) d g 1 G g 1 (t n ) 1 C C A 0 ; 0 B B @ h j (t 1 ) : : h j (t n ) 1 C C A + = 0;
for j = 0; :::g 2 and

* 0 B B @ Y (t 1 ) P g 2 j=0 j h j (t 1 ) d g 1 G g 1 (t 1 ) : : Y (t n ) P g 2 j=0 j h j (t n ) d g 1 G g 1 (t n ) 1 C C A 0 ; 0 B B @ G g 1 (t 1 ) : : G g 1 (t n ) 1 C C A + = 0 Hence X i=1;::;n Y (t i ) g 2 X j=0 j h j (t i ) d g 1 G g 1 (t i ) ! G g 1 (t i ) = 0: (7)
Inserting the orthogonality condition n X i=1 G (t i ) h j (t i ) = 0; for j = 0; :::; g 2;

in (7) we have n X j=1 Y (t j ) G g 1 (t j ) d g 1 n X j=1 G 2 g 1 (t j ) = 0;
and

d g 1 = P n j=1 Y (t j ) G g 1 (t j ) P n j=1 G 2 g 1 (t j )
:

Finally we obtain the explicit form of the estimator of f (c): It holds

Proposition 7

The least square estimator (Gauss Markov) of the extrapolated value f (c) is

d f (c) = ' g 1 (c) d g 1 = ' g 1 (c) P n j=1 Y (t j ) G g 1 (t j ) P n j=1 G 2 g 1 (t j ) :
5 The Optimal extrapolation design for the Chebyshev regression

The support of the optimal design

We determine the support of the optimal design for the extrapolation of f at point c:

Recall that a design is optimal if and only if it produces a Gauss Markov estimator of f (c) with minimal variance among all such estimators built upon other designs.

We note that the variance of d f (c) depends on the x j 's since

var d f (c) = (' g 1 (c)) 2 P n j=1 var (Y (t j )) G 2 g 1 (t j ) P n j=1 G 2 g 1 (t j ) 2 = ( ' g 1 (c)) 2 P n j=1 G 2 g 1 (t j ) :
The design is de…ned through a discrete probability measure 2 M d ([ 1; 1]) with support (x 0 ; ::; x g 1 ) with (x j ) := n j =n and n j equals the number of the t 0 i s which equal x j , for 0 j g 1: We now determine the support of the optimal design denoted .

:= arg min 2M d ([ 1;1]) 1 P g 1 j=0 n j G 2 g 1 (x j ) = arg max 2M X g 1 X i=0 n j G 2 g 1 (x j ) = arg max 2M d ([ 1;1]) g 1 X i=0 n j h g 1 (x i ) g 2 X j=0 j h j (x i ) ! 2 :
The solution can be obtained in a simple way through some analysis of the objective function. By convenience in order to use simple geometric arguments and to simplify the resulting expressions it is more convenient to write the derivation of the optimal design in terms of the t 0 i s:

The function

n X i=1 h g 1 (t i ) g 2 X j=0 j h j (t i ) ! 2 = 0 @ h g 1 (t 1 ) P g 2 j=0 j h j (t 1 ) : h g 1 (t n ) P g 2 j=0 j h j (t n ) 1 A 2
is the distance from the orthogonal projection of the vector h := h g 1 (t 1 ) ::: h g 1 (t n ) 0 on the linear space V generated by the family fh 0 ; :::; h g 2 ; G g 1 g : Therefore by the minimal projection property

n X i=1 h g 1 (t i ) g 2 X j=0 j h j (t i ) ! 2 = min 2V dist (h; ) : Let := 0 ::: g 2 0 :
The optimal design is obtained through a two steps procedure. Fix the frequencies n 0 ; ::; n g 1 with sum n and determine the discrete measure on [ 1; 1] which minimizes var d f (c) among all 0 s with support x := (x 0 ; ::; x g 1 ) and masses (x j ) = n j =n; 0 j g 1. The optimization is performed upon the x 0 j s: The optimal design solves therefore the problem

= arg max 2M d ([ 1;1]) min 2V dist (h; ) = arg max x2[ 1;1] g min 2R g 1 g 1 X i=0 n i h g 1 (x i ) g 2 X j=0 j h j (x i ) ! 2 = arg max 2M d ([ 1;1]) min 2R g 1 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2
(dx) :

The integrand h g 1 (x) P g 2 j=0 j h j (x)

2 is always non negative. Henceforth it is enough to minimize its square root w.r.t. x. This optimization turns therefore to be independent of the n 0 j s: Denote j ; j = 0; :::; g 2; the values which minimize dist (h; ) w.r.t. j . The optimality condition writes

max x2[ 1;1] h g 1 (x) g 2 X j=0 j h j (x) = min 2R g 1 max x2[ 1;1] h g 1 (x) g 2 X j=0 j h j (x) (8) = min p2W max x2[ 1;1] jh g 1 (x) p (x)j
where W := span fh 0 ; ::; h g 2 g :

If we prove that fh 0 ; :::; h g 2 g is a Chebyshev system on [ 1; 1] ; then clearly the support of the optimal measure consists in the points of maximal value in [ 1; 1] for the function

jh g 1 (x) p (x)j
where p is the best uniform approximating polynomial of h g 1 in W: Indeed the support of consists in the set of points where

h g 1 (x) g 2 X j=0 j h j (x)
in (8) attains its maximal value for p = p the best uniform approximation of h g 1 in W:

This is the major argument of the present derivation, which justi…es all of the uniform approximation theory in this context.

De…nition 8

The vector in R g 1 is a Chebyshev vector i¤ it designates the vector of the coe¢cients of p , where p is the best uniform approximating polynomial of h g 1 in W de…ned in (9). It is de…ned through (8). Now writing := 0 ; :::; g 2 0 we de…ne the set of all points e

x in [ 1; 1] where the distance between h g 1 and its best approximation on the h k , 0 k g 2 is minimal. These points are precisely the support of the optimal design . Formally we de…ne We prove that fh 0 ; :::; h g 2 g is a Chebyshev system on [ 1; 1] :

E := min 2R g 1 max x2[ 1;1] h g 1 (x) g 2 X j=0 j h j (x) (10
Proposition 9 (Hoel) The functions h 0 ; :::; h g 2 are a Chebyshev system on [ 1; 1] :

Proof. For any choice of fx 0 ; :::; x g g in [ 1; 1], with x 0 < ::: < x g 1 ; since the family f' 0 ; :::; ' g 1 g is a Chebyshev system on [ 1; 1], we have, by Proposition 2, assuming a positive sign of the determinant, without loss of generality Hence the family f' 0 ; :::; ' g 1 g is a Chebyshev system in C ([ 1; 1]) :

In the same way we can prove that it is a Chebyshev system in [ 1; c] :

The frequencies of the optimal design

Once characterized the points x in supp ; we characterize the values of the (x) 's:The following Proposition produces a su¢cient condition in order that the measure be optimal, which can be phrased as

min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx)
for any in M d ([ 1; 1]) : Uniqueness might not hold.

Proposition 10 (Kiefer-Wolfowitz) Let B ( ) be de…ned as in (11). If is Chebyshev vector and (B ( )) = 1 and if

Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x)
! h i (x) (dx) = 0; for i = 0; :::; g 2;

then is optimal. Proof. Let 2 M d ([ 1; 1]) with (B ( )) = 1: The hypothesis Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! h i (x) (dx) = 0;
for i = 0; :::; g 2; indicates that

h g 1 (x i ) g 2 X j=0 j h j (x i )
is orthogonal to the linear space W generated by fh 0 ; :::; h g 2 g : Thus P g 2 j=0 j h j is the orthogonal projection of h g 1 on W: The inner product is

< v; w >:= Z [ 1;1] v (x) w (x) (dx) :
By the minimal projection property

A ( ) := min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) = Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) = X e x2supp h g 1 (e x) g 2 X j=0 j h j (e x) ! 2 (e x) E 2 X e x2supp (e x) = E 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) ; max v2M d ([ 1;1]) min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) = min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) =: A ( )
The measure v which appears in lines 5 and followings in the above displays are arbitrary measures in M d ([ 1; 1]) :

Since by de…nition

:= arg max v2M d ([ 1;1]) min 2R g 2 Z [ 1;1] h g 1 (x) g 2 X j=0 j h j (x) ! 2 (dx) (12) i.e.
A ( ) A ( ) :

Hence A ( ) = A ( ) :
6 Identi…cation of the optimal design

In this Section we provide an explicit solution for the optimal design and prove its uniqueness.

By the Borel-Chebyshev Theorem 4 there exist at least g points x 0 < ::: < x g 1 in [ 1; 1] on which the best uniform approximation of h g 1 ; namely P g 2 j=0 j h j ; satis…es the following conditions

h g 1 (x i ) g 2 X j=0 j h j (x i ) = ( 1) i E:
We now see that there are exactly g points on which the function h g 1 P g 2 j=0 j h j equals E:

Since fh 0 ; :::; h g 1 g is a Chebyshev system the linear combination which is the absolute value of a linear combination of the Chebyshev system fh 0 ; :::; h g 1 g cannot have more than g roots. Therefore h g 1 P g 2 j=0 j h j cannot have more than g + 1 maximal values.

As seen previously the support of the optimal measure consists in the points of maximal value in [ 1; 1] for the function

h g 1 g 2 X j=0 j h j :
Applying the Borel-Chebyshev Theorem we now determine the support of :

Since E is known the support is the vector x 1 ; :::; x g which solves the linear system h g 1 (x i ) g 2 X j=0 j h j (x i ) = ( 1) i E; i = 1; :::; g:

We apply the su¢cient condition provided by Kiefer and Wolfowitz above, Proposition 10. This condition states that the values (x i ) ; i = 0; :::; g 1; satisfy the system ( P g 1 i=0 h g 1 (x i ) P g 2 j=0 j h j (x i ) h r (x i ) (x i ) = 0 r = 0; :::; g 2 :

In the x i 's it holds

E = h g 1 (x i ) g 2 X j=0 j h j (x i )
and

Z [ 1;1] h g 1 (x)
g 2 X j=0 j h j (x) ! h i (x) (dx) = 0; for i = 0; :::; g 2:

Therefore 0 = g 1 X i=0 h g 1 (x i ) g 2 X j=0 j h j (x i ) ! h r (x i ) (x i ) = E g 1 X i=0
( 1) i h r (x i ) (x i ) , for r = 0; :::; g 2:

The optimal extrapolation design f(x i ; (x i )) : i = 0; :::g 1g thus solves 8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > :

h g 1 (x 0 ) P g 2 j=0 j h j (x 0 ) = +E ::::::::::::::::::::::::::::::::::::::::::::::: h g 1 (x i ) P g 2 j=0 j h j (x i ) = ( 1) i E :::::::::::::::::::::::::::::::::::::::::::::: h g 1 (x g 1 ) P g 2 j=0 j h j (x g 1 ) = ( 1) g 1 E P g 1 i=0 ( 1) i h 0 (x i ) (x i ) = 0 ::::::::::::::::::::::::::::::::::::: P g 1 i=0 ( 1) i h r (x i ) (x i ) = 0 :::::::::::::::::::::::::::::::::: P g 1 i=0 ( 1) i h g 2 (x i ) (x i ) = 0 In practice we …rst evaluate j for 0 j g 2 through (8). Note that E is known by (10). The above system consits in 2g 1 equations in the 2g unknown quantities f(x i ; (x i )) : i = 0; :::g 1g : Add the constraint (x 0 ) + :: + (x g 1 ) = 1 to obtain a linear system with a unique solution.

The …rst g equations determine the nodes, by Borel Chebyshev Theorem. The last g 1 ones determine the values of the n 0 j s by the Proposition of Kiefer and Wolfowitz 10. Hence there is a unique optimal design solving the minimal variance problem for the extrapolation.

' 0 (x 0 ) ' 0 (x 1 ) : ' 0 (x g 2 ) ' 0 (c) ' 1 (x 0 ) ' 1 (x 1 ) : ' 1 (x g 2 ) ' 1 (c) : : : : :

For j = 0; :::; g 1;the operations

do not change the value of the determinant. Hence,

By the Laplace Theorem pertainig to determinants, we get

Therefore the two real numbers ' g 1 (c) and have same sign. Since ' g 1 (c) 6 = 0 we deduce that det 0 B B @ h 0 (x 0 ) h 0 (x 1 ) : h 0 (x g 2 ) h 1 (x 0 ) h 1 (x 1 ) : h 1 (x g 2 ) : : : : h g 2 (x 0 ) h g 2 (x 1 ) : h g 2 (x g 2 ) 1 C C A 6 = 0: