IRSTEA Joachim Huet
email: joachim.huet@irstea.fr

C Druilhe

G Debenest

IRSTEA ; C Druilhe

Study of thermal conductivity in organic solid wastes before composting

. Therefore, it fits perfectly with the fourth priority, recovery, which fosters extraction of useful material or energy from wastes. In this context, interest towards composting has increased continuously during the last few years. During treatment, micro-organisms breakdown organic matter and produce carbon dioxide, water and heat. Heat generated by biological activity modifies moisture content and temperature conditions. These changes result in the appearance of a temperature peak in first days of treatment and most pathogens are killed by the high temperatures reached (around 70-80°C), turning waste into a hygienic product.

Among the various physical parameters taking part in the composting process, thermal conductivity seems to be of major importance, and could be used as an indicator to follow heat transfers within the organic matrix. However, as all physical parameters involved in the process, the initial preparation of the substrate (adjustment of moisture content or C/N ratio, addition and mixing with bulking agent ) has an influence on the physical parameters involved in the process (such as bulk density, Free Air Space or air permeability) and thus, on thermal conductivity. Moreover, difficulties often occur in composting experiments because the effects of compaction on physical properties are ignored, or information about these effects is lacking. As soon as the pile of waste is built, the settlement of the composting matrix begins. This settlement, called primary settlement or physical compressive settlement [START_REF] Gourc | Decoupling MSW settlement into mechanical and biochemical processes -Modelling and validation on large-scale setups[END_REF][START_REF] Yue | Pile settlement and volume reduction measurement during forced-aeration static composting[END_REF] is related to the vertical load and leads to compaction. Despite its major importance, until now little has been written on thermal conductivity in composting, in particular about its link with compaction. That s why this study focus on it and aims to investigate how it evolves with (i) compaction or depth within the pile of waste, (ii) preparation parameters of the substrate and (iii) temperature. This investigation was carried out on mixtures of urban sludge and wooden palettes used as bulking agents.

To understand how preparation parameters of the mixtures would affect thermal conductivity, two moisture content (50 and 65%), two types (fresh and recycled) and two meshes of bulking agent (< and > 20 mm) were tested. The influence of compaction (or depth) was evaluated in two steps. First, a Schaub-Szabo device (strongly inspired by the apparatus described in Schaub- [START_REF] Schaub-Szabo | Characterizing the bulk density of compost[END_REF]) was used to get depth-bulk density profiles in the different substrates. Then, these bulk densities were recreated in a modified air pycnometer where thermal conductivity was measured with a thermal probe directly embedded in the composting sample. Therefore, a link between thermal conductivity and compaction (or depth) could be established. On the other hand, the study of the impact of temperature on thermal conductivity was carried out in 10 liter cells where biological activity was prevented by a nitrogen atmosphere. The cells were filled with the same eight sludge-wooden palettes mixtures as before, and a thermal probe was put directly inside each sample. Then, they were put in different constant-temperature baths with target temperature from 5 to 75°C.

In this study, thermal conductivity increased with depth and a statistical analysis highlighted the fact that it was only significantly impacted by moisture content (among the three preparation parameters cited above). Moreover, the impact of temperature on thermal conductivity was clear and a linear relationship between these two parameters could be established. Each correlation was specific to the substrate but with a similar slope. These results are interesting in two ways: first, until now little has been written on thermal conductivity in organic solid wastes, and in particular about its link with compaction. Besides this originality, the data obtained can now be used in numerical modeling to get a more thorough and accurate way to model heat transfers, an essential part of modeling composting systems.

INTRODUCTION

Nowadays, over 800 million tons of wastes are generated each year in France. To control the rise in the amount of waste produced since several decades now, waste policies have been set up in France as in most developed countries to fight this phenomenon. Current waste policy in the EU follows five main priorities, which in a hierarchical order are: prevention, reuse, recycling, recovery and disposal. Prevention encourages the community to reduce the amount of waste generated and to be more efficient in the use of resource. If waste cannot be prevented, as many of the materials as possible should be recovered, preferably by reuse and recycling. Energy and material recovery aims to foster extraction of useful material or energy from wastes. And disposal, as the least desirable option, must be handled with maximum precautions to minimize negative environmental outcomes. As a consequence, composting fits perfectly in the frame of material recovery and that is why during the last few years, interest has grown towards it as a way to treat organic wastes.

Composting can be defined as the biological organic matter decomposition process which occurs in aerobic conditions. At first sight, it could be seen as a quite simple process: with the proper mixture of water, oxygen and nutrients, micro-organisms breakdown organic matter of waste to produce compost. Yet, at micro scales it involves complex processes. These processes are of three types: biological processes and heat and mass transfers. Biological processes are linked to the activity of the decomposing micro-organisms which, provided with oxygen, turn organic matter into compost and produce heat, water and other gaseous compounds as well (mainly carbon dioxide). Their biological activity (consumption of oxygen and production of water and carbon dioxide) leads to mass and volume reduction of the substrate and organic matter stabilization. Furthermore, most pathogens are killed when the peak of temperature (70-80°C) is reached in the first days of treatment.

In return biological activity generates heat, leading to changes in moisture content and temperature conditions. Regarding heat transfers, thermal conductivity is a key role parameter. It is impacted by numerous parameters and among them moisture content, contact between particles and temperature seems to be of great importance. Contact between particles is link to particle size: the thinner the material, the more contacts there is. As moisture content, particle size is a preparation parameter which is set when the substrate is prepared.

As soon as built, the pile of waste is subject to compaction because its physical structure is unable to bear its own weight. One solution to that issue, widely used, is to add some bulking agents to solidify the pile. Nevertheless, in addition to the fact that these bulking agents added should not exceed a certain amount (to keep an acceptable C/N ratio for instance), this primary settlement or physical compressive settlement , as some authors called it [START_REF] Gourc | Decoupling MSW settlement into mechanical and biochemical processes -Modelling and validation on large-scale setups[END_REF][START_REF] Yue | Pile settlement and volume reduction measurement during forced-aeration static composting[END_REF], cannot be avoided and its impact on initial physical characteristics should not be ignored, as it is usually the case because of the lack of information about it.

The original idea of this paper is to focus on thermal conductivity as it seems to be a physical parameter of major importance in regards to heat transfers in composting and yet, not actually well-considered in the existing literature. Therefore, this study aimed to investigate the links between thermal conductivity and (i) compaction or depth within the pile of waste, (ii) preparation parameters of the substrate and (iii) temperature. All experiments were carried out on urban sludge-bulking agent mixtures, with two different moisture contents, two particle sizes of bulking agent and two types of bulking agent (but quite similar anyway).

In order to perform thermal conductivity measurements at different heights (or under different degrees of compaction) we proceed in two steps, using two sets of apparatus. First a Schaub-Szabo device (named after its designer) was used to simulate compaction occurring in the pile and link compaction, depth and bulk density. Then, the sample was placed in another device, called CPP device, where the previous measured bulk densities were recreated. Here, thermal conductivity was measured using a thermal probe directly embedded in the sample. On the other hand, to evaluate the impact of temperature on thermal conductivity, measurements were made with thermal probes placed inside 10 liter cells, filled with substrate and where biological activity was inhibited with a nitrogen atmosphere.

MATERIALS AND METHODS

Sample preparation and characteristics of materials

The substrate tested in these trials was an urban sludge mixed with bulking agents. Two types of bulking agents were used: recycled and fresh wooden palettes. The urban sludge came from a wastewater treatment plant, the recycled wooden palettes used as bulking agents were collected from a composting platform and had already undergone a dozen composting cycles while the fresh ones were directly bought from a retailer.

For each urban sludge bulking agent mixtures, two moisture contents (50% and 65%) and two meshes of bulking agents (< 20 mm and > 20 mm) were tested, leading to a total of eight different mixtures. As on the industrial site, the volumetric ratio of the sludge/bulking agent mixtures was fixed at 1/3 (which corresponded to a dry mass ratio of 0.147).

To minimize changes in the fresh materials, the sludge was sealed in plastic bags, stored frozen at -20°C and thawed as needed in a refrigerated room at 4°C. Likewise, the prepared mixtures were kept at 4°C between two sets of measurements in the Schaub-Szabo and the CPP devices.

Experimental setup

Schaub-Szabo device

The experimental device used to simulate compaction and to determine the variation of bulk density with depth was inspired by the device designed by [START_REF] Schaub-Szabo | Characterizing the bulk density of compost[END_REF]. It consisted of a cylindrical Plexiglas container for the material and a set of weights with a platform scale to apply vertical loads. The container was approximately 700 mm high with an inside diameter of 388 mm; and a fill line was marked 400 mm from the bottom. Inside the container, a bottom grid allowed potential water loss (leachates) from the sample, and the container was perforated at the bottom so that the water can be collected in a recipient placed below.

The recipient containing the sample was weighed; the sample was placed in the container until it reached the fill line, and the recipient was weighed again. The mass of the sample (M 1 ) was determined by difference and the bulk density of the first layer (BD 1 ) calculated according to:

fl V M BD 1 1 =
where V fl is the volume to fill line (0.04729 m 3 ).

To determine the bulk density in the second layer (BD 2 ), known masses were placed on top of the sample to apply a load equivalent to the mass M 1 to the sample in the container. The masses used were steel discs with a diameter of 380 mm (slightly less than the Plexiglas container) and adjusted by adding a known volume of water on top of them. The loading was applied for 24 hours, until the material stabilized to a constant volume.

At the end of the compression, the settlement of the matrix (h 2 ) was recorded. Then, the loading system was removed and fresh material was added to get a 400 mm (h 0 ) height. The mass added (m 2 ) was calculated as follows:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⋅ ⋅ = 0 2 2 2 1 h h V BD m fl
The above procedure was repeated, layer by layer, to simulate six layers of material, which corresponded to a pile of waste of 2.4 m high. A relationship between bulk density and depth for each mixture was finally obtained.

CPP device and thermal conductivity measurement

The bulk densities previously measured in the Schaub-Szabo device were recreated in a CPP device, which was originally used to measure porosity and permeability of compacted materials ( C stands for compaction, P for porosity and P for permeability in the acronym CPP ). It was the same apparatus as described in Druilhe, Benoist et al. [START_REF] Druilhe | Sludge composting: influence of the waste physical preparation on initial free air space, air permeability and specific surface[END_REF] with a thermal probe directly embedded in the composting material. Thermal conductivity measurements were conducted at each bulk density calculated before in the Schaub-Szabo device. Thus, the relationship between thermal conductivity and depth was brought to light.

To recreate the various degrees of compaction, the sample was placed in a removable basket designed to hold about 40-60 L of mixture. Once the basket loaded and the volume known, the bulk density of the sample was modified by adjusting the depth of a perforated compression plate, manually controlled by an airtight screw. Thermal conductivity was measured with an unsteady state probe method which is a commonly used technique on porous materials [START_REF] Chandrakanthi | Thermal conductivity of leaf compost used in biofilters: An experimental and theoretical investigation[END_REF][START_REF] Iwabuchi | Effect of volumetric water content on thermal properties of dairy cattle feces mixed with sawdust[END_REF][START_REF] Van Ginkel | Physical properties of composting material: Gas permeability oxygen diffusion coefficient and thermal conductivity[END_REF]. The thermal probe consists of a heating wire and a thermocouple which measures the temperature at this source. When supplied with a constant electric power Q, the temperature of the probe increases as a function of time. After a certain period of time, the elevation of temperature reaches an asymptotic regime. Therefore, the graph of the rise in probe temperature versus the logarithm of the time gives a straight line, and the thermal conductivity (W.m -1 .K -1 ) can be calculated from the slope A:

A Q ⋅ = π λ 4
Practically, the probe was embedded directly at middle height into the composting material. The electric power was supplied by a constant power generator and the thermocouple of the probe was connected to a data acquisition terminal. Thermal conductivity measurements were directly made with the dedicated software.

Thermal conductivity cells

In parallel, the impact of temperature on thermal conductivity was also investigated. The experiment was carried out in 10 liter cells filled with substrate. Thermal conductivity was measured with a thermal probe placed inside the sample; and the same protocol described in 2.2.2 was applied to the eight sludge-bulking agent mixtures studied. To ensure that temperature was the only changing parameter during the experiment, biological activity had to be inhibited. Indeed, the physical structure of the matrix could be modified under the effect of micro-organisms breaking down organic matter. Therefore, the cells were supplied with nitrogen to inhibit biological activity during the measurements.

The cells were put in different constant-temperature bath with target temperatures ranging from 20 to 75°C, by steps of 5°C. In addition, thermal conductivity measurements were performed right after the samples were taken out of the refrigerated room (and before they reached the ambient temperature) at 5, 10 and 15°C.

As said earlier, the peak of temperature in composting usually reaches temperature values around 70-75°C, that is why the maximum target value was set at 75°C.

Statistical analysis

A statistical analysis was carried out to identify which factors and their interactions had a significant impact on thermal conductivity. The four factors were considered were depth, moisture content and particle size and type of bulking agent (temperature, although studied in this paper, was not used as a factor in the statistical analysis). Four independent variables were used (namely depth, moisture content, particle size and type of bulking agent) and two normalized values, -1 and +1, were considered for each of them. For depth, these two values corresponded to the first (-1) and sixth layers (+1) from -0.4 to 0 m and from -2.4 to -2 m in depth respectively -simulated with the Schaub-Szabo device, for moisture content to 50 (-1) and 65% (+1), for particle size to < 20 mm (-1) and > 20 mm (+1), and for the type of bulking agent to recycled (-1) and fresh (+1) wooden palettes. Therefore, a negative impact of the variable type of bulking agent on one of the physical responses would mean that this response decreased when switching from recycled palettes to fresh ones.

The significance of variable effects and interactions was determined using a Student test at a confidence level of 95% (P < 0.05). All statistical analyses were performed with the software Statgraphics (Centurion XV, Warrentown, Virginia, USA).

Study of thermal conductivity in organic solid wastes before composting J. HUET, C. Druilhe, G. Debenest ORBIT2012 5

RESULTS AND DISCUSSION

Impact of depth, moisture content, particle size and type of bulking agent on thermal conductivity

The results of the statistical analysis are presented in Table 1 below. It displays the impact of a given parameter or interaction on thermal conductivity, if it is whether positive or negative and its significance. As previously stated by some researchers before [START_REF] References Ahn | Determination of thermal properties of composting bulking materials[END_REF][START_REF] Chandrakanthi | Thermal conductivity of leaf compost used in biofilters: An experimental and theoretical investigation[END_REF], thermal conductivity increased with compaction (or increasing bulk density) due to a reduction in the void space. On Figures 1 and 2 below, thermal conductivity seemed to increase with depth as well but the statistical analysis revealed that depth alone didn t have a significant impact and neither in interaction with another factor (see Table 1). Among the three preparation parameters tested (moisture content, particle size and type of bulking agent), the statistical analysis shows that thermal conductivity was significantly influenced by only one factor, moisture content, and increased with it. This result was expected since in theory, thermal conductivity of water is higher than of air (0.06 W.m -1 .K -1 > 0.026 W.m -1 .K -1 ). Therefore, by increasing moisture content, the pores of the organic matrix were filled with water and then, thermal conductivity increased. These results were consistent with the existing literature, even if the experiments were performed on different substrates such as borage seeds [START_REF] Yang | Determination of thermal conductivity, specific heat and thermal diffusivity of borage seeds[END_REF], grain dusts [START_REF] Chang | Thermal-Conductivity and Specific-Heat of Grain Dust[END_REF], beef manure [START_REF] Houkom | Effect of moisture content on thermal diffusivity of beef manure[END_REF], dairy cattle feces mixed with sawdust [START_REF] Iwabuchi | Effect of volumetric water content on thermal properties of dairy cattle feces mixed with sawdust[END_REF], leaf composts [START_REF] Chandrakanthi | Thermal conductivity of leaf compost used in biofilters: An experimental and theoretical investigation[END_REF] or compost-bulking agent materials [START_REF] References Ahn | Determination of thermal properties of composting bulking materials[END_REF].

The type and particle size of bulking agent, although they did not have significant impacts alone, intervened with significant interactions (as shown in Table 1). The first interaction was negative and between moisture content and particle size. It meant that the increase in thermal conductivity with an increase in moisture content was maximal at low particle size. It also meant that at high particle size (> 20 mm), the impact of moisture content on thermal conductivity became negative. Another interpretation of this interaction is that, at high moisture content (65%), thermal conductivity tended to increase when particle size decreased. These results were in agreement with the study of [START_REF] References Ahn | Determination of thermal properties of composting bulking materials[END_REF] who observed an increase in thermal conductivity for eleven different composting materials when grinding them from 10 cm to 0.5 mm. An explanation to this phenomenon was that decreasing particle size created more thermal contacts between particles which, combined with high moisture content, led to an increase in thermal conductivity.

The second significant interaction brought to light was positive and between particle size and the type of bulking agent (see Table 1). The tendency of thermal conductivity to decrease when switching from recycled palettes to fresh ones was significant with particles < 20 mm: 0.083-0.3 W/m/°C with recycled palettes and 0.063-0.158 W/m/°C with fresh ones.

Evolution of thermal conductivity with temperature

As shown on Figure 3, thermal conductivity ( λ ) clearly increased with temperature (T ). Between 5 and 75°C, thermal conductivity increased between two and five times depending on the mixture. Moreover, a linear relationship of the form b T a + ⋅ = λ was obtained, with correlation coefficients ranging from 0.946 to 0.989. b coefficients were all different, meaning that each correlation was specific to the mixture studied. However, they all had a similar slope a (0.003-0.004 W/m/°C 2 ). These are interesting results because, although thermal properties are very important in composting, information on their values for various composting materials is lacking. In this way, the data provided by this study could be very useful, in particularly regarding numerical modeling.

Furthermore, these results were consistent with the little literature currently available. Van Ginkel [START_REF] Van Ginkel | Physical properties of composting material: Gas permeability oxygen diffusion coefficient and thermal conductivity[END_REF] also showed that thermal conductivity of a mixture of chicken manure and wheat straw increased with temperature; and their relationship showed an almost perfect linearity according to the authors. Unfortunately, they did not give the value of the slopes they obtained so we could not compare. Similarly, Yang et al. [START_REF] Yang | Determination of thermal conductivity, specific heat and thermal diffusivity of borage seeds[END_REF] obtained a mathematical equation linking thermal conductivity of borage seeds to temperature and moisture content.

CONCLUSION

In this paper, the influence of three preparation parameters of the substrate (moisture content, particle size and type of bulking agent), depth (or compaction) and temperature on thermal conductivity in organic solid wastes was investigated. Knowledge about thermal conductivity in such substrates and its link with depth (or compaction) is quite poor at present. That is why this paper aimed to give some answers about what parameters influence thermal conductivity in organic solid wastes before composting.

Among the three preparation parameters tested and depth, thermal conductivity was only significantly influenced by moisture content. Contrary to what could have been expected from theory, depth (or compaction) didn t significantly impact thermal conductivity. On the other hand, temperature had a clear impact and a linear relationship between these two parameters was brought to light. Each correlation was specific to the substrate but with a similar slope. However, it is important to note that such results are only valid in the sludge-wooden palettes mixtures tested, and similar studies should be carried out on other substrates.

The data experimentally obtained can now be used in numerical modeling to improve existing heat transfers models, an essential step to get a more thorough and accurate modeling of whole composting systems.

Nomenclature

Thermal conductivity (W.m -1 .K -1 )

A

Numerical coefficient (-) BD i Bulk density of the layer i (kg.m -3 ) FAS Free Air Space (-) h i Settlement of the layer i (m) h 0 Initial level before settlement (m) M1

Sludge-recycled wooden palettes mixture (particle size < 20mm, 50% moisture content) M2

Sludge-recycled wooden palettes mixture (particle size < 20mm, 65% moisture content) M3

Sludge-recycled wooden palettes mixture (particle size > 20mm, 50% moisture content) M4

Sludge-recycled wooden palettes mixture (particle size > 20mm, 65% moisture content) M'1 Sludge-fresh wooden palettes mixture (particle size < 20mm, 50% moisture content) M'2 Sludge-fresh wooden palettes mixture (particle size < 20mm, 65% moisture content) M'3

Sludge-fresh wooden palettes mixture (particle size > 20mm, 50% moisture content) M'4

Sludge-fresh wooden palettes mixture (particle size > 20mm, 65% moisture content) Q Dissipated electrical power per unit length (W.m -1 ) T Temperature (°C) V fl Volume to fill line (m 3 )

Figure 1 :Figure 2 :

 12 Figure 1: Thermal conductivity as a function of depth at different moisture contents and particle sizes (sludgerecycled palettes mixtures)

Figure 3 :

 3 Figure 3: Thermal conductivity as a function of temperature at different moisture contents and particle sizes


Sylvie Boldo
email: sylvie.boldo@inria.fr

Catherine Lelay
email: catherine.lelay@inria.fr

Guillaume Melquiond
email: guillaume.melquiond@inria.fr

Improving Real Analysis in Coq: a User-Friendly Approach to Integrals and Derivatives ⋆

Verification of numerical analysis programs requires dealing with derivatives and integrals. High confidence in this process can be achieved using a formal proof checker, such as Coq. Its standard library provides an axiomatization of real numbers and various lemmas about real analysis, which may be used for this purpose. Unfortunately, its definitions of derivative and integral are unpractical as they are partial functions that demand a proof term. This proof term makes the handling of mathematical formulas cumbersome and does not conform to traditional analysis. Other proof assistants usually do not suffer from this issue; for instance, they may rely on Hilbert's epsilon to get total operators. In this paper, we propose a way to define total operators for derivative and integral without having to extend Coq's standard axiomatization of real numbers. We proved the compatibility of our definitions with the standard library's in order to leverage existing results. We also greatly improved automation for real analysis proofs that use Coq standard definitions. We exercised our approach on lemmas involving iterated partial derivatives and differentiation under the integral sign, that were missing from the formal proof of a numerical program solving the wave equation.

Introduction

From Newton and Leibniz during the 17th century, many mathematicians have used integrals and derivatives. Their use is both for pure analysis theorems, but also more recently for applied mathematics. For example, numerical analysis aims at solving ordinary differential equations and partial differential equations. When the solutions are not analytic, it provides algorithms to approximate these solutions and bounds to assert their correctness. Typically, it consists in a numerical scheme over a discrete grid and its convergence proof, meaning that the approximation improves when the grid size decreases.

Recent advances in formal proof assistants have shown that they can be applied to various kinds of problems, but analysis and especially numerical analysis was not as much studied as algebra. One reason may be that formalizations of analysis were done years ago and seldom used. This is precisely the case in the standard library of Coq: derivatives and integrals were defined with real numbers a dozen years ago, but the libraries did not evolve with Coq. A more extensive use could have proved the ponderousness of the library. The most blatant example is that derivatives require a proof term to be written. This means that, instead of f ′ (x) we have to handle (f, x, H) where H is a proof that f is differentiable in x. This makes the rewritings clumsy and unpractical. More, this is not the way mathematicians prove their theorems: proofs that functions are regular enough, when present, are side-proofs that do not arise in the main development. Another example is the missing definitions and lemmas about partial derivatives.

As shown in Figure 1, we have extended the standard Coq library with equivalent definitions that are easier to use and with some automations. The developed library can be downloaded and browsed at http://coquelicot.saclay.inria.fr/results.html To validate this approach, we have applied it to an example coming from a numerical analysis program. We had previously proved the correctness of a program solving the wave equation, but we left out the proofs about the existence and regularity of a solution to the partial differential equation.

Section 2 presents the existing formalizations of analysis in proof assistants, and especially in Coq. Section 3 presents the power of the underlying logic of Coq, especially the changes when considering the axiomatic of real numbers. Section 4 presents our design choices for the derivative and the integral. Section 5 presents our application about the wave equation and its required results.

We will present the notions of differentiability and integrability as they are defined in various proof assistant such as Coq1 , the Coq constructive library C-CoRN2 , HOL Light3 , Isabelle/HOL4 , PVS 5 , Mizar 6 , and ACL2(r) 7 . We will first present the definitions of the differentiability and integrability predicates, and then describe the design choices for the derivative and integral functions.

Differentiability and Integrability

The choices made here are to adopt one of the common mathematical definitions. For differentiability, Coq and PVS use ε-δ-definition based on the Landau definition directly or through the limit definition:

∃ℓ ∈ R, ∀ε > 0, ∃δ > 0, ∀h ∈ R, (h = 0 ∧ |h| < δ) ⇒ f (x + h) -f (x) h -ℓ < ε
In HOL Light (and Isabelle/HOL which has inherited that formalization of analysis), the Newton's difference quotient is also used, but the limit is a more generic notion as it is defined in a topological space based on a field [START_REF] Harrison | Theorem Proving with the Real Numbers[END_REF]. Another approach is implemented in additional libraries: the Frechet derivative in a real normed vector space [START_REF] Harrison | A HOL theory of Euclidean space[END_REF]. In C-CoRN, a constructive formalization of real numbers for Coq, the previous definition is modified to get a uniform differentiability [START_REF] Cruz-Filipe | Constructive Real Analysis: a Type-Theoretical Formalization and Applications[END_REF] on a closed interval [a; b], i.e. there is a single δ for all x ∈ [a; b] . In Mizar, differentiability is defined for multi-variable functions [START_REF] Raczkowski | Real function differentiability[END_REF] as usual by the existence of a linear function L such that f

(x + h) -f (x) -L(h) = o(h).
To define Riemann integrability, Coq defines the integral on step functions and then uses the traditional ε-δ-definition: ∀ε > 0, there are two step functions ϕ, ψ :

[a; b] → R, such that (∀t ∈ [a; b], |f (t) -ϕ(t)| ≤ ψ(t)) ∧ | ψ| < ε
The value of the integral is then defined as the limit of ϕ when ε → 0. PVS in [START_REF] Butler | Formalization of the integral calculus in the PVS theorem prover[END_REF] and Mizar in [START_REF] Endou | The definition of the Riemann definite integral and some related lemmas[END_REF] define Riemann integrability as the convergence of Riemann sums. Both definitions are mathematically equivalent. A difference is that the integral value is explicitly given by this second definition while Coq's definition only provides approximations of Riemann integral.

C-CoRN uses a third equivalent definition of Riemann integrability: the convergence of Darboux sequences. As with Riemann sums, the integral value is directly obtained from the definition. HOL Light does not define Riemann integrals, but both gauge and Lebesgue integrals, which are more general but less intuitive. As with derivatives, all the notions are defined for multivariate functions.

In a former library about reals in Isabelle/HOL [START_REF] Fleuriot | On the mechanization of real analysis in Isabelle/HOL[END_REF], and the current library in ACL2(r) [START_REF] Gamboa | Continuity and differentiability in ACL2[END_REF], differentiability was defined from non-standard analysis where the formal notion of "infinitely close", i.e. difference is less than all positive real numbers, replaces the usual notion of "arbitrarily close" corresponding to a common formula stating that ∀ε > 0, . . . , the difference may be made smaller than ε.

Derivative and Integral

In pen-and-paper mathematics, we understand f ′ (x) as the function f is differentiable at x and f ′ (x) is its derivative value. But the corresponding definitions and their uses in formal proof assistants are not straightforward and differ, among others as the underlying logic is different.

PVS is the closest to the mathematical point of view by using Type Correctness Conditions (TCCs). A user may write a statement containing occurrences of f ′ (x) without justification and PVS generates additional goals to prove that f is actually differentiable in x. PVS tries to infer automatically these additional goals from the context. The same approach is used in ACL2(r) [START_REF] Gamboa | Non-standard analysis in ACL2[END_REF].

HOL Light uses its Hilbert's epsilon applied to the differentiability property. The user can then write f ′ (x) without proof, but has to prove that f is differentiable at x before using differentiation rules.

In Coq, the derivative is a partial function taking explicitly a proof term of differentiability. As a derivative cannot be defined or used without this proof term, this is cumbersome to use. One of the goals of the MathClasses project [START_REF] Spitters | Type classes for mathematics in type theory[END_REF] is to provide assistance work with these proof terms by trying to infer them.

Logical Foundations

Rather than introducing our own set of axioms, we have based our formalization on the axioms from the Coq standard library on real numbers. This ensures that our automated proof tools are usable by works based on the standard library, that some results from the standard library could be reused, and that our work is consistent (assuming the standard library is). We have also taken great care to not introduce any other axioms, especially not the excluded middle, which is neither a native concept in Coq logic, nor a consequence of Coq axioms on real numbers.

An Overview of Coq's Logic

The formal system of Coq is an intuitionistic logic called the Calculus of Inductive Constructions. A salient point is that, whenever one wants to prove a property of the form ∃x, P (x), one has to actually build a witness x such that P (x) holds. This is different from classical logic where one could have simply proved ¬∀x, ¬P (x) and be done.

Another peculiarity of Coq's logic is related to its type hierarchy. Logical formulas have type Prop while values and functions have their types in Type. The point of interest is that Prop is non-informative, that is, one can only use the witness of an existential property ∃x, P (x) inside the proof of logical formulas. A witness can never be used inside a logical formula itself, or more importantly to define a value or a function. For instance, being able to prove a formula ∀x : X, ∃y : Y, P (x, y) does not provide a function f : X → Y such that ∀x, P (x, f (x)). Note that strengthening the property so that y exists and is unique does not help either.

The traditional way in Coq to circumvent this issue is in the use of specification types. They are existential types denoted {x : X | P (x)} and they contain dependent pairs (x, p) such that x is an element of type X and p is a proof of the logical formula P (x). The upside is that witnesses are readily available for use inside values and functions. The downside is that these types do not live in Prop and are therefore less natural to manipulate.

For instance, consider the predicate derivable pt from the standard library that states that a function f is differentiable at point x. It is in fact a notation for the specification type that associates a value ℓ with the proof that ℓ is the limit of the slope function of f at point x. As a consequence, knowing that a function is differentiable gives access to the value of its derivative. But it also means that trying to express that both f and g are functions that are differentiable as derivable pt(f, x) ∧ derivable pt(g, x) will be rejected by Coq. Indeed, this formula is ill-typed since neither members are logical formulas. A coercion from specification types to formulas in Prop would avoid this issue as long as one only needs the differentiability property and not the derivative of a function.

In the same way that one can extract a witness from the formula ∃x, P (x) only when performing a proof, the information about the disjunction P ∨ Q cannot be used to make a choice inside a value or a function. Again, types outside Prop have been introduced to offer this possibility; they are denoted {P } + {Q} in Coq and their values can be used to select the branch of an if-then-else.

To end this overview, one should mention that subsets of a type T are usually represented by predicates, that is, functions of type T → Prop. As a consequence, we will indifferently note the containment property by the logical formulas x ∈ S or S(x).

Coq's Standard Axioms for Real Numbers

The formalization of real numbers from the standard library is axiomatic rather than definitional. Instead of building reals as Cauchy sequences or Dedekind cuts of rational numbers and proving their properties, Coq developers have chosen to assume the existence of a set with the usual properties of the real line. In other words, the standard library states that there is a set R, some arithmetic operators -, +, ×, -1 , and a comparison operator <, that have the properties of an ordered field. Except perhaps for the choice for the domain of -1 , the previous axioms are not controversial.

Below are the three axioms that state that R is Archimedean, closed under the supremum bound, and its order is decidable:

archimed: There is a function up : R → Z such that ∀x ∈ R, x < up(x) ≤

x + 1. completeness: As long as one can prove that a subset E of R is not empty (∃x, E(x)) and is bounded (∃M, ∀x, E(x) ⇒ x ≤ M ), one gets a value of type {y : R | y is an upper bound of E and it is the least one}. total order T: Given two real numbers x and y, there is a value of type {x < y}

+ {x = y} + {x > y}.
While all three axioms could have been defined as logical formulas in Prop, they were not. In other words, it is equivalent to having three functions that can compute the integer part of a real number, the supremum of a bounded subset of R, and the order of two numbers. Notice that excluded-middle is not one of the axioms of Coq's standard real numbers. While the standard library sometimes imports this axiom, the CoqTail project8 has shown that it was often unneeded. So we restrict our use of classical reasoning to the proof of goals that are the negation of a logical formula.

Notice also that the completeness axiom has an intuitionistic rather than classical feel to it. Indeed, the axiom is written in such a way that it cannot produce a proof of E(x) for some real x (especially not the real that is the supremum of E), since properties of the form E(x) always appear as premises in the axiom. At best, one can derive the following property which protects the existential quantifier behind a double negation:

∀ε > 0, ¬¬∃x, y -ε ≤ x ≤ y ∧ E(x).

Limited Principle of Omniscience

Let P be a decidable predicate on natural numbers. The limited principle of omniscience (LPO) states that one can decide whether the property never holds, and if it does, return n such that P (n) holds. In Coq syntax, the principle is stated {n | P (n)} + {∀n, ¬P (n)}. It cannot be derived without axioms in Coq, as it would require the ability to test all the values of n at once. Thanks to the axioms on real numbers, it becomes possible.

The way we have proved the LPO is as follows. Since P is decidable, we can build a function f (n) that returns 1/(n+1) if P (n) holds and 0 otherwise. Let us consider the subset of real numbers {f (n) | n ∈ N}. It is nonempty and bounded by 1, thus it has a supremum (completeness). This supremum can be tested against 0 (total order T). If it is zero, we deduce ∀n, ¬P (n). Otherwise we compute its discrete inverse (archimed) which will act as a witness for building a value of type {n | P (n)}. This proof of the LPO, inspired by the CoqTail project, is an improvement over some previous work that depended on the not all ex not consequence of the excluded-middle axiom [START_REF] Kaliszyk | Computing with classical real numbers[END_REF].

Bounds and Limits

Now that we have proved this principle, we can use it to decide whether a subset E of real numbers is bounded, which is a precondition for computing its supremum. First, let us consider the family of subsets E n = {0}∪(E ∩ (-∞; n]). They are nonempty and bounded, so they have a supremum s n . As a consequence, deciding whether E is bounded and computing an upper bound is a matter of applying the LPO to decide the following alternative (Lemma Rbar ub dec of our development):

{M | ∀n, s n ≤ M } + {∀M, ¬(∀n, s n ≤ M )}.
The proof requires that ∀n, s n ≤ M is decidable, which is just a consequence of the LPO applied to the decidable predicate n → M < s n .

Let us define the complete lattice Rbar = R ∪ {-∞, +∞}, inspired by [START_REF] Hölzl | Three chapters of measure theory in Isabelle/HOL[END_REF]. Since we are able to decide whether a set is bounded, we can define supremum and infimum functions for nonempty subset of Rbar. 9 Let us consider a sequence (u n ) n∈N of elements of R (or Rbar). The set of its values is nonempty, thus we can define its superior and inferior limits:

lim sup (u n ) n∈N = inf m∈N (sup n∈N u m+n ) lim inf (u n ) n∈N = sup m∈N ( inf n∈N u m+n )
At this point, we arbitrarily define a function Lim seq from sequences to real numbers as (lim sup(u n ) + lim inf(u n ))/2, or 0 in case of infinities. If (u n ) n is a converging sequence, Lim seq (u n ) n∈N is the actual limit of the sequence, since inferior, superior, and plain limits are then equal. Note that this also gives us a way to decide whether an arbitrary sequence is converging: we just have to compare its inferior and superior limits. With this definition, we also get "limits value" for non-convergent sequences: Lim seq(-1) n = (1 + (-1))/2 = 0.

There were simpler possibilities for defining Lim seq, but this one offers the equality Lim seq(α•u n ) = α•Lim seq(u n ) for any real α without requiring (u n ) n to converge.

We now have an operator able to compute the limit of any converging sequence and otherwise return an undefined value (as far as the user is concerned). We can similarly define the limit of a function f in a point x by

Lim(f, x) = Lim_seq f x + 1 n n∈N .
Again, if the function has a limit at this point, in sense of the usual εδ definition for pointed limit, operator Lim returns it, since lim u→x f (u) is then equal to lim n→∞ f (u n ) for any sequence such that lim(u n ) = x and ∀n ∈ N, u n = x. Otherwise it does return a value, but a meaningless one. We did not gain anything significant by defining the limit with more complicated sequences such as (f (x -1/n) + f (x + 1/n))/2.

Compactness

Limits are the basis for doing real analysis. Another important tool is the property of compactness, which has numerous applications in traditional mathematics. For instance, a function continuous on a compact set is uniformly continuous. Unfortunately, the compactness property is inherently classical, up to the point that constructive mathematics tend to continuity so that it actually means uniform continuity in order to avoid compactness [START_REF] Cruz-Filipe | C-CoRN : the constructive Coq repository at Nijmegen[END_REF]. Our goal is to stay as close as possible to traditional analysis, so dropping compactness is not a solution.

One of the definitions of a compact set is a set such that, from any cover with open sets, one can extract a finite subcover. Yet in most of the proofs we are interested in, we do not need the whole power of this property. Indeed, the extracted sets are useless, only their minimum diameter matters. Moreover, the finiteness property is only useful so that this minimum is nonzero. As a consequence, we can substitute to the traditional definition a property of interval [a, b] related to Cousin covers and gauge functions:

∀δ : R → R + , {δ ′ : R + | ∀x ∈ [a, b], ¬¬∃t ∈ [a, b], |x -t| < δ(t) ∧ δ ′ ≤ δ(t)}.
The proof of this lemma (compactness value) requires us to construct a value δ ′ that satisfies the property. It is defined by the formula

δ ′ = 1 2 sup {d | d ≤ 1 ∧ ∀x ∈ [a, b], ∃t ∈ [a, b], |x -t| < δ(t) ∧ d ≤ δ(t)}.
First, we prove that this supremum is not equal to zero. This is the same idea as with most cover-based proofs: consider the widest interval [a, b ′ ) such that the

property ∃d > 0, ∀x ∈ [a, b ′ ), ∃t ∈ [a, b], |x -t| < δ(t) ∧ d ≤ δ(t) holds.
The bound b ′ is defined as a supremum value and we prove that, if b ′ < b, there is a contradiction. Then we prove that it still holds for the whole interval [a, b], so δ ′ is positive. Second, we prove that δ ′ satisfies the original property. This is immediate, but only because the property contains a double negation. Indeed, as mentioned before, the completeness axiom never provides us with proofs that the property holds for values smaller than the supremum. So the best we can prove is that there would be contraction if δ ′ did not satisfy the property. In practice, this double negation does not matter because we always use the compactness property to exhibit contradictions.

Finally, we have proved the compactness property not just for segments of the real line but for any n-orthotope (Cartesian product of n segments).

Derivative

From the previous definition of a limit function in Section 3.4, we defined a total derivative function:

Derive(f, x) = lim h→0 h =0 f (x + h) -f (x) h
and we proved that, if the derivative number f ′ (x) exists, it is equal to our derivative number using lemmas proved in the standard library. This definition allows to write properties about derivative without any proof term. For example to write "the function f is differentiable on domain D and f ′ satisfies the property P on this domain" in Coq, we previously had to write Lemma pr : forall x, D x → derivable_pt f x. Goal forall x (Dx : D x), P (derive_pt f x (pr x Dx)).

using a dependent pair. With our approach, we can write the same statement as Goal forall x, D x → ex_derive f x. Goal forall x, D x → P (Derive f x). and then we can prove separately differentiability and the property P . Moreover, as our derivative is a limit, some properties such as (α • f ) ′ = α • f ′ do not require the differentiability of f . Therefore, the proof burden is lightened, as theorems have less preconditions. Such a property could also be obtained for a derivative defined on top of Hilbert's epsilon, but not without introducing a cumbersome definition based on conditionals. Our approach also makes it possible to express more easily the n-th derivative:

Derive_n (f, 0, x) = f (x) Derive_n (f, n + 1, x) = Derive(x → Derive_n (f, n, x), x)
and we express the n-th differentiability as

forall i : nat, (i < n)%nat → ex_Derive (Derive_n f i) x

Riemann integral and Riemann sums

We define our integral RInt as the limit of

RInt_val (f, a, b, n) = b -a n + 1 n k=0 f x n k + x n k+1 2 where x n k = a + k • (b -a)/(n +

1). As above, this definition allows to prove

∀α ∈ R, RInt(α • f, a, b) = α • RInt(f, a, b) without hypotheses on f .
The Riemann integrability in the standard library is based on step functions. Unfortunately, step functions from the standard library are hard to use: a step function is built from a function f and two lists lx and ly that must satisfy five conditions. These conditions are difficult both to prove and to use. For example, the last one states that ∀i, ∀x ∈ (lx i ; lx i+1 ), f (x) = ly i . This is impractical as it does not provide any information about the values f (lx i ).

We chose to define new step functions based on Ssreflect sequences [START_REF] Gonthier | A small scale reflection extension for the Coq system[END_REF]:

Record SF_seq {T : Type} := mkSF_seq {SF_h : R ; SF_t : seq (R * T)}.

Using Ssreflect libraries, our step functions were easier to define and use. We where RInt_part a b n is the partition used in this proof and SF_seq_f2 builds the needed step function. Moreover, as the standard library does not provide a global relation between Riemann integrability and Riemann integral, we chose to prove the equivalence between this definition and the following convergence of Riemann sums:

∃If ∈ R, ∀ε > 0, ∃δ > 0, ∀(σ, ξ), max 0≤k≤n |σ k+1 -σ k | < δ ⇒ |S(f, σ, ξ) -If | < ε where S(f, σ, ξ) = n k=0 f (ξ k ) • (σ k+1 -σ k
) is a Riemann sum, n is the length of ξ, n + 1 the length of σ, σ 0 = a, σ n+1 = b and (σ, ξ) is a pointed subdivision, i.e. ∀k ≤ n, σ k ≤ ξ k ≤ σ k+1 . This is the same approach as in PVS and Mizar.

Thanks to this new definition, we can express integrability with the same structure as convergence and differentiability: a function gives the value of the Riemann integral, a predicate states the property of Riemann integrability, and is_RInt (f, a, b, If ) states that If is the Riemann integral of f between a and b.

We can note that our value is the limit of a sequence of specific Riemann sums S(f, σ n , ξ n ) such that ∀i ∈ [[0; n]], σ i+1 -σ i = (b -a)/(n + 1) . The correctness is then immediate for our new definition of Riemann integrability.

To ensure the compatibility with the standard library of Coq, we proved the equivalence between both definitions, so that we can take advantage of existing lemmas.

Application

Case study and d'Alembert's Formula

Our main application is part of a project aiming at proving numerical analysis programs. The case study was a C program that implements a numerical scheme for the resolution of the one-dimensional acoustic wave equation. This corresponds to the oscillation of an attached rope where c is the constant propagation velocity, which depends on the section and density of the string. More precisely, we consider the following initial-boundary value problem: we have the initial values u 0 and u 1 , a source term s and we want to compute an approximation of the exact solution u of

∀t ≥ 0, ∀x ∈ [x min , x max ], ∂ 2 u ∂t 2 (x, t) -c 2 ∂ 2 u ∂x 2 (x, t) = f (x, t), ∀x ∈ [x min , x max ], ∂u ∂t (x, 0) = u 1 (x), ∀x ∈ [x min , x max ], u(x, 0) = u 0 (x), ∀t ≥ 0, u(x min , t) = u(x max , t) = 0
To actually compute an approximation of the solution u, we chose the second order centered finite difference scheme, also known as three-point scheme. The size of the grid is (∆x, ∆t) and the value u k j ≈ u(j∆x, is given by

u k j -2u k-1 j + u k-2 j ∆t 2 -c 2 u k-1 j+1 -2u k-1 j + u k-1 j-1 ∆x 2 = s k-1 j
and similar formulas that depend on u 0 and u 1 for the initializations u 0 j and u 1 j . We proved that rounding errors do not endanger the results of the numerical scheme [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF]. We also formalized in Coq the numerical scheme and proved its convergence, for an infinite rope [START_REF] Boldo | Formal proof of a wave equation resolution scheme: the method error[END_REF]. In that work, the differentiation operator was an arbitrary function. We did not define it nor did we assume any of its properties. We only used the fact that u is a solution of the partial differential equation expressed using this operator. This fully corresponds to the way mathematical proofs are done: we put f ′ (x) or ∂ 2 u ∂t 2 (x, t) and not diff(f, x, H) where H is a proof that f is derivable in x or diff(z → diff(y → u(x, y), z, H 1 ), t, H 2 ) where H 1 and H 2 are adequate proof terms. We also needed the regularity of this solution: it is supposed to be near its Taylor expansion with the usual mathematical bounds (see below).

Later, we proved the full C program for a finite rope [START_REF] Boldo | Wave Equation Numerical Resolution: a Comprehensive Mechanized Proof of a C Program[END_REF]. As we needed to precisely specify what the program was supposed to compute, we defined each derivative as the limit of f (x+h)-f (x) h when h goes to zero. We used the Frama-C platform with the Jessie plugin and the specification of the C program is described in C comments called annotations. As the language of these annotations is first-order logic, we could not define in our specifications a differentiation operator, but had to define each of the four derivatives as a limit (with a ∀ε, ∃δ . . . formula). This fully specifies the derivatives but was impractical and difficult to read. Yet, this meant an equivalence between our previous work and real derivatives. But as our previous work required a differentiation operator, we had to actually provide it. The chosen solution was to define it as a parameter and add an axiom stating that, if the function is differentiable, the result of this operator is the expected derivative. This axiom, similar to a Hilbert ε operator, was not satisfactory. Thanks to the formalization presented in this paper, we can now create this operator as a function and get rid of that axiom.

The other axioms needed by this development are the following ones: the existence of a solution to the partial differential equation and its regularity.

About the existence, the mathematical proof is simpler than expected as this equation has an analytical solution. More precisely, the following d'Alembert's formula

u(x, t) = 1 2 (u 0 (x + ct) + u 0 (x -ct)) α(x,t) + 1 2c x+ct x-ct u 1 (ξ) dξ β(x,t) + 1 2c t 0 x+c(t-τ ) x-c(t-τ ) f (ξ, τ ) dξ dτ γ(x,t)
defines a function that is solution to the previous partial differential equation. We define α, β, and γ, as parts of this formula that will be used below. Just note that they are of increasing difficulty to handle and derive.

Taylor expansions

The regularity of the solution u is the base of the convergence of the chosen numerical scheme. In the scheme statement, since the grid sizes are small, we can recognize discrete derivatives:

u k j -u k-1 j ∆t ≈ ∂u ∂t (j∆x, k∆t) u k j -2u k-1 j + u k-2 j ∆t 2 ≈ ∂ 2 u ∂t 2 (j∆x, k∆t)
and similarly for space derivatives. The discrete equation is the exact discrete analog of the continuous wave equation.

Our definition of the Taylor polynomial is the usual one:

TP n (f, x, t) = n p=0 1 p! p m=0 p m • ∂ p f ∂x m ∂t p-m (x, t) • ∆x m • ∆t p-m
For the main iteration, we need to guarantee that the difference between the function and its order-4 Taylor polynomial is proportional to ( √ ∆x 2 + ∆t 2 ) 4 . For the initializations, we also need this property at level 3.

We first proved the common Taylor-Lagrange theorem and extended it to its two-dimensional version we really needed here. We had to prove the Schwarz theorem to be able to switch derivatives in space and time. Note that the hypotheses required to make this switch possible are strong (existence and continuity of both second-order derivatives).

Automation

As explained, the two unproved properties from the original development are that the solution exists and is sufficiently regular. Existence has already been tackled thanks to a reflexive Coq tactic for proving differentiability [START_REF] Lelay | Différentiabilité et intégrabilité en Coq. Application à la formule de d'Alembert[END_REF]. The tactic was limited: it could handle expressions with one variable only. As a consequence, while it could automatically perform differentiability proofs on the α and β parts of d'Alembert's formula, human intervention was needed for differentiating under the integral sign of γ. The reason for this shortcoming was the need for proof terms in derivatives and integrals. For instance, the term b(x) a(x) ∂f ∂x (x, t)dt contains a proof that ∂f ∂x (x, t) can be integrated for t between a(x) and b(x), while the term ∂f ∂x (x, t) itself contains a proof that f has a first derivative at any point (x, t) of a domain that depends on x. This nesting of values and proof terms ended up being out of reach of our tactic.

For the existence, we only had to consider four partial derivatives of γ. So, despite the absence of automation, we succeeded in formally proving it in Coq. For this work, however, we wanted to prove not only the existence but also the regularity, which means manipulating tens of partial derivatives of γ. This makes it out of reach of a non-automated proof. So we have improved the original tactic now that we have got rid of proof terms in values. The new tactic auto derive still produces side conditions, but the values no longer depend on their proofs, which means the tactic is now able to differentiate below the integral sign. As the former one, this tactic is programmed in Ltac.

The tactic is meant to help proving statements of the form Operator D is then applied to this inductive object. This is a recursive function that differentiates an expression and generates side conditions. For instance, given an object representing x → 2 • f (x), it produces an object representing x → 2 • f ′ (x) and a side condition that f can be differentiated at the considered point. Lemma D correct states that the generated object is the actual derivative when the side conditions hold. The tactic simply applies this lemma to the current goal, thus solving it, assuming the user can prove the side conditions.

Consider the following script that proves that ∂ 2 α ∂x 2 exists and is equal to α 20 .

Definition alpha x t := 1/2 * (u0 (x + c*t) + u0 (x -c*t)). Definition alpha20 x t := 1/2 * (Derive_n u0 2 (x + c*t) + Derive_n u0 2 (x -c*t)).

Lemma alpha_20_lim : forall x t, is_derive_n (fun u ⇒ alpha u t) 2 x (alpha20 x t).

Proof. intros x t. unfold alpha. auto_derive_2.

The auto derive 2 is just an ad-hoc tactic developed for this example. It simply applies auto derive twice in a row, so as to prove properties on iterated derivatives. After executing the tactic, the user is left with several goals to prove. They state that function u 0 can be differentiated with respect to the first variable at points (v + c • t) and (v -c • t) for an arbitrary real v around x. They also state that the first partial derivative of u 0 can be differentiated with respect to the first variable at points (x + c • t) and (x -c • t). Finally, the last goal the user has to prove is the equality between α 20 (x, t) and the expression obtained by automatic differentiation, which is straightforward.

Conclusion

We have presented a Coq development for real analysis that aims at being closer to the traditional way of handling analysis theorems in pen-and-paper proofs. The main idea we have followed is to replace all the partial operators with total operators, so that the user no longer has to manipulate dependent types. Once all the theorems and especially rewriting rules have non-dependent hypotheses, they become much easier to apply, since reasoning is back to being backward: from the goal to the hypotheses.

The standard libraries of HOL Light and Isabelle/HOL also provide such total operators for derivatives and integrals. The main difference with our work is that they are defined thanks to Hilbert ε operator, which is not available in Coq. Instead, we have provided algorithms for these operators. They are not the kind of algorithms that one would find in traditional constructive mathematics, since our model of computation is a bit unusual. It has a decidable order, as in the Real RAM model, but it also has a supremum operator and an integer part.

To obtain this model, we have not added any axiom, we have just reused the axiomatization of real numbers from the Coq standard library. We have also taken great care to never use the axiom of excluded middle, contrarily to what is done in the standard library. Unfortunately, when looking at the assumptions of some of the theorems of our library, there might be occurrences of this axiom. They leak from the standard library through the equivalence lemmas between our definition of integral and the standard one. Work is under way to remove these uses from the standard library and therefore get a formalization that no longer relies on excluded middle.

Our development of more than 500 lemmas provides total operators for limits, derivatives, and integrals, and equivalence lemmas between our constructive definitions and the partial operators from the standard library. We have also extended the theory of real analysis further than what is available in the standard library: iterated partial derivatives, parametric integrals, and so on.

We have applied our formalization to filling the holes in the formal verification of a numerical program: the three-point scheme for solving the one-dimensional wave equation. The original formalization set as an axiom the existence of total operators for partial derivatives; the verification would never have succeeded if it had had to cope with the dependent types needed for expressing fourth derivatives. The work presented in this paper fills this gap.

The original formalization was also assuming the existence and the regularity of a solution to the partial differential equations. When starting from d'Alembert formula, the formal proof of these properties is mostly mechanical. One just has to differentiate the formula as many times as needed (up to 20 times for order-4 regularity). For that purpose, we have also developed a reflexive Coq tactic that is able to perform such repetitive tasks. The strength of our tactic is that it is able to differentiate under the integral sign. As far as we know, no similar strategy has been developed for other provers.

Future Works

For this work, we chose Riemann's definition of integral. Our only motivation for this choice was that it is the integral provided by Coq's standard library and we wanted to check that we were not less expressive than the standard library. If not for this constraint, we would have chosen a different definition, e.g. Lebesgue integral. Indeed, compared to other definitions of integral, Riemann integral does not have much positive points, except for its prestigious name. Lebesgue integral would have been easier for us to define and to manipulate, since it is almost a simple supremum.

In fact, we could presumably go further than Lebesgue integral and define gauge integral. Indeed, our work on compactness has shown that our limited framework (no general excluded middle) was still sufficient to manipulate gauge functions and extract finite subcovers. This is the main property needed to prove Cousin's theorem and therefore to define this integral. It has two main positive points. First, it is no more complicated than Riemann integral. Second, the second fundamental theorem of calculus can now be expressed without any precondition: each differentiable function is the integral of its derivative. This makes its usage in formal proofs straightforward.

For now, we have only defined limits, derivatives, and integrals, as total operators. Our goal is to extend this paradigm to other common operators, e.g. power series and reciprocal, so as to provide all the basic blocks of analysis. We also intend to extend our automated tools beyond just differentiation. An obvious extension is integration, but also automatic proofs of integrability and continuity. Indeed, differentiating under the integral sign tends to generate numerous side conditions about these properties and they would greatly benefit from being automatically discharged by the prover during the symbolic computations.

Fig. 1 .

 1 Fig.1. Sketch of our library: equivalent definitions, additional lemmas, and an automatic tactic, without introducing new axioms.

  define our step functions with a generic type T to use it in the same way with T := R or T := Rbar. For example, to define the step function needed for RInt_val (f, a, b, n), we use Definition SF_val_seq (f : R → R) (a b : R) (n : nat) : SF_seq := SF_seq_f2 (fun x y ⇒ f ((x+y)/2)) (RInt_part a b n) 0.

derivable_pt_lim f x l

  that is, f has a derivative at point x and it is equal to l. A variant of the tactic is able to tackle goal where l is not yet known. The tactic first performs a reification of the function f into an inductive object describing the expression. Variables are encoded using De Bruijn's indexes. For instance, the inductive object for (* addition *) (AppExt 1 g [:: Var 1]) (* application: g(y) *) (Var 0)) (* z *) (Cst 0) (* lower bound of integral: 0 *) (Binary Emult (Cst 2) (Var 0))) (* upper bound: 2y *)

http://coq.inria.fr/stdlib/index.html

http://corn.cs.ru.nl/

http://www.cl.cam.ac.uk/ ~jrh13/hol-light/

http://isabelle.in.tum.de/dist/library/HOL/index.html

http://shemesh.larc.nasa.gov/fm/ftp/larc/PVS-library/

http://www.mizar.org/

http://www.cs.utexas.edu/users/moore/acl2/v5-0/

http://coqtail.sourceforge.net/

Supremum and infimum of a subset P are extended with sup ∅ = -∞ and inf ∅ = +∞ , when one can prove P = ∅ ⇒ ∃x, P (x). This generalization is not needed here.

⋆ This research was supported by the F st project (ANR-08-BLAN-0246-01) of the French national research organization (ANR) and by the Coquelicot project of the Digiteo cluster and the Île-de-France regional council.