Julien Alexandre

Sandretto

Alexandre Chapoutot
email: alexandre.chapoutot@ensta-paristech.fr

Validated Explicit and Implicit Runge-Kutta Methods

A set of validated numerical integration methods based on explicit and implicit Runge-Kutta schemes is presented to solve, in a guaranteed way, initial value problems of ordinary differential equations. Runge-Kutta methods are wellknown to have strong stability properties which make them appealing to be the basis of validated numerical integration methods. A new approach to bound the local truncation error of any Runge-Kutta methods is the main contribution of this article which pushes back the current state of the art. More precisely, an efficient solution to the challenge of making validated Runge-Kutta methods is presented based on the theory of John Butcher. We also present a new interval contractor approach to solve implicit Runge-Kutta methods. A complete experimentation based on Vericomp benchmark is described.

Chapter 4 Validated Runge-Kutta methods

We present in this section our main contribution on the validation of Runge-Kutta methods. In Section 4.1, we present our approach to express the expression of the local truncation error of any (explicit and implicit) Runge-Kutta methods. In Section 4.2, our approach to solve implicit system of equations associated to implicit Runge-Kutta methods is presented. Finally, we show how we can use Runge-Kutta method to define a new interval enclosure of the Picard-Lindelöf operator in Section 4.3.

Chapter 1 Introduction

Many scientific applications such as in mechanics, in robotics, in chemistry or in electronics require the solution of differential equations. In the general case, differential equations can not be integrated formally, and a numerical integration scheme is used to approximate the state of the system. Nevertheless, in many applications, as for example [START_REF] Henzinger | Beyond HYTECH: Hybrid systems analysis using interval numerical methods[END_REF][START_REF] Bouissou | HybridFluctuat: A static analyzer of numerical programs within a continuous environment[END_REF][START_REF] Eggers | Improving SAT modulo ODE for hybrid systems analysis by combining different enclosure methods[END_REF][START_REF] Tucker | A rigorous ode solver and smale's 14th problem[END_REF], an approximation of the solution is not sufficient and a bound of the exact solution is mandatory. A new approach to compute such bounds is presented based on well-known Runge-Kutta methods.

In this article, we are interested in the computation of the solution of interval initial value problem (IIVP ) for autonomous Ordinary Differential Equations (ODEs) defined by ẏ = f (y) with y(0) ∈ [y 0 ] and t ∈ [0, t end ] .

(1.1)

The function f : 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.1) admits a unique solution [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF] for a given initial condition y 0 . Even more, for our purpose, we shall assume that f is continuously differentiable as needed. Note that the initial value is given by an interval, i.e., there is some bounded uncertainties on the initial value. More precisely, we are interested by methods computing the set of solutions y(t; [y 0 ]) of IIVP such that y(t; [y 0 ]) = y(t; y 0 ) : ∀y 0 ∈ [y 0 ] .

Remark 1.0.1 For simplicity only autonomous first order ODE are considered. It is not restrictive since any non-autonomous ODE can be rewritten into an autonomous ODE by increasing the dimension by one such that

ẏ = f (t, y) ⇔ ż = ṫ ẏ = 1 f (t, y) = g(z) .
Moreover, any high order ODE can be rewritten into a system of first order 2 ODEs. For example, a scalar second order problem can written such that

ÿ = f (y, ẏ) ⇔ ẏ1 ẏ2 = y 2 f (y 1 , y 2 )
with y 1 = y and y 2 = ẏ .

The guaranteed or validated solution of IIVP using interval arithmetic is mainly based on two kinds of methods based on: i) Taylor series [START_REF] Moore | Interval Analysis[END_REF][START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[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] ii) Runge-Kutta schemes [START_REF] Bouissou | GRKLib: a Guaranteed Runge Kutta Library[END_REF][START_REF] Gajda | A survey of interval Runge-Kutta and multistep methods for solving the initial value problem[END_REF][START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF]. The former is the oldest method used in interval analysis community because the expression of the bound of a Taylor series is simple to obtain. Nevertheless, the family of Runge-Kutta methods is very important in the field of numerical analysis. Indeed, Runge-Kutta methods have several interesting stability properties which make them suitable for an important class of problems. They are less used in interval analysis community because the expression to bound the approximation error is complex to synthesize.

We present new guaranteed numerical integration schemes based on Runge-Kutta methods. This work is based on [START_REF] Bouissou | GRKLib: a Guaranteed Runge Kutta Library[END_REF] which studied the classical Runge-Kutta method and its extension to any explicit Runge-Kutta methods [START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF]. The main contribution is the extension of these previous work on the definition of a set of guaranteed numerical integration schemes based on implicit Runge-Kutta formulas. Hence, having different guaranteed numerical integration schemes, explicit and implicit Runge-Kutta methods, we can handle more efficiently various kinds of problems.

Outlines In Section 2, we recall the classical algorithm of a validated simulation of an ODE, based on the 2-step Löhner type algorithm. In Section 3, we recall the basics on Runge-Kutta methods and their theory. We present in Section 4 our main contribution on the computation of the bounds of the local truncation error of any Runge-Kutta methods. In Section 5, an algorithm and associated proof for our new approach to compute implicit Runge-Kutta methods is presented. Section 6 contains the results of our method on several examples coming from the Vericomp benchmark. In Section 7, we summarize the main contributions of the paper.

Notations x denotes a real values while x represents a vector of real values.

[x] represents an interval value. An interval [x i ] = [x i , x i ] defines the set of reels x i such that x i ≤ x i ≤ x i . IR denotes the set of all intervals while R denotes the set of real values. The size or width of [

x i ] is w([x i ]) = x i -x i and m([x]) denotes the center of [x]. A vector of intervals, or a box, [x] is the Cartesian product of intervals [x 1 ] × ... × [x i ] × ... × [x n ].
Chapter 2

Technical preliminaries and related work

In Section 2.1, we recall the main steps of the validated method of numerical integration as it can be found in [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF]. We present related work in Section 2.2.

Validated numerical integration: a remainder

We recall the main algorithm used in the context of validated numerical integration and we refer to [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] for a more detailed presentation. The goal of a validated numerical algorithm to solve Equation (1.1) is to compute a sequence of time instants 0 = t 0 < t The main approach in a validated numerical integration method, as presented in [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF], is that each step of a validated integration scheme consists of two phases Phase 1 One computes an a priori enclosure [ỹ j ] of the solution such that • y(t; [y j ]) is guaranteed to exist for all t ∈ [t j , t j+1 ], i.e.n along the current step, and for all y j ∈ [y j ];

• y(t; [y j ]) ⊆ [ỹ j ] for all t ∈ [t j , t j+1 ];
• the step-size h j = t j+1 -t j > 0 is as large as possible in terms of accuracy and existence proof for the IIVP solution.

Phase 2 One computes a tighter enclosure of [y j+1 ] at time t j+1 such that

y(t j+1 , [y j ]) ⊆ [y j+1 ].
The different enclosures computed during one integration step between time t j and t j+1 are shown on Figure 2.1.

4 t y(t) [ỹ j ] [y j ] [y j+1 ] t j t j+1 h j Figure 2
.1: Enclosures appearing during one step Some simple algorithms to perform these two steps are described in the following. We refer to [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] for a description of more advanced algorithms. The main issue in these two phases is to counteract the well known wrapping effect [START_REF] Lohner | On the ubiquity of the wrapping effect in the computation of error bounds[END_REF]. This phenomenon appears when one tries to enclose a set within a box. The reader is refer to [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] to have a clear presentation of the methods to reduce pessimism such as the QR-decomposition. In Section 5, an other approach to reduce pessimism based on affine arithmetic [START_REF] De Figueiredo | Self-Validated Numerical Methods and Applications[END_REF] is presented.

A priori solution enclosure

Phase 1 computes an a priori enclosure of the solution of IIVP over the whole time interval [t j , t j+1 ] based on the application of the Banach fixed point theorem (see Theorem 2.1.1) with the Picard-Lindelöf operator, see Equation (2.1).

Theorem 2.1.1 (Banach fixed-point theorem) Let (K, d) a complete metric space and let g : K → K a contraction that is for all x, y in K there exists c ∈ ]0, 1[ such that d (g(x), g(y)) c • d(x, y) , then g has a unique fixed-point in K.

We consider the space of continuously differentiable functions C 0 ([t j , t j+1 ], R n ) and the Picard-Lindelöf operator

p f (y) = t → y j + t tj f (y(s))ds , (2.1)
with y j the condition at time t j used to solve Equation (1.1). Note that this operator is associated to the integral form of Equation (1.1). As a consequence, if this operator is a contraction then its solution is unique and its solution is the solution of Equation (1.1). One can define an interval counter part of the Picard-Lindelöf operator which can be used to operationally prove the contraction and so the existence and uniqueness of the solution of Equation (1.1). With a first order integration scheme [START_REF] Moore | Interval Analysis[END_REF], that is for f : R n → R n a continuous function and [a] ⊂ IR n , we have

a a f (s)ds ∈ (a -a)f ([a]) = w([a])f ([a]) , (2.2)
we can define a simple enclosure function of Equation (2.1) such that

[p f ]([r]) def = [y j ] + [0, h] • f ([r]) , (2.3)
with h = t j+1 -t j the step-size. In consequence, if one can find [r] such that

[p f ]([r]) ⊆ [r] then [ỹ j ] ⊆ [r]
by the Banach fixed-point theorem.

Once the contraction of [p f ] has been proven, we can also define an interval contractor on [r], in the sense of [START_REF] Chabert | Contractor programming[END_REF], to refine the value of

[r] such that [r] ← [r] ∩ [p f ]([r]) . (2.4)
The operator defined in Equation (2.3) and its associated contractor defined in Equation (2.4) can be improved with more accurate interval enclosure functions for the integral operator. For example, the evaluation of t tj f (s)ds can be improved with any validated integration scheme, such as Taylor polynomial, see [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] for more details.

Tighter enclosure

Once Phase 1 is completed, one has the a priori enclosure [ỹ j ] such that

y(t; [y j ]) ⊂ [ỹ j ] ∀t ∈ [t j , t j+1 ] .
In particular, we have y(t j+1 ; [y j ]) ⊂ [ỹ j ]. The goal of Phase 2 is thus to compute a tighter enclosure of [y j+1 ] at time t j+1 such that

y(t j+1 ; [y j ]) ⊂ [y j+1 ] ⊆ [ỹ j ] .
Each validated numerical integration method is decomposed in two parts

• the approximation part Φ(t, [y j ]) ≈ y(t; [y j ]). This algorithm usually follows a numerical integration method such as Runge-Kutta methods.

•

[y j+1 ] ⊆ [y j ] + h[f ] ([y j ]) + h 2 2 df dt ([ỹ j ]) .
The approximation part is given by

[y j ] + h[f ] ([y j ]) while the local truncation error LTE Euler is given by h 2 2 [ df dt ] ([ỹ j ]).

Related work

Validated numerical integration methods have been intensively developed since the work of R. Moore [START_REF] Moore | Interval Analysis[END_REF] using Taylor series. Indeed, Taylor series became very popular because a simple expression of the local truncation error exists, and because the development of automatic differentiation techniques has offered efficient algorithms to compute high-order derivatives. Several tools based on Taylor series have been developed among AWA [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF], ADIODES [START_REF] Stauning | Automatic Validation of Numerical Solutions[END_REF], Vnode-LP [START_REF] Nedialko | Implementing a rigorous ode solver through literate programming[END_REF], COSY Infinity [START_REF] Makino | COSY INFINITY version 9. Nuclear Instruments and Methods in[END_REF], VSpode [START_REF] Lin | Validated solutions of initial value problems for parametric odes[END_REF], CAPD [START_REF]CAPD -Computer Assisted Proofs in Dynamics group, a C++ package for rigorous numerics[END_REF], Flow * [START_REF] Chen | Taylor model flowpipe construction for non-linear hybrid systems[END_REF]. We propose an other look at validated numerical integration method by using Runge-Kutta methods. In particular, in Taylor series approach implicit schemes [START_REF] Nedialko | An interval hermiteobreschkoff method for computing rigorous bounds on the solution of an initial value problem for an ordinary differential equation[END_REF] are only used in Phase 2 (see Section 2.1) while our contribution shows that it is possible to use implicit Runge-Kutta scheme even for Phase 1. Moreover, except for [START_REF] Lin | Validated solutions of initial value problems for parametric odes[END_REF] and [START_REF] Makino | COSY INFINITY version 9. Nuclear Instruments and Methods in[END_REF] which use Taylor models, Taylor series approaches provide validated result for small uncertainties on the initial values. Taylor models can increase the size of the interval of initial values but it remains costly in term of computation. We provide also an approach based on affine arithmetic [START_REF] De Figueiredo | Self-Validated Numerical Methods and Applications[END_REF], in the same spirit than [START_REF] Kühn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF], in order to increase the width of initial values while keeping low complexity on arithmetic operations. The work of Andrzej Marciniak et al. presented in [START_REF] Marciniak | One-and two-stage implicit interval methods of Runge-Kutta type[END_REF][START_REF] Gajda | Three-and fourstage 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] Gajda | A survey of interval Runge-Kutta and multistep methods for solving the initial value problem[END_REF] is the closest of ours. Indeed, they intensively studied a subclass of implicit Runge-Kutta methods [START_REF] Marciniak | One-and two-stage implicit interval methods of Runge-Kutta type[END_REF][START_REF] Gajda | Three-and fourstage 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] for which they provide insights on how make them guaranteed. The main differences are i) they express "by hand" the local truncation error (i.e., the distance between the exact solution and the numerical one, see Section 3) of a set of particular implicit Runge-Kutta methods. Indeed they claim in [START_REF] Marciniak | On representations of coefficients in implicit interval methods of runge-kutta type[END_REF] that the local truncation error expression ". . . is very complicated and cannot be written in general form for an arbitrary order p". In this article, we provide an algorithm to compute this local truncation error for any Runge-Kutta method, and so we push back the current state of the art, see Section 4.

ii) they only consider fixed step integration methods that is the step-size h is fixed during the simulation. In this article, we present validated Runge-Kutta methods in the standard framework of validated numerical integration as presented in [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] and so our approach benefit of variable step-size techniques, see Section 5.3.

iii) implicit methods require the solution of non-linear system of equations. Marciniak et al. solve this problem using a simple iterative scheme as it can be found in [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF]. In contrary, we use interval contractor approach [START_REF] Chabert | Contractor programming[END_REF] to solve the non-linear system of equations producing a more robust algorithm, see Section 4.2.

Other approaches to define validated numerical integration have been considered. In particular, Valencia-IVP [START_REF] Rauh | A novel interval arithmetic approach for solving differential-algebraic equations with ValEncIA-IVP[END_REF] is based on a defect estimate approach which does not necessitate high order derivative but only f and a approximate trajectory computing by standard numerical algorithms. Moreover, in [START_REF] Fazal | Error bounds for initial value problems by optimization[END_REF], the defect estimate approach is also used but combined with a global optimization approach to bound the solution of Equation (1.1).

Chapter 3

Recall on Numerical Runge-Kutta methods and their theory

When the initial value of IIVP is exactly known that is y(0) = y 0 , an initial value problem (IVP ) is considered. In that case, there are several numerical methods to solve IVPs [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]. Among them, Runge-Kutta methods are very well studied and often used. A Runge-Kutta method, starting from y 0 at time t 0 and a finite time horizon h, produces an approximation y 1 at time t 0 + h of the solution y(t 0 + h; y 0 ). Furthermore, to compute y 1 , a Runge-Kutta method computes s intermediate steps where s is known as the number of stages. More precisely, a Runge-Kutta method, for a non-autonomous system, i.e., ẏ = f (t, y), is defined by

y 1 = y 0 + h s i=1 b i k i , (3.1)
with k i defined by

k i = f t 0 + c i h, y 0 + h s i=1 a i,j k j . (3.2)
In case of autonomous systems as the case considered in Equation (1.1), Equation (3.2) is rewritten as The form of the Butcher tableau, in particular, the form of the matrix A made of the coefficients a ij , determines the Runge-Kutta method, it can be

k i = f y 0 + h s i=1 a i,j k j . ( 3
• explicit, the matrix A is strictly lower triangular so each value k i is only defined from the previous values k j for j < i, e.g., the classical Runge-Kutta method given in Figure 3.1(a);

• diagonally implicit, the matrix A is lower triangular, e.g., the singly diagonally implicit method given in Figure 3.1(b);

• fully implicit, the matrix A is full, e.g., the Runge-Kutta method with a Lobatto quadrature formula given in Figure 3.1(c). Theses different kinds of Runge-Kutta methods are associated to different stability properties. In particular, only implicit Runge-Kutta methods can be A-stable that is they are unconditionally stable1 for any stable linear dynamical systems of the form ẏ = Ay with ρ(A) < 1, where ρ(.) denotes the spectral radius of the matrix A. Nevertheless, implicit methods are more costly than explicit ones. Indeed, for implicit methods, at each integration step a nonlinear system of equations has to be solved to compute the values k i for all i = 1, . . . , s.

A Runge-Kutta has order p if one has

y(t 0 + h; y 0 ) -y 1 ≤ O(h p+1 ) .
We recall in the next section the construction of Runge-Kutta methods of a given order p. The theory of Runge-Kutta methods plays an important role to build a validated version of these methods.

Theory of Runge-Kutta methods: a brief overview

The modern theory of Runge-Kutta methods has been defined by John Butcher with his work presented in [START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF]. Informally, the order of a Runge-Kutta method is defined as the highest order of the first non-zero term of a Taylor series build from the difference between the Taylor series of the exact solution and the Taylor series of the numerical solution. This is known as the order conditions of Runge-Kutta methods, see [20, Chap. III]. One of the major contributions of the work of J. Butcher is to express these two Taylor series on a common basis made of elementary differentials, that are partial derivatives of f given in Equation (1.1). We briefly recall the order condition of Runge-Kutta methods as it plays an important role in our contribution. The presentation and notations follow [20, Chap. III]. In the sequel, we consider IVP problem defined in Equation (1.1) with an exact initial condition y(0) = y 0 .

High order derivative of exact solution It is interested in computing the higher derivatives of the exact solution y(t) at t = 0. In particular, the q-th time derivative of y is defined by y (q) = (f (y)) (q-1) . Using the chain rule and some symmetry of partial derivatives, we get the following first four derivatives

ẏ = f (y) ÿ = f (y) ẏ y (3) = f (y)( ẏ, ẏ) + f (y)ÿ y (4) = f (4) (y)( ẏ, ẏ, ẏ, ẏ) + 3f (y)(ÿ, ẏ) + f (y)y (3) (3.5)
f stands for the first order partial derivatives of f w.r.t. y. In the same, way f stands for the second order partial derivatives of f w.r.t. y and so on. By recursively inserting in the right hand side of Equation (3.5) the definition of ẏ, ÿ, . . . , and removing the argument y, we get

ẏ = f ÿ = f f y (3) = f (f , f ) + f f f y (4) = f (4) (f , f , f , f ) + 3f (f f , f ) + f f f f (3.6)
The expressions appearing in these sums, denoted in the sequel F (τ ), are named elementary differentials. Moreover, it is possible to represent these terms by a rooted tree τ as follows

• each f is a leaf of τ , • each f (k) , k 1, is associated to a node in τ with k branches.
The number of nodes in a rooted tree τ is denoted by |τ |.

F (τ )(y) = f (m) (y) (F (τ 1 )(y), F (τ 2 )(y), . . . , F (τ m )(y)) if τ made of sub-trees τ 1 , τ 2 , . . . , τ m .
The link between rooted trees and elementary differentials is given in Table 3.1. In consequence, Theorem 3.1.1 can be stated to express high order derivative of the exact solution in terms of elementary differentials. Theorem 3.1.1 The q-th derivative of the exact solution at t = 0 is given by

y (q) (0) = |τ |=q α(τ )F (τ )(y 0 ) , α (
τ ) are positive integer values has a combinatorial meaning that is they represents the possible symmetries in rooted trees τ . In particular, a tree τ does not depend on the order of sub-trees τ 1 , τ 2 , . . . , τ m . High order derivative of numerical solution Let hk i = g i hence a Runge-Kutta method can be written as

g i = hf (u i ) , (3.7)
and

u i = y 0 + j a ij g j , y 1 = y 0 + i b i g i . |τ | Trees F (τ ) α(τ ) γ(τ ) ϕ(τ ) 1 f 1 1 i b i 2 f f 1 2 ij b i a ij 3 f (f , f ) 1 3 ijk b i a ij a ik 3 f f f 1 6 ijk b i a ij a jk 4 f (f , f , f ) 1 4 ijkl b i a ij a ik a il 4 f (f f , f ) 3 8 ijkl b i a ij a ik a jl 4 f f (f , f ) 1 12 ijkl b i a ij a jk a jl 4 f f f f 1 24 ijkl b i a ij a jk a kl
Table 3.1: Rooted trees and their associated elementary differentials with their coefficients Note that, u i , g i and y 1 are functions of h. We compute the derivative of Equation (3.7) at h = 0 using g 1) . In consequence, the first three time derivatives are

(q) i = q • (f (u i )) (q-
ġi = 1 • f (y 0 ) gi = 2 • f (y 0 ) ui g (3) i = 3 • (f (y 0 )( ui , ui + f (y 0 )ü i ) (3.8)
where the derivative of g i and u i are evaluated at h = 0. We can remark a similar form than Equation (3.6). Inserting recursively the definition of ui , üi , . . . , in Equation (3.8) and using u

(q) i = j a ij g (q) i , one has ġi = 1 • f ui = 1 •   j a ij   • f (3.9) gi = (1 • 2) •   j a ij   f f üi = (1 • 2) •   jk a ij a jk   • f f (3.10)
and so on. The integer factors 1, (1 • 2), . . . , are denoted by γ(τ ). The factors containing the a ij 's are denoted by g i (τ ) and u i (τ ), then one has

g (q) i | h=0 = |τ |=q γ(τ ) • g i (τ ) • α(τ )F (τ )(y 0 ) (3.11) u (q) i | h=0 = |τ |=q γ(τ ) • u i (τ ) • α(τ )F (τ )(y 0 ) (3.12)
where α(τ ) and F (τ ) have the same meaning as before. Avoiding some more writing rules, see [START_REF] Hairer | Geometric numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]Chap. III] for the details, we get Theorem 3.1.2 which gives the high order derivative of the numerical solution in terms of elementary differentials.

Theorem 3.1.2 The q-th derivative of the numerical solution of a Runge-Kutta method is given by

y (q) 1 | h=0 = |τ |=q γ(τ ) • ϕ(τ ) • α(τ )F (τ )(y 0 ) , (3.13)
with ϕ(τ ) = i b i g i (τ ).

In Figure 3.1 some values of the coefficients γ(τ ) and ϕ(τ ) are given. Note that Theorem 3.1.2 can be applied on explicit and implicit Runge-Kutta methods, once the Butcher tableau of the method is known.

Order condition of Runge-Kutta methods Theorem 3.1.3, associated to the order condition of Runge-Kutta method, can be stated in terms of the coefficients γ(τ ) and ϕ(τ ).

Theorem 3.1.3 A Runge-Kutta method has order p if and only if

ϕ(τ ) = 1 γ(τ ) ∀|τ | p . (3.14)
The main difficulty to build a Runge-Kutta method is that solving Equation (3.14) require to solve a high dimensional under-determined system of polynomial equations.

Bounding the local truncation error

In our purpose to make Runge-Kutta methods validated, we apply Theorem 3.1.3 assuming that the considering Runge-Kutta method has order p. In that case, Theorem 3.1.3 offers an unified approach to express the local truncation error of any Runge-Kutta method.

From Theorem 3.1.1, the Taylor series up to order p + 1, with Lagrange remainder, of the exact solution around t j with y(t j ) = y j is given by

y(t j + h) = p i=0 h i i! |τ |=i α(τ ) • F (τ )(y j ) + h p+1 (p + 1)! |τ |=p+1 α(τ ) • F (τ )(y(ξ)) with ξ ∈]t j , t j+1 [ . (4.1)
Moreover, from Theorem 3.1.2, the Taylor series up to order p + 1, with Lagrange remainder, of the numerical solution around t j with y(t j ) = y j and 15 h = t j+1 -t j is given by

y j+1 = p i=0 h i i! |τ |=i γ(τ ) • ϕ(τ ) • α(τ ) • F (τ )(y j ) + h p+1 (p + 1)! |τ |=p+1 γ(τ ) • ϕ(τ ) • α(τ ) • F (τ )(y(ξ)) with ξ ∈]t j , t j+1 [ . (4.2)
From Equation (2.5), with Equation (4.1) and Equation (4.2), we get

LTE RK (ξ, y) def = h p+1 (p + 1)! |τ |=p+1 α(τ ) [1 -γ(τ )ϕ(τ )] F (τ )(y(ξ)) with ξ ∈]t j , t j+1 [ . (4.3)
Indeed, assuming that the considered Runge-Kutta has order p, we know that for all |τ | p one has ϕ(τ ) γ(τ ) = 1, so all the first p + 1 terms of the Taylor series defined as the difference between Equation (4.1) and Equation (4.2) are zero. We present in Section 5.2, how Equation (4.3) can be computed.

∈]t j , t j+1 [, y(ξ) ∈ [ỹ j ] ⇒ LTE RK (ξ, y) ∈ LTE RK ([t j , t j+1 ], [ỹ j ]) .
As a consequence of Proposition 4.1.1, we are able to validate any explicit and implicit Runge-Kutta methods.

Solving implicit Runge-Kutta methods

Using an implicit Runge-Kutta in an integration scheme needs to solve a system of non-linear equations (see Section 3). In classical numerical methods, it is done with a Newton-like solving procedure which provides generally a good approximation of the k i . While some interval Newton-like procedure exists for solving systems of non-linear interval equations [START_REF] Moore | Interval Analysis[END_REF], we propose a lighter approach described in the following.

First of all, it is interesting to note that each stages of an implicit Runge-Kutta allowing us to compute the intermediate k i can be used as an interval contractor [START_REF] Chabert | Contractor programming[END_REF], see Proposition 4.2.1. Proposition 4.2.1 Each stage of an implicit Runge-Kutta is a natural contractor for k i , i = 1, . . . , s.

Proof: We recall the form of an intermediate stage of an implicit Runge-Kutta method

k i = f (y j + h s n=1 a in k n ) . (4.4)
We also know that for all the Runge-Kutta methods [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]

c i = s n=1 a in ≤ 1, ∀i = 1, . . . , s .
Moreover, by the Picard-Lindelöf operator, we have k i ∈ [ỹ j ], i = 1, . . . , s, because t j + c i h ≤ t j + h. Inserting this inside Equation (4.4) leads to

s n=1 a in k n ⊂ s n=1 a in [ỹ j ] = c i [ỹ j ] .
Then, we can write

y j + h s n=1 a in k n ⊂ y j + h[ỹ j ] .
By Theorem 2.1.1 and property of [ỹ j ] obtained by Picard-Lindelöf operator, f is contracting on y j + h[ỹ j ], and also on y j + h s n=1 a in k n . By using the previous proposition, we write the following contractor scheme:

k i = k i ∩ f   y n + h s j=1 a i,j k j   .
This contractor is used inside a fixpoint presented in Algorithm 1 to form the solver for the implicit Runge-Kutta methods.

Algorithm 1 Solving an implicit Runge-Kutta

Require: [ỹ j ], a in of an implicit RK k i = f ([ỹ j ]), ∀i = 1, . . . , s while at least one k i is contracted do k 1 = k 1 ∩ f (y j + h s n=1 a 1n k n ) . . . k s = k s ∩ f (y j + h s n=1 a sn k n ) end while

A priori enclosure method with Runge-Kutta formula

In this section, we present how Runge-Kutta methods combined with their expression of the local truncation error can be used to define a new a priori enclosure method, i.e., a new algorithm for Phase 1 (see Section 2.1). The challenge to search for new algorithms for Phase 1 is to increase the size of the class of problems that validated numerical methods can be applied on. In particular, solving stiff differential problems remains a challenge. For stiff problem, the integration step-size h is highly related to the stability conditions of numerical methods. Until now, only explicit algorithms, i.e., that only depends on [y j ] to compute [y j+1 ], are used in Phase 1. A novelty of our approach is that we can define new validated numerical integration methods based on implicit Runge-Kutta methods, which are known to have very good stability properties.

Despite that more work have to be done to prove the relevance of our approach on stiff problems, we believe that this opens the way to new researches on the stability of validated numerical methods, in the same spirit than [START_REF] Nedialko | A new perspective on the wrapping effect in interval methods for initial value problems for ordinary differential equations[END_REF].

To define a new a priori enclosure, we consider Runge-Kutta methods for which we clearly introduce the time dependence, i.e.,

k i (t, y j ) = f y j + (t -t j ) s n=1 a in k n , y j+1 (t, y) = y j + (t -t j ) s i=1 b i k i (t, y j ) + LTE RK (ξ, y) with ξ ∈]t j , t[ .
With an integration step-size h = t j+1 -t j , we can define an inclusion function such that

y j+1 ([t j, t j+1 ], [r]) def = [y j ] + [0, h] s i=1 b i k i ([t j , t j+1 ], [y j ]) + LTE RK ([t j , t j+1 ], [r]) . (4.5)
Proving the contraction of such scheme, that is Chapter 5

[r] ⊇ y j+1 ([t j , t j+1 ], [r]) , ( 4

Implementation details

A presentation of the main features of the implementation of the validated numerical integration based Runge-Kutta methods is given.

Affine arithmetic

Example 5.1.1 illustrates the pessimism in numerical integration method introduced by the dependency problem inherited in interval arithmetic. Usually, to fight this problem sharper enclosure functions are used as the centered form.

Example 5.1.1 Consider the ordinary differential equation ẋ(t) = -x solved with the Euler's method with an initial value ranging in the interval [0, 1] and with a step-size of h = 0.5. For one step of integration, we have to compute with interval arithmetic the expression e = x + h × (-x) which produces as a result the interval [-0.5, 1]. Rewriting the expression e such that e = x(1 -h), we obtain the interval [0, 0.5] which is the exact result. Unfortunately, we cannot in general rewrite expressions with only one occurrence of each variable. More generally, it can be shown that for most integration schemes the width of the result can only grow if we interpret sets of values as intervals [START_REF] Neumaier | The wrapping effect, ellipsoid arithmetic, stability and confidence regions[END_REF].

In our work, to avoid or limit pessimism due to the dependency problem, we use an improvement over interval arithmetic named affine arithmetic [START_REF] De Figueiredo | Self-Validated Numerical Methods and Applications[END_REF][START_REF] Rump | Implementation and improvements of affine arithmetic[END_REF] which can track linear correlations between program variables. A set of values in this domain is represented by an affine form x (also called a zonotope), which is a formal expression of the form x = α 0 + n i=1 α i ε i where the coefficients α i are real numbers, α 0 being called the center of the affine form, and the ε i are formal variables ranging over the interval [-1, 1]. Obviously, an interval a = [a 1 , a 2 ] can be seen as the affine form x = α 0 + α 1 ε with α 0 = (a 1 + a 2 )/2 and α 1 = (a 2 -a 1 )/2. Moreover, affine forms encode linear dependencies between variables: if x ∈ [a 1 , a 2 ] and y is such that y = 2x, then x will be represented by the affine form x above and y will be represented as ŷ = 2α 0 + 2α 1 ε.

19

Usual operations on real numbers extend to affine arithmetic in the expected way. For instance, if x = α 0 + n i=1 α i ε i and ŷ = β 0 + n i=1 β i ε i , then with a, b, c ∈ R we have

ax + bŷ + c = (aα 0 + bβ 0 + c) + n i=1 (aα i + bβ i )ε i .
However, unlike the addition, most operations create new noise symbols. Multiplication for example is defined by

x × ŷ = α 0 α 1 + n i=1 (α i β 0 + α 0 β i )ε i + νε n+1 , where ν = ( n i=1 |α i |) × ( n i=1 |β i |)
over-approximates the error between the linear approximation of multiplication and multiplication itself. Example 5.1.2 illustrates the benefit of affine arithmetic.

Example 5.1.2 Consider again e = x + h × (-x) with h = 0.5 and x = [0, 1] which is associated to the affine form x = 0.5 + 0.5ε 1 . Evaluating e with affine arithmetic without rewriting the expression, we obtain [0, 0.5] as a result.

Example 5.1.2 also shows the important role of affine arithmetic when it is combined with numerical integration methods. Most of all, it shows the necessity to keep track the linear dependency between state variables in order to reduce the pessimism.

Other operations, like sin, exp, are evaluated using either the Min-Range method or a Chebychev approximation see [START_REF] De Figueiredo | Self-Validated Numerical Methods and Applications[END_REF][START_REF] Rump | Implementation and improvements of affine arithmetic[END_REF] for more details.

The main challenge by using the affine arithmetic during an integration process is to limit the number of noise symbols. No good solution exist and we have an heuristic to periodically collect all the noise symbols which have a coefficient lower than a given user threshold.

Wrapping effect and affine arithmetic The problem of reducing the wrapping effect has been studied in many different ways. One of the most known and effective is the QR-factorization [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF]. This method improves the stability of the Taylor series in the Vnode-LP tool [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] of the CAPD tool. Nevertheless, affine arithmetic allows to counteract the wrapping effect as shown in Figure 5.1 while keeping a fast computation. Indeed, the geometric interpretation of an affine form x = α 0 + n i=1 α i ε i is a zonotope that is a convex polytope with central symmetry. As a consequence, zonotopes are invariant by rotation so they are well suited for representing sets on which rotation operation may be applied as in numerical integration methods.

with initial values y 1 (0) ∈ [-1, 1], y 2 (0) ∈ [START_REF]CAPD -Computer Assisted Proofs in Dynamics group, a C++ package for rigorous numerics[END_REF][START_REF] Cash | A variable order Runge-Kutta method for ivp with rapidly varying right-hand sides[END_REF]. The exact solution of Equation (5.1) is

y(t) = A(t)y 0 with A(t) = cos(t) sin(t) -sin(t) cos(t)
We compute periodically at t = π 4 n with n = 1, . . . , 4 the solution of Equation (5.1) with different enclosure methods: standard interval, interval with QR-decomposition and affine arithmetic. The evolution of the enclosures are given in Figure 5.1. Implicit Runge-Kutta method and affine arithmetic Algorithm 1, presented in Section 4.2, is light and, according to our tests, as efficient than a Newton-like method. Nevertheless, one major problem occurs. The intersection operator is not available on affine arithmetic because the intersection of two zonotopes is not a zonotope. We then perform the Algorithm 1 with interval arithmetic, and after reaching the fixpoint, we evaluate Equation (4.4) with affine arithmetic. Hence, we can keep tracking linear correlation between state variables and we use the natural contractivity of Equation (4.4) to keep a sharp enclosure of the solution of the non-linear systems of equations. This method is sufficient to counteract the wrapping effect and the dependency problem appearing during the simulation process.

Computing the local truncation error

In Section 4.1, a new expression of the local truncation error for any Runge-Kutta methods has been presented but it remains to show how this formula can be computed. We recall the expression of the local truncation error for any Runge-Kutta methods

LTE RK (ξ, y) def = h p+1 (p + 1)! |τ |=p+1 α(τ ) [1 -γ(τ )ϕ(τ )] F (τ )(y(ξ)) with ξ ∈]t j , t j+1 [ . (5.2)
We remark that many elements in Equations (5.2) can be computed offline. Indeed, for a given Runge-Kutta of order p, we can automatically generate the set of all rooted trees τ such that |τ | = p + 1 following algorithms presented in [START_REF] Bornemann | Runge-kutta methods, trees, and maple -on a simple proof of butcher's theorem and the automatic generation of order condition[END_REF]. Moreover, once the set of rooted trees is given all the coefficients α(τ ), γ(τ ), ϕ(τ ) (the Butcher tableau of the method is also needed for this coefficient) can be also computed off-line, see [START_REF] Lambert | Numerical Methods for Ordinary Differential Systems: The Initial Value Problem[END_REF] for a formal definition of this values based on the structure of a rooted tree τ . In consequence, only F (τ ) is problem dependent and require an online computation.

The elementary differential F (τ ) are associated to Fréchet derivatives. Following the notations of [START_REF] Lambert | Numerical Methods for Ordinary Differential Systems: The Initial Value Problem[END_REF], let z, f (z) ∈ R m , the M -th 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 (5.3) where i f j1j2...j M = ∂ M ∂ j1 z∂ j2 z . . . ∂ j M z K k = [ 1 K 1 , 2 K 2 , . . . , M K M ] ∈ R m , for k = 1, . . . , M .
The notation x stands for the -th component of the vector x, and e i stands for the vector full of zero except for i-th component which is set to one. In the context of the local truncation error, the vectors K k for k = 1, . . . , M are the elementary derivatives associated to the rooted tree τ = [τ 1 , . . . , τ x ] such that |τ | = M .

The advantage of the Fréchet derivative given in Equation (5.3) is that high order derivative of f does not necessitate complex tensor operations because they are computed component by component for the state vector z. From an algorithmic point of view, computing a Fréchet derivative only require the computation of partial derivative of f up to order p + 1 if the Runge-Kutta has order p. Nevertheless, the drawback is that the M -th Fréchet derivative involve an algorithm with M + 1 nested loops such that the complexity of the algorithm is exponential in the dimension n of the function f , i.e., O(n M +1 ). In consequence, only small dimension problems can be efficiently solved.

Algorithm of validated Runge-Kutta methods

In this section, we present the validated numerical integration based on Runge-Kutta method as it is implemented in our tool DynIbex [START_REF] Alexandre | Dynibex library[END_REF]. Algorithm 2 gather all the steps needed for the simulation of IIVP with Runge-Kutta schemes, explicit or implicit. This algorithm uses some functions described below

• RKe: a non guaranteed explicit Runge-Kutta method (RK4 for example)

• RKx: a guaranteed explicit, by an affine evaluation, or implicit, with Algorithm 1, Runge-Kutta method (RK4 or LC3 for examples)

• LTE: the local truncation error associated to RKx (see Section 4.1)

• PL: the Picard-Lindelöf operator based on an integration scheme (rectangular, Taylor or Runge-Kutta, see Section 2.1.1)

The main simulation loop is given between Line 2 and Line 34. The first part of the simulation loop, between Line 6 and Line 19, is the proof of existence and uniqueness of the solution of the IIVP using the interval Picard-Lindelöf operator. First an estimation of the a priori enclosure is computed from the initial condition and the result of a numerical integration (Line 6) using RKe function. Between Line 8 and Line 14, the algorithm computes the post fixedpoint of the interval Picard-Lindelöf operator. Between Line 15 and Line 19, the contractor version of the interval Picard-Lindelöf operator is used to refine the bounds of the a priori enclosure.

Once the a priori enclosure is computed the bounds of the local truncation error is computed at Line 20. After this operation is done, between Line 21 and Line 27 the accuracy of the bounds of the LTE is checked (Line 22). If the accuracy is met then a new step-size h is computed (Line 23). We based our method to vary the step-size on classical formula as found in [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]Chap. II.4]. Otherwise, the step-size is divided by two and we restart the whole integration step. Note that if the interval Picard-Lindelöf operator is unable to compute a post fixed-point (we bound the number of iterations in function of the dimension of the function f , see Line 10) then we also divide the step-size by two an we restart the loop. In case of success, we compute the next enclosure of the solution of IIVP (Line 32), we advance the simulation time and we pursue the simulation.

Algorithm 2 Simulation algorithm Require: RK, f , y 0 , t end , h, atol, rtol 1: t n = t 0 , y n = y 0 , factor = 1 2: while (t n < t end ) do

3: h = h × factor 4: h = min(h, t end -t n ) 5:
Loop:

6:

Initialize ỹ0 = y n ∪ RKe(y n , h) Compute ỹ1 with P L(ỹ 0 )

13:

iter = iter + 1

Experimental results

We used the VERICOMP [START_REF] Auer | Vericomp: a system to compare and assess verified ivp solvers[END_REF] to put to test our approach. All the detailed results can be found in [START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF]. In this paper, we also show the result of our validated Runge-Kutta methods on two examples in order to demonstrate the efficiency of affine arithmetic with an harmonic oscillator, and the ability of Runge-Kutta methods for stiff problem with a chemical reaction problem. We also summarize the results obtained with the VERICOMP benchmark.

Our tool DynIbex [START_REF] Alexandre | Dynibex library[END_REF] is a plug-in of the Ibex library1 which is a C++ library for constraint processing over real numbers. We implemented various validated Runge-Kutta methods2 among Heun method, Midpoint method, Radau-IIA (order 3), classic Runge-Kutta of order 4, Lobatto-IIIA (order 4), Lobatto-IIIC (order 4). In these tests, the Picard-Lindelöf operator is a Taylor one at order 3. It is a good compromise and allows us to compare only the integration scheme and not the a priori enclosure.

Two simple problems 6.1.1 Harmonic oscillator

We start by the simulation of the differential system given by ẏ1 ẏ2 = -y 2 y 1 .

with the interval initial states: y 0 ∈ ([0, 0.1]; [0.95, 1.05]). The result in term of [y j ] is plotted in Figure 6.1. The stability of the boxes size proves that the wrapping effect is well neutralized.

Chemical reaction

The second problem, coming from [START_REF] Cash | A variable order Runge-Kutta method for ivp with rapidly varying right-hand sides[END_REF], is described by the differential system

Vericomp benchmark

This section reports the results of the solution of various problems coming from the VERICOMP3 benchmark [START_REF] Auer | Vericomp: a system to compare and assess verified ivp solvers[END_REF]. For each problem, different validated methods of Runge-Kutta of order 4 are applied among: the classical formula of Runge-Kutta 4 (explicit), the Lobatto-3a formula (implicit) and the Lobatto-3c formula (implicit). Moreover, an homemade version of Taylor series, limited to order 4 and using affine arithmetic, is also applied on each problem.

For each problem, we report the following metrics:

• c5t: user time taken to simulate the problem for 1 second.

• c5w: the final diameter of the solution (infinity norm is used).

• c6t: the time to breakdown the method with a maximal limit of 10 seconds.

• c6w: the diameter of the solution a the breakdown time.

The complete table gathering these results are available in [START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF]. In this paper, we sort the competitors with two filters. Firstly, we count the number of problems for which each method (for each order and each precision) is first in term of solution diameter, second or last. This account is done for the simulation at 1 second and at 10 seconds. The results are summarized in Table 6.1. f course, we are aware that the results are biased by the number of methods we have. Nevertheless, this table allows us to consider that Valencia and Riot are not valid competitors. After this first reduction of competitors, only the best results for our three order-4 Runge-Kutta methods, and for Vnode are kept for comparison. We present in the spider graph on Figure 6.3, respectively on Figure 6.4, the normalized diameters (divided by the median and multiplied by 10) for each problem for a simulation at 1 second, respectively at 10 seconds. The median used to normalize the results is computed with all the methods: Taylor4, RK4, LA3, LC3 Riot, Valencia and Vnode (for all precision and all order). We can easily see on spider Graph 6.3 that the Runge-Kutta methods are more stable, by describing a circle while Vnode results are more in a star shape. Moreover, the implicit methods (LA3 and LC3) provide better results than the explicit RK4 in a majority of problems. This fact is even more clear on Graph 6.4. On this latter graph, we can also see that Vnode fails many times while at least one of our Runge-Kutta methods performs a good simulation for all the considered problems. Finally, if Vnode are the best on many problems, by our stability and our better results for some problems, we can conclude that our tool is a good competitor for Vnode. The last remark but not the least, it is important to remember that we have currently only methods of order 4, when Vnode can use a Taylor expansion at order 25.

Conclusion

We present in this paper a novel approach to create validated numerical integration methods based on explicit and implicit Runge-Kutta schemes. These methods are useful to solve initial value problems of ordinary differential equations, even with interval initial value. We particularly present an elegant manner to bound the local truncation error of any Runge-Kutta methods. Moreover, our contractor based approach allows us to solve implicit RK without Newton-like method. Finally, our affine arithmetic provides us a framework which naturally counteracts wrapping effect. Runge-Kutta methods are already well-known for these stability properties that we can verify through a large benchmark of problems. The current state of the art in term of tool for validated ODE integration is Vnode-LP. If we are not competitive on time computation, the properties provided by our approach can easily balance this weakness for some problems.

Nevertheless, if our local truncation error computation is elegant, it is also too time consuming. An improvement in development is to base this computation on Butcher's series [START_REF] Ferenc | Computing of B-series by automatic differentiation[END_REF] which provides equivalent efficiency than automatic differentiation. This improvement will allow us to exploit higher order methods and be more accurate and fast.

. 3 )c

 3 The coefficients c i , a ij and b i fully characterize a Runge-Kutta method and they are usually given in a Butcher tableau c 1 a 11 a 12 . . . a 1s c 2 a 21 a 22 . . . a 2s . s a s1 a s2 . . . a ss b 1 b 2 . . . b s (3.4)

Figure 3 . 1 :

 31 Figure 3.1: Different kinds of Runge-Kutta methods

Example 3 . 1 . 1

 311 The elementary differentials f (f f , f ) is associated to the following rooted tree f f f f Definition 3.1.1 For a rooted tree τ , the elementary differential is a mapping F (τ ) : R n → R n defined recursively by F (•)(y) = f (y) is |τ | = 1 and

Remark 3 . 1 . 1

 311 The number of rooted trees increases very quickly, for example for |τ | = 11 the number of rooted trees is 1842.

Remark 4 . 1 . 1

 411 In Equation (4.3), if one can find a Runge-Kutta method such that for all τ with |τ | = p + 1 one has 0 < ϕ(τ )γ(τ ) < 1 then the LTE RK is eventually smaller than the LTE a Taylor series. Finding such Runge-Kutta method is an open problem. Proposition 4.1.1 show how a local truncation error of a Runge-Kutta method can be bounded using the a priori enclosure [ỹ j ] at time t j produced by Phase 1 (see Section 2.1). Proposition 4.1.1 Assuming an a priori enclosure [ỹ j ] on the interval [t j , t j+1 ] is given then one has ∀ξ

. 6 )

 6 can prove the existence and the uniqueness of the solution of Equation (1.1) using Theorem 2.1.1. From an algorithmic point of view, solving Equation (4.6), when implemented with implicit Runge-Kutta methods, require to embed Algorithm 1 into an iterative algorithm to compute the post fixed-point [r]. Nevertheless, the operator defined in Equation (4.5) will share many intermediate computations, such as LTE RK , with the algorithm for Phase 2 that the computational cost should remains low.

Figure 5 . 1 :

 51 Figure 5.1: Wrapping effect comparison (black: initial, green: interval, blue: interval from QR, red: zonotope from affine)

  while (ỹ 1 ⊂ ỹ0 ) and (iter < size(f ) + 1) do 11:

Figure 6 . 1 :

 61 Figure 6.1: Simulation of the circle system till t = 100s, with explicit method RK4 (Figure 3.1(a))

  parameter ρ given between 0.1 to 0.0001, and punctual initial states y 0 = (10, 0). The result in term of [ỹ j ] is plotted in Figure6.2, with the maximal stiffness.

Figure 6 . 2 :

 62 Figure 6.2: Simulation of the oil-reservoir system (stiffness=0.0001) till t = 50s, with implicit method LC3 (Figure 3.1(c))

Remark 6 . 3 . 1

 631 For Figure6.4, we truncate the results at diam = 50 for the clarity. It leads to the truncation of Vnode (fifteen times), LC3 (one time) and RK4 (one time).

Figure 6 . 3 :

 63 Figure 6.3: Results in term of normalized diameter gathered in spider graph for a simulation of 1 second, for the methods: RK4, LC3, LA3 and Vnode

Figure 6 . 4 :

 64 Figure 6.4: Results in term of normalized diameter gathered in spider graph for a simulation of 10 seconds, for the methods: RK4, LC3, LA3 and Vnode

Table 6 .

 6 1: Number of times a method produced the sharpest enclosure or the second sharpest enclosure.

	Method	c5w (1st)	c5w (2nd)	c5w (last)	c6w (1st)	c6w (2nd)	c6w (last)
	RK	103	35	8	58	39	8
	Vnode-LP	70	28	9	44	27	8
	Riot	36	11	0	24	12	2
	Valencia	3	3	49	3	2	49

It means that the stability of the numerical method does not depend on the value of the step-size h.

http://www.ibex-lib.org

see https://en.wikipedia.org/wiki/List_of_RungeKutta_methods for a detailed list of Runge-Kutta methods.

http://vericomp.inf.uni-due.de