Sylvie Boldo -François

Clément -Jean-Christophe Filliâtre

Micaela Mayero
email: micaela.mayero@lipn.univ-paris13.fr.

Guillaume Melquiond
email: guillaume.melquiond@inria.fr.

Pierre Weis
email: pierre.weis@inria.fr.

Sylvie Boldo

François Clément
email: pomdapi.francois.clement@inria.fr.

Jean-Christophe Filliâtre
email: christophe.filliatre@inria.fr.

Wave Equation Numerical Resolution: a Comprehensive Mechanized Proof of a C Program

Keywords: Programmation, vérification et preuves Formal proof of numerical program, Convergence of numerical scheme, Proof of C program, Coq formal proof, Acoustic wave equation, Partial differential equation, Rounding error analysis

We formally prove correct a C program that implements a numerical scheme for the resolution of the one-dimensional acoustic wave equation. Such an implementation introduces errors at several levels: the numerical scheme introduces method errors, and floating-point computations lead to round-off errors. We annotate this C program to specify both method error and round-off error. We use Frama-C to generate theorems that guarantee the soundness of the code. We discharge these theorems using SMT solvers, Gappa, and Coq. This involves a large Coq development to prove the adequacy of the C program to the numerical scheme and to bound errors. To our knowledge, this is the first time such a numerical analysis program is fully machine-checked.

Introduction

Ordinary differential equations (ODE) and partial differential equations (PDE) are ubiquitous in engineering and scientific computing. They show up in nuclear simulation, weather forecast, and more generally in numerical simulation, including block diagram modelization. Since solutions to nontrivial problems are non-analytic, they must be approximated by numerical schemes over discrete grids.

Numerical analysis is a part of applied mathematics that is mainly interested in proving the convergence of these schemes [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF], that is, proving that approximation quality increases as the size of discretization steps decreases. The approximation quality represents the distance between the exact continuous solution and the approximated discrete solution; this distance must tend toward zero in order for the numerical scheme to be useful.

A numerical scheme is typically proved to be convergent with pen and paper. This is a difficult, time-consuming, and error-prone task. Then the scheme is implemented as a C/C++ or Fortran program. This introduces new issues. First, we must ensure that the program correctly implements the scheme and is immune from runtime errors such as out-of-bounds accesses or overflows. Second, the program introduces round-off errors due to floating-point computations and we must prove that those errors could not lead to irrelevant results. Typical pen-and-paper proofs do not address floating-point nor runtime errors. Indeed the huge number of proof obligations, and their complexity, make the whole process almost intractable. However, with the help of mechanized program verification, such a proof becomes feasible. In the first place, because automated theorem provers can alleviate the proof burden. More importantly, because the proof is guaranteed to cover all aspects of the verification.

Our case study. We consider the acoustic wave equation in an one-dimensional space domain. The equation describes the propagation of pressure variations (or sound waves) in a fluid medium; it also models the behavior of a vibrating string. Among the wide variety of numerical schemes to approximate the 1D acoustic wave equation, we choose the simplest one: the second order centered finite difference scheme, also known as three-point scheme. To keep it simple, we assume an homogeneous medium (the propagation velocity is constant) and we consider discretization over regular grids with constant discretization steps for time and space. Our goal is to prove the correctness of a C program implementing this scheme.

Method and tools. We use the Jessie plug-in of Frama-C [START_REF] Marché | Jessie: an intermediate language for Java and C verification[END_REF][START_REF]The Frama-C platform for static analysis of C programs[END_REF] to perform the deductive verification of this C program. The source code is augmented with ACSL annotations [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF] to describe its formal specification. When submitted to Frama-C, proof obligations are generated. Once these theorems are proved, the program is guaranteed to satisfy its specification and to be free from runtime errors. Part of the proof obligations are discharged by automated provers, e.g. Alt-Ergo [START_REF] Bobot | The Alt-Ergo automated theorem prover[END_REF], CVC3 [START_REF] Barrett | CVC3[END_REF], Gappa [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF], and Z3 [START_REF] De | Z3, an efficient SMT solver[END_REF]. The more complicated ones, such as the one related to the convergence of the numerical scheme, cannot be proved automatically. These obligations were manually proved with the Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: The Calculus of Inductive Constructions[END_REF][START_REF]The Coq reference manual[END_REF] interactive proof assistant. In the end, we have formally verified all the properties of the C program. To our knowledge, this is the first time this kind of verification is machine-checked. The annotated C program and the Coq sources of the formal development are available from http://fost.saclay.inria.fr/wave_total_error.html

State of the art.

There is an abundant literature about the convergence of numerical schemes, e.g. see [START_REF] William | Numerical Partial Differential Equations: Finite Difference Methods[END_REF][START_REF] Zwillinger | Handbook of Differential Equations[END_REF]. In particular, the convergence of the three-point scheme for the wave equation is well-known and taught relatively early [START_REF] Bécache | Étude de schémas numériques pour la résolution de l'équation des ondes[END_REF]. Unfortunately, no article goes into all the details needed for a formal proof. These mathematical "details" may have been skipped for readability, but some mandatory details may have also been omitted due to oversights.

In the fields of automatic provers and proof assistants, few works have been done for the formalization and mechanical proofs of mathematical analysis, and even fewer works for numerical analysis. The first developments on real numbers and real analysis are from the late 90's, in systems such as ACL2 [START_REF] Gamboa | Nonstandard analysis in ACL2[END_REF], Coq [START_REF] Mayero | Formalisation et automatisation de preuves en analyses réelle et numérique[END_REF], HOL Light [START_REF] Harrison | Theorem Proving with the Real Numbers[END_REF], Isabelle [START_REF] Fleuriot | On the mechanization of real analysis in Isabelle/HOL[END_REF], Mizar [START_REF] Rudnicki | An overview of the MIZAR project[END_REF], and PVS [START_REF] Dutertre | Elements of mathematical analysis in PVS[END_REF]. An extensive work has been done by Harrison regarding R n and the dot product [START_REF] Harrison | A HOL theory of euclidean space[END_REF]. Constructive real analysis [START_REF] Geuvers | Constructive reals in Coq: Axioms and categoricity[END_REF][START_REF] Cruz-Filipe | A Constructive Formalization of the Fundamental Theorem of Calculus[END_REF][START_REF] Krebbers | Type classes for efficient exact real arithmetic in Coq[END_REF]] and further developments in numerical analysis [START_REF] Russell | Certified exact transcendental real number computation in Coq[END_REF][START_REF] Russell | A computer-verified monadic functional implementation of the integral[END_REF] have been carried out at Nijmegen. We can also mention the formal proof of an automatic differentiation algorithm [START_REF] Mayero | Using theorem proving for numerical analysis (correctness proof of an automatic differentiation algorithm)[END_REF].

As explained by Rosinger in 1985, old methods to bound round-off errors were based on "unrealistic linearizing assumptions" [START_REF] Elemer | Propagation of round-off errors and the role of stability in numerical methods for linear and nonlinear PDEs[END_REF]. Further work was done under more realistic assumptions about round-off errors [START_REF] Elemer | Propagation of round-off errors and the role of stability in numerical methods for linear and nonlinear PDEs[END_REF][START_REF] Elemer | L-convergence paradox in numerical methods for PDEs[END_REF], but none of these assumptions were proved correct with respect to the numerical schemes. As Roy and Oberkampf, we believe that round-off errors should not be treated as random variables and that traditional statistical methods should not be used [START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF]. They propose the use of interval arithmetic or increased precision to control accuracy. Note that their example of hypersonic nozzle flow is converging so fast that round-off errors can be neglected [START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF]. Interval arithmetic can also take method error into account [START_REF] Szyszka | An interval method for solving the one-dimensional wave equation[END_REF]. The final interval is then claimed to contain the exact solution. This is not formally proved, though. Additionally, the width of the final interval can be quite large.

There are other tools to bound round-off errors not dedicated to numerical schemes. Some successful approaches are based on abstract interpretation [START_REF] Cousot | The ASTR ÉE analyzer[END_REF][START_REF] Delmas | Towards an industrial use of FLUCTUAT on safety-critical avionics software[END_REF]. In our case, they are of little help, since there is a complex phenomenon of error compensation during the computations. Ignoring this compensation would lead to bounds on round-off errors growing as fast as O(2 k ) (k being the number of time steps). That is why we had to thoroughly study the propagation of round-off errors in this numerical scheme to get tighter bounds. It also means that most of the proofs had to be done by hand to achieve this part of the formal verification.

Outline. Section 2 presents the PDE, the numerical scheme, and their mathematical properties. Section 3 is devoted to the proofs of the convergence of the numerical scheme and the upper bound for the round-off error. Finally, Section 4 describes the formalization, i.e. the tools used, the annotated C program, and the mechanized proofs.

Numerical Scheme for the Wave Equation

A partial differential equation (PDE) modeling an evolution problem is an equation involving partial derivatives of an unknown function of several independent space and time variables. The uniqueness of the solution is obtained by imposing initial conditions, i.e. values of the function and some of its derivatives at initial time. The problem of the vibrating string tied down at both ends, among many other physical problems, is formulated as an initial-boundary value problem where the boundary conditions are additional constraints set on the boundary of the supposedly bounded domain [START_REF] William | Numerical Partial Differential Equations: Finite Difference Methods[END_REF].

This section, as well as the steps taken at Section 3.1 to conduct the convergence proof of the numerical scheme, is inspired by [START_REF] Bécache | Étude de schémas numériques pour la résolution de l'équation des ondes[END_REF].

The Continuous Equation

The chosen PDE models the propagation of waves along an ideal vibrating elastic string that is tied down at both ends, see [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF][START_REF] Brekhovskikh | Mechanics of Continua and Wave Dynamics[END_REF], see also Figure 1. The PDE is obtained from Newton's laws of motion [START_REF] Newton | Axiomata, sive Leges Motus[END_REF].

The gravity is neglected, so the string is supposed rectilinear when at rest. Let x min and x max be the abscissas of the extremities of the string. Let Ω = [x min , x max ] be the bounded space domain. Let p(x, t) be the transverse displacement of the point of the string of abscissa x at time t from its equilibrium position; it is a (signed) scalar. Let c be the constant propagation velocity; it is a positive number that depends on the section and density of the string. Let s(x, t) be the external action on the point of abscissa x at time t; it is a source term, such that t = 0 ⇒ s(x, t) = 0. Finally, let p 0 (x) and p 1 (x) be the initial position and velocity of the point of abscissa x. We

∀t ≥ 0, ∀x ∈ Ω, (L(c) p)(x, t) def = ∂ 2 p ∂t 2 (x, t) + A(c) p(x, t) = s(x, t), (1) ∀x ∈ Ω, (L 1 p)(x, 0) def = ∂p ∂t (x, 0) = p 1 (x), (2) ∀x ∈ Ω, (L 0 p)(x, 0) def = p(x, 0) = p 0 (x), (3) ∀t ≥ 0, p(x min , t) = p(x max , t) = 0 (4)
where the differential operator A(c) is defined by ( 5)

A(c) def = -c 2 ∂ 2 ∂x 2 .
This simple partial derivative equation happens to possess an analytical solution given by the so-called d'Alembert's formula [START_REF] Le Rond D'alembert | Recherches sur la courbe que forme une corde tendue mise en vibrations[END_REF], obtained from the method of characteristics and the image theory [START_REF] John | Partial Differential Equations[END_REF], ∀t ≥ 0, ∀x ∈ Ω, [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF]

p(x, t) = 1 2 (p 0 (x -ct) + p0 (x + ct)) + 1 2c x+ct x-ct p1 (y)dy + 1 2c t 0 x+c(t-σ) x-c(t-σ)
s(y, σ)dy dσ where the quantities p0 , p1 , and s are respectively the successive antisymmetric extensions in space of p 0 , p 1 , and s defined on Ω to the entire real axis R.

We have formally verified d'Alembert's formula as a separate work [START_REF] Lelay | Différentiabilité et intégrabilité en Coq. Application à la formule de d'Alembert[END_REF]. But for the purpose of the current work, we just admit that under reasonable conditions on the Cauchy data p 0 and p 1 and on the source term s, there exists a unique solution p to the initial-boundary value problem (1)-( 4) for each c > 0. Simply supposing the existence of a solution instead of exhibiting it, opens the way to scale our approach to more general cases for which there is no known analytic expression of a solution, e.g. in the case of a nonuniform propagation velocity c.

For such a solution p, it is natural to associate at each time t the positive definite quadratic quantity ( 7)

E(c)(p)(t) def = 1 2 x → ∂p ∂t (x, t) 2 + 1 2 (x → p(x, t)) 2 A(c)
where q, r def = Ω q(x)r(x)dx, q 2 def = q, q and q 2 A(c) def = A(c) q, q . The first term is interpreted as the kinetic energy, and the second term as the potential energy, making E the mechanical energy of the vibrating string.

The Discrete Equations

Let i max be the positive number of intervals of the space discretization. Let the space discretization step ∆x and the discretization function i ∆x be defined as

∆x def = x max -x min i max and i ∆x (x) def = x -x min ∆x .
Let us consider the time interval [0, t max ]. Let ∆t ∈]0, t max [ be the time discretization step. We define

k ∆t (t) def = t ∆t and k max def = k ∆t (t max ).
Now, the compact domain Ω × [0, t max ] is approximated by the regular discrete grid defined by

(8) ∀k ∈ [0..k max ], ∀i ∈ [0..i max ], x k i def = (x i , t k ) def = (x min + i∆x, k∆t).
For a function q defined over Ω × [0, t max ] (resp. Ω), the notation q h denotes any discrete approximation of q at the points of the grid, i.e. a discrete function over [0..i max ] × [0..k max ] (resp. [0..i max ]). By extension, the notation q h is also a shortcut to denote the matrix (q k i ) 0≤i≤imax,0≤k≤kmax (resp. the vector (q i ) 0≤i≤imax ). The notation qh is reserved to the approximation defined on [0..i max ] × [0..k max ] by Let p 0h and p 1h be two discrete functions over [0..i max ]. Let s h be a discrete function over [0..i max ] × [0..k max ]. Then, the discrete function p h over [0..i max ] × [0..k max ] is said to be the solution of the three-point 1 finite difference scheme, as illustrated in Figure 2, when the following set of equations holds:

qk i def = q(x k i ) (resp. qi def = q(x i )). x t x j+1 x j-1 x j t k+1 t k t k-1
(9) ∀k ∈ [2..k max ], ∀i ∈ [1..i max -1], (L h (c) p h ) k i def = p k i -2p k-1 i + p k-2 i ∆t 2 + (A h (c) (i ′ → p k-1 i ′ )) i = s k-1 i , 1
In the sense "three spatial points", for the definition of matrix A h (c).

∀i ∈

[1..i max -1], (L 1h (c) p h ) i def = p 1 i -p 0 i ∆t + ∆t 2 (A h (c) (i ′ → p 0 i ′ )) i = p 1,i , (10) ∀i ∈ [1..i max -1], (L 0h p h ) i def = p 0 i = p 0,i , (11)
∀k ∈ [0..k max ],

p k 0 = p k imax = 0 [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF] where the matrix A h (c), a discrete analog of A(c), is defined for any vector q h , by [START_REF] Boldo | Formal proof of a wave equation resolution scheme: the method error[END_REF] ∀i ∈

[1..i max -1], (A h (c) q h ) i def = -c 2 q i+1 -2q i + q i-1 ∆x 2 .
A discrete analog of the energy is also defined by2

(14) E h (c)(p h ) k+ 1 2 def = 1 2 i → p k+1 i -p k i ∆t 2 ∆x + 1 2 (i → p k i ), (i → p k+1 i ) A h (c)
where, for any vectors q h and r h ,

q h , r h ∆x def = imax i=0 q i r i ∆x, q h 2 ∆x def = q h , q h ∆x , q h , r h A h (c) def = A h (c) q h , r h ∆x , q h 2 A h (c) def = q h , q h A h (c) .
Note that the three-point scheme is parameterized by the discrete Cauchy data p 0h and p 1h , and by the discrete source term s h . Of course, when these discrete inputs are respectively approximations of the continuous functions p 0 , p 1 , and s (e.g. when p 0h = p0h , p 1h = p1h , and s h = sh ), then the discrete solution p h is an approximation of the continuous solution p.

Convergence

Let ξ be in ]0, 1[. The CFL(ξ) condition (for Courant-Friedrichs-Lewy, see [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]) states that the discretization steps satisfy the relation [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF] c∆t ∆x ≤ 1 -ξ.

The convergence error e h measures the distance between the continuous and discrete solutions. It is defined by [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF] ∀k

∈ [0..k max ], ∀i ∈ [0..i max ], e k i def = pk i -p k i .
Note that when p 0h = p0h , then for all i, e 0 i = 0. The truncation error ε h measures at which precision the continuous solution satisfies the numerical scheme. It is defined for k ∈ [2..k max ] and i ∈ [1..i max -1] by

ε k i def = (L h (c) ph ) k i -sk-1 i , (17)
ε 1 i def = (L 1h (c) ph ) i -p1,i , (18)
ε 0 i def = (L 0h ph ) i -p0,i . (19)
Again, note that when p 0h = p0h and p 1h = p1h , then for all i, ε 0 i = 0 and ε 1 i = e 1 i /∆t. Furthermore, when there is also s h = sh , then the convergence error e h is itself solution of the same numerical scheme with inputs defined by, for all i, k,

p 0,i = ε 0 i = 0, p 1,i = ε 1 i = e 1 i ∆t , and
s k i = ε k+1 i .
The numerical scheme is said to be convergent of order 2 if the convergence error tends toward zero at least as fast as ∆x 2 +∆t 2 when both discretization steps tend toward zero. 3 More precisely, the numerical scheme is said to be convergent of order (m,n) uniformly on the interval [0, t max ] if the convergence error satisfies4 [START_REF]The Coq reference manual[END_REF] i → e

k∆t(t) i ∆x = O [0,tmax] (∆x m + ∆t n ).
The numerical scheme is said to be consistent with the continuous problem at order 2 if the truncation error tends toward zero at least as fast as ∆x 2 + ∆t 2 when the discretization steps tend toward 0. More precisely, the numerical scheme is said to be consistent with the continuous problem at order (m, n) uniformly on interval [0, t max ] if the truncation error satisfies [START_REF] Coquand | Inductively defined types[END_REF] i → ε

k∆t(t) i ∆x = O [0,tmax] (∆x m + ∆t n ).
The numerical scheme is said to be stable if the discrete solution of the associated homogeneous problem (i.e. without any source term, s(x, t) = 0) is bounded independently of the discretization steps. More precisely, the numerical scheme is said to be stable uniformly on interval [0, t max ] if the discrete solution of the problem without any source term satisfies

(22) ∃α, C 1 , C 2 > 0, ∀t ∈ [0, t max ], ∀∆x, ∆t > 0, ∆x 2 + ∆t 2 < α ⇒ i → p k∆t(t) i ∆x ≤ (C 1 + C 2 t)( p 0h ∆x + p 0h A h (c) + p 1h ∆x ).
The result to be formally proved at Section 3.1 states that if the continuous solution p is regular enough on Ω × [0, t max ] and if the discretization steps satisfy the CFL(ξ) condition, then the three-point scheme is convergent of order (2, 2) uniformly on interval [0, t max ].

We do not admit (nor prove) the Lax equivalence theorem which stipulates that for a wide variety of problems and numerical schemes, consistency implies the equivalence between stability and convergence. Instead, we establish that consistency and stability implies convergence in the particular case of the one-dimensional acoustic wave equation.

Program

The main part of the C program is listed in Listing 1.

The grid steps ∆x and ∆t are respectively represented by the variables dx and dt, the grid sizes i max and k max by the variables ni and nk, and the propagation velocity c by the variable v. The initial position p 0h is represented by the function p0. The initial velocity p 1h and the source term s h are supposed to be zero and are not represented. The discrete solution p h is represented by the two-dimensional array p of size (i max + 1) × (k max + 1). (This is a simple naive implementation, a more efficient implementation would store only two time steps.)

To assemble the stiffness matrix A h (c), we only have to compute the square of the CFL coefficient c∆t ∆x (lines 1-2). Then, we recognize the space loops for the initial conditions: Equation (11) on lines 6-8, and Equation [START_REF] Bobot | The Alt-Ergo automated theorem prover[END_REF] with p 1h = 0 on lines 14-17. The space-time loop on lines 23-28 for the evolution problem comes from Equation [START_REF] Bertot | Canonical Big Operators[END_REF]. And finally, the boundary conditions of Equation [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF] are spread out on lines 9-10, 18-19, and 29-30.

Bounding Errors

Method Error

We first present the notions necessary to formalize and prove the method error. Then, we detail how the proof is conducted: we establish the consistency, the stability and prove that these two properties imply convergence in the case of the one-dimensional acoustic wave equation.

i<ni ; i++) { dp = p [ i +1][0] -2. * p [ i ] [ 0 ] + p [ i -1 ] [ 0 ] ; p [ i ] [ 1 ] = p [ i ] [ 0 ] + 0.
/ * Time i t e r a t i o n k = space loop . * / f o r ( i =1; i<ni ; i++) { dp = p [ i +1][k ] -2. * p [ i ] [ k ] + p [ i -1][k ] ; p [ i ] [ k+1] = 2. * p [ i ] [ k ] -p [ i ] [ k-1] + a * dp ; } / * Right boundary . * / p [ ni ] [ k+1] = 0 . ; }

Big O, Differentiability, and Regularity

When considering a big O equality a = O(b), one usually assumes that a and b are two expressions defined over the same domain and its interpretation as a quantified formula comes naturally. Here the situation is a bit more complicated. Consider

f (x, ∆x) = O(g(∆x))
when ∆x goes to 0. If one were to assume that the equality holds for any x ∈ R 2 , one would interpret it as

∀x, ∃α > 0, ∃C > 0, ∀∆x, ∆x ≤ α ⇒ |f (x, ∆x)| ≤ C • |g(∆x)|,
which means that constants α and C are in fact functions of x. Such an interpretation happens to be useless, since the infimum of α may well be zero while the supremum of C may be +∞.

A proper interpretation requires the introduction of a uniform big O relation with respect to the additional variable x:

(23) ∃α > 0, ∃C > 0, ∀x ∈ Ω x , ∀∆x ∈ Ω ∆x , ∆x ≤ α ⇒ |f (x, ∆x)| ≤ C • |g(∆x)|.
To emphasize the dependency on the two subsets Ω x and Ω ∆x , uniform big O equalities are now written

f (x, ∆x) = O Ωx,Ω∆x (g(∆x)).
We now precisely define the notion of "sufficiently regular" functions in terms of the full-fledged notation for the big O. The further result on the convergence of the numerical scheme requires that the solution of the continuous equation is actually sufficiently regular. We introduce two operators that, given a real-valued function f defined on the 2D plane and a point in the plane, return the values ∂f ∂x and ∂f ∂t at this point. Given these two operators, we can define the usual 2D Taylor polynomial of order n of a function f :

TP n (f, x) def = (∆x, ∆t) → n p=0 1 p! p m=0 p m • ∂ p f ∂x m ∂t p-m (x) • ∆x m • ∆t p-m .
Let Ω x ⊂ R 2 . We say that the previous Taylor polynomial is a uniform approximation of order n of f on Ω x when the following uniform big O equality holds:

f (x + ∆x) -TP n (f, x)(∆x) = O Ωx,R 2 ∆x n+1 .
A function f is then said to be sufficiently regular of order n uniformly on Ω x when all its Taylor polynomials of order smaller than n are uniform approximations of f on Ω x .

Consistency

The consistency of a numerical scheme expresses that, for ∆x small enough, the continuous solution taken at the points of the grid almost solves the numerical scheme. More precisely, we formally prove that when the continuous solution of the wave equation ( 1)-( 4) is sufficiently regular of order 4 uniformly on [x min , x max ] × [0, t max ], the numerical scheme ( 9)-( 12) is consistent with the continuous problem at order (2, 2) uniformly on interval [0, t max ] (see definition [START_REF] Coquand | Inductively defined types[END_REF] in Section 2.3). This is obtained using the properties of Taylor approximations; the proof is straightforward while involving long and complex expressions.

The key idea is to always manipulate uniform Taylor approximations that will be valid for all points of all grids when the discretization steps goes down to zero.

For instance, to take into account the initialization phase corresponding to Equation [START_REF] Bobot | The Alt-Ergo automated theorem prover[END_REF], we have to derive a uniform Taylor approximation of order 1 for the following continuous function (for any v sufficiently regular of order 3)

((x, t), (∆x, ∆t)) → v(x, t + ∆t) -v(x, t) ∆t - ∆t 2 c 2 v(x + ∆x, t) -2v(x, t) + v(x -∆x, t) ∆x 2 .
Note that the expression of this function involves both x + ∆x and x -∆x, meaning that we need a Taylor approximation which is valid for both positive and negative growths. The proof would have been impossible if we had required 0 < ∆x (as a space grid step) in the definition of the Taylor approximation.

In contrast with the case of an infinite string [START_REF] Boldo | Formal proof of a wave equation resolution scheme: the method error[END_REF], we do not need here a lower bound for c ∆t ∆x .

Stability

The stability of a numerical scheme expresses that the growth of the discrete solution is somehow bounded in terms of the input data (here, the Cauchy data u 0h and u 1h , and the source term s h ).

For the proof of the round-off error (see Section 3.2), we need a statement of the same form as definition ( 22) of Section 2.3. Therefore, we formally prove that, under the CFL(ξ) condition ( 15), the numerical scheme ( 9)-( 12) is stable uniformly on interval [0, t max ].

But, as we choose to prove the convergence of the numerical scheme by using an energetic technique 5 , it is more convenient to formulate the stability in terms of the discrete energy. More precisely, we also formally prove that under the CFL(ξ) condition ( 15), the discrete energy ( 14) satisfies the following overestimation,

E h (c)(p h ) k+ 1 2 ≤ E h (c)(p h ) 1 2 + √ 2 2 2ξ -ξ 2 • ∆t • k k ′ =1 i → s k ′ i ∆x
for all t ∈ [0, t max ] and with k = t ∆t -1. The evolution of the discrete energy between two consecutive time steps is shown to be proportional to the source term. In particular, the energy is constant when the source is inactive. Then, we obtain the following underestimation of the discrete energy, ∀k,

1 2 1 -c ∆t ∆x 2 i → p k+1 i -p k i ∆t ∆x ≤ E h (c)(p h ) k+ 1 2 .
Therefore, the non-negativity of the discrete energy is directly related to the CFL(ξ) condition. Note that this stability result is valid for any input data p 0h , p 1h , and s h .

Convergence

The convergence of a numerical scheme expresses the fact that the discrete solution gets closer to the continuous solution as the discretization steps go down to zero. More precisely, we formally prove that when the continuous solution of the wave equation ( 1)-( 4) is sufficiently regular of order 4 uniformly on [x min , x max ] × [0, t max ], and under the CFL(ξ) condition (15), the numerical scheme ( 9)-( 12) is convergent of order (2, 2) uniformly on interval [0, t max ] (see definition [START_REF]The Coq reference manual[END_REF] in Section 2.3). Firstly, we prove that the convergence error e h is itself the discrete solution of a numerical scheme of the same form but with different input data 6 . In particular, the source term (on the right-hand side) is here the truncation error ε h associated with the initial numerical scheme for p h . Then, the previous stability result holds, and we have an overestimation of the square root of the discrete energy associated with the convergence error E h (c)(e h ) that involves a sum of the corresponding source terms, i.e. the truncation error. Finally, the consistency result also makes this sum a big O of ∆x 2 + ∆t 2 , and a few more technical steps conclude the proof.

Round-off Error

As each operation is done with IEEE-754 floating-point numbers [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], round-off errors will occur and may endanger the accuracy of the final results. On this program, naive forward error analysis gives an error bound that is proportional to 2 k 2 -53 for the computation of a p k i . If this bound was sensible, it would cause the numerical scheme to compute only noise after a few steps. Fortunately, round-off error actually compensate themselves. To take into account the compensations and hence prove a usable error bound, we need a precise statement of the round-off error [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF] to exhibit the cancellations made by the numerical scheme.

Local Round-off Errors

Let δ k i be the (signed) floating-point error made in the two lines computing p k i (lines 26-27 in Listing 1). Floating-point values as computed by the program will be underlined: a, p k i to distinguish them from the discrete values of previous sections. They match the expressions a and p[i][k] in the annotations, while a and p k i can be represented in the annotations by \exact(a) and \exact(p[i][k]), as described in Section 4.1.4.

The δ k i are defined as follow:

δ k+1 i = p k+1 i -(2p k i -p k-1 i + a × (p k i+1 -2p k i + p k i-1 )).
Note that the program explained in Section 2.4 gives us that

p k+1 i = fl 2p k i -p k-1 i + a × (p k i+1 -2p k i + p k i-1 )
where fl(•) means that all the arithmetic operations that appear between the parentheses are actually performed by floating-point arithmetic, hence a bit off.

In order to get a bound on δ k i , we need to have the range of p k i . For this bound to be usable in our correctness proof, we need the range to be [-2, 2]. We have proved this fact by using the bounds on the method error and the round-off error of all the p k and p k-1 .

To prove the bound on δ k i , we perform forward error analysis and then use interval arithmetic to bound each intermediate error. We prove that, for all i and k, we have |δ k i | ≤ 78 × 2 -52 for a reasonable error bound for a, that is to say |a -a| ≤ 2 -49 .

Convolution of Round-off Errors

Note that the global floating-point error ∆ k i = p k i -p k i depends not only on δ k i , but also on all the δ k-l i+j for 0 < l ≤ k and -l ≤ j ≤ l. Indeed round-off errors propagate along floatingpoint computations. Their contributions to ∆ k i , which are independent and linear (due to the structure of the numerical scheme), can be computed by performing a convolution with a function λ : (Z × Z) → R. This function λ represents the results of the numerical scheme when fed with a single unit value:

λ 0 0 = 1 ∀i = 0, λ 0 i = 0 λ 1 -1 = λ 1 1 = a λ 1 0 = 2(1 -a) ∀i ∈ {-1, 0, 1}, λ 1 i = 0 λ k i = a × (λ k-1 i-1 + λ k-1 i+1 ) + 2(1 -a) × λ k-1 i -λ k-2 i Theorem 1. ∆ k i = p k i -p k i = k l=0 l j=-l λ l j δ k-l i+j .
Details of the proof can be found in [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF], but this point of view using convolution is new. The proof mainly amounts to performing numerous tedious transformations of summations until both sides are proved equal.

The previous proof assumes that the double summation is correct for all (i ′ , k ′ ) such that k ′ < k. This would be correct if there was an unbounded set of i where p k i is computed. This is no longer the case for a finite string. Indeed, at the ends of the range (i = 0 or i max ), p k i and p k i are equal to 0, so ∆ k i has to be 0 too.

The solution is to define the successive antisymmetric extension in space (as is done for d'Alembert's formula in Section 2.1) and to use it instead of δ. This ensures that both ∆ k 0 and ∆ k imax are equal to 0. It does not have any consequence on the values of ∆ k i for 0 < i < i max .

Bound on the Global Round-off Error

The analytic expression of ∆ k i can be used to obtain a bound on the round-off error. We will need two lemmas for this purpose.

Lemma 1. +∞ i=-∞ λ k i = k + 1.
Proof. We have

+∞ i=-∞ λ k+1 i = 2ǎ +∞ i=-∞ λ k i + 2(1 -ǎ) +∞ i=-∞ λ k i - +∞ i=-∞ λ k-1 i = 2 +∞ i=-∞ λ k i - +∞ i=-∞ λ k-1 i .
The sum by line verifies a simple linear recurrence. As λ 0 i = 1 and

λ 1 i = 2, we have λ k i = k + 1. Lemma 2. λ k i ≥ 0. Proof.
The demonstration was found out by M. Kauers and V. Pillwein.

If we denote by P the Jacobi polynomial, we have

λ j n = n k=j 2k j + k n + k + 1 2k + 1 (-1) j+k a k = a j n-j k=0 P (2j,0) k (1 -2a)
Now the conjecture follows directly from the inequality of Askey and Gasper [START_REF] Askey | Certain rational functions whose power series have positive coefficients[END_REF], which asserts that n k=0 P (r,0) k (x) > 0 for r > -1 and -1 < x ≤ 1 (see Theorem 7.4.2 in The Red Book [START_REF] Andrews | Special functions[END_REF]).

Theorem 2. ∆ k i = p k i -p k i ≤ 78 × 2 -53 × (k + 1) × (k + 2).
Proof. According to Theorem 1, ∆ k i is equal to k l=0 l j=-l λ l j δ k-l i+j . We know that for all j and l, |δ l j | ≤ 78 × 2 -52 and that λ l i = l + 1. Since the λ k i are nonnegative, the error is easily bounded by 78 × 2 -52 × k l=0 (l + 1).

Total Error

Let E h be the total error. It is the sum of the method error (or convergence error) e h of Sections 2.3 and 3.1.4, and of the round-off error ∆ h of Section 3.2.

From Theorem 2, we can estimate 7 the following upper bound for the spatial norm of the round-off error when ∆x ≤ 1 and ∆t ≤ t max /2: for all t ∈ [0, t max ],

i → ∆ k∆t(t) i ∆x = imax i=0 ∆ k∆t(t) i 2 ∆x ≤ (i max + 1)∆x × 78 × 2 -53 × t max ∆t + 1 × t max ∆t + 2 ≤ √ x max -x min + 1 × 78 × 2 -53 × 3 × t 2 max ∆t 2 . 7 When tmax ∆t ≥ 2, we have tmax ∆t + 1 tmax ∆t + 2 ≤ 3 t 2 max ∆t 2 .

RR n°7826

Thus, from the triangular inequality for the spatial norm, we obtain the following estimation of the total error:

∀t ∈ [0, t max ], ∀∆x, ∆x ≤ min(α e , α ∆ ) ⇒ i → E k∆t(t) i ∆x ≤ C e (∆x 2 + ∆t 2 ) + C ∆ ∆t 2
where the convergence constants α e and C e were extracted from the Coq proof (see Section 3.1.4) and are given in terms of the constants for the Taylor approximation of the exact solution at degree 3 (α 3 and C 3 ), and at degree 4 (α 4 and C 4 ) by

α e = min(1, t max , α 3 , α 4 )
,

C e = 2µt max √ x max -x min C ′ √ 2 + µ(t max + 1)C ′′ with µ = √ 2 √ 2ξ-ξ 2 , C ′ = max(1, C 3 + c 2 C 4 + 1)
, and

C ′′ = max(C ′ , 2(1 + c 2 )C 4 )
, and where the round-off constants α ∆ and C ∆ , as explained above, are given by

α ∆ = min(1, t max /2), C ∆ = 234 × 2 -53 × t 2 max √
x max -x min + 1. To give an idea of the relative importance of both errors, we consider the academic case where the space domain is the interval [0, 1], the velocity of waves is c = 1, and there is no initial velocity (u 1 (x) = 0) nor source term (s(x, t) = 0). We suppose that the initial position is given by u 0 (x) = χ(2(x -x 0 )/l) where x 0 = 0.5, l = 0.25, and χ is the C 4 function defined on [-1, 1] by χ(z) = (cos( π 2 z)) 5 , and with null continuation on the real axis. For this function, we may take

α 3 = α 4 = √ 2/2, C 3 = 5120 √ 2
, and C 4 = 409600/3. The corresponding solution presents two hump-shaped signals that propagate in each direction along the string, see Figure 1.

The upper bound on the total error is represented in Figure 3. Note that everything is in logarithmic scale. Of course, decreasing the size of the grid step decreases the method error, but in the same time, it increases the round-off error. Hence, the existence of a minimum for the upper bound on the total error (about 0.02 in our test case), corresponding to optimal grid step sizes.

Fortunately, the effective total error usually happens to be much smaller than this upper bound (by about a factor of 10 6 in our example).

Even if the effective total error on this example is off by several orders of magnitude with respect to the theoretical bound, this experiment is still reassuring. First, the left side of Figure 3 shows that the optimal choice (the darker part) for choosing ∆x and ∆t is reached near the limit of the CFL condition. This property matches common knowledge from numerical analysis. Second, the right side shows that both the effective error and the theoretical error have the same asymptotic behavior. So the properties we have verified in this work are not intrinsically easier than the best theorems one could state. It is just that the constants of the formulas extracted from the proofs (which we did not tune for this specific purpose) are not optimal for this example.

Mechanization of Proofs

In Sections 3.1 and 3.2, we have mostly described the method and round-off errors introduced when solving the wave equation problem with the given numerical scheme. We do not yet know whether this formalization actually matches the program described in Section 2.4 and fully given in Appendix A. In addition, the program might contain programming errors like out-of-bound accesses, which would possibly be left unattended while comparing the program and its formalization.

To fully verify the program, our process is as follows. First, we annotated the C program with comments specifying its behavioral properties, that is, what the program is supposed to compute. Second, we let Frama-C/Why generate proof obligations that state that the program matches its specification and that its execution is safe. Third, we used automated provers and Coq to prove all of these obligations.

Section 4.1 presents all the tools we have used for verifying the C program. Then Section 4.2 explains how the program was annotated. Finally, Section 4.3 shows how we proved all the obligations, either automatically or with a proof assistant.

Tools

Several software packages are used in this work. The formal proof of the method error has been made in Coq. The formal proof of the round-off error has been made in Coq, and using the Gappa tactic. The certification of the C program has used Frama-C (with the Jessie plug-in), and to prove the produced goals, we used Gappa, SMT provers, and the preceding Coq proofs. This section is devoted to present these tools and necessary libraries.

Coq

Coq8 is a formal system that provides an expressive language to write mathematical definitions, executable algorithms, and theorems, together with an interactive environment for proving them [START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: The Calculus of Inductive Constructions[END_REF]. Coq's formal language is based on the Calculus of Inductive Constructions [START_REF] Coquand | Inductively defined types[END_REF] that combines both a higher-order logic and a richly-typed functional programming language. Coq allows to define functions or predicates, that can be evaluated efficiently, to state mathematical theorems and software specifications, and to interactively develop formal proofs of these theorems. These proofs are machine-checked by a relatively small kernel, and certified programs can be extracted from them to external programming languages like Objective Caml, Haskell, or Scheme [START_REF] Letouzey | A new extraction for Coq[END_REF].

As a proof development system, Coq provides interactive proof methods, decision and semidecision algorithms, and a tactic language for letting the user define its own proof methods. Connection with external computer algebra system or theorem provers is also available.

The Coq library is structured into two parts: the initial library, which contains elementary logical notions and data-types, and the standard library, a general-purpose library containing various developments and axiomatizations about sets, lists, sorting, arithmetic, real numbers, etc.

In this work, we mainly use the Reals standard library [START_REF] Mayero | Formalisation et automatisation de preuves en analyses réelle et numérique[END_REF], that is a classical axiomatization of an Archimedean ordered complete field. We chose Reals to make our numerical proofs because we do not need an intuitionistic formalization.

For floating-point numbers, we use a large Coq library9 initially developed in [START_REF] Daumas | A generic library for floating-point numbers and its application to exact computing[END_REF] and extended with various results afterwards [START_REF] Boldo | Preuves formelles en arithmétiques à virgule flottante[END_REF]. It is a high-level formalization of IEEE-754 with gradual underflow. This is expressed by a formalization where floating-point numbers are pairs (n, e) associated with real values n × β e . The requirements for a number to be in the format (e min , β p ) are |n| < β p and e min ≤ e.

This formalization is convenient for human interactive proofs as testified by the numerous proofs using it. The huge number of lemmas available in the library (about 1400) makes it suitable for a large range of applications. This library has since then been superseded by the Flocq library [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF], but it was not yet available at the time we proved the floating-point results of this work.

Frama-C, Jessie, Why, and the SMT Solvers

We use the Frama-C platform10 to perform formal verification of C programs at the sourcecode level. Frama-C is an extensible framework that combines static analyzers for C programs, written as plug-ins, within a single tool. In this work, we use the Jessie plug-in for deductive verification. C programs are annotated with behavioral contracts written using the ANSI C Specification Language (ACSL for short) [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF]. The Jessie plug-in translates them to the Jessie language [START_REF] Marché | Jessie: an intermediate language for Java and C verification[END_REF], which is part of the Why verification platform [START_REF] Filliâtre | The Why/Krakatoa/Caduceus platform for deductive program verification[END_REF]. This part of the process is responsible for translating the semantics of C into a set of Why logical definitions (to model C types, memory heap, etc.) and Why programs (to model C programs). Finally, the Why platform computes verification conditions from these programs, using traditional techniques of weakest preconditions, and emits them to a wide set of existing theorem provers, ranging from interactive proof assistants to automated theorem provers. In this work, we use the Coq proof assistant (Section 4.1.1), SMT solvers Alt-Ergo [START_REF] Sylvain Conchon | CC(X): Semantical combination of congruence closure with solvable theories[END_REF], CVC3 [START_REF] Barrett | CVC3[END_REF] and Z3 [START_REF] De | Z3, an efficient SMT solver[END_REF], and the automated theorem prover Gappa (Section 4.1.3). Details about automated and interactive proofs can be found in Section 4.3. The dataflow from C source code to theorem provers can be depicted as follows:

ACSL-annotated C program

Frama-C (Jessie plug-in) Why

Coq Alt-Ergo CVC3 Z3

Gappa

More precisely, to run the tools on a C program, we use a graphical interface called gWhy. A screenshot is in Appendix B. In this interface, we may call one prover on one or on many goals. We then get a graphical view of how many goals are proved and by which prover.

In ACSL, annotations are using first-order logic. Following the programming by contract approach, the specifications involve preconditions, postconditions, and loop invariants. Contrary to other approaches focusing on run-time assertion checking, ACSL specifications do not refer to C values and functions, even if pure, but refer instead to purely logical symbols. In the following contract for a function computing the square of an integer x //@ ensures \ r e s u l t == x * x ; i n t square ( i n t x ) ; the postcondition, introduced with ensures, refers to the return value \result and argument x. Both are denoting mathematical integer values, for the corresponding C values of type int. In particular, x * x cannot overflow. Of course, one could give function square a more involved specification that handles overflows, e.g. with a precondition requiring x to be small enough. Simply speaking, we can say that C integers are reflected within specifications as mathematical integers, in an obvious way. The translation of floating-point numbers is more subtle and explained in Section 4.1.4.

Gappa

The Gappa tool11 adapts the interval-arithmetic paradigm to the proof of properties that occur when verifying numerical applications [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF]. The inputs are logical formulas quantified over real numbers whose atoms are usually enclosures of arithmetic expressions inside numeric intervals. Gappa answers whether it succeeded in verifying it. In order to support program verification, one can use rounding functions inside expressions. These unary operators take a real number and return the closest real number in a given direction that is representable in a given binary floating-point format. For instance, assuming that operator • rounds to the nearest binary64 floating-point number, the following formula states that the relative error of the floating-point addition is bounded:

∀x, y ∈ R, ∃ε ∈ R, |ε| ≤ 2 -53 ∧ •(•(x) + •(y)) = (•(x) + •(y)) × (1 + ε).
Converting straight-line numerical programs to Gappa logical formulas is easy and the user can provide additional hints if the tool were to fail to verify a property. The tool is specially designed to handle codes that are performing convoluted manipulations. For instance, it has been successfully used to verify a state-of-the-art library of correctly-rounded elementary functions [START_REF] Florent De Dinechin | Certifying the floatingpoint implementation of an elementary function using Gappa[END_REF]. In the current work, Gappa has been used to check much simpler properties. (In particular, no user hint was needed to discharge a proof automatically.) But the length of their proofs would discourage even the most dedicated users if they were to be manually handled. One of the properties is the round-off error of a local evaluation of the numerical scheme (Section 3.2.1). Other properties mainly deal with proving that no exceptional behavior occurs while executing the program: due to the initial values, all the computed values are sufficiently small to never cause overflow.

The verification of some formulas requires reasonings that are so long and intricate [START_REF] Florent De Dinechin | Certifying the floatingpoint implementation of an elementary function using Gappa[END_REF], that it might cast some doubts on whether an automatic tool actually succeeded in proving them. This is especially true when the tool ends up proving a property stronger than what the user expected. That is why Gappa also generates a formal certificate that can be mechanically checked by a proof assistant. This feature has served as the basis for a Coq tactic for automatically solving goals related to floating-point and real arithmetic [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF]. The tactic reads the current Coq goal, generates a Gappa goal, executes Gappa on it, recovers the certificate, and converts it to a complete proof term that Coq matches against the current goal. At this point, whether Gappa is correct or not no longer matters: the original Coq goal is formally proved by a complete Coq proof.

This tactic has been used whenever a verification condition would have been directly proved by Gappa, if not for some confusing notations or encodings of matrix elements. We just had to apply a few basic Coq tactics to put the goal into the proper form and then call the Gappa tactic to discharge it automatically.

Floating-Point Formalizations

A natural question is the link between the various representations of floating-point numbers. We assume that the execution environment (mostly the processor) complies with the IEEE-754 standard [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], which defines formats, rounding modes, and operations. The C program we consider is compiled in an assembly code that will directly use these formats and operations. We also assume that the compiler optimizations preserve the visible semantics of floating-operations from the original code, e.g. no use of the extended registers. Such optimizations could have been taken into account though, but at a cost [START_REF] Boldo | Proofs of numerical programs when the compiler optimizes[END_REF].

When verifying the C program, the floating-point operations are translated by Frama-C/-Jessie/Why following some previous work by two of the authors [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF]. A floating-point number f is modeled in the logic as a triple of real numbers (r, e, m). Value r simply stands for the real number that is immediately represented by f ; value e stands for the exact value of f , as obtained if no rounding errors had occurred; finally, value m stands for the model of f , which is a placeholder for the value intended to be computed and filled by the user. The two latter values have no existence in the program, but are useful for the specification and the verification. In particular, they help state assertions about the rounding or the model error of a program. In ACSL, the three components of the model of a floating-point number f can be referred to using f, \exact(f), and \model(f), respectively. \round error(f) is a macro for the rounding error, that is, \abs(f -\exact(f)).

For instance, the following excerpt from our C program specifies the error on the content of the dx variable, which represents the grid step ∆x (see Section 2). dx = 1./ ni ; / * @ a s s e r t @ dx > 0 . && dx <= 0.5 && @ \abs (\ exact (dx )dx ) / dx <= 0x1 . p-53; @ * / Note that 0x1.p-53 is a valid ACSL (and C too) literal meaning 2 -53 . Proof obligations are extracted from the annotated C program by computing weakest preconditions and then translated to automated and interactive provers. For SMT provers, the three fields r, e, and m, of floating-point numbers are expressed as real numbers and operations on floatingpoint numbers are uninterpreted relations axiomatized with basic properties such as bounds on the rounding error or monotonicity. For Gappa too, the fields are seen as real numbers. The tool, however, knows about floating-point arithmetic and its relation to real arithmetic. So floating-point operations are translated to the corresponding symbols from Gappa.

For Coq, we use the formalization described in Section 4.1.1 with a limited precision and gradual underflow (so that subnormal numbers are correctly translated). It is based on the real numbers of the standard library, which are also used for the translation of the exact and the model parts of the floating-point number.

While the IEEE-754 standard defines infinities and Not-a-Number as floating-point values, our translation does not take them into account. This does not compromise the correctness of the translation though, as each operation has a precondition that raises a proof obligation to guarantee that no exceptional events occur, such as overflow or division by zero, and therefore no infinities nor Not-a-Number are produced by the program.

To summarize, there is one assumption about the actual arithmetic being executed (IEEE-754 compliant and no overly aggressive optimizations from the compiler) and three formalizations of floating-point arithmetic used to verify the program: one used by Jessie/Why and then sent to the SMT solvers, one used by Gappa, and one used by Coq. The combination of these three different formalizations does not introduce any inconsistency. Indeed, we have formally proved in Coq that Gappa's and Coq's formalizations are equivalent for floating-point formats with limited precision and gradual underflow, that is, IEEE-754 formats. We have also formally proved that the Jessie/Why specifications and the properties for SMT provers are compatible with these formalizations, including the absence of special values (infinity or Not-a-Number) and the possibility to disregard the upper bound on reals representing floating-point numbers.

In fact, there is a fourth formalization of floating-point arithmetic involved, which is the one used internally by the interval computations of Gappa for proving results about real-valued expressions. It is not equivalent to the previous ones, since it is a multi-precision arithmetic, but it has no influence whatsoever on the formalization that Gappa uses for modeling floating-point properties.

Program Annotations

The full annotations are given in Appendix A. We give here hints about how to specify this program.

There are two axiomatics. The first one corresponds to the mathematics: the exact solution of the wave equation and its properties. It defines the needed values (the exact solution p, and its initialization p 0 ). We here assume that s and p 1 are zero functions. It also defines the derivatives of p (psol 1 , first derivative for the first variable of p, and psol 11 , second derivative for the first variable, and psol 2 and psol 22 for the second variable) as functions such that their value is the limit of p(x+∆x)-p(x) ∆x when ∆x → 0. As the ACSL annotations are only first order, these definitions are quite cumbersome: each derivative needs 5 lines to be defined.

We also put as axioms the fact that the solution has the expected properties [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF][START_REF] Andrews | Special functions[END_REF][START_REF] Askey | Certain rational functions whose power series have positive coefficients[END_REF][START_REF] Avigad | A Decision Procedure for Linear "Big O" Equations[END_REF]. The last property needed on the exact solution is its regularity. We require it to be near its Taylor approximations of degrees 3 and 4 on the whole interval [x min , x max ]. For instance, the following annotation states the property for degree 3.

/ * @ axiom p s o l s u f f r e g u l a r 3 : @ 0 < alpha 3 && 0 < C 3 && @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; \ f o r a l l r e a l dx ; \ f o r a l l r e a l dt ;

@ 0 <= x <= 1 == > \ s qr t (dx * dx + dt * dt ) <= alpha 3 == > @ \abs ( ps o l ( x + dx , t + dt ) -ps o l T a y l o r 3 (x , t , dx , dt )) <= @ C 3 * \abs (\pow(\ s qr t (dx * dx + dt * dt ) , 3 ) ) ; @ * /
The second axiomatic corresponds to the properties and loop invariant needed by the program. For example, we require the matrix to be separated: it means that a line of the matrix should not mix with another line (or a modification could alter another point of the matrix). We also state the existence of the loop invariant analytic error that is needed for applying the results of Section 3.2.

The initializations functions are specified, but not stated. This corresponds firstly to the function array2d alloc that initializes the matrix and p zero that produces an approximation of the p 0 function. Our program verification is modular: our proofs are generic with respect to p 0 and its implementation.

The preconditions of the main functions are the following ones:

• i max and k max must be greater than one, but small enough so that i max + 1 and k max + 1 do not overflow;

• the grid sizes ∆x must fulfill some mathematical conditions that are required for the convergence of the scheme;

• the floating-point values computed for the grid sizes must be near their mathematical values;

• to prevent exceptional behavior in the computation of a, the time discretization step must be greater than 2 -1000 and c∆t ∆x must be greater than 2 -500 . There are two postconditions, corresponding to the method and round-off errors. See Sections 3.1 and 3.2 for more details.

Automation and Manual Proofs

This section is devoted to formal specifications and proofs corresponding to the bounds proved in Section 3. We give some key points of the automated proofs.

Big O. In section 3.1.1, we present two interpretations of the big O notation. Usual mathematical pen-and-paper proofs switch from one interpretation to the other depending on which one is the most adapted, without noticing that they may not be equivalent. The formal development was helpful in bringing into light the erroneous reasoning hidden by the usage of big O notations. We introduced the notion of uniform big O in [START_REF] Boldo | Formal proof of a wave equation resolution scheme: the method error[END_REF] in the context of an infinite string. In the present paper, we consider the case of the finite string, hence for compactness reasons, both notions are in fact equivalent. However, we still use the more general uniform big O notion to share most of the proof developments between the finite and the infinite cases. Regarding automation, a decision procedure has been developed in [START_REF] Avigad | A Decision Procedure for Linear "Big O" Equations[END_REF]; unfortunately, those results were not applicable since we needed a more powerful big O.

Differential operators.

As long as we were studying only the method error, we did not have to define the differential operators nor assume anything about them [START_REF] Boldo | Formal proof of a wave equation resolution scheme: the method error[END_REF]. Their only properties appeared through their usage: function p is a solution of the partial differential equation and it is sufficiently regular. This is no longer possible for the annotated C program. Indeed, due to the underlying logic, the annotations have to define p as a solution of the PDE by using first-order formulas stating differentiability, instead of second-order formulas involving differential operators. Since the formalization of Taylor approximations has been left unchanged, the natural way to relate the C annotations with the Coq development is therefore to define the operators as actual differential operators. Note that this has forced us to introduce a small axiom. Indeed, our definition of Taylor approximation depends on differential operators that are total functions, while Coq's standard library defines only partial operators. So we have assumed the existence of some total operators that are equal to the partial ones whenever applied to differentiable functions. The axiom states absolutely nothing about the result of these operators for nondifferentiable functions, so no inconsistencies are introduced this way. This is just a specific instance of Hilbert ε operator [START_REF] Zach | Hilbert's "Verunglueckter Beweis[END_REF], which does not make the logic inconsistent [START_REF] Lee | Proof-irrelevant model of CC with predicative induction and judgmental equality[END_REF].

Method error. The Coq proof of the method error is about 5000-line long. About half of it is dedicated to the wave equation and the other half is re-usable (definition and properties of the dot product, the big O, Taylor expansions. . . ). We formally proved without any axiom that the numerical scheme is convergent of order 2, which is the known mathematical result. An interesting aspect of the formal proof in Coq is that we were able to extract the constants α and C appearing in the big O for the convergence result in order to obtain their precise values. The recursive extraction was fully automatic after making explicit some inlining. The mathematical expressions are given in Section 3.3.

Round-off errors. Except for Lemma 2, all the proofs described in section 3.2 have been done and machined-checked using Coq. In particular, the proof of the bound on δ k i was done automatically by calling Gappa from Coq. Lemma 2 is a technical detail compared to the rest of our work, that is not worth the immense Coq development it would require: keen results on integrals but also definitions and results about the Legendre, Laguerre, Chebychev, and Jacobi polynomials.

The program proof. Given the program code, the Why tool generates 149 verification conditions that have to be proved. While possible, proving all of them in Coq would be rather tedious. Moreover, it would lead to a rather fragile construct: any later modification to the code, however small it is, would cause different proof obligations to be generated, which would then require additional human interaction to adapt the Coq proofs. We prefer to have automated provers (SMT solvers and Gappa) discharge as many of them as possible, so that only the most intricate ones are left to be proven in Coq. The following table shows how many goals are discharged automatically and how many are left to the user. On safety goals (matrix access, loop variant decrease, overflow), automatic provers are helpful: they prove about 90 % of the goals. On behavior goals (loop invariant, assertion, postcondition), automatic provers succeed for half of the goals. As our loop invariant involves an uninterpreted predicate, the automatic provers cannot prove all the behavior goals (they would have been too complicated anyway). That is why we resort to an interactive higher-order theorem prover, namely Coq.

Coq proofs are split into two sets: first, the mathematical proof of convergence and second, the proofs of bounded round-off errors and absence of runtime errors. Appendix C displays the layout of the Coq formalization.

The following tabular gives the compilation times of the Coq files on a 3-GHz dual core machine. Note that most theorem statements regarding round-off and runtime errors are automatically generated (7 321 lines out of 7 737) by the Frama-C/Jessie/Why framework.

The compilation time may seem prohibitive; it is mainly due to the size of the theorems and to calls to the omega decision procedure for Presburger arithmetic. The difficulty does not lie in the arithmetic statement itself, but rather in a large number of useless hypotheses. In order to reduce the compilation time, we could manually massage the hypotheses to speed up the procedure, but this would defeat the point of using an automatic tactic.

Conclusion

In the end, having formally verified the C program means that all of the proof obligations generated by Frama-C/Jessie/Why have been proved, either by automated tools or by Coq formal proofs. These formal proofs depend on some axioms specific to this work: the fact about Jacobi polynomials, the existence and regularity of a solution to the EDP, and the existence of differential operators. The last two have been tackled by subsequent works, which means that the only remaining Coq axiom is the one about Jacobi polynomials.

We succeeded in verifying a C program that implements a numerical scheme for the resolution of the one-dimensional acoustic wave equation. This is comprised of three sets of proofs. First we formalized the wave equation and proved the convergence of a scheme for its numerical resolution. Second we proved that the C program behaves safely: no out-of-bound array accesses and no overflow during floating-point computations. Third we proved that the round-off errors are not causing the numerical results to go astray. This is the first verification of this kind of program that covers all its aspects, both mathematics and implementation.

This work shows a tight synergy between researchers from applied mathematics and logic. Three domains are intertwined here: applied mathematics for an initial proof that was enriched and detailed upon request, computer arithmetic for smart bounds on round-off errors, and formal methods for machine-checking them. This may be the reason why such proofs never appeared before, as this kind of collaboration is uncommon.

Each proof came with its own hurdles. For ensuring the correct behavior of the program, the most tedious point was to prove that setting a result value did not cause other values to change, that is, that all the lines of the matrix are properly separated. In particular, verifying the loop invariant requires checking that, except for the new value, the properties of the memory are preserved. An unexpectedly tedious part was to check that the program actually complies with our mathematical model for the numerical scheme.

Another difficulty lies in the mathematical proof itself. We based our work on proofs found in books, courses, and articles. It appears that pen-and-paper proofs are sometimes sketchy: they may be fuzzy about the needed hypotheses, especially when switching quantifiers. We have also learned that filling the gaps may cause us to go back to the drawing board and to change the basic blocks of our formalization to make them more generic (e.g. devising a big O that needs to be uniform and also generic with respect to a property P ).

An unexpected side effect of having performed this formal verification in Coq is our ability to automatically extract the constants hidden inside the proofs. That way, we are able to explicitly bound the total error rather than just having the usual O(∆x 2 + ∆t 2 ) bound. In particular, we can compare the magnitudes of the method error and round-off error and then decide how to scale the discretization grid.

Coq could have offered us more: it would have been possible to describe and prove the algorithm directly in Coq. The same formalism would have been used all the way long, but we were more interested in proving a real-life program in a real-life language. This has shown us the difficulties lying in the memory handling for matrices. In the end, we have a C code with readable annotations instead of a Coq theorem and that seems more convincing to applied mathematicians.

For this exploratory work, we considered the simple three-point scheme for the one-dimensional wave equation. Further works involve scaling to higher-dimension. The one-dimensional case showed us that summations and finite support functions play a much more important role in the development than we first expected. We are therefore moving to the SSReflect interface and libraries for Coq [START_REF] Bertot | Canonical Big Operators[END_REF], so as to simplify the manipulations of these objects in the higher-dimensional case.

This example also exhibits a major cancellation of rounding errors and it would be interesting to see under which conditions numerical schemes behave so well.

Another perspective is to generalize our approach to other higher-order numerical schemes for the same equation, and to other PDEs. However, the proofs of Section 3.1 are entangled with particulars of the presented problem, and would therefore have to be redone for other problems. So a more fruitful approach would be to prove once and for all the Lax equivalence theorem that states that consistency implies the equivalence between convergence and stability. This would considerably reduce the amount of work needed for tackling other schemes and equations.

A Source Code

/ * @ axiomatic d i r i c h l e t m a t h s { @ @ l o g i c r e a l c ; @ l o g i c r e a l p0 ( r e a l x ) ; @ l o g i c r e a l ps o l ( r e a l x , r e a l t ) ; @ axiom c pos : 0 < c ; @ l o g i c r e a l ps o l 1 ( r e a l x , r e a l t ) ; @ axiom p s o l 1 d e f : @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; @ \ f o r a l l r e a l eps ; \ e x i s t s r e a l C; 0 < C && \ f o r a l l r e a l dx ; @ \abs (dx ) < C == > @ \abs (( ps o l ( x + dx , t )ps o l (x , t )) / dxps o l 1 (x , t )) < eps ; @ l o g i c r e a l psol 11 ( r e a l x , r e a l t ) ; @ axiom ps o l 1 1 de f : @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; @ \ f o r a l l r e a l eps ; \ e x i s t s r e a l C; 0 < C && \ f o r a l l r e a l dx ; @ \abs (dx ) < C == > @ \abs (( ps o l 1 ( x + dx , t )ps o l 1 (x , t )) / dxpsol 11 (x , t )) < eps ; @ l o g i c r e a l ps o l 2 ( r e a l x , r e a l t ) ; @ axiom p s o l 2 d e f : @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; @ \ f o r a l l r e a l eps ; \ e x i s t s r e a l C; 0 < C && \ f o r a l l r e a l dt ; @ \abs ( dt ) < C == > @ \abs (( ps o l (x , t + dt )ps o l (x , t )) / dtps o l 2 (x , t )) < eps ; @ l o g i c r e a l psol 22 ( r e a l x , r e a l t ) ; @ axiom ps o l 2 2 de f : @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; @ \ f o r a l l r e a l eps ; \ e x i s t s r e a l C; 0 < C && \ f o r a l l r e a l dt ; @ \abs ( dt ) < C == > @ \abs (( ps o l 2 (x , t + dt )ps o l 2 (x , t )) / dtpsol 22 (x , t )) < eps ; @ axiom wave eq 0 : \ f o r a l l r e a l x ; 0 <= x <= 1 == > ps o l (x , 0) == p0 ( x ) ; @ axiom wave eq 1 : \ f o r a l l r e a l x ; 0 <= x <= 1 == > ps o l 2 (x , 0) == 0 ; @ axiom wave eq 2 : @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; @ 0 <= x <= 1 == > psol 22 (x , t )c * c * psol 11 (x , t ) == 0 ; @ axiom w a v e e q di r i c hl e t 1 : \ f o r a l l r e a l t ; ps o l (0 , t ) == 0 ; @ axiom w a v e e q di r i c hl e t 2 : \ f o r a l l r e a l t ; ps o l (1 , t ) == 0 ; @ l o g i c r e a l ps o l T a y l o r 3 ( r e a l x , r e a l t , r e a l dx , r e a l dt ) ; @ l o g i c r e a l ps o l T a y l o r 4 ( r e a l x , r e a l t , r e a l dx , r e a l dt ) ; @ l o g i c r e a l alpha 3 ; l o g i c r e a l C 3 ; @ l o g i c r e a l alpha 4 ; l o g i c r e a l C 4 ; @ axiom p s o l s u f f r e g u l a r 3 : @ 0 < alpha 3 && 0 < C 3 && @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; \ f o r a l l r e a l dx ; \ f o r a l l r e a l dt ; @ 0 <= x <= @ l o g i c r e a l T max ; @ axiom T max pos : 0 < T max ; @ l o g i c r e a l C conv ; l o g i c r e a l alpha conv ; @ lemma alpha conv pos : 0 < alpha conv ; @ @ } * / / * @ axiomatic d i r i c h l e t p r o g { @ @ pr e di c a te a n a l y t i c e r r o r {L} @ ( double * * p , i n t e g e r ni , i n t e g e r i , i n t e g e r k , double a , double dt ) @ reads p [ . . ] [ . . ] ; @ @ lemma a n a l y t i c e r r o r l e {L} : @ \ f o r a l l double * * p ; \ f o r a l l i n t e g e r ni ; \ f o r a l l i n t e g e r i ; @ \ f o r a l l i n t e g e r nk ; \ f o r a l l i n t e g e r k ; @ \ f o r a l l double a ; \ f o r a l l double dt ;

@ 0 < ni == > 0 <= i <= ni == > 0 <= k == > @ 0 < \ exact ( dt ) == > @ a n a l y t i c e r r o r (p , ni , i , k , a , dt ) == > @ \ s qr t (1 . / ( ni * ni ) + \ exact ( dt ) * \ exact ( dt )) < alpha conv == > @ k <= nk == > nk <= 7598581 == > nk * \ exact ( dt ) <= T max == > @ \ exact ( dt ) * ni * c <= 1 -0x1 . p-50 == > @
\ f o r a l l i n t e g e r i 1 ; \ f o r a l l i n t e g e r k1 ;

@ 0 <= i 1 <= ni == > 0 <= k1 < k == > @ \abs (p [ i 1 ] [ k1 ]
) <= 2 ; @ @ pr e di c a te separated matrix {L}( double * * p , i n t e g e r l e n i ) = @ \ f o r a l l i n t e g e r i ; \ f o r a l l i n t e g e r j ; @ 0 <= i < l e n i == > 0 <= j < l e n i == > i != j == > @ \base addr (p [ i ] ) != \base addr (p [ j ] ) ; @ @ l o g i c r e a l sqr norm dx conv err {L} @ ( double * * p , r e a l dx , r e a l dt , i n t e g e r ni , i n t e g e r i , i n t e g e r k) @ reads p [ . . ] [ . . ] ; @ l o g i c r e a l sqr ( r e a l x ) = x * x ; @ lemma sqr norm dx conv err 0 {L} : @ \ f o r a l l double * * p ; \ f o r a l l r e a l dx ; \ f o r a l l r e a l dt ; @ \ f o r a l l i n t e g e r ni ; \ f o r a l l i n t e g e r k ; @ sqr norm dx conv err (p , dx , dt , ni , 0 , k) == 0 ; @ lemma sqr norm dx conv err succ {L} : @ \ f o r a l l double * * p ; \ f o r a l l r e a l dx ; \ f o r a l l r e a l dt ; @ \ f o r a l l i n t e g e r ni ; \ f o r a l l i n t e g e r i ; \ f o r a l l i n t e g e r k ; @ 0 <= i == > ) ) ; @ l o g i c r e a l norm dx conv err {L} @ ( double * * p , r e a l dt , i n t e g e r ni , i n t e g e r k) = @ \ s qr t ( sqr norm dx conv err (p , 1 . / ni , dt , ni , ni , k ) ) ; @ @ } * / / * @ r e q u i r e s l e n i >= 1 && l e n j >= 1 ; @ ensures @ \ v a l i d r a ng e (\ r e s u l t , 0 , l e n i -1) && @ (\ f o r a l l i n t e g e r i ; 0 <= i < l e n i == > @ \ v a l i d r a ng e (\ r e s u l t [ i ] , 0 , l e n j -1)) && @ separated matrix (\ r e s u l t , l e n i ) ; @ * / double * * a r r a y 2 d a l l o c ( i n t l e n i , i n t l e n j ) ; / * @ r e q u i r e s ( l != 0 ); @ ensures @ \ r o und e r r o r (\ r e s u l t ) <= 14 * 0x1 . p-52 && @ \ exact (\ r e s u l t ) == p0(\ exact ( x ) ) ; @ * / double p zero ( double xs , double l , double x ) ; / * @ r e q u i r e s @ ni >=

@

B Screenshot

This is a screenshot of gWhy: we have the list of all the verification conditions and if they are proved by the various automatic tools.

Figure 1 :

 1 Figure 1: Space-time representation of the signal propagating along a vibrating string. Each curve represents the string at a different time step.

Figure 2 :

 2 Figure 2: Three-point scheme: p k+1 i

  5 * a * dp ; } / * Right boundary . * / p [ ni ] [ 1 ] = 0 . ; / * Evolution problem and boundary c o ndi ti o ns . * / / * Propagation = time loop . * / f o r (k=1; k<nk ; k++) { / * Le f t boundary . * / p [ 0 ] [ k+1] = 0 . ;

Figure 3 :

 3 Figure 3: Upper bound for the total error in log-scale. Left: for ∆x and ∆t satisfying the CFL condition. The lighter area (in yellow) represents the higher values above 10 4 , whereas the darker area represents the lower values below 10 -1 . Right: for an optimal CFL condition with ∆t = 1-ξ c ∆x. The green crosses represent the effective total error computed by the C program for a few values of the space step.

  double xs , double l ) { / * Output v a r i a b l e . * / double * * p ; / * Local v a r i a b l e s . * / i n t i , k ; double a1 , a , dp , dx ; dx = 1./ ni ; / * @ a s s e r t @ dx > 0 . && dx <= 0.5 && @ \abs (\ exact (dx )dx ) / dx <= 0x1 . p-53; @ * / / * Compute the constant c o e f f i c i e n t of the s t i f f n e s s matrix . * / a1 = dt/dx * v ; a = a1 * a1 ; / * @ a s s e r t @ 0 <= a <= 1 && @ 0 < \ exact (a ) <= 1 && @ \ r o und e r r o r (a ) <= 0x1 . p-49; @ * / / * A l l o c a te space-time v a r i a b l e f o r the d i s c r e t e s o l u t i o n . * / p = a r r a y 2 d a l l o c ( ni +1, nk+1); / * F i r s t i n i t i a l c o ndi ti o n and boundary c o ndi ti o ns . * / / * Le f t boundary . * / p [ 0 ] [ 0 ] = 0 . ; / * Time i t e r a t i o n -1 = space loop . * / / * @ loop i n v a r i a n t @ 1 <= i <= ni && @ a n a l y t i c e r r o r (p , ni , i -1 , 0 , a , dt ) ; @ loop v a r i a nt nii ; * / f o r ( i =1; i<ni ; i++) { p [ i ] [ 0 ] = p zero ( xs , l , i * dx ) ; } / * Right boundary . * / p [ ni ] [ 0 ] = 0 . ; / * @ a s s e r t a n a l y t i c e r r o r (p , ni , ni , 0 , a , dt ) ; * / / * Second i n i t i a l c o ndi ti o n ( with p one=0) and boundary c o ndi ti o ns . * / / * Le f t boundary . * / p [ 0 ] [ 1 ] = 0 . ; / * Time i t e r a t i o n 0 = space loop . * / / *

Listing 1 :

 1 The main part of the C code, without annotations.

	/ * Compute the constant c o e f f i c i e n t of the s t i f f n e s s matrix . * /
	a1 = dt/dx * v ;
	a = a1 * a1 ;
	/ * F i r s t i n i t i a l c o ndi ti o n and boundary c o ndi ti o ns . * /
	/ * Le f t boundary . * /
	p [ 0 ] [ 0 ] = 0 . ;
	/ * Time i t e r a t i o n -1 = space loop . * /
	f o r ( i =1; i<ni ; i++) {
	p [ i ] [ 0 ] = p0 ( i * dx ) ;

} / * Right boundary . * / p [ ni ] [ 0 ] = 0 . ; / * Second i n i t i a l c o ndi ti o n ( with p1=0) and boundary c o ndi ti o ns . * / / * Le f t boundary . * / p [ 0 ] [ 1 ] = 0 . ; / * Time i t e r a t i o n 0 = space loop . * / f o r ( i =1;

  1 == > \ s qr t (dx * dx + dt * dt ) <= alpha 3 == >

	@ axiom p s o l l e :
	@	\ f o r a l l r e a l x ; \ f o r a l l r e a l t ;
	@	0 <= x <= 1 == > 0 <= t == > \abs ( ps o l (x , t )) <= 1 ;
	@	\abs ( ps o l ( x + dx , t + dt ) -ps o l T a y l o r 3 (x , t , dx , dt )) <=
	@	C 3 * \abs (\pow(\ s qr t (dx * dx + dt * dt ) , 3 ) ) ;

@ axiom p s o l s u f f r e g u l a r 4 : @ 0 < alpha 4 && 0 < C 4 && @ \ f o r a l l r e a l x ; \ f o r a l l r e a l t ; \ f o r a l l r e a l dx ; \ f o r a l l r e a l dt ; @ 0 <= x <= 1 == > \ s qr t (dx * dx + dt * dt ) <= alpha 4 == > @ \abs ( ps o l ( x + dx , t + dt )ps o l T a y l o r 4 (x , t , dx , dt )) <= @ C 4 * \abs (\pow(\ s qr t (dx * dx + dt * dt ) , 4 ) ) ;

  sqr norm dx conv err (p , dx , dt , ni , i + 1 , k) == @ sqr norm dx conv err (p , dx , dt , ni , i , k) + @ dx * sqr ( ps o l (1 . * i / ni , k * dt ) -\ exact (p [ i ] [ k ]

  * / double * * forward prop ( i n t ni , i n t nk , double dt , double v ,

		2 && nk >= 2 && l != 0 &&
	@	dt > 0 . && \ exact ( dt ) > 0 . &&
	@	\ exact ( v ) == c && \exact ( v ) == v &&
	@	0x1 . p-1000 <= \ exact ( dt ) &&
	@	ni <= 2147483646 && nk <= 7598581 &&
	@	nk

* \ exact ( dt ) <= T max && @ \abs (\ exact ( dt )dt ) / dt <= 0x1 . p-51 && @ 0x1 . p-500 <= \exact ( dt ) * ni * c <= 1 -0x1 . p-50 && @ \ s qr t (1 . / ( ni * ni ) + \exact ( dt ) * \ exact ( dt )) < alpha conv ; @ @ ensures @ \ f o r a l l i n t e g e r i ; \ f o r a l l i n t e g e r k ; @ 0 <= i <= ni == > 0 <= k <= nk == > @ \ r o und e r r o r (\ r e s u l t [ i ] [ k ] ) <= 78. / 2 * 0x1 . p-52 * (k + 1) * (k + 2 ); @ @ ensures @ \ f o r a l l i n t e g e r k ; 0 <= k <= nk == > @ norm dx conv err (\ r e s u l t , \ exact ( dt ) , ni , k) <= @ C conv * (1 . / ( ni * ni ) + \ exact ( dt ) * \ exact ( dt ) ) ; @

  } / * Right boundary . * / p [ ni ] [ 1 ] = 0 . ; / * @ a s s e r t a n a l y t i c e r r o r (p , ni , ni , 1 , a , dt ) ; * / / * Evolution problem and boundary c o ndi ti o ns . * / / * Propagation = time loop . * /

	@ loop i n v a r i a n t
	@	1 <= i <= ni &&
	@	a n a l y t i c e r r o r (p , ni , i -1 , 1 , a , dt ) ;
	@ loop v a r i a nt ni -i ; * /
	f o r ( i =1; i<ni ; i++) {
	/ * @ loop i n v a r i a n t
	@	1 <= k <= nk &&
	@	a n a l y t i c e r r o r (p , ni , ni , k , a , dt ) ;
	@ loop v a r i a nt nk -k ; * /
	f o r (k=1; k<nk ; k++) {
	/ * Le f t boundary . * /
	p [ 0 ] [ k+1] = 0 . ;
		dp = p [ i +1][k ] -2. * p [ i ] [ k ] + p [ i -1][k ] ;
		p [ i ] [ k+1] = 2. * p [ i ] [ k ] -p [ i ] [ k-1] + a * dp ;
	}
	/ * Right boundary . * /
	p [ ni ] [ k+1] = 0 . ;
	/ * @ a s s e r t a n a l y t i c e r r o r (p , ni , ni , k + 1 , a , dt ) ; * /
	}
	return p ;
	}

/ * @ a s s e r t \abs (p [ i -1][0]) <= 2 ; * / / * @ a s s e r t \abs (p [ i ] [ 0 ] ) <= 2 ; * / / * @ a s s e r t \abs (p

[ i +1][0]) <= 2 ; * / dp = p [ i +1][0] -2. * p [ i ] [ 0 ] + p [ i -1 ] [ 0 ] ; p [ i ] [ 1 ] = p [ i ] [ 0 ] + 0.5 *

a * dp ; / * Time i t e r a t i o n k = space loop . * / / * @ loop i n v a r i a n t @ 1 <= i <= ni && @ a n a l y t i c e r r o r (p , ni , i -1 , k + 1 , a , dt ) ; @ loop v a r i a nt nii ; * / f o r ( i =1; i<ni ; i++) { / * @ a s s e r t \abs (p [ i -1][k ] ) <= 2 ; * / / * @ a s s e r t \abs (p [ i ] [ k ] ) <= 2 ; * / / * @ a s s e r t \abs (p [ i +1][k ] ) <= 2 ; * / / * @ a s s e r t \abs (p [ i ] [ k-1]) <= 2 ; * /

RR n°7826

By convention, the energy is defined between steps k and k + 1, hence the notation k + 1 2 .RR n°7826

Note that ∆x tending toward 0 actually means that imax goes to infinity.

See Section 3.1.1 for the precise definition of the big O notation. INRIA

The popular alternative, using the Fourier transform, would have required huge additional Coq developments.

Of course, there is no associated continuous problem.RR n°7826

INRIA

http://coq.inria.fr/ RR n°7826

http://lipforge.ens-lyon.fr/www/pff/

http://www.frama-c.cea.fr/

http://gappa.gforge.inria.fr/ RR n°7826

Note that verification conditions might be discharged by one or several automated provers.

This research was supported by the ANR projects CerPAN (ANR-05-BLAN-0281-04) and F st (ANR-08-BLAN-0246-01).

C Dependency Graph

In the following graph, the ellipse nodes are Coq files formalizing the wave equation and the convergence of its numerical scheme. The octagon nodes are Coq files that deal with proof obligations generated from the dirichlet.c program file, that is, propagation of round-off errors and error-free execution. Arrows represent dependencies between the Coq files.