Olivier Mullier
email: mullier@ensta.fr

Alexandre Chapoutot
email: chapoutot@ensta.fr

Julien Alexandre
email: alexandre@ensta.fr

Dit Sandretto

Validated Computation of the Local Truncation Error of Runge-Kutta Methods with Automatic Differentiation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Differential equations are used in many domains to describe the temporal evolution or the dynamics of systems (e.g., in chemistry, biology or physics). Most of the time, simpler form of differential equations that is Ordinary Differential Equation (ODE) is considered and has to be solved. A closed form of the solution of an ODE is usually no computable, except for some particular classes of ODEs, e.g., linear ordinary differential equations, then this is usually done by numerical integration.

More precisely, numerical integration methods are able to solve initial value problem (IVP ) of non-autonomous ODE defined by ẏ = f (t, y) with y(0) = y 0 and t ∈ [0,

t end ] . (1)
The function f : R × R n → R n is the flow, y ∈ R n is the vector of state variables, and ẏ is the derivative of y with respect to time t. We shall always assume at least that f is globally Lipschitz in y, so Equation (1) admits a unique solution [START_REF] Hairer | Solving Ordinary Differential Equations I Springer[END_REF] for a given initial condition y 0 . Even more, for our purpose, we shall assume that f is continuously differentiable as needed. The exact solution of Equation ( 1) is denoted by y(t; y 0 ). The goal of a numerical integration is to compute a sequence of time instants 0 = t 0 < t 1 < • • • < t N = t end (not necessarily equidistant) and a sequence of states y 0 , . . . , y N produced by the numerical integration method such that ∀ ∈ [0, N ], y ≈ y(t , y -1 ).

To solve IVP-ODE problems with numerical integration, several methods have been studied and defined for a long time [START_REF] Hairer | Solving Ordinary Differential Equations I Springer[END_REF]. These methods are gathered in different classes following different criteria such as single-step or multi-step methods; as fixed step-size or variable step-size methods; or also as explicit or implicit methods. This work is focused on single-step, variable step-size, explicit and implicit methods belonging to the familly of Runge-Kutta methods.

Validated numerical integration is an appealing approach to produce rigorous results on the solution of IVP-ODEs. Most of the rigorous techniques defined so far, since Ramon Moore's seminal work [START_REF] Moore | Interval Analysis[END_REF], are based on Taylor series approach, see for example [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF][START_REF] Makino | COSY INFINITY version 9[END_REF][START_REF] Lin | Validated solutions of initial value problems for parametric odes[END_REF][START_REF] Dzetkulič | Rigorous integration of non-linear ordinary differential equations in Chebyshev basis[END_REF] and the references therein. Nevertheless, it is unlikely that only one kind of methods is adapted to all various classes of ODEs. So, more recent work [START_REF] Gajda | Three-and four-stage implicit interval methods of Runge-Kutta type[END_REF][START_REF] Marciniak | On representations of coefficients in implicit interval methods of runge-kutta type[END_REF][START_REF] Marciniak | Implicit interval methods for solving the initial value problem[END_REF][START_REF] Bouissou | GRKLib: a Guaranteed Runge Kutta Library[END_REF][START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF] deals with the adaptation of Runge-Kutta methods for ODEs and DAEs (Differential-Algebraic Equations), see [START_REF] Alexandre | Validated simulation of differential algebraic equations with runge-kutta methods[END_REF]. The motivation having validated Runge-Kutta methods is to benefit the strong stability properties of these methods such as the A-stability which is a relevant requirement when dealing with stiff problems. A new guaranteed numerical integration scheme based on Runge-Kutta methods has already been introduced in [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF] by the authors. The challenge is in the computation of an enclosure of the method error. This error is defined by the Local Truncation Error (LTE) which involve high order temporal derivatives of f . This previous work was not based on Automatic Differentiation [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] but on symbolic differentiation of the flow f . The symbolic approach introduces some limitations on the size of the IVP-ODEs which has to be solved.

In this article, we propose an enhancement of this guaranteed numerical integration method by computing the local truncation error using B-series which stand for Butcher series. B-series were introduced by E. Hairer and G. Wanner in [START_REF] Hairer | On the butcher group and general multi-value methods[END_REF] after the seminal work of J. Butcher [START_REF] Butcher | The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods[END_REF]. They follow the results from Cayley [START_REF] Cayley | On the theory of the analytical forms called trees[END_REF] and Merson [START_REF] Rh Merson | An operational method for the study of integration processes[END_REF] on the expressions arising when studying the derivatives of the solution of an ODE. Since then, a lot of attention has been made in the study of B-series (see, e.g., [START_REF] Munthe-Kaas | Lie-butcher theory for runge-kutta methods[END_REF], [START_REF] Chartier | Algebraic structures of B-series[END_REF] or [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF]). More recently, results on the computation of B-series were published in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] on which this work strongly rely. Nevertheless, enhancements of this algorithm have been defined and are the main results of this article which are

• an instance of the approach based on automatic differentiation defined in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] is given, to efficiently compute the local truncation error of explicit and implicit Runge-Kutta methods;

• a set-membership method to compute an enclosure of the local truncation error producing rigorous bounds of the error method is presented;

• an experimental evaluation of the conjectured time complexity of the approach defined in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] is performed.

This article is presented as follows: Section 2 recalls results of the literature on Runge-Kutta theory which are mandatory to the comprehension of this work. In Section 3, the previous results on the computation of B-series are briefly introduced before the main result of this work, validated numerical integration based on Runge-Kutta methods, is presented in Section 4. Experiments are provided in Section 5 to study the time complexity of our method using automatic differentiation compared to a previous one using symbolic differentiation before the conclusion.

A Recall on Modern Runge-Kutta Theory

This section is dedicated to the introduction of the Runge-Kutta methods and the relevant work on it for the purpose of this article.

Runge-Kutta Methods

Runge-Kutta methods [START_REF] Runge | Über die numerische auflösung von differentialgleichungen[END_REF][START_REF] Wilhelm | Beitrag zur naherungsweisen integration totaler differential gleichungen[END_REF] are numerical methods to the computation of the solution of IVP-ODE defined in Equation [START_REF] Hairer | Solving Ordinary Differential Equations I Springer[END_REF].

Definition 1 (Runge-Kutta method). Let f : R × R n → R n be a flow and the point y 0 = y(t 0 ) at time t 0 . y 1 ≈ y(t 0 + h; y 0 ) is computed using a given step-size h and following the process

     k i = f (t 0 + c i h, y 0 + Y i ), i = 1, . . . , s Y i = h s j=1 a ij k j , i = 1, . . . , s y 1 = y 0 + h s j=1 b j k j (2) with k i , Y i ∈ R n and c i , a ij , b j ∈ R.
Note that c i = s j=1 a ij and the real coefficients c i , a ij and b i fully characterise a Runge-Kutta method, see [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF]. They can be represented by a partitioned tableau, known as the Butcher tableau [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF]:

c A b T (3) with A = (a ij ) 1 i,j s , b = (b i ) 1 i s and c = (c i ) 1 i s .
The form of the square matrix A determines the Runge-Kutta method: explicit for A strictly lower triangular, diagonally implicit when A is lower triangular and fully implicit for a full matrix A.

In order to simply the notation, autonomous systems of ODEs will be considered. This restriction is not severe as it is always possible to transform non-autonomous systems of ODEs as defined in Equation ( 1) into an autonomous one using the rewriting rule

ẏ = f (t, y) with y(0) = y 0 ⇐⇒ ż = ẏ ẋ = f (x, y) 1 = g(z) with z(0) = y 0 0 .

Order Theory of Runge-Kutta Methods

The accuracy of numerical integration methods is related to the notion of order. More precisely, the order of a method is p if and only if the local truncation error (LTE) is such that

||y(t ; y -1 ) -y || = O(h p+1 )
with || • || some norm in R n and h = t -t -1 . In the context of Runge-Kutta methods, the order p means that the Taylor expansion of the exact solution y(t ; y -1 ) and the one of the Runge-Kutta method y share the same p first Taylor components. From the work of John Butcher in [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF], these coefficients for the exact and the numerical solution can be expressed on the same basis of elementary differentials. If we write the Taylor expansion of the exact solution of Equation (1), i.e., y(t ; y -1 ), one has

y(t ; y -1 ) = y(t -1 ) + h ẏ(t -1 ) + 1 2 h 2 ÿ(t -1 ) + 1 3! h 3 y (3) (t -1 ) + 1 4! h 4 y (4) (t -1 ) + . . . , (4)
the components can be rewritten in terms of these elementary differentials

ẏ = f (y) = f ÿ = f (1) (y) ẏ = f (1) f y (3) = f (2) (y)( ẏ, ẏ) + f (1) (y)ÿ = f (2) (f , f ) + f (1) f (1) f y (4) = f (1) (y)(f (1) (y)(f (1) (y)( ẏ))) + 3f (2) (y)( ẏ, f (1) (y) ẏ) + f (1) (y)(f (2) (y)( ẏ, ẏ)) +f (3) (y)( ẏ, ẏ, ẏ) = f (1) f (1) f (1) f + 3f (2) (f , f (1) f ) + f (1) (f (2) (f , f )) + f (3) (f , f , f ) . . .
with f (k) (y) the k-th derivative of the vector field f w.r.t. y, which is a multilinear map. Equation (4) can then be rewritten as a sum of elementary differentials. From work by Merson [START_REF] Rh Merson | An operational method for the study of integration processes[END_REF] (and developed in [START_REF] John | Numerical methods for ordinary differential equations[END_REF]), elementary differentials can be in one-to-one correspondence with rooted trees which will be a useful data structure to compute these elementary differentials. Rooted trees are connected acyclic graphs with a vertex designated as the root. For any rooted tree τ , we denote F(τ ) its corresponding elementary differential.

Example 1. The elementary differentials f (2) (f , f ) is associated to the tree f f f (2) and f (1) (f (1) (f )) corresponds to the tree described by

f (1)
f (1) f A rooted tree τ is associated to properties given by (see [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF] for more details): the order denoted r(τ ) is the number of vertices in τ , the symmetry σ(τ ) is order of automorphism group, the density denoted γ(τ ) is the product over all vertices of the order of the subtree rooted to that vertex, and α(τ ) is the number of ways of labeling with ordered set. These functions will play an essential part in the computation of the elementary differentials.

A rooted tree τ can be represented by the list of its corresponding subtree(s), i.e., the list of trees appearing when the root is removed.

Example 2. The representation of the tree in term of list of subtrees is the list , .

Furthermore isomorphic subtrees can be factorised meaning that a tree τ can be represented by the list of its subtrees [τ m1 1 τ m2 2 . . . τ m k k ] with m i the number of occurrences of τ i as a subtree of τ . Using this representation, we now recall the formal statement of recursions for the computation of r(τ ), σ(τ ) and γ(τ ).

Theorem 2.1 (Theorem 301A in [START_REF] John | Numerical methods for ordinary differential equations[END_REF], p. 140). Let τ = [τ m1 1 τ m2 2 . . . τ m k k ] be a rooted tree where τ 1 , . . . , τ k are distinct. Then, for the tree consisting in only a root,

r( ) = σ( ) = γ( ) = 1 and r(τ ) = 1 + k i=1 m i r(τ i ), σ(τ ) = k i=1 m i !σ(τ i ) mi , γ(τ ) = r(τ ) k i=1 γ(τ i ) mi .
The elementary differentials are made of sums of partial derivatives of f with respect to the components of y. An other achievement in [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF] was to relate these partial derivatives of order q to a combinatorial problem to enumerate all the rooted trees τ with exactly q nodes (r(τ ) = q).

It is then possible to express the Taylor components of the exact solution of an IVP-ODE term of elementary differentials. This is expressed in the following theorem.

Theorem 2.2. The q-th derivative w.r.t. time of the exact solution is given by

y (q) = r(τ )=q α(τ )F(τ )(y 0 ) with α(τ ) = r(τ )! σ(τ )γ(τ )
. It is obvious that expressing the Taylor expansion of the numerical methods is a tedious work. The clever idea of John Butcher in [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF] is to express them also in terms of elementary differentials, i.e., on the same basis of functions that the Taylor expansion of the exact solution. Moreover, a relation between the coefficients c i , a ij and b i of the Butcher tableau with these elementary differentials are then possible. In consequence, the rooted tree structure can be used to compute also elementary weights which are particular coefficients of elementary differentials expressed in terms of coefficients of the Butcher tableau. We denote by φ(τ ) the elementary weight of τ based based on the coefficients defining a Runge-Kutta method. Following the structure of τ , φ(τ ) are computed as defined in Definition 2.

Definition 2 (Elementary weight [START_REF] Denholm | Numerical methods for ordinary differential systems: the initial value problem[END_REF]). Let τ = [τ m1 1 τ m2 2 . . . τ m k k ] be a rooted tree where τ 1 , . . . , τ k are distinct trees, is the tree with only a root. For i = 1, 2, . . . , s, s + 1, we define

φ i ( ) = s j=1 a ij = c i φ i (τ ) = s j=1 a ij k =1 φ j (τ ) m .
We then define φ( ) = i b i and φ(τ ) = φ s+1 (τ ).

The next theorem gives the definition of the q-nth derivative of the numerical solution y 1 with respects to time t.

Theorem 2.3. The q-th derivative w.r.t. time of the numerical solution is given by

y (q) 1 = r(τ )=q γ(τ )φ(τ )α(τ )F(τ )(y 0 ) with α(τ ) = r(τ )! σ(τ )γ(τ )
. Finally, using theorems 2.2 and 2.3, the order condition of a Runge-Kutta can be exhibited.

Theorem 2.4 (Order condition). A Runge-Kutta method has order p iff

φ(τ ) = 1 γ(τ ) ∀τ, r(τ ) p .
The consequence of Theorem 2.4 is that the definition of new Runge-Kutta methods can be stated as the solution of a non-linear systems of equations made of all the conditions of the form φ(τ ) = 1 γ(τ ) for all τ such that r(τ ) p. Note that this system of nonlinear equations is generally under-constrained and so the computation of the coefficients c i , a ij and b i of the Butcher tableau is usually hard.

In context of validated numerical integration based on Runge-Kutta methods, the problem is much simpler than building a new method because the starting point is the coefficients of the Butcher tableau. The challenge is then an efficient computation of the higher order derivatives expressed in terms of elementary differentials.

Local Truncation Error of Runge-Kutta Methods

As previously mentioned, and following the order condition given in Theorem 2.4, a Runge-Kutta method is of order p if the p first terms of the Taylor form associated to the numerical solution y are equal to the terms of the exact solution of (1) that is y(t ; y -1 ) assuming the same initial condition. In this case the LTE corresponds to the difference of the two Taylor remainders. The challenge is to compute these Taylor remainders. In consequence of Theorems 2.2-2.4, if the order p of the considered Runge-Kutta method is known, its LTE is then defined by

y(t j ; y(t ; y -1 ) -y = h p+1 (p + 1)! r(τ )=p+1 α(τ ) 1 -γ(τ )φ(τ ) F(τ )(y(ξ)) ξ ∈]t j , t [ (5)
with F(τ )(y(ξ)) the elementary differential associated to τ computed at point y(ξ). Remark 1. Equation ( 5) is valid for explicit and implicit Runge-Kutta methods. Indeed, elementary weights can be computed following Definition 2 and this computation only depends of the structure of rooted trees and the coefficients of the Butcher tableau.

Computation Methods of B-Series

Following Section 2.3, the computation of the LTE of an explicit or implicit Runge-Kutta method involves the computation of high order derivatives of f . For that goal, two different ways of computing Equation ( 5) can be considered: by symbolic differentiation or by automatic differentiation approach. In other terms, automatic differentiation (see, e.g., [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] and references therein) produces the derivative of a function on a particular point of its domain while symbolic differentiation produces an analytical expression of the derivative. The main drawback of the symbolic differentiation approach is its lack of scaling according to the dimension of the problem. B-series are sums of elementary differentials (such as the one occurring in Equation ( 4)). They are a fundamental tool to design Runge-Kutta methods. Remaining that they allow to exhibit the necessary conditions a Runge-Kutta method must fulfil to be of a specific order (see [START_REF] Munthe-Kaas | Lie-butcher theory for runge-kutta methods[END_REF]). B-series play also a central role in validated numerical integration method as the LTE is based on such mathematical object.

The previous work [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] defines an algorithm to compute B-series. The elementary differentials are vectorfields and have to be computed at point x ∈ R n . We have F( )(x) = f (x) and, from Equation (2.7) in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF], for the tree τ

= [τ m1 1 τ m2 2 . . . τ m k k ], F(τ )(x) = ∂ i mi f (x + z 1 F(τ 1 )(x) + z 2 F(τ 2 )(x) + • • • + z p F(τ k )(x)) ∂z m1 1 ∂z m2 2 . . . ∂z m k k z=0 . (6)
In Section 3.1, we briefly recall our previous work on computing the LTE of Runge-Kutta methods [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF] based on symbolic differentiation. In Section 3.2, we also recall the basis of the algorithm defined in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] based on automatic differentiation.

A Symbolic Approach

A symbolic approach for the computation of these B-series was discussed in [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF]. The elementary differentials, as expressed in [START_REF] Dzetkulič | Rigorous integration of non-linear ordinary differential equations in Chebyshev basis[END_REF], can be be associated to Fréchet derivatives: Definition 3 (M -th Fréchet derivatives [START_REF] Denholm | Numerical methods for ordinary differential systems: the initial value problem[END_REF]). Let f : R n → R m be a function and z ∈ R n . The M -Fréchet derivative of f is defined by

f (M ) (z)(K 1 , K 2 , . . . , K M ) = m i=1 m j1=1 m j2=1 • • • m j M =1 i f j1j2...j M j1 K 1 j2 K 2 . . . j M K M e i (7)
with

f j1j2...j M = ∂ M ∂ j1 z∂ j2 z . . . ∂ j M z and K k = [ 1 K 1 , 2 K 2 , . . . , M K M ] ∈ R m , for k = 1, .
. . , M . The notation x is the -th component of x and e i is the vector of 0 except for its i-th component which is 1.

Note that in Definition 3, the vector made of K i element can be associated to the structure of a rooted tree τ with r(τ ) = M . In consequence, the main idea of the approach in [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF] is to enumerate all the partial derivatives of f up to order p and then use the structure of the Fréchet derivative to combine them. The drawback of this approach is that the number of partial derivatives grows exponentially with the order of the method then it requires an algorithm with a complexity in O(n M +1 ) to produce these M + 1 sums in [START_REF] Gajda | Three-and four-stage implicit interval methods of Runge-Kutta type[END_REF]. While efficient it is then only suitable for small dimension problems.

An Automatic Differentiation Approach

Following [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] it is possible to make the correspondence between the computation of elementary differentials and the one of Taylor coefficients of a particular univariate function. This is resumed in the following proposition Proposition 3.1 (Taylor to tensor conversion (Proposition 10.2 in [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF])). Let f : R n → R m be at least d-times continuously differentiable at some point x ∈ R n and denote by F r (x, s) the rth Taylor coefficient of the curve f (x + ts) at t = 0 for some direction s ∈ R n . Then we have for any seed matrix S = (s j ) 1 j p ∈ R n×p and any multi-index i = (i 1 , . . . , i p ) ∈ N p with |i| = p r=1 i r d the identity

∂ |i| f (x + z 1 s 1 + z 2 s 2 + • • • + z p s p ) ∂z i1 1 ∂z i2 2 . . . ∂z ip p z=0 = |j|=d ϕ(i, j)F |i| (x, Sj), (8)
where the constant coefficients ϕ(i, j) are given by the finite sums of product of multinomials

ϕ(i, j) = 0<k i (-1) |i-k| i k dk/|k| j |k| d |i| (9)
Using Equation ( 6) and Proposition 3.1, the computation of the elementary differential F(t) is then the computation of a sum of Taylor coefficients of an univariate function with the seed matrix S = [F(τ 1 ), . . . , F(τ k )] the matrix for which columns are elementary differentials of the subtrees

τ i of τ = [τ m1 1 τ m2 2 . . . τ m k k ]
. This is done with Proposition 3.1 using automatic differentiation with the techniques described in [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF]. The next section focuses on the adaptation of these techniques in the case of set-membership computation.

Validated Numerical Integration Based on Runge-Kutta Methods

This section is dedicated to the adaptation of the results given in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] to the computation of bounds on the local truncation error of Runge-Kutta methods. Firstly, the IVP-ODE is stated when uncertainties occur and a brief introduction to set-membership computation is provided. It then follows the main results of this article.

Problem Statement

When dealing with validated computation, mathematical representation of an IVP-ODE is as follows:

ẏ = f (t, y) y(0) ∈ Y 0 ⊆ R n . ( 10
)
The set Y 0 of initial conditions is used to model some (bounded) uncertainties. For a given initial condition y 0 ∈ Y 0 , the solution at time t > 0 when it exists is denoted y(t; y 0 ). The goal, for validated (or rigorous) numerical integration methods, is then to compute the set of solutions of (10), i.e., the set of possible solutions at time t given the initial condition in the set of initial conditions Y 0 :

y(t; Y 0 ) = {y(t; y 0 ) | y 0 ∈ Y 0 }. (11)
Validated numerical integration schemes using set-membership framework aims at producing the solution of the IVP-ODE that is the set defined in [START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF]. It results in the computation of an outer approximation of y(t; Y 0 ).

Set-Membership Framework

Interval Arithmetic

Interval analysis [START_REF] Moore | Introduction to interval analysis[END_REF] is a method designed to produce outer-approximation of the set of possible values for variables occurring in some computations in a sound manner. Hereafter, an interval is denoted

[x] = [x, x] with x
x and the set of intervals is

IR = {[x] = [x, x] | x, x ∈ R, x x}.
The Cartesian product of intervals [x] ∈ IR n is a box. The main result of interval analysis is its fundamental theorem stating that the evaluation of an expression using intervals leads to an outer-approximation of the resulting set of values for this expression whatever the values considered in the intervals. In order to deal with interval functions, an interval inclusion function also known as interval extension of a function can be defined.

Definition 4 (Interval extension). Consider a function

f : R n → R m . The interval function [f ] : IR n → IR m is an interval extension of f if the evaluation of [f ] ([x]) gives an outer approximation of the image of [x] by the function f , noted f ([x]) = {f (x)|x ∈ R n }: ∀ [x] ∈ IR n , f ([x]) = {f (x)|x ∈ R n } ⊆ [f ] ([x]). ( 12
)
Many interval extensions of functions can be defined thenceforth they verify [START_REF] Alexandre | Validated simulation of differential algebraic equations with runge-kutta methods[END_REF] (see [START_REF] Moore | Interval Analysis[END_REF]). We can cite the natural extension [START_REF] Moore | Introduction to interval analysis[END_REF] which replaces the operations on reals by their interval counterparts using interval arithmetic.

Example 3. Considering a scalar initial value problem defined by ẏ = -y with y(0) ∈ [0, 1]. Applying the Euler numerical integration method, one get

y n+1 = y n + h(-y n ) (13)
with h the integration step-size. Evaluating this expression into intervals, one get

y n+1 = [0, 1] + h × (-[0, 1]) = [0, 1] + [-h, 0] = [-h, 1] .
Another interval extension is the mean value extension [START_REF] Moore | Introduction to interval analysis[END_REF] which linearizes the function around its mean value.

→ R m be a C p function, p 1, [x] ∈ IR n an interval, x ∈ [x] and J [x] f ∈ IR m×n an interval matrix such that    ∂f i ∂x j 1 i m 1 j n (x) | x ∈ [x]    ⊆ J [x] f
The function [f c ] : IR n × R n → IR m defined as follows:

[f c ] ([x] , x) = f (x) + J [x] f ([x] -x) (14)
is an interval extension of f . Example 4. Let denote g(y) = y -hy as defined in Example 3. The computation of the centered form requires an outer approximation of the derivative of g along y. We consider this outer approximation as the natural inclusion of the function g (y) = 1 -h. Then the computation of the centered form, using ỹ as m([y]), the midpoint of [y] gives

[g c ] ([y] , m([y])) = g(m([y])) + (1 -h)([y] -m([y]) = 1 2 - 1 2 h + (1 -h)([-0.5, 0.5]) = 1 2 (1 -h) + - 1 2 (1 -h), 1 2 (1 -h) = 1 2 (1 -h) - 1 2 (1 -h), 1 2 (1 -h) + 1 2 (1 -h) = [0, 1 -h]
which is here a better result than the one computed with natural extension in Example 3.

Affine Arithmetic

As interval arithmetic, affine arithmetic is a model to produce outer-approximation of ranges of set-valued expression which aims is to improve interval arithmetic in certain situations. It is designed to keep track of the linear dependencies between variables occurring in a computation. An affine form x representing the set of values taken by a variable x is denoted by:

x = α x 0 + n i=1 α x i ε i , (15)
with α x i ∈ R for all i. Each noise symbol ε i ∈ [-1, 1] is unknown. It represents an independent component of the global uncertainty on x. An interval x = [x, x] can easily be converted to an affine form x:

x = x + x 2 + x -x 2 ε 1 . (16)
From an affine form, the associated interval can be computed by replacing each noise symbol by the interval [-1, 1] and use interval arithmetic:

[x] = α x 0 + n i=1 α x i [-1, 1] (17)
which is an outer approximation of the values x can take. When nonlinear operations occur on affine forms, the nonlinear dependencies are represented by a new noise symbol η with its associated partial deviation. They represent an outer approximation of the associated nonlinear dependency. For example, a possible way to compute the multiplication between two affine forms

x = α x 0 + n i=1 α x i ε i and ŷ = α y 0 + n i=1 α y i ε i is x × ŷ = α x 0 α y 0 + 1 2 n i=1 α x i α y i + n i=1 (α x 0 α y i + α y 0 α x i )ε i +     1 2 n i=1 α x i α y i + n i=1 n j=1 j =i α x i α y j     η (18)
Example 5. The function defined in Example 3 is once more considered. The affine form associated to y is then ŷ = 1 2 + 1 2 ε 1 . The computation of g(y) using affine arithmetic then gives

g(ŷ) = 1 2 + 1 2 ε 1 + h - 1 2 + 1 2 ε 1 = 1 2 - 1 2 h α0 + 1 2 - 1 2 h α1 ε 1
which is associated to the interval value, using Equation ( 17),

1 2 - 1 2 h + 1 2 - 1 2 h [-1, 1] = [0, 1 -h]
The choice between using interval arithmetic or affine arithmetic is essentially based on the trade-of between precision and time computation. However some intrinsic properties of affine arithmetic makes them more suitable for our work. Indeed, the well-known wrapping effect which appears in interval arithmetic when some geometric rotation appears, typically when matrix-vector multiplication appears, need special treatment to reduce over-approximation. For example, see [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] for a more detailed explanation. In contrary, affine arithmetic which as a geometric representation based on zonotopes benefits of a rotational stability property. Here, affine arithmetic will be used when such precision is mandatory to keep sound approximation without introducing too much pessimism.

Validated Numerical Integration Method

The use of set-membership computation for the problem described in Section 4.1 makes possible the design of an inclusion function for the computation of [y] (t; [y 0 ]) which is an outer approximation of y(t; [y 0 ]) defined in [START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF]. To do so, a sequence of time instants t 1 , . . . , t n such that t 1 < • • • < t n = t and a sequences of boxes

[y 1 ] , . . . , [y n ] such that y(t i+1 ; [y i ]) ⊆ [y i+1 ], ∀i ∈ [0, n -1] are computed. From [y i ], computing the box [y i+1 ] is a classical 2-step method (see [3]): Phase 1 compute an a priori enclosure [y(ξ)] of the set {y(t k ; y i ) | t k ∈ [t i , t i+1 ] , y i ∈ [y i ]} such that y(t k ; [y i ])
is guaranteed to exist, Phase 2 compute a tight enclosure of the solution [y i+1 ] at time t i+1 .

Validated Computation of the Local Truncation Error

In summary, our approach to design a validated Runge-Kutta method (see [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF]) is defined by

     [k i ] = f ([y 0 ] + [Y i ]), i = 1, . . . , s [Y i ] = h s j=1 a ij [k j ] , i = 1, . . . , s [y 1 ] = [y 0 ] + h s j=1 b j [k j ] + LTE(t, [y(ξ)]) (19) with [y(ξ)] ⊇ {y(t k ; y i ) | t k ∈ [t i , t i+1 ] , y i ∈ [y i ]} computed in Phase 1.
it is bounded following classical approach in validated numerical integration methods. More precisely, a variant of the Picard operator, see [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF], is used in combination with affine arithmetic.

Note that in case of implicit Runge-Kutta methods, the equations [k i ], for i = 1, . . . , s, form a contracting system of equations. In consequence, we can easily build an interval contractor from the system of k i and so we can solve it easily.

Independently to the given IVP-ODE, computing the LTE of Runge-Kutta methods always requires the same amount of elementary differentials which depends on the order of the used method. As described in [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF], we use a structure of weighted directed acyclic graphs (wDAG) of rooted trees to represent the elementary differentials. Since a method of order p requires to produce the rooted trees of order p + 1 (see Eq. ( 5)), the wDAG representing these rooted trees is produced.

As an example, we represent in Figure 1 the DAG of the rooted trees of order 4, mandatory to the computation of the LTE for methods of order 3. Each node of the wDAG will represent a rooted tree and each edge from it is a subtree of the one represented by the node. The weighting represents the number of identical subtrees of a rooted tree. Once the dynamical system is known, the LTE is computed using the wDAG in a depth first search manner for the computation of the elementary differentials. The structure of wDAG makes possible the memoization of the result of the intermediate elementary differentials in order to not recompute it for others rooted trees.

For example, in Figure 1, both first and second rooted trees of order 4 require the computation of the one of order 2. The same for the rooted tree of order 1 which is required for any elementary differential computation. In this case, the corresponding elementary differential is only computed once.

Each node of the wDAG representing a tree τ will also contain all the required values for the computation of the LTE (using Equation ( 5)), namely r(τ ), σ(τ ) and γ(τ ) (see Theorem 2.1) as well as the elementary weight φ(τ ) (see Definition 2). The computation of these associated values do not depend on the considered system so they can be computed offline.

Eventually, a validated computation of the LTE on the interval [y(ξ)] uses the computation of For a given rooted tree, the computation of its associated elementary differential strongly rely on the computation of the one for its subtree and therefore generate a lot of redundant variables which is a well known trouble for interval analysis to produce tights outer approximations. It is why affine arithmetic is preferred to deal with this issue.

Bounding Taylor coefficients of univariate functions

For each node of the wDAG, its corresponding elementary differential must be computed. As seen in Eq. ( 8) it corresponds to the computation of the Taylor coefficients of a unary function. To achieve this, we use the results in [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF] using automatic differentiation. The first step is to produce the computational graph of the function decomposing the function to its atomic operations (such as addition, multiplication, sinus, exponential, etc.). The method described in [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF] to compute the Taylor coefficients of these operations are extended for the computation using affine arithmetic. It results in the computation of an outer approximation of the Taylor coefficients of these considered operations. Example 6. Let f (x) = x sin(x). We want to compute the Taylor coefficients up to degree d for the variable x for which its set of possible values is represented by the affine form x. The corresponding computational graph of f is then the one described in figure 2. We extend the formulae provided by Table 10.2 in [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF] for the computation of the Taylor coefficients of sin(x) and cos(x) using the affine form x. The Taylor coefficients of the leaves of the tree which are x are then (x, 1, 0, 0, . . . ). For the sinus function, its k-th Taylor coefficient is

ŝk = k j=1 j xj ĉk-j k with ĉk = k j=1 -j xj ŝk-j k (20)
and for the multiplication of two affine forms û and ŵ, the k-th Taylor coefficient is k j=0 ûj ŵk-j with ûi and ŵi the i-th coefficient of û and ŵ respectively.

Algorithmic Complexity Analysis

The complexity of the algorithm is hard to expose but, thanks to the previous works [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF][START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF], a hint of it can be produced. We focus on the computation of the elementary differential of a tree given by [START_REF] Marciniak | On representations of coefficients in implicit interval methods of runge-kutta type[END_REF]. The knowledge of the time complexity to compute this can be separated in the computation of its different parts. Firstly it requires to produce the seed matrix S that is the computation of the elementary differentials of all the sub-tree of the considered tree. If we consider this computation for the rooted trees of order d, as for example the one represented in Figure 1 for d = 4, the elementary differential at each node is computed only once and recorded in the structure of the wDAG thanks to memoization. Then, the total number of computation of elementary differentials equals the number of non-isomorphic rooted trees τ with r(τ ) d, that is, one for each vertex of the wDAG. In [START_REF] Otter | The number of trees[END_REF], An approximation of the number of non-isomorphic rooted trees τ with r(τ ) = d is given:

ρ(d) = d m=1 βα 3/2 2 √ π α -m m 3/2 (21)
with α ≈ 0.3383219 and β ≈ 7.924780. The number of vertices in a wDAG can then be approximated by d i=1 ρ(i). Secondly, A. Griewank in [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] provides the time complexity of computing a d-th univariate Taylor component which is in O(nd 2 ). Thirdly, this computation has to be done each time a coefficient ϕ(i, j) is different than zero. As described in [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate taylor series[END_REF], this number for i, j ∈ N n 0 is less or equal than

p(d, n) = d m=1 n m d m m + d -1 d . (22)
Finally, it results that the complexity for computing (8) is in 23given dimensions n and orders d. Bold values mean a worse time complexity than the symbolic approach and a better one otherwise.

τ 1<|τ | d p(|τ |, e(τ ))O(n|τ | 2 ) (23
In Table 1 is shown the time complexity for different values of the dimension of the system and the order of the approach. These computation directly uses the generation of the corresponding wDAG to evaluate precisely the number of non-isomorphic trees of a given order and the number of subtrees e(τ ) of a given tree τ but uses the outer approximation p(d, n) for the number of non vanishing coefficients γ(i, j). These result are compared to the time complexity of the symbolic approach which is in O(n d ). From these results one can deduce that the time complexity is between O(ne d ) and O(n3 d ). The time complexity is experimentally examined in the next section to assess the improvement of the automatic differentiation approach on time complexity compared to the symbolic one.

Experiments

Our automatic differentiation approach has been implemented in the library DynIbex: a tool for Constraint Reasoning with Differential Equations using the library IBEX as it was the case for the symbolic approach1 . In this Section, this implementation is used to compare the symbolic approach described in [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF] with the new one. Two experiments were designed, one to show the impact of the dimension of the state space on the computation time and another one on the impact of the order of the Runge-Kutta method also on the computation time. Both experiments use the same problem that is the numerical integration of a water tank system as introduced in [START_REF] Althoff | Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization[END_REF]. This system can be modeled as following using Toricelli's law:

                                                               ẏ(t) =             0.1 + κ(4 -y n (t)) -k 1 2gy 1 (t) k 1 2gy 1 (t) -k 2 2gy 2 (t) . . . k i-1 2gy i-1 (t) -k i 2gy i (t) . . . k n-1 2gy n-1 (t) -k n 2gy n (t)             y(0) ∈                 [1.9,
                . ( 24
)
The parameter n represents the number of tanks that is considered, y i (t) is the water level of tank i at time t; v ∈ [-0.005, 0.005], κ = 0.01; k i ∈ [0.0149, 0.015] are tank specific parameters and g is the gravity constant. AD approach Symbolic approach Figure 3: Time computation for different state space dimensions for example in [START_REF] Wilhelm | Beitrag zur naherungsweisen integration totaler differential gleichungen[END_REF]. n = 6, . . . , 50.

Computation time against state space dimension

For this first experimentation, example in ( 24) is taken for different state space dimensions n = 6, . . . , 50. The Runge-Kutta method which is in use is fixed as an explicit one: the RK4 method (Cf. Appendix A.1 for the Butcher tableau of the method). The results are presented in Figure 3 for the total time computation of the solution of y(t) at time t = 400s for each dimension n. This illustrates the tremendous improvement of the automatic differentiation approach compared to the symbolic one when the dimension of the state space is sufficiently large. On our example, the benefit of using the automatic differentiation approach starts at n = 6.

Computation time against method order

We now exhibit the relation between time computation and the order of the Runge-Kutta method that is used. The example that is taken is still the one described in [START_REF] Wilhelm | Beitrag zur naherungsweisen integration totaler differential gleichungen[END_REF] for a state dimension n = 3. Figure 4 shows the results for explicit methods and Figure 5 for implicit methods. For explicit methods, this experiment shows the improvement in computation time for the automatic differentiation approach compared to the symbolic one. In the case of implicit methods, the use of a contractor based method to solve the contracting system of equations (see Section 4.3) explains the difference in computation time against explicit methods. This shows an even better improvement compared to the symbolic approach. The Butcher tableau of the different implicit and explicit methods used for these experiments are described in Appendix A.

Conclusion

In this article is presented an automatic differentiation approach to produce validated bounds of the local truncation errors of explicit and implicit Runge-Kutta methods. It is based on the work in [START_REF] Bartha | Computing of B-series by automatic differentiation[END_REF] on the computation of B-series. This makes possible the design of a validated numerical integration of the initial value problem of ordinary differential equations using automatic differentiation and set-membership computation. It has been implemented in a C++ library used to show the benefits of this automatic differentiation approach compared to a symbolic approach previously published by the authors [START_REF] Alexandre | Validated explicit and implicit runge-kutta methods[END_REF]. This was illustrated for space dimensions of the ODE problem as well as the considered Runge-Kutta method order with several experiments showing the decrease of time complexity. AD approach Symbolic approach Figure 5: Time computation for different implicit methods for example in [START_REF] Wilhelm | Beitrag zur naherungsweisen integration totaler differential gleichungen[END_REF]. n = 3.

Theorem 4 . 1 (

 41 Mean value extension (a.k.a. centered form[START_REF] Moore | Introduction to interval analysis[END_REF])). Let the function f : R n

3 Figure 1 :

 31 Figure 1: Representation of the rooted trees of order 4 as a wDAG.

  |j|=d γ(i, j) [F ] |i| ([y(ξ)] ; [S] j) with [F ] |i| an inclusion function for F |i| (y(ξ); Sj) the |i|-th Taylor coefficient of f (y(ξ)+Sj) and the computation [F ] |i| ([y(ξ)] ; [Sj]) an outer approximation bounding the set {F |i| (y(ξ); Sj) | y(ξ) ∈ [y(ξ)] , S ∈ [S]} which is the set of |i|-th Taylor coefficients of f (y(ξ) + Sj), y(ξ) ∈ [y(ξ)] and S ∈ [S]. The interval matrix [S] ⊇ {[F(τ 1 )(y(ξ)), . . . , F(τ m )(y(ξ))] | y(ξ) ∈ [y(ξ)]}.

Figure 2 :

 2 Figure 2: Computational graph for the function f (x) = x sin(x).

Figure 4 :

 4 Figure4: Time computation for different explicit methods for example in[START_REF] Wilhelm | Beitrag zur naherungsweisen integration totaler differential gleichungen[END_REF]. n = 3.

Table 1 :

 1 Time complexity described in Equation

	d\n	2	4	6	8	10	20	30	40
	1	6	12	18	24	30	60	90	120
	2	18	36	54	72	90

) with e(τ ) the number of resulting sub-trees if the root of τ is removed. It correspond to the sum over each rooted tree τ of the number of non vanishing coefficients ϕ(i, j) multiplied by the time complexity of producing the d-th univariate Taylor component.

The library is available at the website http://perso.ensta-paristech.fr/~chapoutot/dynibex/

A Runge-Kutta methods used in the experiments

In this appendix are provided the Butcher tableaux of the Runge-Kutta methods that have been used in the experiments in Section 5.

A.1 Explicit methods


Shlomo Barak
email: shlomo@tagapro.com

Electromagnetism as the Geometrodynamics of Space

Keywords: Electric charge, Space lattice, Electromagnetism

. This, surprisingly enough, yields the Maxwell theory of electrostatics, with no phenomenology. Adopting the idea that space is a lattice (cellular structure), Space Density  is defined as the number of space cells per unit volume (denoted  0 for space with no deformations). Based on this we define Electric Charge Density as: q=1/4π • (ρ-ρ 0 )/ρ. This charge density is positive if  >  0 and negative if  <  0 . This is all we need to derive electrostatics. In this paper, by adding the Lorentz Transformation to the definition of q, we derive the entire Maxwell Electromagnetic theory with no phenomenology.

Introduction

The essence of electric charge has been a mystery. Recently it has been resolved in [START_REF] Barak | On the Essence of Electric Charge, Part 1: Charge as Deformed Space[END_REF] and [START_REF] Barak | On the Essence of Electric Charge, Part 2: How Charge Curves Space[END_REF]. Our definition of charge density, in [START_REF] Barak | On the Essence of Electric Charge, Part 1: Charge as Deformed Space[END_REF], yields the Maxwell theory of electrostatics with no phenomenology. It also realizes Einstein's vision that in "future physics" there will be no physical separation between a charge and its field. In this paper we derive the entire theory of Electromagnetism by adding the Lorentz Transformation to the definition of.

By relating the field energy density to charge density our theory becomes non-linear. This is the specificity of our theory, which is applicable both for an elementary charge as it is for an ensemble of elementary charges.

Note that reading our papers [START_REF] Barak | On the Essence of Electric Charge, Part 1: Charge as Deformed Space[END_REF] [2] [3] [START_REF] Barak | On the Essence of Gravitation and Inertia, Part 2: The Curving of Space by an Elementary Particle[END_REF] is helpful for a full understanding of this paper.

The Electric Field of a Moving Charge

This section is based on [START_REF] Purcell | [END_REF]. Let a point charge Q reside at the origin of a frame k, Fig. [START_REF] Barak | On the Essence of Electric Charge, Part 1: Charge as Deformed Space[END_REF].

The electric field E has the magnitude HQ/r 2 , and for a positive charge it is directed radially outward.

In the xz plane its components at any point (x, z) are: The frame k' moves in the negative x direction, with speed v. According to the Lorentz transformation, the relations between the coordinates of an event (a spacetime point) in the two frames are:

    2
                  c ' x ' t t ' z z ' y y ' ct ' x x (2) x x z z E ' E and E ' E   
. We express the field components x z ' E and ' E in terms of the coordinates in k' using equations ( 1) and [START_REF] Barak | On the Essence of Electric Charge, Part 2: How Charge Curves Space[END_REF]. At the instant t' = 0, when Q passes the origin in k':

z ' z k ' k Q E Z E E X ( x , z ) r  x x ' v ( a ) z ' k ' Q E ' Z E ' E ' X r ' ' x ' v a t t ' = 0 ( b )
        2 3 2 2 z z 2 3 2 2 x x z' x' HQz' E E' z' x' HQx' E E'            (3)
Note first that E' z /E' x = z'/x'. This tells us that the vector E' makes the same angle with the x' axis as does the radius vector r'. Hence E' points radially outward along a line drawn from the instantaneous position of Q, as shown in Fig. (1b).

To find the strength of the field, we calculate

 

            3 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 4 2 2 2 2 3 2 2 2 2 2 2 2 2 z 2 x 2 z' x' z' 1 z' x' 1 Q H z' z' x' z' x' Q H z' x' z' x' Q H E' E' E'                             (4)
Let r' denote the distance from the charge Q, which is momentarily at the origin, to the point (x',z') where the field is measured:

  2 1 2 2 z' x' r'  
. Let ' denote the angle between this radius vector and the velocity of the charge Q, which is moving in the positive x' direction in the frame k'. Then since z' = r' sin ', the magnitude of the field can be written as:

  2 3 2 2 2 2 ' sin 1 1 r' HQ E'       (5)

Fig. (2) Representation of the Field of a Uniformly Moving Charge

For low speeds the field reduces simply to E'  HQ/r' 2 . But if  2 is not negligible, at the same distance from the charge, the field is stronger at right angles to the motion than in the direction of the motion, A simple representation of the field is shown in Fig. [START_REF] Barak | On the Essence of Electric Charge, Part 2: How Charge Curves Space[END_REF]. A crosssection through the field with some field lines in the x' z' plane is indicated. For the field in the x'y' plane we get an identical representation. [START_REF] Purcell | [END_REF]: This is a remarkable electric field. It is not spherically symmetrical, which is not surprising because in this frame there is a preferred direction, the direction of motion of the charge. Also, it is a field that no stationary charge distribution, whatever its form, could produce. For in this field the line integral of E' is not zero around every closed path. Consider, for example, the closed path ABCD in Fig. [START_REF] Barak | On the Essence of Electric Charge, Part 2: How Charge Curves Space[END_REF]. The circular arcs contribute nothing to the line integral, being perpendicular to the field; on the radial sections, the field is stronger along BC than along DA, so the circulation of E' on this path is not zero. 2) is a natural result of the charge being a longitudinal wavepacket. As such it becomes contracted in the direction of motion, [START_REF] Barak | On the Essence of Gravitation and Inertia, Part 1: Inertia and Free Fall of an Elementary Particle[END_REF] [6] [START_REF] Pierce | Almost All About Waves[END_REF]. Charge density is related to space curvature [START_REF] Barak | On the Essence of Electric Charge, Part 2: How Charge Curves Space[END_REF], hence the strong field perpendicular to the direction of motion.

3

Forces between Moving Charges

We consider, [START_REF] Resnick | Introduction to Special Relativity[END_REF], two particles of equal charge, Q, moving with equal uniform velocity v. Let v be along the x-axis of frame k, and let the particles have the same x-coordinate, and their separation be r, see Fig. (3a). We expect the charges to exert forces upon one another, a repulsive electric force and an attractive magnetic force, as we show.

Let k' move with uniform velocity v relative to k along common x-x' axis, so that the charges are at rest in k', see Fig.

Fig. (3) The Magnetic Field of a Moving Point Charge

v v K y x 0 Q Q (a) v r K' y' x' 0' Q Q (b) r
Here, there is no magnetic force at all and the electric force is repulsive: the charges would tend to move apart along the y'direction, each exerting a force on the other, of magnitude:

2 2 y r HQ F'  (6)
The force on the upper charge is + F y ' and that on the lower charge is -F' y . Note that charge invariance is assumed and that the separation r is unchanged by the transformation from k to k'. We now use the general force transformation equations (Equations. 3-33 in [START_REF] Resnick | Introduction to Special Relativity[END_REF]) for a force with components F x ', F y ' and F z ', in the S' frame, to obtain the force components F x , F y and F z acting on the particle in the k-frame.

F x = F x ' = 0 F z = F z '/ = 0 and: 2 2 2 2 2 2 2 2 y y y r Q H /c v 1 r HQ /c v 1 ' F / ' F F         (7)
Hence, the net force is in the positive y-direction, but smaller than in [START_REF] Lorentz | Electromagnetic Phenomena In A System Moving With Any Velocity Less Than That Of Light[END_REF]. The charged particles repel one another and the coulomb electric force of repulsion must exceed the force of attraction which we call the magnetic force.

For v < c the electric force is always greater than the magnetic force. Only when v = c does the net force become zero.

When v  0, we return to the static result wherein only an electric force exists.

The magnetic force, which exists only when v  0, is a second-order effect compared to the electric force, that is, it enters as (v/c) 2 .

The Magnetic Field B, the Vector Potential A and the Lorentz Force

The Magnetic Field 4.1

The electric field of a moving charge at velocity v, perpendicular to v is: This field applies a repulsive force ' y F on another charge Q that runs in parallel, which is:

2 2 y r HQ F   
On the other hand, equation ( 7), shows that the net force is smaller:

2 2 y r Q H F   (7)
where:

y y y F F F     
and we want to find the attractive force y F  :

                                             QE β r HQ γβ y F 2 2 2 2      
and for the net force:

  2 2 2 2 2 2 2 β 1 QE E Qβ QE r HQ γβ r HQ γ y F      
This result is for v  E, but we can show that the general expression for the attractive force, F  , is:

  2 c Q E v v F      (8)
Thus the net force F is:

            2 c Q E v v E F (9)
For different velocities v 1 and v 2 of the charged particles, we can show that in general:

            2 1 2 c Q E v v E F (10)
where v 1 is the velocity of the particle that creates the field E and v 2 is the velocity of the other particle in this field, and vice versa. The expression:

  E v B    c 1 Magnetic Field (11)
is defined as the magnetic field created by a particle moving at speed v ' .

The dimensions of B and E are the same: [E] = [B] = LT -2 , see [START_REF] Barak | On the Essence of Electric Charge, Part 1: Charge as Deformed Space[END_REF], and both express the elastic displacement.

5

The Lorentz Force

Using [START_REF] Konopinski | [END_REF] and [START_REF] Spiegel | Vector Analysis, Shaum's Series[END_REF] gives the force on a charge Q moving at speed v:

           B v E F c 1 Q Lorentz Force (12)
6

The Vector Potential A

The magnetic field is merely the displacement vector u created by a moving charge. In this case, u is not considered as a gradient of a scalar function , namely a polar vector, but an axial vector which is a rotor of a vector potential A. In the static case: E = - whereas, in the general case, we have to add the axial component:

B =  x A [A] = [] = L 2 T -2 Note that: u = (1/H) (E+B) (13)
According to [START_REF] Spiegel | Vector Analysis, Shaum's Series[END_REF]:

  E v A     c 1  (14)
and hence:

r c H r Q c H i v A    . ( 15
) v i   Q [i] = L 4 T -1 and v j   q [j] = LT -1 ( 16
)
are defined as the electric current and the current density respectively. Hence:

        ' d c H r r j A (17)     dσ j i (18)
If the electric current flows in a closed circle the displacement B creates a torsional space deformation and, hence, space contraction. This phenomenon is related to gravitation [START_REF] Barak | On the Essence of Gravitation and Inertia, Part 1: Inertia and Free Fall of an Elementary Particle[END_REF] and [START_REF] Barak | On the Essence of Gravitation and Inertia, Part 2: The Curving of Space by an Elementary Particle[END_REF]. Note the relation of the current to the helical motion of a positron (or that of an electron moving in the opposite direction).This helical motion creates B, which is a component of the displacement vector u.

Examples of B and A Created by Constant Currents

Fig. (5) The Vector Potential A

The Current Equation of Continuity

We have seen, (18), that:

    dσ j i (18)
For a given zone , with a surface , in which there is a charge Q we get:

      d t q dt d dt dQ , but:       d j dt dQ n ( 19
)
where n is the unit vector perpendicular to the surface  at a given point.

According to Gauss Lemma:

A j                   d t q dτ d j n j
This implies the equality of the integrands, namely:

0 t q      j  Current Equation of Continuity (20)
In the GDM, both j and q include the contribution of the field energy.

8

Ampère's Law

      ' ' d ' c H r' r j A (17)
therefore:

         ' ' d ' c H r' r j A B   (   B = 0 ) (21)         ' d ' c H 2 r' r r' r j B (22)
but:

       r' r j B ' d ' c H    (23)
The identity x(xA) = -(A) - 2

A for an arbitrary vector A, transforms (23) into:

                       ' d ' 1 ' H ' d ' 1 ' H r r r j r r r j B      ( 24
)
If we use, see [START_REF] Jackson | Classical Electrodynamics[END_REF], the relations:

                    ' 1 ' 1 ' r r r r   and:   ' 4 ' 1 r r r r              
the integrals in (24) can be written:

    r j r' r 1 r j                 4 ' d ' ' H    
Integration by parts yields:

           ' d ' ' H H 4 r' r r j j    
But for steady-state magnetic phenomena   j = 0 , hence:

j B H 4    Maxwell equation (25)
From ( 25), by applying Stokes' theorem, we obtain the integral equation:

           dr j dσ Β H 4  ( 26
)
and transforming (26

) into:        c H 4  d j dl B
gives:

    c H 4 i dl B Ampère's law (27)
Maxwell "fixed up" (25) by adding to it the displacement current

t c 1  E : t c 1 c H 4       E j B 

Maxwell equation

(28)

9

Comments on the Electric Field E, Magnetic Field B, Potential Vector A and Lorentz Force

Equation ( 13), notated here (29), is:

u = (1/H) (E+B) (29)
E expresses the polar displacement and B the axial displacement. This is compatible with E and B being the components of the four-field anti-symmetric EM tensor:

k ijk ij i oi B F E F     (30)
where ε is the anti-symmetric Levi-Civita symbol , see [START_REF] Jackson | Classical Electrodynamics[END_REF].

                      0 B B E B 0 B E B B 0 E E E E 0 F 1 2 3 1 3 2 2 3 1 3 2 1 (31)
Space has an interesting feature: the vector potential A expresses a kind of potential momentum, that creates an axial displacement

A B   
. This situation is similar to the wellknown behavior of a spring: load it and it twists (the radius increases), see Fig. [START_REF] Lorentz | Electromagnetic Phenomena In A System Moving With Any Velocity Less Than That Of Light[END_REF].

Fig. (6) The Vector Potential as a Potential Momentum

Note that equation [START_REF] Barak | The Geometrodynamics of Space hal-01435685v2[END_REF], which here is notated (32), can be written as (33):

         B v E F c Q Lorentz force (32)
F              A v   c Q (33) 10

Faraday's Law of Induction

When currents change, there is also a change in A and B. But a change in B means a change in the displacement u. In this case, a charge is not displaced because of a direct force, but because of a displacement of the zone of space in which it is immersed. In other words, the background in which the charge is immersed moves.

In this case, we should look upon qA as the momentum potential density in space of the M charge density, q, immersed in it, and on the time derivative of qA as the force required to move it. A good guess of the force density is, therefore, the expression:

t A c q q         (34)
and for the force:

t c Q Q       A F  (35)
F = QE, and therefore:

t A c 1 E        (36)
By applying the operator  to both sides, we get the:

t c 1      B E  Maxwell equation (37)
From the Maxwell differential equation, we get the integral relation:

  c σ      E d E dl  (38)
and on the other hand:

                d B d E t c 1 ( 39
)
and therefore:

∮ E•dl= - 1 c c ∂ ∂t ∫ B•dσ σ Faraday's law of induction (40)      Bd  (41)
where  is the magnetic flux through the surface .

The Vector Potential as the Momentum Potential in Space

We show that A c Q can be considered as the potential momentum of the electromagnetic field, i.e., the momentum stored in space, momentum available for conversion to kinetic energy, see [START_REF] Konopinski | [END_REF].

By substituting:

A B    and: t c 1       A E 
in the Lorentz law of force:

B v E F    c Q Q we obtain: A v A F            c Q t c Q Q . ( 42
)
Note identity [START_REF] Resnick | Introduction to Special Relativity[END_REF] in the book [START_REF] Spiegel | Vector Analysis, Shaum's Series[END_REF]:

          v A A A v v A A v A v                (43)
Assuming that Q moves at a constant velocity, v, in a field with a uniform A, the second, third and fourth terms in (43) equal zero. We therefore get:

    A v A v       (44)
Substituting (44) in (42) enables the reorganization of (42) such that:

                  A v A v c 1 Q c Q M dt d  (45)
We can relate to (45) as follows:

The force is the time derivative of the generalized momentum:

        A v P Q c 1 M
that is equal to the gradient of the potential energy:

          A v c 1 Q U hence: U t    

P

Let us adopt the known expressions:

A v c Q M  Conjugate Momentum (46)          A v c 1 Q Interaction Energy (47)
To the Potential Energy:

Qφ ( 48
)
we can add Potential Momentum:

A c Q (49)
Which, in quantum theory, is related to the operator  i  .

12

The Electromagnetic Equations Classical Electromagnetism 12.1 Electromagnetic theory can be fully expressed by the four Maxwell's equations:

I q 4   E  Coulomb's Law see [1] II t c 1      B E  Faraday's Law (37) III B   = 0 No magnetic monopoles IV t c 1 c 4       E j B  Ampère's Law (28)
These equations can also be expressed in terms of differential equations for the scalar and vector potentials:

 q 4    (48)  j A c 4   (49)
Using the potential 4-vector

  A , A    and the current density 4-vector     j j cq, , j j 0    we get:       j c H 4 A
as the potential 4-vector equation. Since:

F = A -A      
the covariant Maxwell's equation can also be written as:

       j c 4 F (50)

Gauge Transformation 12.2

The correspondence between F μ and A μ is not unique. Maxwell's equations are invariant under the gauge transformation:

A μ → A μ ' = A μ -∂ μ ψ ( 51
)
where ψ is an arbitrary spacetime dependent scalar function (the gauge function). A change of A μ does not change the EM representation by E and B fields. The relevant electric field E, in this case, is t

A c 1     E
and the magnetic field is simply B =  x A.

13

The GDM ElectroMagnetic Equations

The Equivalence of Energy Density and Charge Density 13.1

In the GDM charge density is q=1/4π • (ρ-ρ 0 )/ρ. This charge density is positive if  >  0 and negative if  <  0 . But space density  depends on the space energy density ϵ , [START_REF] Barak | The Geometrodynamics of Space hal-01435685v2[END_REF].

GDM EM equations 13.2

To obtain the GDM EM non-linear equations we use Maxwell's equations and modify the expression     v j q , cq cq, j    by adding the contribution of the field energy (this, however, is out of the scope of this paper).

  2 2 source s H 8 q q      
, where H = 1 [H] = T -2 see [START_REF] Barak | On the Essence of Electric Charge, Part 1: Charge as Deformed Space[END_REF], [   ] = LT -2 , s = 1 , [s] = LT -1 , see [START_REF] Barak | On the Essence of Electric Charge, Part 2: How Charge Curves Space[END_REF] and [q] = 1.

Thus both gauge invariance and charge conjugation are lost.

On Charge Conjugation and Gauge Invariance 13.3

The GDM-modified classical potential field equation includes the contributions of the energy density and the current of energy density of the field, which are equivalent to the charge density and current density, respectively. Maxwell's equations are the approximated equations of GDM Electromagnetism where the energy density of the field is neglected. The electric charge, Q, of an elementary particle is not a constant. As the distance from the particle increases, the absolute value of the charge decreases. E is, therefore, proportional not only to1/r 2 but also to the changing Q.

In the field there is a change in space density that affects c. Thus, c depends on the potential as in the case of gravitation, and we can therefore expect light bending. However, this bending is different for positive and negative charges. Thus, there is no charge conjugation, nor gauge invariance close to a charged particle.

Charge conjugation: the results of experiments are independent of alternating the sign of all charges.

Gauge invariance: the results of experiments are independent of the choice of the gauge for the potentials.

Fig. ( 1 )

 1 Fig. (1) The Electric Field of a Point Charge, Q

  the square of the magnitude of the field, E' 2 .

Fig. (

  Fig. (2) is a natural result of the charge being a longitudinal wavepacket. As such it becomes

  (3b).

Fig. ( 4 )Fig. ( 5

 45 Fig. (4) The Vector Potential A

Acknowledgment

We would like to thank Mr. Roger M. Kaye for his linguistic contribution and technical assistance.

14

Summary

In the classical Maxwell EM, space does not play any role, whereas in our theory EM is merely the Geometrodynamics of space. This new understanding enables the turning of the form of the equations from linear into non-linear, which is also the case for QED. This is the specificity of the theory. The new understanding adds tangibility and enables the derivation of the theory from merely the definition of charge density and the Lorentz Transformation.


domain: spi scco

text:
Pascal Sire

Gilles Haeffelé

Sebastien Dubois

TRIZ as a tool to develop a TRIZ educational method by learning it

Keywords: TRIZ, education system, teaching method development, innovative evaluation

Since last year, in France, an official education program includes more than 50 hours of innovation methodology and training mainly based on TRIZ methods! But the problem is that teachers of the classes (sciences and technology for engineers and sustainable development) do not have any training on TRIZ, or will only receive it in the next years... So the contradiction is that teachers have to be taught about TRIZ methods to be able to teach it, but teachers could not be taught on TRIZ methods as there are no time resources to do it. With the aid of an external TRIZ expert, it has been proposed to use as a main resource a team of 6 teachers and the time of course preparation to both teach teachers and develop a new, innovative teaching method. This experiment started this year for 2 classes to inventively teach TRIZ to teachers "in real-time" while students learn it by practicing the innovation methodology and creativity tools through a project. The paper describes the 3-level approach to "develop a new educational method", to "teach the teachers" and to enable students to learn through a "creativity challenge for innovation". The identified key problems related to this new educational method and the results of this new approach will be presented. As results, one can notice a new way of evaluating where everybody wins (education system, teachers, students), based on contradictory parameters in a "2D radar" highlighting students' inventive abilities. From initial state to Ideal Final result, the paper explores the key success factors, the next steps and the set of problems linked to the implied actors: teachers and students, and also the educational system. The originality of the approach is that typical solutions used in learning with industries have been applied in the educational system, whereas it is traditionally the opposite…

Introduction

Since last year, in France, an official education program includes more than 50 hours of innovation methodology and training mainly based on TRIZ methods! But the problem is that teachers of the classes (sciences and technology for engineers and sustainable development) do not have any training on TRIZ, or will only receive it in the next years... Thus school teachers faced a contradiction: they had to teach methods for which they did not have any knowledge. This initial problem has in fact been recognized as a set of problems teachers were facing: teaching without knowing, defining a program without having the required knowledge, and so on. Fortunately, the staff of teachers met a TRIZ expert, Pascal Sire, glad to help with the transmission of TRIZ knowledge to young students. Based on the pedagogical competences of the teachers, and based on TRIZ expertise, it was then necessary to define the educational method to transmit TRIZ know-how on problem resolution to students. This educational method has to take into account students' specificities, knowing that they are young students, without the normal technical background that can be expected from an engineering student.

How to enable teachers to be taught at the same time that they had to teach students?

Teachers do not have time to learn, and moreover they do not have any resources to obtain knowledge from, nor TRIZ-based methods prior to teaching them to their students. So their only useful resource was the expertise of the TRIZ expert. This totally new situation led to the definition of an educational method at two levels: a new way to teach teachers and a new way to transfer know-how to students.

The teachers' training had to enable teachers to always keep a step ahead in terms of acquiring methods in comparison with the students. This was partially solved by the fact that teachers already had a technical background, already knew about systemic thinking, and about function description. So students had to familiarize themselves with these new approaches and moreover to acquire know-how on problemsolving methods in the context of design. But teachers also had to take into account that all these differences (methods, models of description, and way of thinking) could lead to too big a cognitive gap between them and the students. So three levels of expertise were identified:

1. Students, with no particular knowledge in design, no technical knowledge, and not involved in problem solving methods training sessions; 2. Teachers as technical experts also involved in the process of the problem-solving educational methods; 3. TRIZ expert, only focused on transmission of problem-solving methods. In this particular situation, one favorable condition was the fact that teachers had to practice and learn new knowledge in the same context, so at least they did not face one of teaching difficulty. They generally acquire knowledge in a particular passive context and had to transmit the knowledge in a frontfacing active context.

What to propose to students?

For the first time, students were also facing a new situation; they would not have to learn knowledge, as they were used to, but to develop skills, know-how, moreover in a context they are not familiar with: design situation. Designing requires technical knowledge and the development of intelligent thinking. Intelligent thinking is recognized by [START_REF] Perkins | Thinking Frames[END_REF] as a sum of tactics, power and content. "The Power component relates to those elements of intelligence that depend on neurological efficiency. The Content component of intelligence expresses the need for a strong knowledge base in the discussed domain, such as mathematics, physics, or technology. The Tactics component consists of strategies, tactics, and techniques used in solving a problem." [START_REF] Perkins | Thinking Frames[END_REF] in [START_REF] Barak | Idea focusing versus idea generating: a course for teachers on inventive problem solving[END_REF] The main objective of the sessions defined with students was not to increase the content aspect, as it was neither an objective to feed students with technical knowledge, nor to change their neurological efficiency, so it was limited to give them efficient strategies, tactics to face and to solve problems. To enable them to solve problems, without spending time on technical knowledge, it was decided to process in project conditions, letting students choose their problems, choose the object they aimed to design, and thus letting them look by themselves for the knowledge they required. Apart from time management, this project approach also had a second benefit: to motivate the students, as they totally appropriated their projects and this let them be proactive in the education process.

How to evaluate when no predefined solution exists?

A last challenge in this project was to find a way to motivate students, and enable them not to be blocked or afraid by the evaluation. So the evaluation process had to enable the promotion of each student's personality.

Another aspect of the evaluation process is linked to the staff of teachers. For them it is also difficult to accept that the problem they proposed to their students does not have one a priori known solution, and to accept that the evaluation mode has to be redefined.

So a double constraint had to be taken into account concerning the evaluation process: giving enough confidence to the students to face a new situation when for the first time they did not have to reproduce something they were taught to do; and also building confidence for the teachers to be able to evaluate a situation for which there were no predefined and no better solution.

Boden, in [START_REF] Boden | Dimensions of creativity[END_REF], proposes to differentiate H-creative (historical) and P-creative (psychological) ideas. "The H-creative idea brings a creative product as a result related to market or society represented by patent; whereas the P-creative idea brings a creative product as a result to process related to whom solve the problem". As the students were not experts in the domain they were analyzing problems, it was obvious to choose to evaluate creativity as P-creative ideas. But then a trick had to be found not to enable students to propose trivial solutions and also not to be blocked by the lack of expert knowledge.

Known partial solutions

Training without being trained

So teachers were facing a real contradiction, they had to be trained, taught TRIZ-based methods to be able to teach them to their students, but it was not possible for them to be trained and taught as they had no time resources available. With the help of the TRIZ expert, it was decided that the development of a new, innovative pedagogy would be the first application of a TRIZ-based method; and thus that the definition of this pedagogy would be considered as a problem to solve and would be the case of training teachers. So the teachers were trained by the TRIZ expert using the direct application of the TRIZ-based method to their own problem. To be put in the context of learning application is a common and wellknown pedagogy: Problem-Based Learning (PBL), which has been applied successfully in different conditions [START_REF] Murphy | Merging Problem-Based Learning and Simulation as an Innovative Pedagogy in Nurse Education[END_REF], [START_REF] Nakagawa | Education and training of creative problem solving thinking with TRIZ/USIT[END_REF]. This approach is in accordance with the way Dreyfuss described the expert knowledge. In [START_REF] Dreyfuss | La portée philosophique du connexionnisme[END_REF] was described that one has to abandon the traditional idea that a beginner starts with specific cases and then, as he progresses, will bring out by abstraction and internalize increasingly sophisticated rules. On the contrary, the acquisition of know-how would be to move from abstract rules to particular cases.

Another aspect of the pedagogy that had to be developed was to build a trusting relationship between the teachers and the students. For that the students had to feel that teachers had enough knowledge to help them in the project, but at the same time the teachers had to be able to explain their expectations a way that could be understood by the students.

Then a second contradiction was faced: the cognitive distance between students and teachers had to be high to build a trusting relationship, but the cognitive distance between students and teachers had to be low to create a good communication and a good transmission of knowledge.

A third aspect, which is also very important for teachers, is that generally they are trained in a learning context, totally different from the context in which they will have to teach. And so they have to acquire both knowledge and understanding of the studied subject, but also specific know-how and good behavior to convey to the students.

Problem-based learning

PBL is based on the principles of adult learning. Knowles [START_REF] Knowles | The modern practice of adult education[END_REF], the father of adult learning theory, proposed that a learning environment which is characterized by physical comfort, mutual respect and freedom of expression is conducive for adult learning. Differences are accepted, the learners perceive learning goals as their own and accept partial responsibility for planning and conducting the learning sessions and their active participation in the learning process is encouraged. PBL has been applied in many domains, principally in the sphere of medicine [START_REF] Murphy | Merging Problem-Based Learning and Simulation as an Innovative Pedagogy in Nurse Education[END_REF], [START_REF] Shankar | Problem-based learning: A Review[END_REF], [START_REF] Epstein | Learning from the problems of problem-based learning[END_REF], but also in engineering domains [START_REF] Brodeur | Problem-Based Learning in Aerospace Engineering Education[END_REF]. As noted in [START_REF] Murphy | Merging Problem-Based Learning and Simulation as an Innovative Pedagogy in Nurse Education[END_REF], "The combination of PBL and simulation offers students a variety of ways of learning and transforms knowledge in different yet controlled and safe environments. The ultimate goal is to challenge students to think for themselves. This merger cultivates an awareness of creative thinking, critical analysis, and decision-making abilities from extrapolating and relating the theoretical and practical knowledge presented in each cycle."

The main goal of problem-based learning is to provide students with opportunities to apply knowledge, not just acquire it. PBL focuses on problem formulation as well as problem solving. It seeks to simulate real-world engineering research and development. Barrows [11] describes the main features of PBL in this way:

• "Learning is student centered, i.e., students make choices about how and what they want to learn. • Learning occurs in small student groups and promotes collaborative learning.

• Teachers are facilitators or guides or coaches.

• Problems form the organizing focus and stimulus for learning.

• Problems are a vehicle for the development of authentic problem-solving skills.

• New information is acquired through self-directed learning."

The PBL could thus solve problems to develop the pedagogy both for teachers training and for students transfer: with the TRIZ expert playing the role of coach for the teachers and teachers being coaches for the students.

Project-based learning

Another well spread approach to learn to develop professional competencies is the project-based learning (also quoted PBL). Project-based learning forces learners to explore "important and meaningful questions through a process of investigation and collaboration" [START_REF] Frank | Implementing the Project-Based Learning Approach in an Academic Engineering Course[END_REF]. Among the many benefits of project-based learning, the following ones are cited in [START_REF] Thomas | A REVIEW OF RESEARCH ON PROJECT-BASED LEARNING[END_REF] and can be recognized as meaningful for our situation:

• Project-based learning projects are the curriculum, "students learn the central concepts of the discipline via the project" • The projects can be based on an ill-defined problem, for which no a priori solution exists.

• The projects imply a "goal-directed process that involves inquiry, knowledge building, and resolution". • The projects are not scripted, teacher-led, the students build their own approach to the curriculum. Project-based learning and problem-based learning are quite similar in their outcomes and in their multidisciplinary orientation. But some differences are noted in [START_REF] Mills | Engineering Education -Is problem-based or project-based learning the answer?[END_REF], among them:

• Project tasks are closer to professional reality and therefore take a longer period of time than problem-based learning problems • Project work is more directed to the application of knowledge, whereas problem-based learning is more directed to the acquisition of knowledge • Self-direction is stronger in project work, compared with problem-based learning, since the learning process is less directed by the problem As in the situation, the objective is more to enable students to acquire basic knowledge and develop competencies, as the allocated time is quite short and as young students are not able to confront professional situation, the problem-based learning approach seems more adequate.

Developing creativity

As cited in [START_REF] Barak | Teaching methods for inventive problem-solving in junior high school[END_REF] "Creativity is often defined as the ability to produce work that is both novel (original, unexpected, imaginative) and appropriate (useful, adaptive regarding task constraints)". Then the authors explained that to enhance creativity it is necessary to have a set of suitable heuristics that narrow the scope of a problem. In their experiments, they propose to make students work in groups of two to three on designing and constructing an original technological system. The conclusion was that the use of a systematic way to analyze the problem was helpful for students to direct their attention and discovery.

In [START_REF] Johari | Systematic Infusion of Creativity in Engineering Design Courses[END_REF] the authors also demonstrated that training students on inventive techniques such as Synectics, Fermi Approach, TRIZ and so on, helped in increasing their creativity and especially if this training takes place as early as possible.

In [START_REF] Belski | TRIZ course enhances thinking and problem solving skills of engineering students[END_REF] the author also tried to evaluate the benefits of TRIZ learning for students. In an experiment the students were trained on four thinking tools of TRIZ: Situation Analysis, Method of the Ideal Result, Contemporary Substance-Field Analysis, 40 Innovative Principles with the Contradiction Table . The result of analyzing of the impact of such training for students was the increase in confidence for students to solve any problem they will face and especially in tackling unfamiliar problems. In [START_REF] Ogot | Systematic Creativity Methods in Engineering Education : A Learning Styles Perspective[END_REF], "the authors were able to conclude that the introduction of TRIZ does lead to an increase in creative output".

Then, Toru Nakagawa presented in [START_REF] Nakagawa | Education and training of creative problem solving thinking with TRIZ/USIT[END_REF] a way to introduce TRIZ-based learning in education of undergraduate students by the definition of a 3-step program. But this program is dedicated to people with a minimal technical and design knowledge.

Proposed solution

To learn TRIZ, it was decided to start at the beginning in the situation of a concrete project, with the willingness to venture beyond the seemingly intractable problem "teach TRIZ without knowing," and to learn how to solve the problem "teach TRIZ" by solving the initial training "know TRIZ."

In the working group, several teachers intuitively formulated the contradiction in terms of an ideal, which facilitates the change of perspective necessary for their learning. Moreover, these teachers themselves willingly accepted being in situations of the TRIZ project with the expert, just as students would be in the next level with them, under the TRIZ principle known as "Matrioska". As stated by Sokol and al. in [START_REF] Sokol | The development of inventive thinking skills in the upper secondary language classroom[END_REF], "The TA (Thinking Approach) is based on the idea of a non-linear nature of learning and thus non-linear organisation of the learning/teaching process. a linear or cyclic curriculum model the TA offers a modular course based on a number of learning technologies."

The expert in innovation methodology was identified in the close circle of teachers, by using networks of colleagues of teachers with representatives of engineers' associations working regionally with high schools (linkage of the worlds of business and education). He was immediately attracted by this opportunity to "plant a few seeds" and test his concept of the "knowledge spiral" from the business world.

He first met the team of 6 teachers 4 times from October 2011 to January 2012 to teach them basic TRIZ notions and methods. After that introduction, the teachers started to work with the students on the projects 4 hours a week for 5 weeks. During this period, the expert just came one time to help the students and the teachers in the use of TRIZ.

The 60 students worked by groups of 4 to solve problems related to housing and they had to defend their results in front of all the groups and a board of examiners (teachers and experts in innovation).

These students are experimenting a new educational program including creativity for the first time. They are 16-17 years old and will be graduated in 2 years.

"Everybody wins", ideal philosophy of TRIZ!

In the first meetings with the TRIZ expert, once a month during the first semester, teachers are confronted with the major concepts illustrated by specific examples drawn from everyday life to help them to appropriate more easily, which corresponds more to the original request for "evidence of a concrete practice of TRIZ by an industry expert known in the region". As noted by Clarissa S. Thompson [START_REF] Thompson | Powerful pedagogy: Learning from and about teaching in an elementary literacy course[END_REF]

, "In learning in and from practice novices use small pieces of the broader terrain of teaching and learning (examples of practice) to understand more deeply what happens in teaching and learning."

The transfer to the student level was then done in the second half by an introductory lecture (same general principles as those previously learned by the teachers, traditional path to education [Fig. 1]) as a prelude to a "challenge of creativity for innovation" to stimulate learning and motivate students themselves (groups of four students from two classes and with different specialties), rather than a contest that would have discouraged them by an individual competition too common in schools.

As outlined by King and Kovacs in [START_REF] King | New Imagery for Schools and Schooling-Challenging, Creating and Connecting[END_REF], the quality of learning varies from extrinsic to intrinsic motivation. "If we wish to improve the quality of learning, then we need to develop a relationship among teachers and students that is towards the right end of this continuum. In other words, a key feature of improving the quality of learning is to increase the extent to which students are responsible for their own learning and intrinsically motivated to learn. The teacher is a skilled facilitator of learning drawing upon their pedagogical knowledge, knowledge of ways to teach 'content' and their profound knowledge of theory for improvement. To enable students to take on increased responsibility for learning, we must equip them with the skills and capacities necessary for them to do so."

The educational establishment was careful from the beginning of this experiment, through the presence of the regional pedagogical inspector (as early as during the second meeting with teachers), then during the development of a new evaluation method, where he actively encouraged it (please note that the risk was low because the rating of the project creativity will be needed in the formal curriculum of the students only the following year).

"Knowledge spiral", universal attractor?

To channel the whirlwind of innovative ideas, their efficient transmission to students and the broad dissemination of creative methods in the high school, the problem of setting in motion these linked actions was resolved by the close relationship established immediately between the TRIZ expert and the group of teachers and between teachers and their students.

The coupling effect is maximal if the between successive spirals and the time of implementation are adapted in real time [Fig. 2] according to three criteria for self linkage: first "skate and rub", then "join and adhere", finally "follow and empower". The TRIZ evolution laws underlying these criteria are harmonization, conductivity of the energy, and transmission to the super-system.

If the purpose of education is to make students freer to make their choices and independent in their learning, then teachers are placed more involved in the strategies to adopt, obtaining resources, and the necessary impetus for the conservation movement of the whole class of students. An additional benefit: teachers as evaluators are also themselves TRIZ learners!

Dimension change for a positive evaluation

Traditionally, students scoring is done by a single score in the range 0 to 20, leaving little flexibility for teachers to include non-measurable parameters in assessing their students. The TRIZ expert proposed to move from this one-dimensional "linear" dimension to a "surface" measurement, by analogy with the "radar charts" used in the industry to incorporate different characteristics in the evaluation of complex situations.

To assess measurable outcomes of students but also their personal and group qualities, this evaluation "in square meters" can not only take into account several factors, but may also enhance the contradictory qualities that an ideal innovator should possess. The "2D Radar" [Fig. 3] therefore includes opposed pairs of criteria, encouraging each students group to express these qualities present in one of them or their group. As noted by Ferreira and Andrade in [START_REF] Ferreira | Conception of a management tool of technology enhanced learning environments[END_REF], "This type of chart, because of its versatility in knowledge representation, is often used in the analysis of organizational development and measurement of quality. In an academic environment still in transition, where there is an increasing penetration of technology associated with pedagogical changes, radar charts are suitable and adaptable."

As innovation is "an idea with market success", the criteria used after discussion with teachers therefore highlight technical and commercial considerations, internal or external communication, methodology or chance, cost and simplicity. By adding a surprise indicator (the famous "WOW factor") and a "like" (want to buy or follow, now well understood attitude in social networks), this new evaluation system quantifies all aspects of the proposed and evaluated innovation.

Best result of the students

The board of examiners had difficulty choosing the winner due to the quality and variety of the teams' delivery (half of the 14 teams could have been in the first 3 places on the podium), so it was decided to recognize the biggest well-shaped "2D radar" with the highest rate of "WOW" and "like" if a difference was needed! The winner is... "the door in the door" as the team itself summarized their project (in a successful marketing way) in front of their 60+ classmates, after a clear and simple presentation of their study, the problem to solve, the TRIZ process they followed, and the concept of solution they arrived at.

Problem to solve: how to change the principle of a door to always open the door by pushing it (whatever the direction of arrival at the door) without modification of the environment (ie. keeping existing frame of the door to avoid replacement cost)?

Model of problem: the physical contradiction [Fig. 4] is to open the door in both directions, without changing the frame.

Model of solution: apply the usual principle to solve physical contradictions and allocate a property P to the system and a contradictory property -P to a subpart of the system, ie. primary property P is "the door opens to the right" and secondary property -P is "the door opens to the left".

Solution concept: Evolution law No. 7 suggests that the work is performed in a system, first on the object itself, then on a specific part of this object, when the system evolves [Fig. 5].

Fig. 4: how to open a door by always pushing it? Fig. 5: the "door in the door"! Last but not least, the winning team not only received an official certificate and a few articles in the regional newspapers, with the recognition of their peers and visibility at the ETRIA conference, but also the protection of their idea (in France, a "Soleau" envelope ensures the date of an idea) offered by the High School, so that the team could claim royalties in the future if an industry partner develop their idea.

Discussion

To evaluate the success of our approach, we noted several key points at all levels, as well as the usual concluding remarks by authors like Moshe Barak, in [START_REF] Barak | Teaching methods for inventive problem-solving in junior high school[END_REF], "The most important aspect of teaching pupils an inventive problem-solving course of the type discussed in this study is creating a climate of thinking and problem-solving in the class. Beyond the direct teaching of terms, thinking schemes and heuristics related to inventive thinking and problem-solving, it is important to give the pupils the time and opportunity to develop their own thinking methods and explain or justify their ideas."

The teachers managed to both train themselves and train students during a school year despite starting with strong contradictions at the beginning. A traditional approach (train teachers and students thereafter), would not have been able to train students (at least not during the current year). In addition to the basics of TRIZ tools, teachers also had practical experience of using these tools in different projects, and a beneficial teaching experience of these tools which remains the heart of their business. As an example of other positive output induced by the challenge of creativity for innovation, the mini company "ErgoTechJr" of the High School was awarded the prize for "excellence and innovation" in the regional championship of mini-companies.

The results of the projects submitted by students were rated at a good level. For example, to assign the 2nd and 3rd places in the challenge, the board of examiners faced six teams (from a total of 14) who all deserved to be on the podium. This particular situation was observed by Ogot and Okudan in [START_REF] Ogot | Systematic Creativity Methods in Engineering Education : A Learning Styles Perspective[END_REF], "From a learning styles perspective it appears that methods such as TRIZ that match the learning styles preferences of a broader proportion of the student population during their use, as compared to the sole use of traditional methods such as brainstorming, should result in an increase in the creative potential and output of the students."

The teachers noted also a stimulation of all stakeholders (experts, teachers, students): given the close links between each of them, the progress of some of these stakeholders was immediately beneficial to others, whereas in a sequential mode of training, the teachers do not take advantage of the progress of their students.

This positive dynamism also spread outside the walls of the school and made teachers reach an upper level by pushing them to formalize the process that they undertook, and to present it to ETRIA 2012.

Given the speed of technological change and increasingly rapid obsolescence of knowledge, this approach of continuous teacher training, with its advantages in terms of resource savings in time, is promises to be deployed widely.

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig. 1: Traditional education curves (1=students, 2=teachers, 3=expert)

Acknowledgements

This project was made possible due to J.J. Ostermeier, Regional pedagogical inspector, C. Buttner, High School Headmaster, D. Gigant, Technical education manager and A. De-Carvalho, O. Genty and J.M. Leurent, teachers, also Mrs C. Spindler and S. Ludmann, English teachers, and Mrs A. Martin-Boillat, Senior advisor for education, who enjoyed being part of the board of examiners, and the students who showed enthusiasm, methodology and inventivity in the "creativity challenge for innovation".