A Verified ODE Solver and the Lorenz Attractor

A Verified ODE Solver and the Lorenz Attractor J Autom Reasoning (2018) 61:73–111 https://doi.org/10.1007/s10817-017-9448-y A Verified ODE Solver and the Lorenz Attractor Fabian Immler Received: 19 April 2017 / Accepted: 27 December 2017 / Published online: 22 January 2018 © The Author(s) 2018. This article is an open access publication Abstract A rigorous numerical algorithm, formally verified with Isabelle/HOL, is used to certify the computations that Tucker used to prove chaos for the Lorenz attractor. The ver- ification is based on a formalization of a diverse variety of mathematics and algorithms. Formalized mathematics include ordinary differential equations and Poincaré maps. Algo- rithms include low level approximation schemes based on Runge–Kutta methods and affine arithmetic. On a high level, reachability analysis is guided by static hybridization and adap- tive step-size control and splitting. The algorithms are systematically refined towards an implementation that can be executed on Tucker’s original input data. Keywords Isabelle/HOL · Ordinary differential equation · Rigorous numerics · Poincaré map · Lorenz attractor 1 Introduction Computer assisted proofs, i.e., mathematical proofs that rely on the output of a computer program, depend crucially on the correctness of the program. Important computer assisted proofs for e.g., the Kepler conjecture or the Four Color Theorem, have therefore been formally verified. In this article, we consider the Lorenz attractor—perhaps one of the most prominent examples of deterministic chaos—and its computer assisted proof by Warwick Tucker. The proof relies on a rigorous numerical ordinary differential equation (ODE) solver. In this article, we describe the long-term project of formally verifying (in Isabelle/HOL [39]) an ODE solver that is capable of certifying Tucker’s computations. B Fabian Immler immler@in.tum.de Institut für Informatik,Technische Universität München, Munich, Germany 123 74 F. Immler 1.1 History In 1963, meteorologist Edward Lorenz [31] introduced a system of ODEs as a simplified model for atmospheric dynamics. He observed that even the smallest perturbation in initial values would lead to completely different long-term behavior of the system. Referring to the original motivation, he asked: “Does the Flap of a Butterfly’s Wings in Brazil Set Off a Tornado in Texas?” and the term butterfly effect entered popular culture. The Lorenz system tends to evolve to a complicated structure (Fig. 2), which became an iconic example of deterministic chaos: According to Sparrow [45] “the number of man, woman, and computer hours spent on [the Lorenz equations …] must be truly immense”. Despite its popularity and the amount of effort put into its study, nobody managed to prove that the Lorenz attractor is chaotic in a rigorous mathematical sense. The problem of rigorously proving chaos in the Lorenz attractor even made it into a list of 18 important problems for the twenty-first century that Field’s medalist Stephen Smale composed in 1998 [43]. Shortly after, Warwick Tucker managed to give an affirmative answer by presenting a computer-assisted proof [47,48]. Tucker’s programs were written in C++ and are not for- mally verified. Tucker even discovered (and fixed) some bugs in it [46,49]. Formal verification of the numerical results needed for the proof is therefore a worthwhile goal. 1.2 The Lorenz Attractor We start with describing the Lorenz attractor and some of the properties that were conjectured from numerical simulations. In his proof, Tucker considers the following three dimensional ODE for fixed parameters k ,λ : 1,2,3 1,2,3 x˙ = λ x − k (x + y)z 1 1 y ˙ = λ y + k (x + y)z 2 1 z ˙ = λ z + (x + y)(k x + k y) 3 2 3 As intuition, an ODE describes the velocity vector (x˙ , y ˙, z ˙) in which a particle at a point (x , y, z) moves. The evolution of a particle subject to the ODE is described by the so-called flow φ. A particle x ∈ R will be at position φ(x , t ) after time t ∈ R. 0 0 Figure 1 depicts the numerical simulation of the evolution of a particle starting at (0.1, 0, 0): It moves to right (x ≈ 15) and up (z ≈ 50) at time t ≈ 0.5, then down to about z = 27 and oscillates with increasingly larger amplitude around z = 27. Figure 2 depicts the trace of a long-term evolution in the three dimensional phase space, it indicates Property 1: Property 1 Solutions remain in a bounded region of the phase space. Particles that approach the origin (0, 0, 0) from above exhibit a very sensitive dependence on initial conditions: a slight perturbation can make the particle flow to either the left or right branch of the Lorenz attractor, which we call Property 2: Property 2 Solutions exhibit sensitive dependence on initial conditions. This dependence is such that arbitrarily small initial sets will eventually spread over the whole attractor. 1 8 That is, Lorenz’ original equations for the classical parameters β = ,σ = 10,ρ = 28 in Jordan normal σ σ −1+τ σ −1−τ −σ −1+τ form using τ:= (σ + 1) + 4σ(ρ − 1), k := , k := , k := , λ := , 1 2 3 1 τ 2σ 2σ 2 −σ −1−τ λ := ,and λ :=−β. 2 3 123 A Verified ODE Solver and the Lorenz Attractor 75 Fig. 1 Temporal evolution of t → φ((0.1, 0, 0), t ) Fig. 2 Simulation of a part of the Lorenz attractor (φ) and Poincaré section (Σ) 1.3 Tucker’s Proof How does Tucker go about proving those properties? First of all, he uses a standard technique: he introduces a so-called Poincaré section. This is a distinguished set in the phase space, in this case a square on the plane z = 27, namely Σ =[−6, 6]×[6, 6]×{27}. Compare also Fig. 2. On a Poincaré section Σ, one defines the so-called Poincaré map P: For a particle x ∈ Σ, the Poincaré map P(x ) is the point where the flow first returns to Σ. This reduces the three- dimensional, continuous dynamics φ to (discrete) iterations of the two-dimensional map P. Tucker then analyzes the dynamics of P. Trapping Region Regarding Property 1, Tucker proves that there is a (compact) trapping region N ⊆ Σ, such that solutions starting in N will remain in N. He does so by subdividing N into a large number of small rectangles. For every small rectangle, Tucker’s program computes safe numeric bounds for all solutions evolving from the small rectangle. In a number of time-discretization steps, the evolution is followed until it eventually returns to Σ. Upon return, the program checks that the returned enclosure is contained in N.Ifthis process succeeds for every small rectangle, one can conclude the following theorem. 123 76 F. Immler Theorem 1 (Trapping Region) ∀x ∈ N − Γ. P(x ) ∈ N Note that there exists a set Γ on which P is not defined: Γ is the set of points, from which solutions tend to the origin in infinite time. Γ is therefore explicitly excluded in the above theorem. Sensitive Dependence Regarding Property 2, sensitive dependence on initial conditions can be quantified with the help of the derivative: A deviation in the direction of a vector v ∈ R is propagated (in linear approximation) like the derivative at x, i.e., P(x + v) ≈ x + DP| · v. Here DP| · v is the matrix of partial derivatives of P (the Jacobian matrix) at the point x, multiplied with the vector v. A mathematically precise notion of chaos is given by the class of singular hyperbolic sys- tems [36]. The term singular denotes the special case where the system contains a hyperbolic fixed point (which renders P undefined on Γ ). A hyperbolic system contracts deviations in stable directions and expands deviations in unstable directions. Both are relevant for the dynamics of the attractor: Stable directions make solutions tend to the attractor, whereas unstable directions lead to sensitive dependence on initial conditions. Tucker proves that the Lorenz attractor is hyperbolic (in fact singular hyperbolic, we discuss how the hyperbolic fixed point is addressed with normal form theory in the next paragraph) by providing safe overapproximations for the unstable direction: Every x ∈ N is equipped with a cone C(x ) (compare Fig. 3), which contains the unstable direction. This is also verified by Tucker’s computer program: In addition to the Poincaré map, the program keeps bounds on its matrix of partial derivatives. The program tracks how initial deviations (inside the cone associated to an initial rectangle) are propagated by the derivative DP.The cone field needs to be forward invariant (otherwise it would not contain the unstable direction) and the expansion needs to be large enough that the enclosed directions are actually expand- −1 ing. Tucker’s program establishes factors E (x ) and E (x ), which quantify the expansion properties of P: Theorem 2 (Derivatives, Cones, and Expansion) 1. ∀x ∈ N − Γ. ∀v ∈ C(x ). DP| · v ∈ C(P(x )) 2. ∀x ∈ N − Γ. ∀v ∈ C(x ). DP| · v ≥ E (x ) |v| −1 3. ∀x ∈ N − Γ. ∀v ∈ C(x ). DP| · v ≥ E (P(x )) v This theorem states that 1., the cone field C is forward invariant under the action of the derivative of P: the image of every cone is slimmer than the cones onto which they are mapped. 2., the vectors v satisfy lower bounds on how much they are expanded: the length DP| · v of the return of the deviation vector v is lower bounded by its length v times an −1 expansion factor E (x ). They also satisfy a pre-expansion bound E (x ) (this does not denote ) for the pre-image of x, which is required for technical reasons in Tucker’s proof. Normal Form Theory For the Lorenz equations, the origin (0, 0, 0) is a hyperbolic fixed point. The origin is a fixed point, because the ODE evaluates to 0. It is hyperbolic, because solutions tend to it in the two stable directions given by the y and z axis and expands in the unstable direction given by the x-axis. This hyperbolic fixed point poses problems for the aforementioned approach of using rigorous numerical methods: there are solutions that tend to the origin as time goes to infin- ity. In such a situation, a time-discretization algorithm is at a loss, because it would need infinitely many steps. To remedy this problem, Tucker’s program interrupts computations in 123 A Verified ODE Solver and the Lorenz Attractor 77 Fig. 3 Enclosures for the flow and cones evolving from X =[4.375, 4.4]×[2.77, 2.79]×{27} with a ◦ ◦ representation of a cone between 1.5 and 11.5 (detail on the right) a small cube L =[−0.1, 0.1]×[−0.1, 0.1]×[−0.1, 0.1] around the origin. Inside the cube, where the numerical methods would fail, the evolution of solutions can be described with classical, analytical means: more than half of Tucker’s thesis is devoted to accurate analytical expressions for the flow inside the cube L. These expressions can be used to provide explicit bounds on how solutions exit the cube L and continue with numerical computations. Informal Theorem 3 There is an explicit form that bounds the dynamics inside the cube L =[−0.1, 0.1]×[−0.1, 0.1]×[−0.1, 0.1]. 1.4 Outline of This Article This article is about the formalization of a rigorous ODE solver that computes the Poincaré map P and its derivative DP. It is sufficiently efficient and precise to certify Tucker’s numer- ical results. In particular, computing with the verified algorithm proves Theorems 1 and 2. In fact, also Fig. 3 was created from the output of this verified algorithm. As a matter of course, since the ODE solver computes Poincaré maps and derivatives thereof, it is proved correct with regard to a formalization of these concepts in Isabelle/HOL. This formalization, as well as the verified algorithms are generic: they are independent of the underlying ODE or dimension. This article summarizes previous work of the author (and contributions by Johannes Hölzl and Christoph Traut) [18–26] and describes the state of the formalization to which it evolved over time. It shows how the various parts fit together for the final result of certifying Tucker’s computations. The following explains how this article is structured and details on the relation to earlier publications. The abstract mathematics needed for the formalization is a theory of ODEs and Poincaré maps together with their (differentiable) dependence on initial conditions. This is presented in Sect. 2, which summarizes a more comprehensive journal article [25], in which the author extends the conference paper [26] with a formalization of Poincaré maps. The foundations, in particular existence and uniqueness of solutions was proved in [24]. 123 78 F. Immler Tucker used the library Profil/BIAS for an implementation of the Euler method in interval arithmetic. Our approach to rigorous numerics is at first agnostic about the concrete type of enclosures (Sect. 3). The main instantiation, affine arithmetic, is presented in Sect. 3.4.Parts of this section were part of earlier publications [18,21]. When working with affine arithmetic, the enclosures are zonotopes (centrally symmetric polytopes) instead of intervals. Because zonotopes are a more complex than intervals, geo- metric operations like intersections with the Poincaré section Σ are more challenging. The verification leads to a digression into computational geometry in Sect. 4, which is based on an earlier publication [19]. Tucker’s program adaptively splits reachable sets and therefore maintains a collection of sets. Section 5 describes how we formalize generic data structures to maintain such collec- tions and use stepwise refinement from nondeterministic abstract specifications to a concrete deterministic implementation. This has not been published before. A more ad-hoc formalization of similar algorithms has been described in earlier work [21], which forms the basis of Sect. 6. Here the presentation is w.r.t. the new, generic framework and extended for derivatives of Poincaré maps. Section 7 is a novel contribution: it describes how the trapping region N, the cone field C −1 and the expansion estimates E , E are defined formally and how the verified ODE solver is set up to certify the results of all of Tucker’s computations. Earlier work [20] was not capable of handling derivatives, and had no formalization of Poincaré map. All of the development described here is available for Isabelle2017 in the Archive of Formal Proof [22,23]. In particular everything displayed as Theorem possesses a formal counterpart. 2Mathematics The required mathematical background consists of mostly standard results which can be found in every textbook on ODEs and dynamical systems. Thanks to a sufficient background theory, the formalization can mostly follow the presentations in such textbooks. We therefore focus on peculiarities of our formalization which are due to Isabelle/HOL and its type system, in particular the use of type classes and the lack of dependent types. Because of the type class based formalization of topological structures, type definitions are used to formalize function spaces, where the Transfer and Lifting tools [15] provide excellent support. Moreover, there are no dependent types in Isabelle/HOL. In situations where this would be more natural, an encoding (e.g., ext-cont in the following Sect. 2.1.2) is necessary. 2.1 Type Classes for Mathematics in Isabelle/HOL In Isabelle/HOL, many of the mathematical concepts (in particular spaces with a certain structure) are formalized using type classes. Isabelle/HOL features axiomatic type classes [11, 38]. The purpose of an axiomatic type class is to specify operations which satisfy given properties for a class of types. The advantage of type class based reasoning is that most of the reasoning is generic: formalizations are carried out in the context of type classes and can then be used for all types inhabiting that type class. For generic formalizations, we use Greek letters α, β, γ and name their type class con- straints in prose. I.e., when we write that we “consider a topological space” α, then this result is formalized generically for every type α that fulfills the properties of a topological space. 123 A Verified ODE Solver and the Lorenz Attractor 79 The spaces we consider are topological spaces with open sets, (real) vector spaces with addition +: α → α → α and scalar multiplication (_)(_) : R → α → α. Normed vector spaces come with a norm |(_)|: α → R. Complete normed vector spaces are called Banach spaces. Much of the theory has been ported from Harrison’s theory of Euclidean space [12] and has been generalized to the hierarchy of type classes for mathematics in Isablle/HOL [14]. 2.1.1 Vectors in Euclidean Space Because of Isabelle/HOL’s restrictive type system (no dependent types), the abstract concept of vectors is notorious for demanding workarounds. In Isabelle/HOL, one tends to use a type class based encoding. We work with a type class for Euclidean space that fixes an order on the Basis elements and therefore enables operations eucl-of-list : R list → α and list-of-eucl : α → R list if α is a Euclidean space. All (finite) vectors of real numbers are instances of the class Euclidean space. This includes real numbers R, complex numbers C, tuples α × β for Euclidean spaces α, β, and Harrison- 2 ι style [12] vectors α for a finite type ι. 2.1.2 Bounded Continuous Function We motivate bounded continuous functions with the Picard–Lindelöf theorem, which guar- antees the existence of a unique solution to an initial value problem. For an ODE f with initial value x at time t , a unique solution on the time interval [t , t ] is constructed by 0 0 0 1 considering iterations of the following operator for continuous functions φ :[t , t ]→ R : 0 1 P(φ):= λt. x + f (τ, φ (t )) dτ From a mathematician’s point of view, P operates on the Banach space of continuous func- tions on the compact domain [t , t ] and therefore the Banach fixed point theorem guarantees 0 1 the existence of a unique fixed point (which is by construction the unique solution). In order to formalize this in Isabelle/HOL, there are two obstructions to overcome: First, the concept of Banach space is a type class in Isabelle/HOL, so we need to introduce a type for the mappings φ :[t , t ]→ R from above. But this poses the second problem: functions 0 1 in Isabelle/HOL are total and types must not depend on term parameters like t and t . 0 1 We work around these restrictions by introducing a type of bounded continuous func- tions, which is a Banach space and comprises (with a suitable choice of representations) all continuous functions on all compact domains. n m typedef α → β:={ f : R → R | f continuous on α ∧ (∃B. ∀t. ft≤ B)} bc In order to define operations on type α → β, the Lifting and Transfer package [15]is bc an essential tool: operations on the plain function type α → β are automatically lifted to definitions on the type α → β when supplied with a proof that functions in the result are bc bounded continuous under the assumption that argument functions are bounded continuous. We write application $ of a bounded continuous function f : α → β to an element bc bc x : α as follows. Vectors of length n are represented by a type of functions ι → α,where n equals the cardinality of the finite type ι. 123 80 F. Immler Definition 1 (Application of Bounded Continuous Functions) ( f $ x ) : β bc Bounded continuous functions form a normed vector space. The norm on α → β is the bc supremum of the range and the vector space operations +, · are defined pointwise. Definition 2 (Normed Vector Space of Bounded Continuous Functions) f : = sup { f $ x| x ∈ α} bc ( f + g) $ x := f $ x + g $ x bc bc bc (a · f ) $ x := a · ( f $ x )) bc bc The type → with the above operations forms a complete normed vector space (a Banach bc space). This allows us to use the Banach fixed point theorem for operators on this type. In order to be able to use this for the operator P from above, we represent functions on a compact interval [a, b] as an element of type → by extending the function continuously bc outside the domain with the help of clamp : clamp x := if x ≤ a then a else (if x ≥ b then b else x ) [a,b] (ext-cont f ) $ x := f (clamp x )) [a,b] bc [a,b] n n With the help of ext-cont we can apply P : (R → R ) → (R → R ) to a continuous bc bc function φ : R → R (that is assumed to be continuous on an interval [a, b]) by writing P(ext-cont φ). According to the Banach fixed point theorem there exists a unique fixed [t ,t ] 0 1 point φ : R → R where P(φ ) = φ and the unique solution of the initial value bc bc bc bc problem is the function λt.φ $ t of type R → R . bc bc The usage of the type → caused minor technical obstructions, but otherwise enabled a bc natural and abstract proof. 2.1.3 Bounded Linear Functions Similar to the type of bounded continuous functions, we also introduce a type of bounded linear functions (also known as continuous linear functions) For vector spaces α and β, a linear function is a function f : α → β that is compatible with addition and scalar multiplication. linear f :=∀xy c. f (c · x + y) = c · f (x ) + f (y) Let us assume normed vector spaces α and β. Linear functions are continuous if the norm of the result is linearly bounded by the norm of the argument. We cast bounded linear functions α → β as a type α → β in order to make it an instance of Banach space. bl { } typedef α → β:= f : α → β | linear f ∧∃K . ∀x .  f (x )≤ K x bl The construction is very similar to bounded continuous functions and we write bounded linear function application ( f · x ). Vector space operations are also analogous to → .The bl bc usual choice of a norm for bounded linear functions is the operator norm: the maximum of the image of the bounded linear function on the unit ball. With this norm, α → β forms a bl normed vector space and we prove that it is Banach if α is a normed vector space and β is Banach. Definition 3 (Norm in Banach Space → )For f : α → β, bl bl f : = onorm (λy. f · y) = max { f · y|y≤ 1} bl bl 123 A Verified ODE Solver and the Lorenz Attractor 81 Having (bounded) linear functions as a separate type makes many formulations eas- ier. For example, consider Harrison’s formalization of multivariate analysis (from which Isabelle/HOL’s analysis descended). In Harrison’s formalization continuity is formalized for n m functions f of type R → R . (continuous f (at x )) = (∀e > 0. ∃d > 0. ∀y. x − y < d ⇒ fx − fy < e) Most of Harrison’s formalization is geared towards viewing derivatives as linear functions n m of type R → R . For continuously differentiable functions, one therefore needs to reason n n m about functions f : R → (R → R ),where f x is the derivative of f at a point x. Continuity of f is written in an explicit ε-δ form and involves the operator norm onorm : n m (R → R ) → R, which is quite verbose: (∀e > 0. ∃d > 0. ∀y. |x − y| < d ⇒ onorm (λv. f x v − f y v) < e)) The ε-δ form could of course be captured in a separate definition, but this would be very similar to the definition of continuity and would introduce redundancy. In the Isabelle/HOL formalization, continuous is defined for functions f : α → β for topological spaces α and β.If α and β are normed vector spaces, the above equality for continuous holds in Isabelle/HOL, too. And indeed, the norm of bounded linear functions is defined using onorm such that onorm (λv. ( f x ) · v − ( f y) · v) = f x − f y bl bl holds. Then, continuity of a derivative f : α → (α → β) can simply be written as bl (continuous f (at x )), which is a better abstraction to work with and also avoids redundant formalizations for different kinds of continuity. 2.2 Dynamical Systems An ODE induces a continuous dynamical system via the notion of flow. A standard technique to reason about such systems is its Poincaré map. We keep the presentation at a high level, since details can be found in the publications [25,26]. 2.2.1 The Flow We consider an autonomous ODE with right hand side f . Under mild assumptions, there exists a solution φ(t ), which is unique for an initial condition x (0) = x and satisfies the differential equation: x˙(t ) = f (x (t )) To emphasize the dependence on the initial condition, one writes φ(x , t ) for the solution. This solution depending on initial conditions is called the flow of the differential equation: Definition 4 (Flow)The flow φ(x , t ) is the (unique) solution of the ODE x˙(t ) = f (x (t )) with initial condition φ(0) = x The flow is only well-defined on the so-called existence interval of the solution, which depends on the initial value. Definition 5 (Existence Interval) t ∈ ex-ivl (x ) ⇒ φ(x , t ) = f (φ (x , t )) 0 0 0 The flow φ and the existence interval ex-ivl provide a clean interface to talk about solu- tions of ODEs. The property of the generic notion of flow makes it possible to easily state composition of solutions and to algebraically reason about them. Flowing from x for time s + t is equivalent to first flowing for time s, and from there flowing for time t: 123 82 F. Immler Theorem 4 (Flow Property) {s, t, s + t}⊆ ex-ivl (x ) ⇒ φ(x , s + t ) = φ(φ(x , s), t ) 0 0 0 For Tucker’s proof, one needs to study how sensitive the flow depends on perturbations of the initial value. We use two main results: One, the flow depends continuously on initial values. Two, if the ODE f is continuously differentiable, then so is the flow. We first take a look at the domain Ω = {(x , t ) | t ∈ ex-ivl (x )} ⊆ X × T of the flow. (t, x ) ∈ Ω means that we can flow a point starting at x for at least time t. Intuitively, solutions starting close to x can be followed for times that are close to t. In topological parlance, the state space is open. Theorem 5 (Open State Space) open Ω One can show that solutions deviate at most exponentially fast: ∃K . φ(x , t ) −φ(y, t ) < K |t | x −ye (using Grönwall’s lemma). Therefore, by choosing x and y close enough, one can make the distance of the solutions arbitrarily small. In other words, the flow is a continuous function on the state space: Theorem 6 (Continuity of Flow) continuous-on Ωφ Continuity states that small deviations in the initial values result in small deviations of the flow. One can be more precise on the way initial deviations propagate. The propagation of initial deviations through the flow (φ := λx.φ(x , t )) can be approximated by a linear function, the derivative Dφ| · v ≈ φ(x , t ) − φ(x + v, t ). We formalize the fact that the derivative of the flow is the solution of a differential equation in the space of bounded linear mappings, the so-called variational equation. Theorem 7 (Variational Equation) W (t ) = D f | · W (t ) φ(x ,t ) bl W (0) = 1 bl Solving this ODE numerically gives a means to obtain numerical bounds on the derivative, which is the approach that we pursue in our algorithms. 2.2.2 The Poincaré Map The Poincaré map is an important tool for studying dynamical systems. Whereas the flow describes the evolution of a continuous system with respect to time, it is the Poincaré map that allows us to study the evolution with respect to some space variables. A Poincaré section is a subset Σ of the state space, which is in general given as an implicit surface Σ ={x | s(x ) = c} with continuously differentiable s. For Tucker’s proof, one chooses s(x , y, z) = z and c = 27. The Poincaré map P(x ) is defined as the point where the flow starting from x first hits the Poincaré section Σ. It is defined with the help of the first return time τ(x ). τ depends on the flow φ (and therefore on the ODE f ) and the Poincaré section Σ, but we keep those dependencies implicit. Definition 6 (First Return Time) τ(x ) is the least t > 0 such that φ(x , t ) ∈ Σ. Obviously, τ is only well-defined for values that actually return to Σ, which we encode in the predicate returns-to : 123 A Verified ODE Solver and the Lorenz Attractor 83 Definition 7 returns-to (Σ, x ):=∃t > 0.φ(x , t ) ∈ Σ The return time can then be used to define the Poincaré map as follows: Definition 8 (Poincaré map) P(x ):= φ(x,τ(x )) It is interesting to note that this way of defining the return time and Poincaré map differs from the usual approach in textbooks. Textbooks study Poincaré maps in a neighborhood around a periodic point x ∈ Σ, i.e., P(x ) = x. This makes it easy to directly apply the implicit function theorem and transfer continuity and differentiability from the flow to the Poincaré map while guaranteeing that τ and P are well-defined. Also, one views P as a mapping on Σ, i.e. P : Σ → Σ. Tucker’s proof, however, requires a more flexible notion of Poincaré map and our notion of τ is more flexible: it is well-defined also for values outside of Σ. This enables reasoning about intermediate sections: Tucker, e.g., composes a sequence of local Poincaré maps between intermediate sections Σ, Σ , ··· ,Σ ,Σ in order to get bounds on the global Poincaré map 1 n Σ → Σ. The goal of Tucker’s computations is a sensitivity analysis of the flow of the Lorenz system and of its Poincaré map. Its derivative can be given in terms of the derivative of the flow and the function s defining the implicit surface for Σ ={x | s(x ) = c}. Theorem 8 (Derivative of Poincaré map) Ds| · (Dφ| · h) P(x ) (x ,τ (x )) DP| · h = Dφ| · h − f (P(x )) x (x ,τ (x )) Ds| · ( f (P(x ))) P(x ) For a rough intuition, the derivative DP| ·h of the Poincaré map is related to the derivative of the flow Dφ| · h. But it needs to corrected in the direction f (P(x )) in which the (x ,τ (x )) flow passes through Σ, because P varies only on Σ and not through it. This correction factor also depends on the tangent space Ds| of the section Σ at P(x ). P(x )) 3 Rigorous Numerics Rigorous (or guaranteed) numerics means computing with sets that are guaranteed to enclose the real quantities of interest. Enclosures can in principle be any data structure that represents sets of real values. Popular choices are intervals, zonotopes, or Taylor models. For formal- izations, it is useful to have a deep embedding of arithmetic expressions as done, e.g., by Dorel and Melquiond [34]aswellasHölzl [13]. This work builds on Hölzl’s language of arithmetic expressions given in excerpts in Fig. 4. Independent of the choice of representation of enclosures, a rigorous approximation scheme approx needs to satisfy that for all lists of real values xs in an enclosure XS,the interpretation [[e]] of the expression e is contained in the result of approximation scheme xs evaluated for the enclosure XS. We could call this the fundamental property of rigorous numerics. x ∈ XS ⇒ [[e]] ∈ approx eXS xs 123 84 F. Immler Fig. 4 Data type of arithmetic expressions and interpretation function (λexs. [[e]] ) : aexp → R list → R xs The key approach in our formalization is to remain agnostic about the concrete approx- imation scheme as long as possible and formalize results on the level of deeply embedded aexp expressions. The central result is an implementation of a second order Runge–Kutta method on the level of aexp expressions. As concrete instance for an approximation scheme approx , we use affine arithmetic [6], an improvement over interval arithmetic that tracks linear dependencies between program variables. 3.1 Expressions for Vectors To represent vectors, we use lists of expressions. A list of expressions es : aexp list is interpreted with (λes vs. [[es ]] ) : aexp list → R list → α componentwise as Euclidean vs space α: [[es ]] = eucl-of-list (map (λe. [[e]] ) es) vs vs In contrast to the interpretation, approximation of a list of expressions should not be componentwise: an approximation function for lists of expressions should be of type approx : aexp list → R list set → R list set , which allows approx to keep track of dependencies between the components of the result. If the type were e.g., aexp list → R set list → R set list , this could only represent the Cartesian product of the component enclosures. n m For a function f : R → R , a deep embedding f is a list of expressions (of length m), that is interpreted over a list of n variables. [[ f ]] = f (eucl-of-list xs) e xs The derivatives with respect to one variable can be computed symbolically from the struc- ture of the expression. This can also be used to compute partial derivatives on the level of expressions. In the multivariate setting, the derivative D f | of f at x is the matrix of its partial derivatives. In general, we can represent matrices as a flat list (according to eucl-of-list /list-of-eucl which are also defined for matrices). For computing derivatives, how- ever, we directly produce an expression that is interpreted as the product of the derivative matrix with a vector: [[D (n, f ,v )]] := D f | · [[v ]] e e e xs eucl-of-list xs bl e xs D takes the derivative with respect to the first n variables, and xs is assumed to be of length at least n. This way, we can produce expressions for higher derivatives: D (n, f ,v ):= f e e e i +1 i D (n, f ,v ):= D (n, D (n, f ,v ), v ) e e e e e e e e Note that the proper interpretation can only be expressed in Isabelle’s type system for fixed values of i: the resulting object is an i-linear function, so the resulting type depends on 123 A Verified ODE Solver and the Lorenz Attractor 85 a term argument. This could also be encoded as functions taking lists as arguments, but fixed values of i suffice for our purposes and we find the interpretation as curried linear mappings more natural. E.g., [[D (n, f ,v )]] =[[ f ]] e e x e x [[D (n, f ,v )]] =D f | · [[v ]] e e x eucl-of-list x bl e [[D (n, f ,v )]] =D(λy. D f | )| · [[v ]] · [[v ]] e e x y eucl-of-list x bl e bl e ··· 3.2 A Runge–Kutta Method On the level of expressions, we verified a two-stage Runge–Kutta method rk (x ) = x + 1 1 h · ψ (x ), with ψ (x ) = (1 − ) f (x ) + f (x + hpf (x )). This Runge–Kutta method h h 2 p 2 p rk (x ) approximates the solution up to third order: |φ(x , h) − rk (x )|∈ O(h ).The h 0 0 h 0 third order term stems from (multivariate) Taylor series expansions of the solution φ and the approximation scheme rk .Ifweset f := λx . D f | and f := λx . D f | , then the h x x remainder is contained in the convex hull of any set that contains rk-remainder (s , s ) for h 1 2 all s , s ∈[0, 1]. 1 2 h 1 rk-remainder (s , s ):= · f x (hs + t ) h 1 2 1 2 3 · f (x (hs + t )) · f (x (hs + t )) bl 1 bl 1 + f x (hs + t ) · f (x (hs + t )) · ( f (x (hs + t ))) 1 bl 1 bl 1 − f x (t ) + hps f (x (t )) · f (x (t )) · f (x (t )) 2 bl bl Theorem 9 (Runge–Kutta Method with Remainder Term) φ(x , h) ∈ rk (x ) + convex-hull (rk-remainder ([0, 1], [0, 1])) 0 h 0 h In order to use this theorem for rigorous numerical computations, we produce a deep embedding of the expressions for rk and rk-remainder . We do so for an arbitrary ODE h h n n f : R → R with deep embedding [[ f ]] = f (eucl-of-list x ). Expressions for f and f e xs 1,2 are computed symbolically from f via D from the previous Sect. 3.1. 3.3 Straight Line Programs The expression for rk-remainder from the previous Sect. 3.2 contains common subexpres- sions. This is not desirable because one would need to perform redundant computations. We therefore follow Dorel and Melquiond’s [34] approach and use straight-line programs with static single assignment instead of plain expressions. For us, a straight line program is just a list of arithmetic expressions, which is interpreted according to function slp : aexp list → R list → R list : slp [] x = x slp (e::es) x = slp es ([[e]] ::x ) The idea is that a straight line program only contains unary or binary operations, although this is not required by the definition. The result of the operation is put on top of the evaluation stack. 123 86 F. Immler The following example illustrates sharing the term x + y in the expression (x + y)(x + y): slp [Add (Var 0)(Var 1), Mult (Var 0)(Var 0)][x , y]=[(x + y)(x + y), x + y, x , y] We provide a function slp-of , which eliminates common subexpressions by traversing an expression bottom-up and saving subexpressions in a map that gives the index of the subex- pression in the resulting straight line program. At run-time (this is important to be able to use the ODE solver as a stand-alone tool), in an initialization phase, the ODE solver computes symbolically the derivatives in the expression for rk and rk-remainder , does constant propagation (as derivatives can produce 0 constants, this is beneficial) and then compiles the resulting expression with slp-of into a straight-line program, which is then used in the course of approximating the ODE in a series of steps. 3.4 Affine Arithmetic Up to now, we have kept the discussion on the level of expressions, let us now motivate affine arithmetic as a concrete approximation scheme. The most basic data structure to represent sets is closed intervals [a, b]={x | a ≤ x ≤ b}, but those suffer from the wrapping effect: rotated boxes cannot be represented without large overapproximations. Moreover dependencies between variables are lost, e.g. for an enclosure x ∈[0, 2], the expression x − x evaluates to [−2, 2] in interval arithmetic whereas the exact result would be representable as the interval [0, 0]. Affine arithmetic [6] improves over interval arithmetic by tracking linear dependencies. An affine form A is a function where only finitely many arguments map to nonzero values. It is interpreted for a valuation ε : N → R : affine ε A:= A + ε A 0 i i Looking at the interpretation, one often calls the terms ε noise symbols and A generators. i i The idea is that noise symbols are shared between affine forms and that they are treated symbolically, as formal parameters: the sum of two affine forms is given by the pointwise sum of their generators, and multiplication with a constant factor is also done componentwise. affine ε(A + B):= (A + B ) + ε (A + B ) 0 0 i i i affine ε(cA):= cA + ε (cA ) 0 i i The range of an affine form is the set of all affine evaluations where the noise symbols range over the closed interval [−1, 1]. For the range of a list of affine forms, those are evaluated jointly for the same valuation of the noise symbols, reflecting the intuition that those are shared. range A:= {affine ε A |∀i. − 1 ≤ ε ≤ 1} { } joint-range AS:= map (affine ε) AS |∀i. − 1 ≤ ε ≤ 1 As a concrete example, let us examine how affine arithmetic handles the dependency problem in the introductory example x − x for x ∈[0, 2]. The interval [0, 2] is represented by the affine form 1 + 1 · ε . This is the affine form given by the function X:= (λi. if i = 0∨i = 1 then 1 else 0). For this function, range X =[0, 2] holds. Then, in affine arithmetic, 123 A Verified ODE Solver and the Lorenz Attractor 87 (1 + 1 · ε ) − (1 + 1 · ε ) = 0 + 0 · ε , which corresponds to the constant zero function. 1 1 1 Therefore range (X − X ) ={0}. In general, with the help of range and joint-range , we can express correctness of a binary operation like addition e.g., as follows: [a, b]∈ joint-range [A, B]⇒ (a + b) ∈ range (A + B) Nonlinear operations like multiplication or division are linearized, adding the linearization error as a fresh noise symbol. Consider e.g., multiplication: (affine ε A) ∗ (affine ε B) = A B + ε (A B + A B ) + ε A ε B 0 0 i 0 i i 0 i i i i i i >0 i >0 For a proper valuation with ε ∈[−1, 1], the last summand on the right can be bounded by ( |A |)( |B |). Therefore, if k is fresh in A and B, one can set i i i >0 i >0 affine ε(A ∗ B):= A B + ε (A B + A B ) + ε |A | |B | 0 0 i 0 i i 0 k i i i i >0 i >0 and the k-th generator bounds the linearization error such that multiplication of affine forms is conservative: [a, b]∈ joint-range [A, B]⇒ a ∗ b ∈ range (A ∗ B) Similar to the additional noise symbol for a linearization error, also round-off errors can be included as additional noise symbols. We provide affine approximations for the primitive functions listed in Fig. 4. Expressions (and straight line programs) involving these functions can then be approximated by recursively keeping track of the next fresh noise symbols. During longer computations more and more noise symbols will be added to the affine form, impairing performance in the long run. The number of noise symbols can be reduced by summarizing (or condensing) several noise symbols into a new one. This process discards the correlation mediated by the summarized noise symbols, so a trade-off needs to be found between precision and efficiency. We consider a list of affine forms AS and use the notation AS := map (λA. A ) AS. We call total deviation |AS|: = map (λA. |A |) AS the i i i componentwise sum of absolute values. We summarize all symbols i with |AS |≤ r |AS| for a given summarization threshold r. We found that it is important to perform the above comparison componentwise and not take (like proposed for other implementations of affine arithmetic [6,42]) the infinity norm on both sides. This is of particular importance when components differ a lot in magnitude. Apart from looking at affine forms as a formal sum, the joint-range of a list of affine forms can also be interpreted geometrically as zonotopes: centrally symmetric, convex polytopes. A zonotope can be visualized as the Minkowski sum (the set of all possible sums of elements of two sets X ⊕ Y ={x + y. x ∈ X ∧ y ∈ Y }) of the line segments defined by the generators. For example, Fig. 6 depicts a two-dimensional zonotope with three generators, [−1, 1]a ⊕[−1, 1]a ⊕[−1, 1]a . Figure 5 contains a three-dimensional zonotope with 1 2 3 three generators (a parallelotope), namely the zonotope Z defined as follows. ⎛ ⎞ [1ε + 0ε + 1ε , 1 2 3 ⎝ ⎠ 0ε + 2ε + 5ε , Z:= joint-range 1 2 3 0ε + 0ε + 20ε ] 1 2 3 123 88 F. Immler Z G 20 G -5 -10 -15 -20 -2 -4 x1 -5 -6 -10 x2 -8 -15 Fig. 5 Three dimensional zonotope Z and intersection with hyperplane G 4 Computational Geometry An important step for Tucker’s proof is the reduction to a Poincaré map: intersecting the flow of the ODE with a plane in the state space. In our algorithms, the flow is approximated with affine arithmetic expressions, therefore enclosed by zonotopes. In order to compute where the flow intersects the hyperplane, one needs to compute the intersection of the enclosing zonotope with the hyperplane (see Fig. 5). This is an interesting geometric problem and we verified an approximative algorithm due to Girard and Le Guernic [7]. At its core, the algorithm is similar to convex hull compu- tations. We can build on a nice abstraction to reason about it, namely Knuth’s theory of counterclockwise (ccw) systems [27]. We needed, however, to extend Knuth’s theory from discrete to continuous sets of points. 4.1 Girard and Le Guernic’s Algorithm The complexity for computing the exact intersection of a zonotope with a hyperplane grows exponentially with the number of generators. An overapproximation of the zonotope before computing the intersection is possible but can lead to overly coarse approximations. Therefore Girard and Le Guernic [7] proposed a way to directly compute overapproximations to the intersection. The first idea is to overapproximate a given set X tightly from a set D of directions, which can be chosen arbitrarily. For every direction d ∈ D ⊆ R , the infimum m and supremum M of the sets {x , d. x ∈ X } needs to be determined (_, _ denotes the inner product, also known as the dot product). Geometrically speaking, m and M give the position d d of two hyperplanes with normal vector d. The two hyperplanes bound X from below and above, respectively. An overapproximation P is then given by the points between all of these hyperplanes: X ⊆ P ={x ∈ R . ∀d ∈ D. m ≤x , d≤ M } d d The second observation of Girard and Le Guernic is that when the set X is the intersection of some set Z with a hyperplane G ={x . x , g= c}, then the computation of the overap- proximation P can be reduced to a two-dimensional problem with the linear transformation n 2 Π : R → R , Π (x ) = (x , g, x , d). g,d g,d x3 A Verified ODE Solver and the Lorenz Attractor 89 Fig. 6 Corners c and edges of a zonotope { ε · a |−1 ≤ ε ≤ 1}, generators a , a , a , intersecting i i i i 1 2 3 line L Lemma 1 (Reduction to Dimension Two) {x , d. x ∈ Z ∩ G}={y.(c, y) ∈ Π (Z )} g,d This lemma is an easy consequence of the definitions of G and Π . For every direction g,d d, the theorem allows to reduce the computation of the intersection Z ∩ G on the left-hand side to the intersection of the projected two-dimensional zonotope Π (Z ) with the vertical g,d line L ={(x , y). x = c}. Computing the intersection of a two-dimensional zonotope like the one given in Fig. 6 and a vertical line L can by done by computing bounds on the intersection of the vertical line L with every edge. This is easy and intuitive. The more challenging part is to compute the set of edges of a two-dimensional zonotope, which we sketch in the following. 4.2 Computing the Set of Edges. First of all, one can assume that all generators point upwards. One then starts at the lowest corner (c in Fig. 6) and appends to it the “rightmost” generator a (twice) to reach c .One 0 1 1 then continues with the “rightmost” of the remaining generators, a and is in the process essentially “wrapping up” the hull of the zonotope. In order to verify such a process, we need a way to reason about “rightmost” vectors (a total order). Similar ideas of “wrapping up” a set of points also occur for convex hull algorithms. 4.3 Knuth’s CCW System In order to verify geometric algorithms, one needs a formal notion of the geometric concepts involved. For convex hull algorithms, Knuth [27] has given a small theory that axiomatizes the notion of orientation of points. The intuition is that for three points p, q, r in the plane, visiting them in order requires either a counterclockwise (ccw) turn (written pqr) or clockwise (¬pqr) turn. Knuth observed that already few of the properties fulfilled by the ccw predicate pqr suffice to define a theory rich enough to formalize many concepts in computational geometry. The notion of ccw system is a set of points together with a ccw predicate written pqr for points p, q, r. The ccw predicate needs to satisfy the following properties, inspired by 123 90 F. Immler (a) (b) (c) Fig. 7 ccw axioms: dashed predicates are implied by solid ones. a Cyclic symmetry, b interiority, c transitivity the relations satisfied by points in the plane. For all axioms in the following, there is the additional implicit assumption that the involved points are pairwise distinct. For three points, only simple axioms need to be fulfilled: – Cyclic Symmetry: pqr ⇒ qrp – Antisymmetry: pqr ⇒ ¬prq – Nondegeneracy: pqr ∨ prq Cyclic symmetry and the more interesting case of interiority, which involves four points, are illustrated in Fig. 7. Interiority states that if one point t is left of three lines pq qr r p, then the three other points are oriented in a triangle according to pqr. – Interiority: tpq ∧ tqr ∧ trp ⇒ pqr The most important tool for reasoning is transitivity, which involves five points and works if three points p, q, r lie in the half-plane left of the line ts, i.e., tsp ∧ tsq ∧ tsr. Then, fixing t as first element for the ccw relation, we have transitivity in the second and third element: tpq ∧ tqr ⇒ tpr (see Fig. 7c). – Transitivity: tsp ∧ tsq ∧ tsr ∧ tpq ∧ tqr ⇒ tpr The same intuition also holds for the other side of the half-plane: – Dual Transitivity: stp ∧ stq ∧ str ∧ tpq ∧ tqr ⇒ tpr Knuth shows that under the assumptions of Cyclic Symmetry, Antisymmetry, and Non- degeracy, Transitivity holds if and only if Dual Transitivity holds. Knuth requires more than half a page of low level reasoning, but as this reasoning is carried out abstractly in a small first order theory, sledgehammer (Isabelle’s interface to various automatic theorem provers) is able to find a proof that consists of just one single invocation of an automated prover. 4.3.1 Total Order from CCW As sketched earlier, in order to compute the edges of a zonotope, we need to be able to select a “rightmost” element of a set of vectors. With the transitivity relation presented before, we can obtain a total order on vectors which allows us to do just that: Given a center t and another point s, the orientation predicate tpq can be used to define a total order on points p, q in the half-plane left of ts, i.e., p < q iff tpq. From Antisymmetry and Nondegeneracy of the ccw 123 A Verified ODE Solver and the Lorenz Attractor 91 system, we get antisymmetry and totality for the order <. Transitivity of the order < follows from the axiom Transitivity of the ccw system and the assumption that all points are in the half-plane left of ts. This ordering is then used to sort the list of generators such that they actually “wrap up” the zonotope. 4.3.2 Instantiation for Points in the Plane Up to now, our reasoning was based abstractly on ccw systems, but of course we also want to use the results for a concrete ccw predicate. Well known from analytic geometry is the fact that ccw orientation is given by the sign of the following determinant |pqr |: p p 1 x y p ∗ q + p ∗ r + q ∗ r − x y y x x y |pqr |: = q q 1 = x y r ∗ q + r ∗ p + q ∗ p x y y x x y r r 1 x y Points are collinear iff |pqr|= 0. Under the assumption that one works with a finite set of points where no three points are collinear, the following predicate pqr satisfies the axioms of a ccw system. pqr :=|pqr | > 0 Most axioms can easily be proved using Isabelle/HOL’s rewriting for algebraic structures. Transitivity is slightly more complicated, but can also be solved automatically after a proper instantiation of Cramer’s rule, which is easily proved automatically: |tqr ||stp|+|tpq||str | |tpr|= , if |stq| = 0 |stq| 4.3.3 CCW on a Vector Space Knuth presented his axioms with a finite set of discrete points in mind, in our case we need to talk about orientation of arbitrary points in a continuous set. We therefore require consistency of the orientation predicate when vector space operations are involved. One obvious requirement is that orientation is invariant under translation (Fig. 8a). With translation invariance, we can reduce every ccw triple to a triple with 0 as origin, and from there it is easy to state consistency with respect to scaling: If at q, there is a ccw turn to r,then every point on the ray from 0 through q will induce a ccw turn to r (Fig. 8b). Negative scalars can be treated by requiring that reflecting one point at the origin inverts the ccw predicate (Reflection). Furthermore, the addition of vectors q and r, which are both ccw of a line p needs to be ccw of p as well. > > – Translation: (p + s)(q + s)(r + s) = pqr > > – Scaling: α> 0 ⇒ 0(α · q)r = 0qr – Reflection: 0(−p)q = 0qp – Addition: 0pq ⇒ 0 pr ⇒ 0 p(q + r ) The predicate pqr simplifies much of the reasoning, because it satisfies the axioms of a ccw system. It does, however, ignore collinear points and therefore all the points of the zonotope that lie on its edges. In order to also include those into the reasoning, we define the slightly relaxed ccw predicate pqr , which holds for all points on the line through pq and for all points on the half-plane left of pq. pqr :=|pqr|≥ 0 123 92 F. Immler (a) (b) Fig. 8 ccw axioms on a vector space. a Invariance under translation. b Invariance under scaling The situation is as follows: pqr is the actual specification that we care about. But it does not satisfy the axioms of a ccw system, which makes reasoning very convenient. We therefore first prove the corresponding properties for the ccw system pqr . With a simple argument > ≥ on continuity, the results about pqr carry over to pqr and therefore the whole zonotope. 5 Program and Data Refinement We use two different approaches to turn abstract formalizations into executable constructs: In situations where abstract operations directly correspond to concrete, executable ones, we use light-weight data refinement via the code generator. In more demanding situations, we employ a dedicated framework for nondeterministic specifications and stepwise program refinement, namely the Autoref tool. 5.1 Light-Weight Data Refinement For light-weight data refinement [10] via the code generator, abstract operations need to be mapped directly to concrete, executable ones. Examples of such abstract types are affine forms and real numbers. Consider the type of real numbers R. We call it abstract, because as an uncountable set, R is obviously not computable. But one can restrict oneself to working on a computable subset of the real numbers. In our case, we use software floating point numbers F ={m ·2 | m, e ∈ Z} for (unbounded) integers m, e ∈ Z. We instruct the code generator to use an uninterpreted constructor Real-of-Float : F → R to represent real numbers. Operations on real numbers are then computed by pattern matching on that constructor and executing the corresponding concrete operation, e.g., for addition: (Real-of-Float f ) + (Real-of-Float g) = (Real-of-Float ( f + g)) Since this implementing equation is proved as a theorem, such a setup does not change the trusted code base. All one has to ensure is that all abstract operations that occur in the code can be executed in the concrete representation. Affine forms are abstractly a subtype of functions N → R. They are implemented using a type of association lists (N × R) list that are (reverse) strictly sorted according to the keys. This sparse representation is useful because the largest index of a non-zero generator can be directly read off by inspecting only the first element. Adding a fresh generator can be 123 A Verified ODE Solver and the Lorenz Attractor 93 done by simply prepending the new element. Binary operations are efficiently and easily implemented by merging the two lists of generators. 5.2 Autoref We use Lammich’s framework Autoref for (automatic) refinement [29,30] in Isabelle/HOL. Autoref allows the user to specify algorithms on an abstract level and provides tool support for stepwise refinement [1] to concrete, executable implementations. In this section we present a setup of Autoref for rigorous numerical algorithms: We provide abstract specifications for elementary operations of common rigorous numerical algorithms as well as suitable implementations. 5.2.1 Nondeterministic Specifications An important insight when verifying algorithms that use rigorous numerical enclosures is the fact that, for correctness, any enclosure of the real quantity suffices. We model this with appropriate nondeterministic specifications. Autoref is based on a nondeterminism monad α nres , where programs can either fail or yield a set of values as result. datatype α nres = FAIL | RES (α set ) The refinement relation ≤ on nres has FAIL as top element and RES S ≤ RES T iff S ⊆ T . For deterministic results, we write return x:= RES {s}. We write for specifications spec P:= RES {x . Px } the result of all values satisfying the predicate P. This allows one to specify correctness, e.g, of a program f whose inputs x satisfy the precondition P and every possible value y in its nondeterministic result satisfies the post- condition Q: ∀x . Px ⇒ fx ≤ spec (λy. Qy). In this setting, we specify a set of operations that are useful in the context of verifying rigorous numerical algorithms, i.e., algorithms that manipulate enclosures. These operations are best modeled nondeterministically, because one is often only interested in some safe result. Subdivisions are a means to maintain precision, we therefore have the following abstract specifications for splitting a set (with or without the possibility to perform overapproxima- tions): split-spec X:= spec (λ(A, B). X ⊆ A ∪ B) split-spec X:= spec (λ(A, B). X = A ∪ B) The following specifications yield some lower/upper bound on the set, not necessarily exact: Inf-spec X:= spec (λi. ∀x ∈ X. i ≤ x ) Sup-spec X:= spec (λs. ∀x ∈ X. x ≤ s) Depending on the concrete representation of sets, one might not be able to decide certain properties, but only give a positive answer if the precision is sufficient. We therefore have a specification that may guarantee disjointness. disjoint-spec XY:= spec (λb. b ⇒ X ∩ Y =∅) As seen in the previous Sect. 4, depending on the data structure, one can not (or does not want to) compute an exact representation for the intersection of sets. These specifications 123 94 F. Immler allow one to overapproximate an intersection, while guaranteeing that the result does not exceed one of the arguments: inter-spec XY:= spec (λR. X ∩ Y ⊆ R ∧ R ⊆ X ) inter-spec XY:= spec (λR. X ∩ Y ⊆ R ∧ R ⊆ Y ) To bridge the gap to concrete numerical computations and the results from Sect. 3,weuse a specification for overapproximations of evaluating straight line programs: approx-slp-spec slp X:= spec (λR. ∀x ∈ X. [[slp]] ∈ R) 5.2.2 Refinement Relations In Autoref, specifications in the nres monad are transferred to executable constructs in a structured way. Autoref is centered around a collection of so-called transfer rules. Transfer rules relate abstract with concrete operations. A transfer rule involves a transfer relation R::(γ × α) set , which relates a concrete implementation c::γ with an abstract element a::α and is of the following form. (c::γ, a::α) ∈ R Transfer rules are used to structurally synthesize concrete algorithms from abstract ones. Relations and relators (which combine relations) are used to express the relationship between concrete and abstract types. br is used to build a relation from an abstraction function a::γ → α andaninvariant I on the concrete type. br aI:={(c,ac) | Ic} Natural Relators For the types of functions, products, sets, or data types like lists and nres , one uses the natural relators A → B, A × B, Aset , Alist-rel , Anres with relations r r r r A, B for the argument types. ( f, f ) ∈ A → B ⇐⇒ ∀(x , y) ∈ A.(fx , f x ) ∈ B ((a, b), (a , b )) ∈ A × B ⇐⇒ (a, a ) ∈ A ∧ (b, b ) ∈ B (X, X ) ∈Aset ⇐⇒ (∀x ∈ X. ∃x ∈ X .(x , x ) ∈ A) ∧ (∀x ∈ X . ∃x ∈ X.(x , x ) ∈ A) (xs, xs ) ∈Alist-rel ⇐⇒ length xs = length xs ∧ (∀i < length xs.(xs , xs ) ∈ A) (RES X, RES X ) ∈Anres ⇐⇒ (X, X ) ∈Aset r r Representing Vectors We represent vectors (an arbitrary type α of class Euclidean space) as lists of real numbers where the length matches the dimension of the Euclidean space. lv-rel := br eucl-of-list (λxs. len xs = DI M (α)) This way, the concrete algorithm is monomorphic, which has the advantage that it can be generated once and for all and can therefore be used as a stand-alone tool. 123 A Verified ODE Solver and the Lorenz Attractor 95 Representing Enclosures We provide several implementations for the sets that can be used as enclosures. Intervals are represented by pairs of element types (which, in turn are implemented via some relation A): Aivl :={((a , b ), [a, b]) | (a , a) ∈ A ∧ (b , b) ∈ A} Zonotopes are represented using the joint range joint-range of affine forms affine := br (λA. eucl-of-list (joint-range A)) (λ_. True ) We use a symbolic representation of planes using the data type constructor Sctn that keeps normal vector n and translation c of a hyperplane. It is interpreted using plane-of (Sctn nc):={x |x , n= c} halfspace (Sctn nc):={x |x , n≤ c} for the hyperplane itself or for the halfspace below the hyperplane, where x , n is the inner product (also called dot product). Asctn is the natural relator that allows one to change the representation of the normal vector. With this, we can give a concrete implementation for hyperplanes and half-spaces. Aplane :=Asctn ◦ br plane-of (λ_. True ) r r Ahalfspace :=Asctn ◦ br halfspace (λ_. True ) r r For those relations, plane-of , halfspace,and ∩ are easily implemented with the identity function or as pair. On the abstract level, they describe useful objects that are convenient to reason about. ((λx . x ), plane-of ) ∈Asctn → Aplane r r r ((λx . x ), halfspace ) ∈Asctn → Ahalfspace r r r ((λxy.(x , y)), ∩) ∈ A → B → A, Binter r r r We will see that in some algorithms, one maintains a collection of enclosures, but abstractly one likes to see them as just one enclosure. For a relation A : (β × α set ) set that implements single enclosures for sets of type α with some concrete representation of type β, and a relation S : (σ × β set ) set that implements sets of concrete elements β, we define a relation that represents the union of all those elements as follows: S, AUnion : (σ × α set ) set S, AUnion := S ◦Aset ◦ br λX. x (λ_. True ) r r x ∈X Currently, we use lists to implement the set of concrete representations S, for which we write AUnion :=list-set , AUnion , and operations like union or extracting one element lr r r (with the specification split-spec ) can be implemented with the respective operations on lists/sets: (λxs ys. return (xs@ys), ∪) ∈AUnion → AUnion → AUnion lr r lr r lr (λx . return (hd x , tl x ), split-spec ) ∈AUnion → A × AUnion nres = lr r r lr r 123 96 F. Immler Relations to Guide Heuristics Often, in particular to guide heuristics, an algorithm needs to carry around information, which does not influence correctness proofs. An ODE solver for example, modifies its step size, also based on previous values. An implementation needs to carry this information around, but for verifying the algorithm, this only introduces unneces- sary clutter. We therefore introduce a relation that carries more information (implemented via A) in the implementation, but keeps the abstract semantics (implemented via B): A, Binfo :={((a , b ), b) |∃a.(a , a) ∈ A ∧ (b , b) ∈ B} Adding information is simply done by using a pair in the implementation side. Semantically, this information is simply discarded (put-info ab:= b). Information can be extracted with get-info , which is semantically just an arbitrary element (get-info b:= spec (λ_. True )). The implementations are straightforward: (λab.(a, b), put-info ) ∈ A → B → A, Binfo r r r ((λ(a, b).return a), get-info ) ∈A, Binfo →Anres r r An example of its usage is illustrated later on in Algorithm 1. 6 Reachability Analysis Overall, we design an algorithm that computes a Poincaré map with a list of intermediate Poincaré sections. The global idea (illustrated in Fig. 9) is as follows: starting from a set X , perform a series of single time discretization steps. If reachable sets grow above a given threshold, subdivide them (Sects. 6.3 and 6.4). Stop before an intermediate (or the final) section would be hit, then resolve the Poincaré map at that section (Sect. 6.5). For Tucker’s proof, it is important to also track the matrix of partial derivatives together with the solution. To this end, one can encode the derivative as a higher-dimensional ODE and use essentially the same algorithms as before. This instrumentation is presented in Sect. 6.7. 6.1 The Framework We use the high-level constructors and abstract specifications from the previous Sect. 5. We remain agnostic about the type of enclosures, for which we assume a relation encl and implementations for the abstract operations that are needed for the reachability analysis algo- rithms: an approximation scheme for expressions approx-slp-spec , enclosures from intervals using an implementation encl-of-ivl , lower and upper bounds with Inf-spec , Sup-spec,inter- sections with a plane inter-spec (note that the relation fixes the second argument to represent a plane, abstractly inter-spec is just intersection on sets): – (approx-slp , approx-slp-spec ) ∈ slp → encl → encl option nres r r r r r r r – (λxy. encl-of-ivl xy,λxy. [x , y]) ∈ lv-rel → lv-rel → encl r r r – (inf-encl , Inf-spec ) ∈ encl → lv-rel nres r r r – (sup-encl , Sup-spec ) ∈ encl → lv-rel nres r r r – (split-encl , split-spec ) ∈ real → nat el → encl → encl × encl nres ⊆ r r r r r r r r r r – (inter-encl-plane , inter-spec ) ∈ encl → lv-rel plane → encl nres 2 r r r r r r Currently, the only instantiation of this scheme is with affine arithmetic (in this case we set encl to affine ). Nevertheless, this structure keeps the formalization modular and one r r can imagine to add further instantiations—with e.g., Taylor models or centered forms—in the future. 123 A Verified ODE Solver and the Lorenz Attractor 97 Fig. 9 Continuous reachability and intermediate Poincaré sections 6.2 The Specification Our algorithms are supposed to compute enclosures for solutions of the ODE. We formalize the enclosure of an evolution from an initial set X to some other set Y with the ternary predicate ,where X  Y holds if the evolution flows every point of X ⊆ R to some point in Y ⊆ R and does not leave the set C in the meantime. We call C the flowpipe from X to Y . Definition 9 (Flows-to Predicate) X  Y:=∀x ∈ X. ∃t ≥ 0.φ(x , t ) ∈ Y ∧ (∀0 ≤ s ≤ t.φ(x , s) ∈ C ) C 0 0 X 6.3 Single Step In order to compute enclosures for a single step, one needs to first certify that a solution exists, which is the case for an initial value x and stepsize h if the iteration given by the Picard iteration from Sect. 2.1.2 has a unique fixed point. This is the standard approach from Bouissou et al. [4], also described in the setting of Isabelle [18]. The idea is that the expression Q (X ) = X +[0, h]· f (X ), which we can evaluate using approx-slp-spec , h 0 is an overapproximation of the Picard iteration and a post-fixed point certifies existence and a crude enclosure for solutions up to time h. This crude enclosure can be used as an overapproximation for the terms x (hs + t ) in the Runge–Kutta approximation scheme from Sect. 3.2. The function rk-step implements this and actually evaluates the Runge–Kutta approximation scheme twice: once for time h and once for the time interval [0, h], because this gives a much better enclosure for the flowpipe up to time h than the crude overapproximation from the Picard iteration. We prove the following specification. Theorem 10 rk-step Xh ≤ spec (λ(ε, C, Y ). X  Y ) The returned value  is an estimate for the approximation error. This is used for an adaptive step size control. Algorithm 1 shows an example how to use this heuristic (and another 123 98 F. Immler heuristic to split large sets), while (almost trivially, because the additional operations are either vacuous, the identity or overapproximations) satisfying the same specification. Theorem 11 single-step Xh ≤ spec (λ(C, Y ).X  Y ) The information on the last (and next) step size is only reflected in the refinement relation for the implementation of single-step : (single-step , single-step ) ∈ impl real , encl info → real , encl info Union , encl Union nres r r r r r r lr r lr r But this information does not clutter the verification of single-step or the statement of The- orem 11, which is very convenient. Algorithm 1 Single Step n n n 1: function single-step (X)  single-step : R set → (R set × R set ) 2: w ← width-spec X  semantically, this is spec (λ_. True ) 3: h ← get-stepsize X 4: if w ≤ max-width then  global parameter max-width 5: (ε, Y, C ) ← rk-step Xh 6: h ← adapt-step-size h ε  spec (λ_. True ) 7: return (put-info h Y, C ) 8: else 9: (Y, Z ) ← split-spec X 10: return (put-info h (Y ∪ Z ), (Y ∪ Z )) 6.4 Continuous Reachability Algorithm 2 Continuous Reachability Loop 1: function reach-cont (sctns, X ) 2: X ← X ; C ←∅; I ←∅ 3: while λ(X, C, I ). X =∅ do 4: (X , X ) ← split-spec X 1 2 5: (Y , C ) ← single-step X 1 1 1 6: d ← disjoint-spec C sctns 7: if d then 8: X ← X ∪ Y 2 1 9: C ← C ∪ C 10: I ← I 11: else 12: X ← X 13: C ← C 14: I ← I ∪ X 15: return (C, I ) Note that single-step returns (from an implementation point of view) a collection of enclosures, so we need some sort of work-list algorithm to resolve all currently reachable sets. Algorithm 2 does so. It maintains three kinds of sets (see also Fig. 9): X is the collection of currently “live” sets. C is the collection of all flowpipes explored so far. I is the collection 123 A Verified ODE Solver and the Lorenz Attractor 99 of sets where reachability analysis has stopped because of an intersection with a Poincaré section from sctns. The algorithm takes one element out of the “work-list” X by splitting the collection of enclosures using split-spec , performs a single step, checks for an intersection with one of the Poincaré sections and updates X, C, and I accordingly. The loop invariant of reach-cont is roughly the following: Elements from X flow via C to X and I , while avoiding sctns. X  (X ∪ I ) ∧ C ∩ sctns =∅ 0 C The specification of reach-cont follows immediately: Theorem 12 reach-cont (sctns, X ) ≤ spec (λ(C, I ). X  I ∧ C ∩ sctns =∅) 0 0 C It is worth noticing that the simplicity of the statement of this correctness theorem is due to the fact that the work-list and heuristic info is hidden via refinement relations in the implementation. If this were represented on the specification level (e.g., by using sets of enclosures paired with their current step size), the specification would have to be much more cluttered: ⎛ ⎞ ⎛ ⎛ ⎞ ⎞ ⎝ ⎠ ⎝ ⎝ ⎠ ⎠ x  x ∪ I (h,x )∈X (h,x )∈X Such a specification distracts the user and also automatic proof tools—we are therefore happy to hide this in the abstraction. 6.5 Resolve Intersection Algorithm 2 performs reachability analysis until each enclosure is just about to intersect an hyperplane. We compute the intersection from an enclosure X with a hyperplane H again in an iteration: continuous reachability steps are repeated as long as flowpipes intersect the hyperplane. The intersection of the individual flowpipes is computed with the geometric algorithm from Sect. 4. When the intersection is computed by flowing the reachable set through the hyperplane step by step, we get a set I consisting of individual intersections I . Many of the sets I i i usually overlap, in order to avoid redundant enclosures, the overlap is resolved with an overapproximative operation: all of the sets I are covered with an interval, which is repeatedly subdivided and shrinked to a given precision (see also section 3.6 in [21] for a more detailed discussion). 6.6 Intermediate Poincaré Maps The overall algorithm repeats the alternation of continuous reachability and resolution of Poincaré sections. When there is a sequence H ,..., H of intermediate Poincaré maps to 1 n be computed, it is important to ensure that, while flowing towards or resolving section H , one must not intersect with any of the H with j > i. Otherwise the later computation of the Poincaré map H might be incorrect because the actual first return time was reached before. 6.7 Derivatives For Tucker’s proof, it is necessary to compute not only the Poincaré map, but also its derivative. The derivative of the flow can be encoded as a higher dimensional ODE according to the variational equation (Theorem 7). 123 100 F. Immler n n n n∗n For an ODE with right hand side f : R → R , a new ODE of type R × R with right hand side (x , W ) → (fx , D f | · W ) is constructed. Here the first component contains the solution, and the second component its matrix of derivatives. We first extend the flows-to predicate  to a predicate  which also takes derivatives into account. Definition 10 (Flows-to Predicate Extended with Derivative) X  Y:=∀(x , d) ∈ X. ∃t ≥ 0.(φ(x , t ), Dφ| · d) ∈ Y ∧ 0 (x ,t ) bl C 0 (∀0 ≤ s ≤ t.(φ(x , s), Dφ| · d) ∈ C ) 0 (x ,s) bl X With this extended predicate for reachability, we can show that reach-cont , i.e., reach-cont for the extended ODE satisfies the specification reach-cont sctns X ≤ (λ(C, Y ). X  Y ). 0 C The Poincaré map, however requires extra care, because we cannot simply intersect the derivative of the flow with the Poincaré section: the derivative of the Poincaré map is given according to the expression in Theorem 8. For a hyperplane H ={x |x , n= c},the derivative is given as follows (for x ∈{x |x , n= c}): Dφ| · d, n (x ,τ (x )) DP| · d = Dφ| · d − f (P(x )) ϕ(t ) (x ,τ (x )) f (P(x )), n We can evaluate this expression using affine arithmetic. But we need to be able to enclose all quantities that occur on the right hand side, in particular P(x ) = ϕ(x,τ(x )) and Dφ| ·d. (x ,τ (x )) But we can enclose those: assume a step in computing an intersection, i.e., X  Y.Let us assume for simplicity that (X ∪ Y ) ∩ H =∅ and X and Y are on opposite sides of the hyperplane. Then the intersection of the flowpipe C with the section H encloses the Poincaré map: P(X ) ={ϕ(x,τ(x )) | x ∈ X}⊆ C ∩ H. For an extended flow X   Y , n∗n this means {(φ (x,τ(x )), Dφ| · d) | (x , d) ∈ X }⊆ C ∩ H × R . Therefore (x ,τ (x )) bl both P(x ) = φ(x,τ(x )) and Dφ| · d are enclosed by the result of the intersection (x ,τ (x )) n∗n C ∩ H × R for which we can use the regular intersection algorithm from Sect. 4. 6.8 Correctness Theorem We call the main algorithm that we outlined in the beginning of this Sect. 6 poincare:It resolves a sequence of intermediate Poincaré maps (together with their derivative). It is veri- fied to compute guaranteed enclosures for Poincaré maps and their derivative. The algorithm poincare takes as arguments an initial set X : R set , and initial matrix of partial derivatives n∗n n DX : R set and a target Poincaré section Σ : R set . It is further parameterized by a list of intermediate Poincaré sections, but they are irrelevant for the final correctness theorem. We formally verify partial correctness: If the algorithm returns a result, then this result encloses n∗n the Poincaré map P(x ) and its derivative DP| · DX for every x ∈ X and DX : R set . Theorem 13 (Correctness of ODE solver with Poincaré maps) poincare XDX Σ = R ⇒ ∀x ∈ X.(P(x ), DP| · DX ) ∈ R 7 Application to Lorenz Attractor In this section, we present how the verified algorithm poincare of the previous section is used to certify Tucker’s computations. We show in particular how we formally prove the Theorems 1 and 2. It helps to recall the roles of the forward invariant set N, the cone field C and the expansion estimates E in Tucker’s proof, as outlined in Sect. 1.3. 123 A Verified ODE Solver and the Lorenz Attractor 101 7.1 The Input Data and its Interpretation It is not necessary to verify precisely the set N that Tucker used, but coming up with a forward invariant set is slightly more involved than certifying one. We therefore use the output of Tucker’s program as a starting point to set up the input for our ODE solver. Since any other forward invariant with suitable cone field and expansion estimates would do just, we are free to modify Tucker’s data slightly. The output of Tucker’s program is available online as a file containing 7258 lines. We preprocessed this file by merging the information of some of the lines and slightly coarsening some of the numerical bounds. The coarsening accounts for slight differences between Tucker’s and our approximations. This results in a list of 400 elements, which we call input-data and will be the basis for all further interpretations: Definition 11 (Input Data) input-data ::result list is a list of 400 elements of type result . datatype result = Result (invoke-nf : B) − + (angle : R)(angle : R) (expansion : R)(preexpansion : R) − + − + (gridx : Z)(gridx : Z)(gridy : Z)(gridy : Z) − − + + (retx : Z)(rety : Z)(retx : Z)(rety : Z) Elements res of type result are interpreted as initial rectangles as follows. The properties − + − + gridx , gridx , gridy ,and gridy encode a rectangle on the Poincaré section Σ (recall Fig. 2), which we denote by N (res). The union of all elements of input-data represents the upper branch N of the forward invariant set N. It is plotted in Fig. 10. Definition 12 − −8 − −8 N (res):=[((gridx res − 1) · 2 ,(gridy res − 1) · 2 , 27), + −8 + −8 ((gridx res + 1) · 2 ,(gridy res + 1) · 2 , 27)] N := N (res) res∈input-data − + N :={(−x , −y, z) | (x , y, z) ∈ (N )} + − N:= N ∪ N The input data also contains information on the image of an initial rectangle. It is encoded − − + + in retx , rety , retx , rety : We select the elements within those bounds with return-of : return-of res:={res ∈ input-data | −  − + gridx res ∈[retx res, retx res]∧ −  − + gridy res ∈[rety res, rety res]} − + angle and angle define the cone C associated with the rectangle: the conic hull of the line segment between the boundary vectors. http://www2.math.uu.se/~warwick/main/rodes/ResultFile. 123 102 F. Immler + − + + Fig. 10 N in gray (N the upper and N = S(N ) the lower branch) and enclosure of P(N ) in black. This is a subset of the Poincaré section Σ (as in Fig. 2) Definition 13 C res =cone hull (segment − − (cos(rad (angle res), sin(rad (angle res), 0) + + (cos(rad (angle res)), sin(rad (angle res), 0))))) x ·π There rad x = is the radian of the angle given in degrees, segment xy is the line segment {(1 − u) · a + u · b | u ∈[0, 1]},and cone hull S ={c · x | 0 ≤ c ∧ x ∈ S} the conic hull of a set S. The elements in input-data also encode a conefield C and expansion estimates as follows. results-at (x ) yields the set of result elements that cover a point x (the rectangles overlap at −1 the boundary). We need to respect this to ensure that C, E,and E are well defined. Definition 14 results-at (x ):={res ∈ input-data | x ∈ N (res)} C(x ):= C(res) res∈results-at (x ) E (x ):= min expansion (res) res∈results-at (x ) −1 E (x ):= min preexpansion (res) res∈results-at (x ) One last property is invoke-nf , which encodes if the numerical computations need to be interrupted and the results of the normal form need to be invoked. First, we define abstractly when this is necessary, namely on the stable manifold of the origin. That is, the set of all points, which tend to the origin in infinite time. We restrict our attention to the part of the stable manifold whose trajectories do not intersect Σ for positive time. 123 A Verified ODE Solver and the Lorenz Attractor 103 Definition 15 Γ := {x |[0, ∞] ⊆ ex-ivl x ∧ (∀t > 0.φ(x , t)/ ∈ Σ) ∧ (φ (x , t ) −→ (0, 0, 0))} t →∞ When invoke-nf is true, the computations will be interrupted once the reachable sets arrive at the small cube L =[−0.1, 0.1]×[−0.1, 0.1]×[−0.1, 0.1] inside which the normal form estimates are valid. In our computations, solutions are guaranteed to enter the cube L through a rectangle T and the tangent vectors are in the cone that contains DT : ⎛ ⎞ [0.8, 1.7] ⎝ ⎠ T:= ([−0.1, 0.1], [−0.00015, 0.00015], 0.1) × [0.0005, 0.002] That is, sets are very slim in the y-direction, and the expanding direction is closely around the x axis. From Tucker’s analysis ([48, Proposition 3.1]), we devised the following bounds for the sets E , E2 (and corresponding cones in DE , DE ) through which solutions emanating 1 1 2 from T exit the cube L: ⎛ ⎞ ⎝ ⎠ E := ([−0.12, −0.088], [−0.024, 0.024], [−0.012, 0.13]) × [−0.56, 0.56] [−0.6, −0.08] ⎛ ⎞ ⎝ ⎠ E := ([0.088, 0.12], [−0.024, 0.024], [−0.012, 0.13]) × [−0.56, 0.56] [0.08, 0.6] When we interrupt computations close to L, we check that the sets entering L do so within T and continuous computations from E ∪ E . Since we have not verified Tucker’s normal 1 2 form theory, we need to trust the following assumption: Assumption 14 (Normal Form Theory Bounds) T  (E ∪ E ) 1 2 7.2 Checking the Input Data In the previous section, we only defined what the input-data encodes. Now we check if the numerical bounds prescribed by the input-data are actually correct. This involves three steps: First, we need to find a suitable setup to be able to use the algorithm poincare,which computes derivatives and not cones. Second, we set up the check that a single element of the input-data is correct. Third, we check all elements of the input-data , from which we conclude the formal counterparts of Theorems 1 and 2. 7.2.1 Representation of Cones Concerning the checking of cone conditions, first note that C res is an infinite cone, i.e., an unbounded set of vectors. In contrast to that, all of our numerical algorithms are tailored towards bounded enclosures. We therefore perform the computations with the line segment connecting the two tangent vectors with unit length. matrix-segment x y x y e encodes 1 1 2 2 a line segment (parameterized by e) in a matrix (such that it can be used as matrix initial condition DX of poincare , compare Theorem 13). mat-seg-of-deg uses this to define the line segment between the endpoints of unit vectors with given angles u,v to the x axis. A cone can therefore represented with the help of mat-seg-of-deg : 123 104 F. Immler Lemma 2 (Matrix Representation of Cone) ⎧ ⎛ ⎞ ⎫ ⎨ m ⎬ (1,1) − + ⎝ ⎠ C(res) = cone hull m m ∈ mat-seg-of-deg (angle res)(angle res) (2,1) ⎩ ⎭ with ⎛ ⎞ x + e · (x − x ) 00 1 2 1 ⎝ ⎠ matrix-segment x y x y2 e:= y + e · (y − y ) 00 1 1 2 1 2 1 mat-seg-of-deg u v:= matrix-segment (cos (rad u))(sin (rad u))(cos (rad v))(sin (rad v))[0, 1] 7.2.2 Checking a Single Result Element Algorithm 3 Check Result 1: function check-line-c1 (res) 2: X ← N (res) − + 3: DX ← mat-seg-of-deg (angle res)(angle res) 4: RES ← poincare X DX Σ 0 0 5: (P × DP ) ← split-along NRES i i i − − + + 6: RET ← get-results (retx res, rety res)(retx res, rety res) 7: return ∀i. ∃ret ∈ RET . returns-within res X DX ret i i Algorithm 3 outlines how to check that a single result element res ∈ input-data represents correct numerical bounds. It works as follows: X is the initial rectangle, DX the initial 0 0 data for the derivatives, which encodes the associated cone with angles angle res and angle res. Then the ODE solver returns with a union of return images RES, which are split along the boundaries of the individual rectangles making up N. This splitting ensures that each individual element (X , DX ) resulting from the splitting is contained in exactly i i on individual element of N. We write singleton parts of the result of this splitting X , DX . i i In RET , there are all elements of the input-data within which res is specified to return. The final check makes sure that every part X , DX of the splitting returns within one element i i ret of the collection RET . It is defined as follows and precisely formulates that X and DX, which emanate from a result res and hit the result ret, satisfy the prescribed bounds on cones and expansion. returns-within res X DX ret:= X ⊆ N (ret ) ∧ − + check-cone-bounds (angle res)(angle res) XDX ∧ −1 DX≥ E (res) ∧DX≥ E (ret ) check-cone-bounds is checked using affine arithmetic: It assumes that u and u are on the x y line segment encoding a cone according to mat-seg-of-deg , therefore checks that u = 0and ignores the other entries of the argument matrix. It further checks that the segment is on the right side (0 < u ) and that the boundary angles L and U (given in degrees) also represent a 123 A Verified ODE Solver and the Lorenz Attractor 105 cone pointing to the right side. The main purpose is in the last line, the check that the angle of the vector (u , u ) with the horizontal axis is between L and U. x y ⎛ ⎞ ⎛ ⎞ x u v w x x x ⎝ ⎠ ⎝ ⎠ check-cone-bounds LU y u v w := y y y z u v w z z z − 90 < L ∧ L ≤ U ∧ U < 90 ∧ 0 < u ∧ u = 0 ∧ x z u u y y tan(rad L) ≤ ∧ ≤ tan(rad U ) u u x x Correctness of check-line-c1 states that the set N (res) is mapped into the part return-of res of the forward invariant set. Vectors in the cone C(res) are mapped by the derivative DP into the cone field with the prescribed expansion estimates. The theorem states that the derivative exists and is defined when approaching x within Σ ={(x , y, z) | z ≤ 27}. Theorem 15 (Correctness of check-line-c1) check-line-c1 (res) = return True ⇒ ∀x ∈ N (res) − Γ. ∀dx ∈ C(res). returns-to Σ x ∧ P(x ) ∈ N (return-of res) ∧ (∃DP.(P has-derivative DP)(at x within Σ ) (DP(dx )≥ E (res) ·dx ) ∧ (∃ret ∈ return-of res. −1 P(x ) ∈ N (ret ) ∧ DP(dx ) ∈ C(ret ) ∧DP(dx )≥ E (ret ) ·dx )) The theorem follows rather directly from the definition of algorithm 3 and the specifications and definitions of the occurring functions. 7.2.3 Checking All Results We have indeed the theorem that all input-data is correct: Theorem 16 (Global Numerical Results) ∀res ∈ input-data . check-line-c1 res = return True We prove formally that under the Assumption 14, Theorem 16 implies Theorems 1 and 2, which is the main result of this article. It follows from combining the individual instances of Theorem 15 in a suitable way. Theorem 16 is proved by computing check-line-c1 (res) for every res ∈ input-data.The computations are carried out using by evaluating the statement Parallel.forall (λres. check-line-c1 res) input-data with Isabelle/HOL’s evaluation engine eval. Parallel.forall results in parallel processing of the 400 individual elements of input-data . Further parallelism is introduced when enclosures are split during reachability analysis. Split sets can be processed in parallel until they reach the next (intermediate) Poincaré section, where they might be (partially merged) upon resolving the intersection (Sect. 6.5). Figure 11 shows the plot of an enclosure for the Lorenz attractor resulting from the verified computation. The plot hints at the intermediate Poincaré sections that were manually set up 123 106 F. Immler (for some initial rectangles) at about z = 27, z = 30, x =±5, x =±1.5, x =±1, x =±0.75, x =±0.1, and z = 0.1. The black part of Fig. 10 is an enclosure for P(N ) resulting from these computations, and it is as expected and verified contained in N. The timing results of a computation on a machine with 22 cores are given below: – Threads: 22 – Elapsed Time: 6 h 33 min 9 s – CPU Time: 131 h 52 min 40 s – Parallelization Factor: 20.13 – Garbage Collection Time: 42 min 36 s To compare this with Tucker’s C++ program, I compiled Tucker’s program in a Vir- tual Machine running Ubuntu 4.20 and gcc version 3.3.4 on a machine with a 2,6 GHz Intel® Core™ i7 CPU and 16 GB RAM. Tucker’s program finished after a total CPU time of 30h and 24min. The algorithms and data structures are very different, so a direct comparison is hard. But with regard to the total CPU time (131 h) of our algorithm, a factor of less than five compared to a C++ program signifies reasonable performance for a verified algorithm. In earlier developments [20], an enclosure for the Lorenz attractor was computed with neither derivative nor cones. This earlier version verified an enclosure for the Lorenz attractor in about 7000 CPU hours. With the present version algorithms, such a computation (without derivatives and cones) can be performed in about 3 CPU hours. The speedup compared to the earlier work is mostly due to less aggressive splitting of reachable sets, and a smaller number of intermediate Poincaré sections: In the earlier work [20], intermediate Poincaré sections were introduced heuristically on-the-fly, and in the present work only where they are really effective. This is beneficial, because resolving the intersection incurs overapproximations. 8 Conclusion This article presented an overview of the diversity of algorithms, abstract results and tech- niques for the formal verification of a general-purpose ODE solver and its application to the rigorous numerical part of Tucker’s proof. 8.1 Future Work Future work would be to lift the assumption 14 on the normal form theory to formal founda- tions. This involves in particular multivariate formal power series and a number of analytic estimations, proving existence and convergence of specifically constructed formal power series. This involves more than 24 pages of Tucker’s article and is therefore a larger formal- ization effort. It also includes computer assisted parts: smalldiv.cc and coeffs.cc help devising the normal form, but neither their specification nor their implementation is particularly challenging; they essentially only evaluate a large number of fixed arithmetic expressions. Another direction for future work would be to formally conclude from the numerical results that the Lorenz equations support a robust singular hyperbolic attractor. One part would be to prove (similar to Tucker’s expansion.cc) that the expansion estimates given by E are sufficiently large to guarantee the locally eventually onto property. We validated this with a non-verified computation for the expansion estimates prescribed by our input-data.Further Intel Xeon CPU E5-2699 v4 @ 2.20 GHz. 123 A Verified ODE Solver and the Lorenz Attractor 107 + + Fig. 11 Enclosure of φ{(x , [0; τ(x )]) | x ∈ N }, during the computation of P(N ) mathematical foundations are required in order to conclude from the computed forward invariant conefield C and the expansion bounds that there exists a stable foliation, and that this foliation can be used to reduce the two-dimensional dynamics on Σ to a one-dimensional map. This requires in particular the formalization of differentiable manifolds, and theorems like the existence of differentiable invariant sections for fiber contractions [41]. 8.2 Size of Code Table 1 shows some statistics on the size in terms of lines of code of several programs related to this verification. RODES is the rigorous ODE solver used by Tucker, it consists of 3800 lines of C++ code and builds on a library for interval arithmetic (Profil/BIAS) of about twice the size. Similar to the sum of those two is the size of the generated SML code. The verification required more effort, but the largest part is generic: the part specific to the Lorenz attractor makes up less than 10% of the total number of lines of code. 8.3 Trust Base We use the evaluation oracle eval in Isabelle/HOL. This is common practice to speed up rewriting. Isabelle/HOL equations are mapped to function definitions in Standard ML. These are compiled and evaluated in PolyML. We also use the common extension http://polyml.org/. 123 108 F. Immler Table 1 Size of code and formalization Sections Language Lines of code/proof RODES C++ 3800 Profil/BIAS C++ 8852 generated ODE solver SML 13200 Flow, Poincaré map 2 Isabelle/HOL theory 12000–16500 Affine Arithmetic 3 8500 Intersection 4 5000 Refinement/Enclosures 3 5000 Reachability Analysis 6 10000 Lorenz Attractor 7 3000 (HOL-Library.Code_Target_Numeral in Isabelle2017) of the code generator that maps Isabelle/HOL integers to the integer type IntInf.int of PolyML, which can be based on either PolyML’s own implementation of arbitrary precision integers or GMP [8]. This setup means that the trusted code base is increased: The translation of Isabelle/HOL terms to SML code is not verified. One needs to trust PolyML and its compiler, but PolyML is Isabelle’s implementation language and therefore anyhow part of the trusted code base. Reducing the trusted code base is an orthogonal issue: there is ongoing work [16,17]for verified code generation to CakeML [28], a verified implementation of ML. Isabelle/HOL’s eval speeds up evaluation by translating terms to the implementation language of the proof checker (PolyML). In view of this, it is more similar to Coq’s native_compute [2], which evaluates terms after translation to Coq’s implementation language OCaml, than to Coq’s virtual machine [9]. 8.4 Related Work Integrals and Differential Equations in Proof Assistants Spitters and Makarov [33]imple- ment Picard iteration to calculate solutions of ODEs in the interactive theorem prover Coq, but restricted to relatively short existence intervals. Boldo et al. [3] approximate the solution of one particular partial differential equation with a C-program and verify its correctness in Coq. Mahboubi and Sibut-Pinote [32] compute rigorous approximations of integrals with Taylor models. Rigorous Numerics in Proof Assistants Rigorous numerical approximation of arithmetic expressions has been done in Coq [34] for different types of enclosures (Taylor models, inter- vals, centered forms). Muñoz and Lester [37] use rational interval arithmetic in PVS. Rigorous numerics with first order interval approximations has been implemented by Solovyev for the Flyspeck project [44] in HOL Light. This work is remarkable in that it is not relying on code generation but uses only primitive inference rules of HOL Light’s kernel. Computational Geometry Pichardie and Bertot [40] were the first to formalize the ccw system and verify a functional convex hull algorithm in Coq. Meikle and Fleuriot [35] formalized an imperative algorithm and verified it using Hoare logic in Isabelle/HOL. Brun et al. [5]verify an algorithm based on hypermaps to compute the convex hull. 123 A Verified ODE Solver and the Lorenz Attractor 109 Acknowledgements I would like to thank Jeremy Avigad for drawing my attention to this particular appli- cation of rigorous ODE solving. I am very grateful to Tobias Nipkow for encouraging and supporting me to pursue this multi-faceted project. The anonymous reviewers provided helpful and constructive feedback. This research was financially supported by DFG RTG 1480 (PUMA) and DFG Koselleck Grant NI 491/16-1. I would like to thank Johannes Hölzl for supervising my student project on the Picard–Lindelöf theorem. I thankfully acknowledge Christoph Traut’s work for his student project on formalizing the variational equation. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. Back, R.J., Wright, J.: Refinement Calculus: A Systematic Introduction. Springer, New York (2012) 2. Boespflug, M., Dénès, M., Grégoire, B.: Full reduction at full throttle. In: International Conference on Certified Programs and Proofs, pp. 362–377. Springer, New York (2011) 3. Boldo, S., Clément, F., Filliâtre, J.C., Mayero, M., Melquiond, G., Weis, P.: Wave equation numerical resolution: a comprehensive mechanized proof of a C program. Journal of Automated Reasoning 50(4), 423–456 (2013). https://doi.org/10.1007/s10817-012-9255-4 4. Bouissou, O., Chapoutot, A., Djoudi, A.: Enclosing temporal evolution of dynamical systems using numerical methods. In: Brat, G., Rungta, N., Venet, A. (eds.) NASA Formal Methods, LNCS, vol. 7871, pp. 108–123. Springer, New York (2013). https://doi.org/10.1007/978-3-642-38088-4_8 5. Brun, C., Dufourd, J.F., Magaud, N.: Designing and proving correct a convex hull algorithm with hyper- maps in Coq. Computational Geometry 45(8), 436–457 (2012). https://doi.org/10.1016/j.comgeo.2010. 06.006.(Geometric Constraints and Reasoning) 6. de Figueiredo, L., Stolfi, J.: Affine arithmetic: concepts and applications. Numerical Algorithms 37(1–4), 147–158 (2004) 7. Girard, A., Le Guernic, C.: Zonotope/hyperplane intersection for hybrid systems reachability analysis. In: Egerstedt, M., Mishra, B. (eds.) Hybrid Systems: Computation and Control, LNCS, vol. 4981, pp. 215–228. Springer, New York (2008). https://doi.org/10.1007/978-3-540-78929-1_16 8. Granlund, T., et al.: GNU MP 6.0 Multiple Precision Arithmetic Library. Samurai Media Limited, Hong Kong (2015) 9. Grégoire, B., Leroy, X.: A compiled implementation of strong reduction. ACM SIGPLAN Not. 37(9), 235–246 (2002) 10. Haftmann, F., Krauss, A., Kuncar ˇ , O., Nipkow, T.: Data Refinement in Isabelle/HOL, pp. 100–115. Springer, Berlin (2013). https://doi.org/10.1007/978-3-642-39634-2_10 11. Haftmann, F., Wenzel, M.: Constructive type classes in Isabelle. In: International Workshop on Types for Proofs and Programs, pp. 160–174. Springer, New York (2006) 12. Harrison, J.: A HOL theory of Euclidean space. In: Hurd, J., Melham, T. (eds.) 18th International Confer- ence on Theorem Proving in Higher Order Logics (TPHOLs 2005), Lecture Notes in Computer Science, vol. 3603, pp. 114–129 (2005) 13. Hölzl, J.: Proving inequalities over reals with computation in Isabelle/HOL. In: Reis, G.D. Théry L. (eds.) Programming Languages for Mechanized Mathematics Systems (ACM SIGSAM PLMMS’09), pp. 38–45 (2009) 14. Hölzl, J., Immler, F., Huffman, B.: Type classes and filters for mathematical analysis in Isabelle/HOL. In: Blazy, S., Paulin-Mohring, C., Pichardie, D. (eds.) ITP 2013, pp. 279–294. Springer, New York (2013) 15. Huffman, B., Kuncar ˇ , O.: Lifting and transfer: a modular design for quotients in Isabelle/HOL. In: Inter- national Conference on Certified Programs and Proofs, pp. 131–146. Springer, New York (2013) 16. Hupel, L.: Dictionary Construction. Archive of Formal Proofs (2017). http://isa-afp.org/entries/Dict_ Construction.html, Formal proof development 17. Hupel, L., Nipkow, T.: A verified compiler from Isabelle/HOL to CakeML. Springer LNCS (2018). https:// lars.hupel.info/pub/isabelle-cakeml.pdf 18. Immler, F.: Formally verified computation of enclosures of solutions of ordinary differential equations. In: Badger, J.M., Rozier, K.Y. (eds.) NFM 2014, pp. 113–127. Springer, New York (2014) 123 110 F. Immler 19. Immler, F.: A verified algorithm for geometric zonotope/hyperplane intersection. In: CPP 2015, pp. 129– 136. ACM (2015) 20. Immler, F.: A verified enclosure for the Lorenz attractor (rough diamond). In: Urban, C., Zhang, X. (eds.) ITP 2015, pp. 221–226. Springer, Berlin (2015) 21. Immler, F.: Verified reachability analysis of continuous systems. In: Baier, C., Tinelli, C. (eds.) TACAS 2015, pp. 37–51. Springer, Berlin (2015) 22. Immler, F.: Affine Arithmetic. Archive of Formal Proofs (2017). http://isa-afp.org/entries/Affine_ Arithmetic.shtml, Formal proof development 23. Immler, F., Hölzl, J.: Ordinary Differential Equations. Archive of Formal Proofs (2017). http://isa-afp. org/entries/Ordinary_Differential_Equations.shtml, Formal proof development 24. Immler, F., Hölzl, J.: Numerical analysis of ordinary differential equations in Isabelle/HOL. In: Interna- tional Conference on Interactive Theorem Proving, pp. 377–392. Springer, New York (2012) 25. Immler, F., Traut, C.: The flow of ODEs: Formalization of variational equation and Poincaré map. J. Autom. Reason. (2018). https://doi.org/10.1007/s10817-018-9449-5 26. Immler, F., Traut, C.: The flow of ODEs. In: Blanchette, J.C., Merz, S. (eds.) ITP 2016, pp. 184–199. Springer, Cham (2016) 27. Knuth, D.: Axioms and Hulls. Springer, Berlin (1992). Number 606 in Lecture Notes in Computer Science 28. Kumar, R., Myreen, M.O., Norrish, M., Owens, S.: CakeML: a verified implementation of ML. In: ACM SIGPLAN Notices, vol. 49, pp. 179–191. ACM (2014) 29. Lammich, P.: Refinement for monadic programs. Archive of Formal Proofs (2012). http://afp.sf.net/ entries/Refine_Monadic.shtml, Formal proof development 30. Lammich, P.: Automatic data refinement. In: International Conference on Interactive Theorem Proving, pp. 84–99. Springer, Heidelberg (2013) 31. Lorenz, E.N.: Deterministic nonperiodic flow. J. Atmos. Sci. 20(2), 130–141 (1963). https://doi.org/10. 1175/1520-0469(1963)020<0130:DNF>2.0.CO;2 32. Mahboubi, A., Melquiond, G., Sibut-Pinote, T.: Formally Verified Approximations of Definite Integrals, pp. 274–289. Springer, Cham (2016). https://doi.org/10.1007/978-3-319-43144-4_17 33. Makarov, E., Spitters, B.: The Picard algorithm for ordinary differential equations in Coq. In: Blazy, S., Paulin-Mohring, C., Pichardie, D. (eds.) Interactive Theorem Proving, LNCS, vol. 7998, pp. 463–468. Springer, Berlin (2013) 34. Martin-Dorel, E., Melquiond, G.: Proving tight bounds on univariate expressions with elementary func- tions in Coq. J. Autom. Reason. 57(3), 187–217 (2016). https://doi.org/10.1007/s10817-015-9350-4 35. Meikle, L., Fleuriot, J.: Mechanical theorem proving in computational geometry. In: Hong, H., Wang, D. (eds.) Automated Deduction in Geometry, Lecture Notes in Computer Science, pp. 1–18. Springer, Berlin (2006). https://doi.org/10.1007/11615798_1 36. Morales, C., Pacifico, M., Pujals, E.: Singular hyperbolic systems. Proc. Am. Math. Soc. 127(11), 3393– 3401 (1999) 37. Muñoz, C., Lester, D.: Real number calculations and theorem proving. In: Hurd, J., Melham, T. (eds.) Theorem Proving in Higher Order Logics (TPHOLs 2005), LNCS, pp. 195–210 (2005) 38. Nipkow, T.: Order-sorted polymorphism in Isabelle. In: Logical Environments, pp. 164–188. Cambridge University Press (1993) 39. Nipkow, T., Paulson, L.C., Wenzel, M.: Isabelle/HOL: A Proof Assistant for Higher-Order Logic, LNCS. Springer, Heidelberg (2002) 40. Pichardie, D., Bertot, Y.: Formalizing convex hull algorithms. In: Boulton, R., Jackson, P. (eds.) Theorem Proving in Higher Order Logics, Lecture Notes in Computer Science, vol. 2152, pp. 346–361. Springer, Berlin (2001). https://doi.org/10.1007/3-540-44755-5_24 41. Robinson, C.: Dynamical Systems—Stability, Symbolic Dynamics, and Chaos. CRC Press, Boca Raton (1999). https://doi.org/10.1007/978-1-4613-0003-8 42. Rump, S.M., Kashiwagi, M.: Implementation and improvements of affine arithmetic. Nonlinear Theory Appl. IEICE 6(3), 341–359 (2015). https://doi.org/10.1587/nolta.6.341 43. Smale, S.: Mathematical problems for the next century. Math. Intell. 20(2), 7–15 (1998). https://doi.org/ 10.1007/BF03025291 44. Solovyev, A., Hales, T.C.: Formal verification of nonlinear inequalities with Taylor interval approxima- tions. In: NASA Formal Methods Symposium, pp. 383–397. Springer, New York (2013) 45. Sparrow, C.: The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors. No. 41 in Applied Mathematical Sciences. Springer. New York. https://doi.org/10.1007/978-1-4612-5767-7 46. Tucker, W.: My thesis: The lorenz attractor exists. http://www2.math.uu.se/~warwick/main/pre_thesis. html 47. Tucker, W.: The Lorenz attractor exists. Comptes Rendus de l’Académie des Sciences-Series I- Mathematics 328(12), 1197–1202 (1999) 123 A Verified ODE Solver and the Lorenz Attractor 111 48. Tucker, W.: A rigorous ODE solver and Smale’s 14th problem. Found. Comput. Math. 2(1), 53–117 (2002) 49. Viana, M.: What’s new on Lorenz strange attractors? Math. Intell. 22(3), 6–19 (2000) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Automated Reasoning Springer Journals

A Verified ODE Solver and the Lorenz Attractor

Free
39 pages
Loading next page...
 
/lp/springer_journal/a-verified-ode-solver-and-the-lorenz-attractor-Y0f8SSKOFs
Publisher
Springer Netherlands
Copyright
Copyright © 2018 by The Author(s)
Subject
Computer Science; Mathematical Logic and Formal Languages; Artificial Intelligence (incl. Robotics); Mathematical Logic and Foundations; Symbolic and Algebraic Manipulation
ISSN
0168-7433
eISSN
1573-0670
D.O.I.
10.1007/s10817-017-9448-y
Publisher site
See Article on Publisher Site

Abstract

J Autom Reasoning (2018) 61:73–111 https://doi.org/10.1007/s10817-017-9448-y A Verified ODE Solver and the Lorenz Attractor Fabian Immler Received: 19 April 2017 / Accepted: 27 December 2017 / Published online: 22 January 2018 © The Author(s) 2018. This article is an open access publication Abstract A rigorous numerical algorithm, formally verified with Isabelle/HOL, is used to certify the computations that Tucker used to prove chaos for the Lorenz attractor. The ver- ification is based on a formalization of a diverse variety of mathematics and algorithms. Formalized mathematics include ordinary differential equations and Poincaré maps. Algo- rithms include low level approximation schemes based on Runge–Kutta methods and affine arithmetic. On a high level, reachability analysis is guided by static hybridization and adap- tive step-size control and splitting. The algorithms are systematically refined towards an implementation that can be executed on Tucker’s original input data. Keywords Isabelle/HOL · Ordinary differential equation · Rigorous numerics · Poincaré map · Lorenz attractor 1 Introduction Computer assisted proofs, i.e., mathematical proofs that rely on the output of a computer program, depend crucially on the correctness of the program. Important computer assisted proofs for e.g., the Kepler conjecture or the Four Color Theorem, have therefore been formally verified. In this article, we consider the Lorenz attractor—perhaps one of the most prominent examples of deterministic chaos—and its computer assisted proof by Warwick Tucker. The proof relies on a rigorous numerical ordinary differential equation (ODE) solver. In this article, we describe the long-term project of formally verifying (in Isabelle/HOL [39]) an ODE solver that is capable of certifying Tucker’s computations. B Fabian Immler immler@in.tum.de Institut für Informatik,Technische Universität München, Munich, Germany 123 74 F. Immler 1.1 History In 1963, meteorologist Edward Lorenz [31] introduced a system of ODEs as a simplified model for atmospheric dynamics. He observed that even the smallest perturbation in initial values would lead to completely different long-term behavior of the system. Referring to the original motivation, he asked: “Does the Flap of a Butterfly’s Wings in Brazil Set Off a Tornado in Texas?” and the term butterfly effect entered popular culture. The Lorenz system tends to evolve to a complicated structure (Fig. 2), which became an iconic example of deterministic chaos: According to Sparrow [45] “the number of man, woman, and computer hours spent on [the Lorenz equations …] must be truly immense”. Despite its popularity and the amount of effort put into its study, nobody managed to prove that the Lorenz attractor is chaotic in a rigorous mathematical sense. The problem of rigorously proving chaos in the Lorenz attractor even made it into a list of 18 important problems for the twenty-first century that Field’s medalist Stephen Smale composed in 1998 [43]. Shortly after, Warwick Tucker managed to give an affirmative answer by presenting a computer-assisted proof [47,48]. Tucker’s programs were written in C++ and are not for- mally verified. Tucker even discovered (and fixed) some bugs in it [46,49]. Formal verification of the numerical results needed for the proof is therefore a worthwhile goal. 1.2 The Lorenz Attractor We start with describing the Lorenz attractor and some of the properties that were conjectured from numerical simulations. In his proof, Tucker considers the following three dimensional ODE for fixed parameters k ,λ : 1,2,3 1,2,3 x˙ = λ x − k (x + y)z 1 1 y ˙ = λ y + k (x + y)z 2 1 z ˙ = λ z + (x + y)(k x + k y) 3 2 3 As intuition, an ODE describes the velocity vector (x˙ , y ˙, z ˙) in which a particle at a point (x , y, z) moves. The evolution of a particle subject to the ODE is described by the so-called flow φ. A particle x ∈ R will be at position φ(x , t ) after time t ∈ R. 0 0 Figure 1 depicts the numerical simulation of the evolution of a particle starting at (0.1, 0, 0): It moves to right (x ≈ 15) and up (z ≈ 50) at time t ≈ 0.5, then down to about z = 27 and oscillates with increasingly larger amplitude around z = 27. Figure 2 depicts the trace of a long-term evolution in the three dimensional phase space, it indicates Property 1: Property 1 Solutions remain in a bounded region of the phase space. Particles that approach the origin (0, 0, 0) from above exhibit a very sensitive dependence on initial conditions: a slight perturbation can make the particle flow to either the left or right branch of the Lorenz attractor, which we call Property 2: Property 2 Solutions exhibit sensitive dependence on initial conditions. This dependence is such that arbitrarily small initial sets will eventually spread over the whole attractor. 1 8 That is, Lorenz’ original equations for the classical parameters β = ,σ = 10,ρ = 28 in Jordan normal σ σ −1+τ σ −1−τ −σ −1+τ form using τ:= (σ + 1) + 4σ(ρ − 1), k := , k := , k := , λ := , 1 2 3 1 τ 2σ 2σ 2 −σ −1−τ λ := ,and λ :=−β. 2 3 123 A Verified ODE Solver and the Lorenz Attractor 75 Fig. 1 Temporal evolution of t → φ((0.1, 0, 0), t ) Fig. 2 Simulation of a part of the Lorenz attractor (φ) and Poincaré section (Σ) 1.3 Tucker’s Proof How does Tucker go about proving those properties? First of all, he uses a standard technique: he introduces a so-called Poincaré section. This is a distinguished set in the phase space, in this case a square on the plane z = 27, namely Σ =[−6, 6]×[6, 6]×{27}. Compare also Fig. 2. On a Poincaré section Σ, one defines the so-called Poincaré map P: For a particle x ∈ Σ, the Poincaré map P(x ) is the point where the flow first returns to Σ. This reduces the three- dimensional, continuous dynamics φ to (discrete) iterations of the two-dimensional map P. Tucker then analyzes the dynamics of P. Trapping Region Regarding Property 1, Tucker proves that there is a (compact) trapping region N ⊆ Σ, such that solutions starting in N will remain in N. He does so by subdividing N into a large number of small rectangles. For every small rectangle, Tucker’s program computes safe numeric bounds for all solutions evolving from the small rectangle. In a number of time-discretization steps, the evolution is followed until it eventually returns to Σ. Upon return, the program checks that the returned enclosure is contained in N.Ifthis process succeeds for every small rectangle, one can conclude the following theorem. 123 76 F. Immler Theorem 1 (Trapping Region) ∀x ∈ N − Γ. P(x ) ∈ N Note that there exists a set Γ on which P is not defined: Γ is the set of points, from which solutions tend to the origin in infinite time. Γ is therefore explicitly excluded in the above theorem. Sensitive Dependence Regarding Property 2, sensitive dependence on initial conditions can be quantified with the help of the derivative: A deviation in the direction of a vector v ∈ R is propagated (in linear approximation) like the derivative at x, i.e., P(x + v) ≈ x + DP| · v. Here DP| · v is the matrix of partial derivatives of P (the Jacobian matrix) at the point x, multiplied with the vector v. A mathematically precise notion of chaos is given by the class of singular hyperbolic sys- tems [36]. The term singular denotes the special case where the system contains a hyperbolic fixed point (which renders P undefined on Γ ). A hyperbolic system contracts deviations in stable directions and expands deviations in unstable directions. Both are relevant for the dynamics of the attractor: Stable directions make solutions tend to the attractor, whereas unstable directions lead to sensitive dependence on initial conditions. Tucker proves that the Lorenz attractor is hyperbolic (in fact singular hyperbolic, we discuss how the hyperbolic fixed point is addressed with normal form theory in the next paragraph) by providing safe overapproximations for the unstable direction: Every x ∈ N is equipped with a cone C(x ) (compare Fig. 3), which contains the unstable direction. This is also verified by Tucker’s computer program: In addition to the Poincaré map, the program keeps bounds on its matrix of partial derivatives. The program tracks how initial deviations (inside the cone associated to an initial rectangle) are propagated by the derivative DP.The cone field needs to be forward invariant (otherwise it would not contain the unstable direction) and the expansion needs to be large enough that the enclosed directions are actually expand- −1 ing. Tucker’s program establishes factors E (x ) and E (x ), which quantify the expansion properties of P: Theorem 2 (Derivatives, Cones, and Expansion) 1. ∀x ∈ N − Γ. ∀v ∈ C(x ). DP| · v ∈ C(P(x )) 2. ∀x ∈ N − Γ. ∀v ∈ C(x ). DP| · v ≥ E (x ) |v| −1 3. ∀x ∈ N − Γ. ∀v ∈ C(x ). DP| · v ≥ E (P(x )) v This theorem states that 1., the cone field C is forward invariant under the action of the derivative of P: the image of every cone is slimmer than the cones onto which they are mapped. 2., the vectors v satisfy lower bounds on how much they are expanded: the length DP| · v of the return of the deviation vector v is lower bounded by its length v times an −1 expansion factor E (x ). They also satisfy a pre-expansion bound E (x ) (this does not denote ) for the pre-image of x, which is required for technical reasons in Tucker’s proof. Normal Form Theory For the Lorenz equations, the origin (0, 0, 0) is a hyperbolic fixed point. The origin is a fixed point, because the ODE evaluates to 0. It is hyperbolic, because solutions tend to it in the two stable directions given by the y and z axis and expands in the unstable direction given by the x-axis. This hyperbolic fixed point poses problems for the aforementioned approach of using rigorous numerical methods: there are solutions that tend to the origin as time goes to infin- ity. In such a situation, a time-discretization algorithm is at a loss, because it would need infinitely many steps. To remedy this problem, Tucker’s program interrupts computations in 123 A Verified ODE Solver and the Lorenz Attractor 77 Fig. 3 Enclosures for the flow and cones evolving from X =[4.375, 4.4]×[2.77, 2.79]×{27} with a ◦ ◦ representation of a cone between 1.5 and 11.5 (detail on the right) a small cube L =[−0.1, 0.1]×[−0.1, 0.1]×[−0.1, 0.1] around the origin. Inside the cube, where the numerical methods would fail, the evolution of solutions can be described with classical, analytical means: more than half of Tucker’s thesis is devoted to accurate analytical expressions for the flow inside the cube L. These expressions can be used to provide explicit bounds on how solutions exit the cube L and continue with numerical computations. Informal Theorem 3 There is an explicit form that bounds the dynamics inside the cube L =[−0.1, 0.1]×[−0.1, 0.1]×[−0.1, 0.1]. 1.4 Outline of This Article This article is about the formalization of a rigorous ODE solver that computes the Poincaré map P and its derivative DP. It is sufficiently efficient and precise to certify Tucker’s numer- ical results. In particular, computing with the verified algorithm proves Theorems 1 and 2. In fact, also Fig. 3 was created from the output of this verified algorithm. As a matter of course, since the ODE solver computes Poincaré maps and derivatives thereof, it is proved correct with regard to a formalization of these concepts in Isabelle/HOL. This formalization, as well as the verified algorithms are generic: they are independent of the underlying ODE or dimension. This article summarizes previous work of the author (and contributions by Johannes Hölzl and Christoph Traut) [18–26] and describes the state of the formalization to which it evolved over time. It shows how the various parts fit together for the final result of certifying Tucker’s computations. The following explains how this article is structured and details on the relation to earlier publications. The abstract mathematics needed for the formalization is a theory of ODEs and Poincaré maps together with their (differentiable) dependence on initial conditions. This is presented in Sect. 2, which summarizes a more comprehensive journal article [25], in which the author extends the conference paper [26] with a formalization of Poincaré maps. The foundations, in particular existence and uniqueness of solutions was proved in [24]. 123 78 F. Immler Tucker used the library Profil/BIAS for an implementation of the Euler method in interval arithmetic. Our approach to rigorous numerics is at first agnostic about the concrete type of enclosures (Sect. 3). The main instantiation, affine arithmetic, is presented in Sect. 3.4.Parts of this section were part of earlier publications [18,21]. When working with affine arithmetic, the enclosures are zonotopes (centrally symmetric polytopes) instead of intervals. Because zonotopes are a more complex than intervals, geo- metric operations like intersections with the Poincaré section Σ are more challenging. The verification leads to a digression into computational geometry in Sect. 4, which is based on an earlier publication [19]. Tucker’s program adaptively splits reachable sets and therefore maintains a collection of sets. Section 5 describes how we formalize generic data structures to maintain such collec- tions and use stepwise refinement from nondeterministic abstract specifications to a concrete deterministic implementation. This has not been published before. A more ad-hoc formalization of similar algorithms has been described in earlier work [21], which forms the basis of Sect. 6. Here the presentation is w.r.t. the new, generic framework and extended for derivatives of Poincaré maps. Section 7 is a novel contribution: it describes how the trapping region N, the cone field C −1 and the expansion estimates E , E are defined formally and how the verified ODE solver is set up to certify the results of all of Tucker’s computations. Earlier work [20] was not capable of handling derivatives, and had no formalization of Poincaré map. All of the development described here is available for Isabelle2017 in the Archive of Formal Proof [22,23]. In particular everything displayed as Theorem possesses a formal counterpart. 2Mathematics The required mathematical background consists of mostly standard results which can be found in every textbook on ODEs and dynamical systems. Thanks to a sufficient background theory, the formalization can mostly follow the presentations in such textbooks. We therefore focus on peculiarities of our formalization which are due to Isabelle/HOL and its type system, in particular the use of type classes and the lack of dependent types. Because of the type class based formalization of topological structures, type definitions are used to formalize function spaces, where the Transfer and Lifting tools [15] provide excellent support. Moreover, there are no dependent types in Isabelle/HOL. In situations where this would be more natural, an encoding (e.g., ext-cont in the following Sect. 2.1.2) is necessary. 2.1 Type Classes for Mathematics in Isabelle/HOL In Isabelle/HOL, many of the mathematical concepts (in particular spaces with a certain structure) are formalized using type classes. Isabelle/HOL features axiomatic type classes [11, 38]. The purpose of an axiomatic type class is to specify operations which satisfy given properties for a class of types. The advantage of type class based reasoning is that most of the reasoning is generic: formalizations are carried out in the context of type classes and can then be used for all types inhabiting that type class. For generic formalizations, we use Greek letters α, β, γ and name their type class con- straints in prose. I.e., when we write that we “consider a topological space” α, then this result is formalized generically for every type α that fulfills the properties of a topological space. 123 A Verified ODE Solver and the Lorenz Attractor 79 The spaces we consider are topological spaces with open sets, (real) vector spaces with addition +: α → α → α and scalar multiplication (_)(_) : R → α → α. Normed vector spaces come with a norm |(_)|: α → R. Complete normed vector spaces are called Banach spaces. Much of the theory has been ported from Harrison’s theory of Euclidean space [12] and has been generalized to the hierarchy of type classes for mathematics in Isablle/HOL [14]. 2.1.1 Vectors in Euclidean Space Because of Isabelle/HOL’s restrictive type system (no dependent types), the abstract concept of vectors is notorious for demanding workarounds. In Isabelle/HOL, one tends to use a type class based encoding. We work with a type class for Euclidean space that fixes an order on the Basis elements and therefore enables operations eucl-of-list : R list → α and list-of-eucl : α → R list if α is a Euclidean space. All (finite) vectors of real numbers are instances of the class Euclidean space. This includes real numbers R, complex numbers C, tuples α × β for Euclidean spaces α, β, and Harrison- 2 ι style [12] vectors α for a finite type ι. 2.1.2 Bounded Continuous Function We motivate bounded continuous functions with the Picard–Lindelöf theorem, which guar- antees the existence of a unique solution to an initial value problem. For an ODE f with initial value x at time t , a unique solution on the time interval [t , t ] is constructed by 0 0 0 1 considering iterations of the following operator for continuous functions φ :[t , t ]→ R : 0 1 P(φ):= λt. x + f (τ, φ (t )) dτ From a mathematician’s point of view, P operates on the Banach space of continuous func- tions on the compact domain [t , t ] and therefore the Banach fixed point theorem guarantees 0 1 the existence of a unique fixed point (which is by construction the unique solution). In order to formalize this in Isabelle/HOL, there are two obstructions to overcome: First, the concept of Banach space is a type class in Isabelle/HOL, so we need to introduce a type for the mappings φ :[t , t ]→ R from above. But this poses the second problem: functions 0 1 in Isabelle/HOL are total and types must not depend on term parameters like t and t . 0 1 We work around these restrictions by introducing a type of bounded continuous func- tions, which is a Banach space and comprises (with a suitable choice of representations) all continuous functions on all compact domains. n m typedef α → β:={ f : R → R | f continuous on α ∧ (∃B. ∀t. ft≤ B)} bc In order to define operations on type α → β, the Lifting and Transfer package [15]is bc an essential tool: operations on the plain function type α → β are automatically lifted to definitions on the type α → β when supplied with a proof that functions in the result are bc bounded continuous under the assumption that argument functions are bounded continuous. We write application $ of a bounded continuous function f : α → β to an element bc bc x : α as follows. Vectors of length n are represented by a type of functions ι → α,where n equals the cardinality of the finite type ι. 123 80 F. Immler Definition 1 (Application of Bounded Continuous Functions) ( f $ x ) : β bc Bounded continuous functions form a normed vector space. The norm on α → β is the bc supremum of the range and the vector space operations +, · are defined pointwise. Definition 2 (Normed Vector Space of Bounded Continuous Functions) f : = sup { f $ x| x ∈ α} bc ( f + g) $ x := f $ x + g $ x bc bc bc (a · f ) $ x := a · ( f $ x )) bc bc The type → with the above operations forms a complete normed vector space (a Banach bc space). This allows us to use the Banach fixed point theorem for operators on this type. In order to be able to use this for the operator P from above, we represent functions on a compact interval [a, b] as an element of type → by extending the function continuously bc outside the domain with the help of clamp : clamp x := if x ≤ a then a else (if x ≥ b then b else x ) [a,b] (ext-cont f ) $ x := f (clamp x )) [a,b] bc [a,b] n n With the help of ext-cont we can apply P : (R → R ) → (R → R ) to a continuous bc bc function φ : R → R (that is assumed to be continuous on an interval [a, b]) by writing P(ext-cont φ). According to the Banach fixed point theorem there exists a unique fixed [t ,t ] 0 1 point φ : R → R where P(φ ) = φ and the unique solution of the initial value bc bc bc bc problem is the function λt.φ $ t of type R → R . bc bc The usage of the type → caused minor technical obstructions, but otherwise enabled a bc natural and abstract proof. 2.1.3 Bounded Linear Functions Similar to the type of bounded continuous functions, we also introduce a type of bounded linear functions (also known as continuous linear functions) For vector spaces α and β, a linear function is a function f : α → β that is compatible with addition and scalar multiplication. linear f :=∀xy c. f (c · x + y) = c · f (x ) + f (y) Let us assume normed vector spaces α and β. Linear functions are continuous if the norm of the result is linearly bounded by the norm of the argument. We cast bounded linear functions α → β as a type α → β in order to make it an instance of Banach space. bl { } typedef α → β:= f : α → β | linear f ∧∃K . ∀x .  f (x )≤ K x bl The construction is very similar to bounded continuous functions and we write bounded linear function application ( f · x ). Vector space operations are also analogous to → .The bl bc usual choice of a norm for bounded linear functions is the operator norm: the maximum of the image of the bounded linear function on the unit ball. With this norm, α → β forms a bl normed vector space and we prove that it is Banach if α is a normed vector space and β is Banach. Definition 3 (Norm in Banach Space → )For f : α → β, bl bl f : = onorm (λy. f · y) = max { f · y|y≤ 1} bl bl 123 A Verified ODE Solver and the Lorenz Attractor 81 Having (bounded) linear functions as a separate type makes many formulations eas- ier. For example, consider Harrison’s formalization of multivariate analysis (from which Isabelle/HOL’s analysis descended). In Harrison’s formalization continuity is formalized for n m functions f of type R → R . (continuous f (at x )) = (∀e > 0. ∃d > 0. ∀y. x − y < d ⇒ fx − fy < e) Most of Harrison’s formalization is geared towards viewing derivatives as linear functions n m of type R → R . For continuously differentiable functions, one therefore needs to reason n n m about functions f : R → (R → R ),where f x is the derivative of f at a point x. Continuity of f is written in an explicit ε-δ form and involves the operator norm onorm : n m (R → R ) → R, which is quite verbose: (∀e > 0. ∃d > 0. ∀y. |x − y| < d ⇒ onorm (λv. f x v − f y v) < e)) The ε-δ form could of course be captured in a separate definition, but this would be very similar to the definition of continuity and would introduce redundancy. In the Isabelle/HOL formalization, continuous is defined for functions f : α → β for topological spaces α and β.If α and β are normed vector spaces, the above equality for continuous holds in Isabelle/HOL, too. And indeed, the norm of bounded linear functions is defined using onorm such that onorm (λv. ( f x ) · v − ( f y) · v) = f x − f y bl bl holds. Then, continuity of a derivative f : α → (α → β) can simply be written as bl (continuous f (at x )), which is a better abstraction to work with and also avoids redundant formalizations for different kinds of continuity. 2.2 Dynamical Systems An ODE induces a continuous dynamical system via the notion of flow. A standard technique to reason about such systems is its Poincaré map. We keep the presentation at a high level, since details can be found in the publications [25,26]. 2.2.1 The Flow We consider an autonomous ODE with right hand side f . Under mild assumptions, there exists a solution φ(t ), which is unique for an initial condition x (0) = x and satisfies the differential equation: x˙(t ) = f (x (t )) To emphasize the dependence on the initial condition, one writes φ(x , t ) for the solution. This solution depending on initial conditions is called the flow of the differential equation: Definition 4 (Flow)The flow φ(x , t ) is the (unique) solution of the ODE x˙(t ) = f (x (t )) with initial condition φ(0) = x The flow is only well-defined on the so-called existence interval of the solution, which depends on the initial value. Definition 5 (Existence Interval) t ∈ ex-ivl (x ) ⇒ φ(x , t ) = f (φ (x , t )) 0 0 0 The flow φ and the existence interval ex-ivl provide a clean interface to talk about solu- tions of ODEs. The property of the generic notion of flow makes it possible to easily state composition of solutions and to algebraically reason about them. Flowing from x for time s + t is equivalent to first flowing for time s, and from there flowing for time t: 123 82 F. Immler Theorem 4 (Flow Property) {s, t, s + t}⊆ ex-ivl (x ) ⇒ φ(x , s + t ) = φ(φ(x , s), t ) 0 0 0 For Tucker’s proof, one needs to study how sensitive the flow depends on perturbations of the initial value. We use two main results: One, the flow depends continuously on initial values. Two, if the ODE f is continuously differentiable, then so is the flow. We first take a look at the domain Ω = {(x , t ) | t ∈ ex-ivl (x )} ⊆ X × T of the flow. (t, x ) ∈ Ω means that we can flow a point starting at x for at least time t. Intuitively, solutions starting close to x can be followed for times that are close to t. In topological parlance, the state space is open. Theorem 5 (Open State Space) open Ω One can show that solutions deviate at most exponentially fast: ∃K . φ(x , t ) −φ(y, t ) < K |t | x −ye (using Grönwall’s lemma). Therefore, by choosing x and y close enough, one can make the distance of the solutions arbitrarily small. In other words, the flow is a continuous function on the state space: Theorem 6 (Continuity of Flow) continuous-on Ωφ Continuity states that small deviations in the initial values result in small deviations of the flow. One can be more precise on the way initial deviations propagate. The propagation of initial deviations through the flow (φ := λx.φ(x , t )) can be approximated by a linear function, the derivative Dφ| · v ≈ φ(x , t ) − φ(x + v, t ). We formalize the fact that the derivative of the flow is the solution of a differential equation in the space of bounded linear mappings, the so-called variational equation. Theorem 7 (Variational Equation) W (t ) = D f | · W (t ) φ(x ,t ) bl W (0) = 1 bl Solving this ODE numerically gives a means to obtain numerical bounds on the derivative, which is the approach that we pursue in our algorithms. 2.2.2 The Poincaré Map The Poincaré map is an important tool for studying dynamical systems. Whereas the flow describes the evolution of a continuous system with respect to time, it is the Poincaré map that allows us to study the evolution with respect to some space variables. A Poincaré section is a subset Σ of the state space, which is in general given as an implicit surface Σ ={x | s(x ) = c} with continuously differentiable s. For Tucker’s proof, one chooses s(x , y, z) = z and c = 27. The Poincaré map P(x ) is defined as the point where the flow starting from x first hits the Poincaré section Σ. It is defined with the help of the first return time τ(x ). τ depends on the flow φ (and therefore on the ODE f ) and the Poincaré section Σ, but we keep those dependencies implicit. Definition 6 (First Return Time) τ(x ) is the least t > 0 such that φ(x , t ) ∈ Σ. Obviously, τ is only well-defined for values that actually return to Σ, which we encode in the predicate returns-to : 123 A Verified ODE Solver and the Lorenz Attractor 83 Definition 7 returns-to (Σ, x ):=∃t > 0.φ(x , t ) ∈ Σ The return time can then be used to define the Poincaré map as follows: Definition 8 (Poincaré map) P(x ):= φ(x,τ(x )) It is interesting to note that this way of defining the return time and Poincaré map differs from the usual approach in textbooks. Textbooks study Poincaré maps in a neighborhood around a periodic point x ∈ Σ, i.e., P(x ) = x. This makes it easy to directly apply the implicit function theorem and transfer continuity and differentiability from the flow to the Poincaré map while guaranteeing that τ and P are well-defined. Also, one views P as a mapping on Σ, i.e. P : Σ → Σ. Tucker’s proof, however, requires a more flexible notion of Poincaré map and our notion of τ is more flexible: it is well-defined also for values outside of Σ. This enables reasoning about intermediate sections: Tucker, e.g., composes a sequence of local Poincaré maps between intermediate sections Σ, Σ , ··· ,Σ ,Σ in order to get bounds on the global Poincaré map 1 n Σ → Σ. The goal of Tucker’s computations is a sensitivity analysis of the flow of the Lorenz system and of its Poincaré map. Its derivative can be given in terms of the derivative of the flow and the function s defining the implicit surface for Σ ={x | s(x ) = c}. Theorem 8 (Derivative of Poincaré map) Ds| · (Dφ| · h) P(x ) (x ,τ (x )) DP| · h = Dφ| · h − f (P(x )) x (x ,τ (x )) Ds| · ( f (P(x ))) P(x ) For a rough intuition, the derivative DP| ·h of the Poincaré map is related to the derivative of the flow Dφ| · h. But it needs to corrected in the direction f (P(x )) in which the (x ,τ (x )) flow passes through Σ, because P varies only on Σ and not through it. This correction factor also depends on the tangent space Ds| of the section Σ at P(x ). P(x )) 3 Rigorous Numerics Rigorous (or guaranteed) numerics means computing with sets that are guaranteed to enclose the real quantities of interest. Enclosures can in principle be any data structure that represents sets of real values. Popular choices are intervals, zonotopes, or Taylor models. For formal- izations, it is useful to have a deep embedding of arithmetic expressions as done, e.g., by Dorel and Melquiond [34]aswellasHölzl [13]. This work builds on Hölzl’s language of arithmetic expressions given in excerpts in Fig. 4. Independent of the choice of representation of enclosures, a rigorous approximation scheme approx needs to satisfy that for all lists of real values xs in an enclosure XS,the interpretation [[e]] of the expression e is contained in the result of approximation scheme xs evaluated for the enclosure XS. We could call this the fundamental property of rigorous numerics. x ∈ XS ⇒ [[e]] ∈ approx eXS xs 123 84 F. Immler Fig. 4 Data type of arithmetic expressions and interpretation function (λexs. [[e]] ) : aexp → R list → R xs The key approach in our formalization is to remain agnostic about the concrete approx- imation scheme as long as possible and formalize results on the level of deeply embedded aexp expressions. The central result is an implementation of a second order Runge–Kutta method on the level of aexp expressions. As concrete instance for an approximation scheme approx , we use affine arithmetic [6], an improvement over interval arithmetic that tracks linear dependencies between program variables. 3.1 Expressions for Vectors To represent vectors, we use lists of expressions. A list of expressions es : aexp list is interpreted with (λes vs. [[es ]] ) : aexp list → R list → α componentwise as Euclidean vs space α: [[es ]] = eucl-of-list (map (λe. [[e]] ) es) vs vs In contrast to the interpretation, approximation of a list of expressions should not be componentwise: an approximation function for lists of expressions should be of type approx : aexp list → R list set → R list set , which allows approx to keep track of dependencies between the components of the result. If the type were e.g., aexp list → R set list → R set list , this could only represent the Cartesian product of the component enclosures. n m For a function f : R → R , a deep embedding f is a list of expressions (of length m), that is interpreted over a list of n variables. [[ f ]] = f (eucl-of-list xs) e xs The derivatives with respect to one variable can be computed symbolically from the struc- ture of the expression. This can also be used to compute partial derivatives on the level of expressions. In the multivariate setting, the derivative D f | of f at x is the matrix of its partial derivatives. In general, we can represent matrices as a flat list (according to eucl-of-list /list-of-eucl which are also defined for matrices). For computing derivatives, how- ever, we directly produce an expression that is interpreted as the product of the derivative matrix with a vector: [[D (n, f ,v )]] := D f | · [[v ]] e e e xs eucl-of-list xs bl e xs D takes the derivative with respect to the first n variables, and xs is assumed to be of length at least n. This way, we can produce expressions for higher derivatives: D (n, f ,v ):= f e e e i +1 i D (n, f ,v ):= D (n, D (n, f ,v ), v ) e e e e e e e e Note that the proper interpretation can only be expressed in Isabelle’s type system for fixed values of i: the resulting object is an i-linear function, so the resulting type depends on 123 A Verified ODE Solver and the Lorenz Attractor 85 a term argument. This could also be encoded as functions taking lists as arguments, but fixed values of i suffice for our purposes and we find the interpretation as curried linear mappings more natural. E.g., [[D (n, f ,v )]] =[[ f ]] e e x e x [[D (n, f ,v )]] =D f | · [[v ]] e e x eucl-of-list x bl e [[D (n, f ,v )]] =D(λy. D f | )| · [[v ]] · [[v ]] e e x y eucl-of-list x bl e bl e ··· 3.2 A Runge–Kutta Method On the level of expressions, we verified a two-stage Runge–Kutta method rk (x ) = x + 1 1 h · ψ (x ), with ψ (x ) = (1 − ) f (x ) + f (x + hpf (x )). This Runge–Kutta method h h 2 p 2 p rk (x ) approximates the solution up to third order: |φ(x , h) − rk (x )|∈ O(h ).The h 0 0 h 0 third order term stems from (multivariate) Taylor series expansions of the solution φ and the approximation scheme rk .Ifweset f := λx . D f | and f := λx . D f | , then the h x x remainder is contained in the convex hull of any set that contains rk-remainder (s , s ) for h 1 2 all s , s ∈[0, 1]. 1 2 h 1 rk-remainder (s , s ):= · f x (hs + t ) h 1 2 1 2 3 · f (x (hs + t )) · f (x (hs + t )) bl 1 bl 1 + f x (hs + t ) · f (x (hs + t )) · ( f (x (hs + t ))) 1 bl 1 bl 1 − f x (t ) + hps f (x (t )) · f (x (t )) · f (x (t )) 2 bl bl Theorem 9 (Runge–Kutta Method with Remainder Term) φ(x , h) ∈ rk (x ) + convex-hull (rk-remainder ([0, 1], [0, 1])) 0 h 0 h In order to use this theorem for rigorous numerical computations, we produce a deep embedding of the expressions for rk and rk-remainder . We do so for an arbitrary ODE h h n n f : R → R with deep embedding [[ f ]] = f (eucl-of-list x ). Expressions for f and f e xs 1,2 are computed symbolically from f via D from the previous Sect. 3.1. 3.3 Straight Line Programs The expression for rk-remainder from the previous Sect. 3.2 contains common subexpres- sions. This is not desirable because one would need to perform redundant computations. We therefore follow Dorel and Melquiond’s [34] approach and use straight-line programs with static single assignment instead of plain expressions. For us, a straight line program is just a list of arithmetic expressions, which is interpreted according to function slp : aexp list → R list → R list : slp [] x = x slp (e::es) x = slp es ([[e]] ::x ) The idea is that a straight line program only contains unary or binary operations, although this is not required by the definition. The result of the operation is put on top of the evaluation stack. 123 86 F. Immler The following example illustrates sharing the term x + y in the expression (x + y)(x + y): slp [Add (Var 0)(Var 1), Mult (Var 0)(Var 0)][x , y]=[(x + y)(x + y), x + y, x , y] We provide a function slp-of , which eliminates common subexpressions by traversing an expression bottom-up and saving subexpressions in a map that gives the index of the subex- pression in the resulting straight line program. At run-time (this is important to be able to use the ODE solver as a stand-alone tool), in an initialization phase, the ODE solver computes symbolically the derivatives in the expression for rk and rk-remainder , does constant propagation (as derivatives can produce 0 constants, this is beneficial) and then compiles the resulting expression with slp-of into a straight-line program, which is then used in the course of approximating the ODE in a series of steps. 3.4 Affine Arithmetic Up to now, we have kept the discussion on the level of expressions, let us now motivate affine arithmetic as a concrete approximation scheme. The most basic data structure to represent sets is closed intervals [a, b]={x | a ≤ x ≤ b}, but those suffer from the wrapping effect: rotated boxes cannot be represented without large overapproximations. Moreover dependencies between variables are lost, e.g. for an enclosure x ∈[0, 2], the expression x − x evaluates to [−2, 2] in interval arithmetic whereas the exact result would be representable as the interval [0, 0]. Affine arithmetic [6] improves over interval arithmetic by tracking linear dependencies. An affine form A is a function where only finitely many arguments map to nonzero values. It is interpreted for a valuation ε : N → R : affine ε A:= A + ε A 0 i i Looking at the interpretation, one often calls the terms ε noise symbols and A generators. i i The idea is that noise symbols are shared between affine forms and that they are treated symbolically, as formal parameters: the sum of two affine forms is given by the pointwise sum of their generators, and multiplication with a constant factor is also done componentwise. affine ε(A + B):= (A + B ) + ε (A + B ) 0 0 i i i affine ε(cA):= cA + ε (cA ) 0 i i The range of an affine form is the set of all affine evaluations where the noise symbols range over the closed interval [−1, 1]. For the range of a list of affine forms, those are evaluated jointly for the same valuation of the noise symbols, reflecting the intuition that those are shared. range A:= {affine ε A |∀i. − 1 ≤ ε ≤ 1} { } joint-range AS:= map (affine ε) AS |∀i. − 1 ≤ ε ≤ 1 As a concrete example, let us examine how affine arithmetic handles the dependency problem in the introductory example x − x for x ∈[0, 2]. The interval [0, 2] is represented by the affine form 1 + 1 · ε . This is the affine form given by the function X:= (λi. if i = 0∨i = 1 then 1 else 0). For this function, range X =[0, 2] holds. Then, in affine arithmetic, 123 A Verified ODE Solver and the Lorenz Attractor 87 (1 + 1 · ε ) − (1 + 1 · ε ) = 0 + 0 · ε , which corresponds to the constant zero function. 1 1 1 Therefore range (X − X ) ={0}. In general, with the help of range and joint-range , we can express correctness of a binary operation like addition e.g., as follows: [a, b]∈ joint-range [A, B]⇒ (a + b) ∈ range (A + B) Nonlinear operations like multiplication or division are linearized, adding the linearization error as a fresh noise symbol. Consider e.g., multiplication: (affine ε A) ∗ (affine ε B) = A B + ε (A B + A B ) + ε A ε B 0 0 i 0 i i 0 i i i i i i >0 i >0 For a proper valuation with ε ∈[−1, 1], the last summand on the right can be bounded by ( |A |)( |B |). Therefore, if k is fresh in A and B, one can set i i i >0 i >0 affine ε(A ∗ B):= A B + ε (A B + A B ) + ε |A | |B | 0 0 i 0 i i 0 k i i i i >0 i >0 and the k-th generator bounds the linearization error such that multiplication of affine forms is conservative: [a, b]∈ joint-range [A, B]⇒ a ∗ b ∈ range (A ∗ B) Similar to the additional noise symbol for a linearization error, also round-off errors can be included as additional noise symbols. We provide affine approximations for the primitive functions listed in Fig. 4. Expressions (and straight line programs) involving these functions can then be approximated by recursively keeping track of the next fresh noise symbols. During longer computations more and more noise symbols will be added to the affine form, impairing performance in the long run. The number of noise symbols can be reduced by summarizing (or condensing) several noise symbols into a new one. This process discards the correlation mediated by the summarized noise symbols, so a trade-off needs to be found between precision and efficiency. We consider a list of affine forms AS and use the notation AS := map (λA. A ) AS. We call total deviation |AS|: = map (λA. |A |) AS the i i i componentwise sum of absolute values. We summarize all symbols i with |AS |≤ r |AS| for a given summarization threshold r. We found that it is important to perform the above comparison componentwise and not take (like proposed for other implementations of affine arithmetic [6,42]) the infinity norm on both sides. This is of particular importance when components differ a lot in magnitude. Apart from looking at affine forms as a formal sum, the joint-range of a list of affine forms can also be interpreted geometrically as zonotopes: centrally symmetric, convex polytopes. A zonotope can be visualized as the Minkowski sum (the set of all possible sums of elements of two sets X ⊕ Y ={x + y. x ∈ X ∧ y ∈ Y }) of the line segments defined by the generators. For example, Fig. 6 depicts a two-dimensional zonotope with three generators, [−1, 1]a ⊕[−1, 1]a ⊕[−1, 1]a . Figure 5 contains a three-dimensional zonotope with 1 2 3 three generators (a parallelotope), namely the zonotope Z defined as follows. ⎛ ⎞ [1ε + 0ε + 1ε , 1 2 3 ⎝ ⎠ 0ε + 2ε + 5ε , Z:= joint-range 1 2 3 0ε + 0ε + 20ε ] 1 2 3 123 88 F. Immler Z G 20 G -5 -10 -15 -20 -2 -4 x1 -5 -6 -10 x2 -8 -15 Fig. 5 Three dimensional zonotope Z and intersection with hyperplane G 4 Computational Geometry An important step for Tucker’s proof is the reduction to a Poincaré map: intersecting the flow of the ODE with a plane in the state space. In our algorithms, the flow is approximated with affine arithmetic expressions, therefore enclosed by zonotopes. In order to compute where the flow intersects the hyperplane, one needs to compute the intersection of the enclosing zonotope with the hyperplane (see Fig. 5). This is an interesting geometric problem and we verified an approximative algorithm due to Girard and Le Guernic [7]. At its core, the algorithm is similar to convex hull compu- tations. We can build on a nice abstraction to reason about it, namely Knuth’s theory of counterclockwise (ccw) systems [27]. We needed, however, to extend Knuth’s theory from discrete to continuous sets of points. 4.1 Girard and Le Guernic’s Algorithm The complexity for computing the exact intersection of a zonotope with a hyperplane grows exponentially with the number of generators. An overapproximation of the zonotope before computing the intersection is possible but can lead to overly coarse approximations. Therefore Girard and Le Guernic [7] proposed a way to directly compute overapproximations to the intersection. The first idea is to overapproximate a given set X tightly from a set D of directions, which can be chosen arbitrarily. For every direction d ∈ D ⊆ R , the infimum m and supremum M of the sets {x , d. x ∈ X } needs to be determined (_, _ denotes the inner product, also known as the dot product). Geometrically speaking, m and M give the position d d of two hyperplanes with normal vector d. The two hyperplanes bound X from below and above, respectively. An overapproximation P is then given by the points between all of these hyperplanes: X ⊆ P ={x ∈ R . ∀d ∈ D. m ≤x , d≤ M } d d The second observation of Girard and Le Guernic is that when the set X is the intersection of some set Z with a hyperplane G ={x . x , g= c}, then the computation of the overap- proximation P can be reduced to a two-dimensional problem with the linear transformation n 2 Π : R → R , Π (x ) = (x , g, x , d). g,d g,d x3 A Verified ODE Solver and the Lorenz Attractor 89 Fig. 6 Corners c and edges of a zonotope { ε · a |−1 ≤ ε ≤ 1}, generators a , a , a , intersecting i i i i 1 2 3 line L Lemma 1 (Reduction to Dimension Two) {x , d. x ∈ Z ∩ G}={y.(c, y) ∈ Π (Z )} g,d This lemma is an easy consequence of the definitions of G and Π . For every direction g,d d, the theorem allows to reduce the computation of the intersection Z ∩ G on the left-hand side to the intersection of the projected two-dimensional zonotope Π (Z ) with the vertical g,d line L ={(x , y). x = c}. Computing the intersection of a two-dimensional zonotope like the one given in Fig. 6 and a vertical line L can by done by computing bounds on the intersection of the vertical line L with every edge. This is easy and intuitive. The more challenging part is to compute the set of edges of a two-dimensional zonotope, which we sketch in the following. 4.2 Computing the Set of Edges. First of all, one can assume that all generators point upwards. One then starts at the lowest corner (c in Fig. 6) and appends to it the “rightmost” generator a (twice) to reach c .One 0 1 1 then continues with the “rightmost” of the remaining generators, a and is in the process essentially “wrapping up” the hull of the zonotope. In order to verify such a process, we need a way to reason about “rightmost” vectors (a total order). Similar ideas of “wrapping up” a set of points also occur for convex hull algorithms. 4.3 Knuth’s CCW System In order to verify geometric algorithms, one needs a formal notion of the geometric concepts involved. For convex hull algorithms, Knuth [27] has given a small theory that axiomatizes the notion of orientation of points. The intuition is that for three points p, q, r in the plane, visiting them in order requires either a counterclockwise (ccw) turn (written pqr) or clockwise (¬pqr) turn. Knuth observed that already few of the properties fulfilled by the ccw predicate pqr suffice to define a theory rich enough to formalize many concepts in computational geometry. The notion of ccw system is a set of points together with a ccw predicate written pqr for points p, q, r. The ccw predicate needs to satisfy the following properties, inspired by 123 90 F. Immler (a) (b) (c) Fig. 7 ccw axioms: dashed predicates are implied by solid ones. a Cyclic symmetry, b interiority, c transitivity the relations satisfied by points in the plane. For all axioms in the following, there is the additional implicit assumption that the involved points are pairwise distinct. For three points, only simple axioms need to be fulfilled: – Cyclic Symmetry: pqr ⇒ qrp – Antisymmetry: pqr ⇒ ¬prq – Nondegeneracy: pqr ∨ prq Cyclic symmetry and the more interesting case of interiority, which involves four points, are illustrated in Fig. 7. Interiority states that if one point t is left of three lines pq qr r p, then the three other points are oriented in a triangle according to pqr. – Interiority: tpq ∧ tqr ∧ trp ⇒ pqr The most important tool for reasoning is transitivity, which involves five points and works if three points p, q, r lie in the half-plane left of the line ts, i.e., tsp ∧ tsq ∧ tsr. Then, fixing t as first element for the ccw relation, we have transitivity in the second and third element: tpq ∧ tqr ⇒ tpr (see Fig. 7c). – Transitivity: tsp ∧ tsq ∧ tsr ∧ tpq ∧ tqr ⇒ tpr The same intuition also holds for the other side of the half-plane: – Dual Transitivity: stp ∧ stq ∧ str ∧ tpq ∧ tqr ⇒ tpr Knuth shows that under the assumptions of Cyclic Symmetry, Antisymmetry, and Non- degeracy, Transitivity holds if and only if Dual Transitivity holds. Knuth requires more than half a page of low level reasoning, but as this reasoning is carried out abstractly in a small first order theory, sledgehammer (Isabelle’s interface to various automatic theorem provers) is able to find a proof that consists of just one single invocation of an automated prover. 4.3.1 Total Order from CCW As sketched earlier, in order to compute the edges of a zonotope, we need to be able to select a “rightmost” element of a set of vectors. With the transitivity relation presented before, we can obtain a total order on vectors which allows us to do just that: Given a center t and another point s, the orientation predicate tpq can be used to define a total order on points p, q in the half-plane left of ts, i.e., p < q iff tpq. From Antisymmetry and Nondegeneracy of the ccw 123 A Verified ODE Solver and the Lorenz Attractor 91 system, we get antisymmetry and totality for the order <. Transitivity of the order < follows from the axiom Transitivity of the ccw system and the assumption that all points are in the half-plane left of ts. This ordering is then used to sort the list of generators such that they actually “wrap up” the zonotope. 4.3.2 Instantiation for Points in the Plane Up to now, our reasoning was based abstractly on ccw systems, but of course we also want to use the results for a concrete ccw predicate. Well known from analytic geometry is the fact that ccw orientation is given by the sign of the following determinant |pqr |: p p 1 x y p ∗ q + p ∗ r + q ∗ r − x y y x x y |pqr |: = q q 1 = x y r ∗ q + r ∗ p + q ∗ p x y y x x y r r 1 x y Points are collinear iff |pqr|= 0. Under the assumption that one works with a finite set of points where no three points are collinear, the following predicate pqr satisfies the axioms of a ccw system. pqr :=|pqr | > 0 Most axioms can easily be proved using Isabelle/HOL’s rewriting for algebraic structures. Transitivity is slightly more complicated, but can also be solved automatically after a proper instantiation of Cramer’s rule, which is easily proved automatically: |tqr ||stp|+|tpq||str | |tpr|= , if |stq| = 0 |stq| 4.3.3 CCW on a Vector Space Knuth presented his axioms with a finite set of discrete points in mind, in our case we need to talk about orientation of arbitrary points in a continuous set. We therefore require consistency of the orientation predicate when vector space operations are involved. One obvious requirement is that orientation is invariant under translation (Fig. 8a). With translation invariance, we can reduce every ccw triple to a triple with 0 as origin, and from there it is easy to state consistency with respect to scaling: If at q, there is a ccw turn to r,then every point on the ray from 0 through q will induce a ccw turn to r (Fig. 8b). Negative scalars can be treated by requiring that reflecting one point at the origin inverts the ccw predicate (Reflection). Furthermore, the addition of vectors q and r, which are both ccw of a line p needs to be ccw of p as well. > > – Translation: (p + s)(q + s)(r + s) = pqr > > – Scaling: α> 0 ⇒ 0(α · q)r = 0qr – Reflection: 0(−p)q = 0qp – Addition: 0pq ⇒ 0 pr ⇒ 0 p(q + r ) The predicate pqr simplifies much of the reasoning, because it satisfies the axioms of a ccw system. It does, however, ignore collinear points and therefore all the points of the zonotope that lie on its edges. In order to also include those into the reasoning, we define the slightly relaxed ccw predicate pqr , which holds for all points on the line through pq and for all points on the half-plane left of pq. pqr :=|pqr|≥ 0 123 92 F. Immler (a) (b) Fig. 8 ccw axioms on a vector space. a Invariance under translation. b Invariance under scaling The situation is as follows: pqr is the actual specification that we care about. But it does not satisfy the axioms of a ccw system, which makes reasoning very convenient. We therefore first prove the corresponding properties for the ccw system pqr . With a simple argument > ≥ on continuity, the results about pqr carry over to pqr and therefore the whole zonotope. 5 Program and Data Refinement We use two different approaches to turn abstract formalizations into executable constructs: In situations where abstract operations directly correspond to concrete, executable ones, we use light-weight data refinement via the code generator. In more demanding situations, we employ a dedicated framework for nondeterministic specifications and stepwise program refinement, namely the Autoref tool. 5.1 Light-Weight Data Refinement For light-weight data refinement [10] via the code generator, abstract operations need to be mapped directly to concrete, executable ones. Examples of such abstract types are affine forms and real numbers. Consider the type of real numbers R. We call it abstract, because as an uncountable set, R is obviously not computable. But one can restrict oneself to working on a computable subset of the real numbers. In our case, we use software floating point numbers F ={m ·2 | m, e ∈ Z} for (unbounded) integers m, e ∈ Z. We instruct the code generator to use an uninterpreted constructor Real-of-Float : F → R to represent real numbers. Operations on real numbers are then computed by pattern matching on that constructor and executing the corresponding concrete operation, e.g., for addition: (Real-of-Float f ) + (Real-of-Float g) = (Real-of-Float ( f + g)) Since this implementing equation is proved as a theorem, such a setup does not change the trusted code base. All one has to ensure is that all abstract operations that occur in the code can be executed in the concrete representation. Affine forms are abstractly a subtype of functions N → R. They are implemented using a type of association lists (N × R) list that are (reverse) strictly sorted according to the keys. This sparse representation is useful because the largest index of a non-zero generator can be directly read off by inspecting only the first element. Adding a fresh generator can be 123 A Verified ODE Solver and the Lorenz Attractor 93 done by simply prepending the new element. Binary operations are efficiently and easily implemented by merging the two lists of generators. 5.2 Autoref We use Lammich’s framework Autoref for (automatic) refinement [29,30] in Isabelle/HOL. Autoref allows the user to specify algorithms on an abstract level and provides tool support for stepwise refinement [1] to concrete, executable implementations. In this section we present a setup of Autoref for rigorous numerical algorithms: We provide abstract specifications for elementary operations of common rigorous numerical algorithms as well as suitable implementations. 5.2.1 Nondeterministic Specifications An important insight when verifying algorithms that use rigorous numerical enclosures is the fact that, for correctness, any enclosure of the real quantity suffices. We model this with appropriate nondeterministic specifications. Autoref is based on a nondeterminism monad α nres , where programs can either fail or yield a set of values as result. datatype α nres = FAIL | RES (α set ) The refinement relation ≤ on nres has FAIL as top element and RES S ≤ RES T iff S ⊆ T . For deterministic results, we write return x:= RES {s}. We write for specifications spec P:= RES {x . Px } the result of all values satisfying the predicate P. This allows one to specify correctness, e.g, of a program f whose inputs x satisfy the precondition P and every possible value y in its nondeterministic result satisfies the post- condition Q: ∀x . Px ⇒ fx ≤ spec (λy. Qy). In this setting, we specify a set of operations that are useful in the context of verifying rigorous numerical algorithms, i.e., algorithms that manipulate enclosures. These operations are best modeled nondeterministically, because one is often only interested in some safe result. Subdivisions are a means to maintain precision, we therefore have the following abstract specifications for splitting a set (with or without the possibility to perform overapproxima- tions): split-spec X:= spec (λ(A, B). X ⊆ A ∪ B) split-spec X:= spec (λ(A, B). X = A ∪ B) The following specifications yield some lower/upper bound on the set, not necessarily exact: Inf-spec X:= spec (λi. ∀x ∈ X. i ≤ x ) Sup-spec X:= spec (λs. ∀x ∈ X. x ≤ s) Depending on the concrete representation of sets, one might not be able to decide certain properties, but only give a positive answer if the precision is sufficient. We therefore have a specification that may guarantee disjointness. disjoint-spec XY:= spec (λb. b ⇒ X ∩ Y =∅) As seen in the previous Sect. 4, depending on the data structure, one can not (or does not want to) compute an exact representation for the intersection of sets. These specifications 123 94 F. Immler allow one to overapproximate an intersection, while guaranteeing that the result does not exceed one of the arguments: inter-spec XY:= spec (λR. X ∩ Y ⊆ R ∧ R ⊆ X ) inter-spec XY:= spec (λR. X ∩ Y ⊆ R ∧ R ⊆ Y ) To bridge the gap to concrete numerical computations and the results from Sect. 3,weuse a specification for overapproximations of evaluating straight line programs: approx-slp-spec slp X:= spec (λR. ∀x ∈ X. [[slp]] ∈ R) 5.2.2 Refinement Relations In Autoref, specifications in the nres monad are transferred to executable constructs in a structured way. Autoref is centered around a collection of so-called transfer rules. Transfer rules relate abstract with concrete operations. A transfer rule involves a transfer relation R::(γ × α) set , which relates a concrete implementation c::γ with an abstract element a::α and is of the following form. (c::γ, a::α) ∈ R Transfer rules are used to structurally synthesize concrete algorithms from abstract ones. Relations and relators (which combine relations) are used to express the relationship between concrete and abstract types. br is used to build a relation from an abstraction function a::γ → α andaninvariant I on the concrete type. br aI:={(c,ac) | Ic} Natural Relators For the types of functions, products, sets, or data types like lists and nres , one uses the natural relators A → B, A × B, Aset , Alist-rel , Anres with relations r r r r A, B for the argument types. ( f, f ) ∈ A → B ⇐⇒ ∀(x , y) ∈ A.(fx , f x ) ∈ B ((a, b), (a , b )) ∈ A × B ⇐⇒ (a, a ) ∈ A ∧ (b, b ) ∈ B (X, X ) ∈Aset ⇐⇒ (∀x ∈ X. ∃x ∈ X .(x , x ) ∈ A) ∧ (∀x ∈ X . ∃x ∈ X.(x , x ) ∈ A) (xs, xs ) ∈Alist-rel ⇐⇒ length xs = length xs ∧ (∀i < length xs.(xs , xs ) ∈ A) (RES X, RES X ) ∈Anres ⇐⇒ (X, X ) ∈Aset r r Representing Vectors We represent vectors (an arbitrary type α of class Euclidean space) as lists of real numbers where the length matches the dimension of the Euclidean space. lv-rel := br eucl-of-list (λxs. len xs = DI M (α)) This way, the concrete algorithm is monomorphic, which has the advantage that it can be generated once and for all and can therefore be used as a stand-alone tool. 123 A Verified ODE Solver and the Lorenz Attractor 95 Representing Enclosures We provide several implementations for the sets that can be used as enclosures. Intervals are represented by pairs of element types (which, in turn are implemented via some relation A): Aivl :={((a , b ), [a, b]) | (a , a) ∈ A ∧ (b , b) ∈ A} Zonotopes are represented using the joint range joint-range of affine forms affine := br (λA. eucl-of-list (joint-range A)) (λ_. True ) We use a symbolic representation of planes using the data type constructor Sctn that keeps normal vector n and translation c of a hyperplane. It is interpreted using plane-of (Sctn nc):={x |x , n= c} halfspace (Sctn nc):={x |x , n≤ c} for the hyperplane itself or for the halfspace below the hyperplane, where x , n is the inner product (also called dot product). Asctn is the natural relator that allows one to change the representation of the normal vector. With this, we can give a concrete implementation for hyperplanes and half-spaces. Aplane :=Asctn ◦ br plane-of (λ_. True ) r r Ahalfspace :=Asctn ◦ br halfspace (λ_. True ) r r For those relations, plane-of , halfspace,and ∩ are easily implemented with the identity function or as pair. On the abstract level, they describe useful objects that are convenient to reason about. ((λx . x ), plane-of ) ∈Asctn → Aplane r r r ((λx . x ), halfspace ) ∈Asctn → Ahalfspace r r r ((λxy.(x , y)), ∩) ∈ A → B → A, Binter r r r We will see that in some algorithms, one maintains a collection of enclosures, but abstractly one likes to see them as just one enclosure. For a relation A : (β × α set ) set that implements single enclosures for sets of type α with some concrete representation of type β, and a relation S : (σ × β set ) set that implements sets of concrete elements β, we define a relation that represents the union of all those elements as follows: S, AUnion : (σ × α set ) set S, AUnion := S ◦Aset ◦ br λX. x (λ_. True ) r r x ∈X Currently, we use lists to implement the set of concrete representations S, for which we write AUnion :=list-set , AUnion , and operations like union or extracting one element lr r r (with the specification split-spec ) can be implemented with the respective operations on lists/sets: (λxs ys. return (xs@ys), ∪) ∈AUnion → AUnion → AUnion lr r lr r lr (λx . return (hd x , tl x ), split-spec ) ∈AUnion → A × AUnion nres = lr r r lr r 123 96 F. Immler Relations to Guide Heuristics Often, in particular to guide heuristics, an algorithm needs to carry around information, which does not influence correctness proofs. An ODE solver for example, modifies its step size, also based on previous values. An implementation needs to carry this information around, but for verifying the algorithm, this only introduces unneces- sary clutter. We therefore introduce a relation that carries more information (implemented via A) in the implementation, but keeps the abstract semantics (implemented via B): A, Binfo :={((a , b ), b) |∃a.(a , a) ∈ A ∧ (b , b) ∈ B} Adding information is simply done by using a pair in the implementation side. Semantically, this information is simply discarded (put-info ab:= b). Information can be extracted with get-info , which is semantically just an arbitrary element (get-info b:= spec (λ_. True )). The implementations are straightforward: (λab.(a, b), put-info ) ∈ A → B → A, Binfo r r r ((λ(a, b).return a), get-info ) ∈A, Binfo →Anres r r An example of its usage is illustrated later on in Algorithm 1. 6 Reachability Analysis Overall, we design an algorithm that computes a Poincaré map with a list of intermediate Poincaré sections. The global idea (illustrated in Fig. 9) is as follows: starting from a set X , perform a series of single time discretization steps. If reachable sets grow above a given threshold, subdivide them (Sects. 6.3 and 6.4). Stop before an intermediate (or the final) section would be hit, then resolve the Poincaré map at that section (Sect. 6.5). For Tucker’s proof, it is important to also track the matrix of partial derivatives together with the solution. To this end, one can encode the derivative as a higher-dimensional ODE and use essentially the same algorithms as before. This instrumentation is presented in Sect. 6.7. 6.1 The Framework We use the high-level constructors and abstract specifications from the previous Sect. 5. We remain agnostic about the type of enclosures, for which we assume a relation encl and implementations for the abstract operations that are needed for the reachability analysis algo- rithms: an approximation scheme for expressions approx-slp-spec , enclosures from intervals using an implementation encl-of-ivl , lower and upper bounds with Inf-spec , Sup-spec,inter- sections with a plane inter-spec (note that the relation fixes the second argument to represent a plane, abstractly inter-spec is just intersection on sets): – (approx-slp , approx-slp-spec ) ∈ slp → encl → encl option nres r r r r r r r – (λxy. encl-of-ivl xy,λxy. [x , y]) ∈ lv-rel → lv-rel → encl r r r – (inf-encl , Inf-spec ) ∈ encl → lv-rel nres r r r – (sup-encl , Sup-spec ) ∈ encl → lv-rel nres r r r – (split-encl , split-spec ) ∈ real → nat el → encl → encl × encl nres ⊆ r r r r r r r r r r – (inter-encl-plane , inter-spec ) ∈ encl → lv-rel plane → encl nres 2 r r r r r r Currently, the only instantiation of this scheme is with affine arithmetic (in this case we set encl to affine ). Nevertheless, this structure keeps the formalization modular and one r r can imagine to add further instantiations—with e.g., Taylor models or centered forms—in the future. 123 A Verified ODE Solver and the Lorenz Attractor 97 Fig. 9 Continuous reachability and intermediate Poincaré sections 6.2 The Specification Our algorithms are supposed to compute enclosures for solutions of the ODE. We formalize the enclosure of an evolution from an initial set X to some other set Y with the ternary predicate ,where X  Y holds if the evolution flows every point of X ⊆ R to some point in Y ⊆ R and does not leave the set C in the meantime. We call C the flowpipe from X to Y . Definition 9 (Flows-to Predicate) X  Y:=∀x ∈ X. ∃t ≥ 0.φ(x , t ) ∈ Y ∧ (∀0 ≤ s ≤ t.φ(x , s) ∈ C ) C 0 0 X 6.3 Single Step In order to compute enclosures for a single step, one needs to first certify that a solution exists, which is the case for an initial value x and stepsize h if the iteration given by the Picard iteration from Sect. 2.1.2 has a unique fixed point. This is the standard approach from Bouissou et al. [4], also described in the setting of Isabelle [18]. The idea is that the expression Q (X ) = X +[0, h]· f (X ), which we can evaluate using approx-slp-spec , h 0 is an overapproximation of the Picard iteration and a post-fixed point certifies existence and a crude enclosure for solutions up to time h. This crude enclosure can be used as an overapproximation for the terms x (hs + t ) in the Runge–Kutta approximation scheme from Sect. 3.2. The function rk-step implements this and actually evaluates the Runge–Kutta approximation scheme twice: once for time h and once for the time interval [0, h], because this gives a much better enclosure for the flowpipe up to time h than the crude overapproximation from the Picard iteration. We prove the following specification. Theorem 10 rk-step Xh ≤ spec (λ(ε, C, Y ). X  Y ) The returned value  is an estimate for the approximation error. This is used for an adaptive step size control. Algorithm 1 shows an example how to use this heuristic (and another 123 98 F. Immler heuristic to split large sets), while (almost trivially, because the additional operations are either vacuous, the identity or overapproximations) satisfying the same specification. Theorem 11 single-step Xh ≤ spec (λ(C, Y ).X  Y ) The information on the last (and next) step size is only reflected in the refinement relation for the implementation of single-step : (single-step , single-step ) ∈ impl real , encl info → real , encl info Union , encl Union nres r r r r r r lr r lr r But this information does not clutter the verification of single-step or the statement of The- orem 11, which is very convenient. Algorithm 1 Single Step n n n 1: function single-step (X)  single-step : R set → (R set × R set ) 2: w ← width-spec X  semantically, this is spec (λ_. True ) 3: h ← get-stepsize X 4: if w ≤ max-width then  global parameter max-width 5: (ε, Y, C ) ← rk-step Xh 6: h ← adapt-step-size h ε  spec (λ_. True ) 7: return (put-info h Y, C ) 8: else 9: (Y, Z ) ← split-spec X 10: return (put-info h (Y ∪ Z ), (Y ∪ Z )) 6.4 Continuous Reachability Algorithm 2 Continuous Reachability Loop 1: function reach-cont (sctns, X ) 2: X ← X ; C ←∅; I ←∅ 3: while λ(X, C, I ). X =∅ do 4: (X , X ) ← split-spec X 1 2 5: (Y , C ) ← single-step X 1 1 1 6: d ← disjoint-spec C sctns 7: if d then 8: X ← X ∪ Y 2 1 9: C ← C ∪ C 10: I ← I 11: else 12: X ← X 13: C ← C 14: I ← I ∪ X 15: return (C, I ) Note that single-step returns (from an implementation point of view) a collection of enclosures, so we need some sort of work-list algorithm to resolve all currently reachable sets. Algorithm 2 does so. It maintains three kinds of sets (see also Fig. 9): X is the collection of currently “live” sets. C is the collection of all flowpipes explored so far. I is the collection 123 A Verified ODE Solver and the Lorenz Attractor 99 of sets where reachability analysis has stopped because of an intersection with a Poincaré section from sctns. The algorithm takes one element out of the “work-list” X by splitting the collection of enclosures using split-spec , performs a single step, checks for an intersection with one of the Poincaré sections and updates X, C, and I accordingly. The loop invariant of reach-cont is roughly the following: Elements from X flow via C to X and I , while avoiding sctns. X  (X ∪ I ) ∧ C ∩ sctns =∅ 0 C The specification of reach-cont follows immediately: Theorem 12 reach-cont (sctns, X ) ≤ spec (λ(C, I ). X  I ∧ C ∩ sctns =∅) 0 0 C It is worth noticing that the simplicity of the statement of this correctness theorem is due to the fact that the work-list and heuristic info is hidden via refinement relations in the implementation. If this were represented on the specification level (e.g., by using sets of enclosures paired with their current step size), the specification would have to be much more cluttered: ⎛ ⎞ ⎛ ⎛ ⎞ ⎞ ⎝ ⎠ ⎝ ⎝ ⎠ ⎠ x  x ∪ I (h,x )∈X (h,x )∈X Such a specification distracts the user and also automatic proof tools—we are therefore happy to hide this in the abstraction. 6.5 Resolve Intersection Algorithm 2 performs reachability analysis until each enclosure is just about to intersect an hyperplane. We compute the intersection from an enclosure X with a hyperplane H again in an iteration: continuous reachability steps are repeated as long as flowpipes intersect the hyperplane. The intersection of the individual flowpipes is computed with the geometric algorithm from Sect. 4. When the intersection is computed by flowing the reachable set through the hyperplane step by step, we get a set I consisting of individual intersections I . Many of the sets I i i usually overlap, in order to avoid redundant enclosures, the overlap is resolved with an overapproximative operation: all of the sets I are covered with an interval, which is repeatedly subdivided and shrinked to a given precision (see also section 3.6 in [21] for a more detailed discussion). 6.6 Intermediate Poincaré Maps The overall algorithm repeats the alternation of continuous reachability and resolution of Poincaré sections. When there is a sequence H ,..., H of intermediate Poincaré maps to 1 n be computed, it is important to ensure that, while flowing towards or resolving section H , one must not intersect with any of the H with j > i. Otherwise the later computation of the Poincaré map H might be incorrect because the actual first return time was reached before. 6.7 Derivatives For Tucker’s proof, it is necessary to compute not only the Poincaré map, but also its derivative. The derivative of the flow can be encoded as a higher dimensional ODE according to the variational equation (Theorem 7). 123 100 F. Immler n n n n∗n For an ODE with right hand side f : R → R , a new ODE of type R × R with right hand side (x , W ) → (fx , D f | · W ) is constructed. Here the first component contains the solution, and the second component its matrix of derivatives. We first extend the flows-to predicate  to a predicate  which also takes derivatives into account. Definition 10 (Flows-to Predicate Extended with Derivative) X  Y:=∀(x , d) ∈ X. ∃t ≥ 0.(φ(x , t ), Dφ| · d) ∈ Y ∧ 0 (x ,t ) bl C 0 (∀0 ≤ s ≤ t.(φ(x , s), Dφ| · d) ∈ C ) 0 (x ,s) bl X With this extended predicate for reachability, we can show that reach-cont , i.e., reach-cont for the extended ODE satisfies the specification reach-cont sctns X ≤ (λ(C, Y ). X  Y ). 0 C The Poincaré map, however requires extra care, because we cannot simply intersect the derivative of the flow with the Poincaré section: the derivative of the Poincaré map is given according to the expression in Theorem 8. For a hyperplane H ={x |x , n= c},the derivative is given as follows (for x ∈{x |x , n= c}): Dφ| · d, n (x ,τ (x )) DP| · d = Dφ| · d − f (P(x )) ϕ(t ) (x ,τ (x )) f (P(x )), n We can evaluate this expression using affine arithmetic. But we need to be able to enclose all quantities that occur on the right hand side, in particular P(x ) = ϕ(x,τ(x )) and Dφ| ·d. (x ,τ (x )) But we can enclose those: assume a step in computing an intersection, i.e., X  Y.Let us assume for simplicity that (X ∪ Y ) ∩ H =∅ and X and Y are on opposite sides of the hyperplane. Then the intersection of the flowpipe C with the section H encloses the Poincaré map: P(X ) ={ϕ(x,τ(x )) | x ∈ X}⊆ C ∩ H. For an extended flow X   Y , n∗n this means {(φ (x,τ(x )), Dφ| · d) | (x , d) ∈ X }⊆ C ∩ H × R . Therefore (x ,τ (x )) bl both P(x ) = φ(x,τ(x )) and Dφ| · d are enclosed by the result of the intersection (x ,τ (x )) n∗n C ∩ H × R for which we can use the regular intersection algorithm from Sect. 4. 6.8 Correctness Theorem We call the main algorithm that we outlined in the beginning of this Sect. 6 poincare:It resolves a sequence of intermediate Poincaré maps (together with their derivative). It is veri- fied to compute guaranteed enclosures for Poincaré maps and their derivative. The algorithm poincare takes as arguments an initial set X : R set , and initial matrix of partial derivatives n∗n n DX : R set and a target Poincaré section Σ : R set . It is further parameterized by a list of intermediate Poincaré sections, but they are irrelevant for the final correctness theorem. We formally verify partial correctness: If the algorithm returns a result, then this result encloses n∗n the Poincaré map P(x ) and its derivative DP| · DX for every x ∈ X and DX : R set . Theorem 13 (Correctness of ODE solver with Poincaré maps) poincare XDX Σ = R ⇒ ∀x ∈ X.(P(x ), DP| · DX ) ∈ R 7 Application to Lorenz Attractor In this section, we present how the verified algorithm poincare of the previous section is used to certify Tucker’s computations. We show in particular how we formally prove the Theorems 1 and 2. It helps to recall the roles of the forward invariant set N, the cone field C and the expansion estimates E in Tucker’s proof, as outlined in Sect. 1.3. 123 A Verified ODE Solver and the Lorenz Attractor 101 7.1 The Input Data and its Interpretation It is not necessary to verify precisely the set N that Tucker used, but coming up with a forward invariant set is slightly more involved than certifying one. We therefore use the output of Tucker’s program as a starting point to set up the input for our ODE solver. Since any other forward invariant with suitable cone field and expansion estimates would do just, we are free to modify Tucker’s data slightly. The output of Tucker’s program is available online as a file containing 7258 lines. We preprocessed this file by merging the information of some of the lines and slightly coarsening some of the numerical bounds. The coarsening accounts for slight differences between Tucker’s and our approximations. This results in a list of 400 elements, which we call input-data and will be the basis for all further interpretations: Definition 11 (Input Data) input-data ::result list is a list of 400 elements of type result . datatype result = Result (invoke-nf : B) − + (angle : R)(angle : R) (expansion : R)(preexpansion : R) − + − + (gridx : Z)(gridx : Z)(gridy : Z)(gridy : Z) − − + + (retx : Z)(rety : Z)(retx : Z)(rety : Z) Elements res of type result are interpreted as initial rectangles as follows. The properties − + − + gridx , gridx , gridy ,and gridy encode a rectangle on the Poincaré section Σ (recall Fig. 2), which we denote by N (res). The union of all elements of input-data represents the upper branch N of the forward invariant set N. It is plotted in Fig. 10. Definition 12 − −8 − −8 N (res):=[((gridx res − 1) · 2 ,(gridy res − 1) · 2 , 27), + −8 + −8 ((gridx res + 1) · 2 ,(gridy res + 1) · 2 , 27)] N := N (res) res∈input-data − + N :={(−x , −y, z) | (x , y, z) ∈ (N )} + − N:= N ∪ N The input data also contains information on the image of an initial rectangle. It is encoded − − + + in retx , rety , retx , rety : We select the elements within those bounds with return-of : return-of res:={res ∈ input-data | −  − + gridx res ∈[retx res, retx res]∧ −  − + gridy res ∈[rety res, rety res]} − + angle and angle define the cone C associated with the rectangle: the conic hull of the line segment between the boundary vectors. http://www2.math.uu.se/~warwick/main/rodes/ResultFile. 123 102 F. Immler + − + + Fig. 10 N in gray (N the upper and N = S(N ) the lower branch) and enclosure of P(N ) in black. This is a subset of the Poincaré section Σ (as in Fig. 2) Definition 13 C res =cone hull (segment − − (cos(rad (angle res), sin(rad (angle res), 0) + + (cos(rad (angle res)), sin(rad (angle res), 0))))) x ·π There rad x = is the radian of the angle given in degrees, segment xy is the line segment {(1 − u) · a + u · b | u ∈[0, 1]},and cone hull S ={c · x | 0 ≤ c ∧ x ∈ S} the conic hull of a set S. The elements in input-data also encode a conefield C and expansion estimates as follows. results-at (x ) yields the set of result elements that cover a point x (the rectangles overlap at −1 the boundary). We need to respect this to ensure that C, E,and E are well defined. Definition 14 results-at (x ):={res ∈ input-data | x ∈ N (res)} C(x ):= C(res) res∈results-at (x ) E (x ):= min expansion (res) res∈results-at (x ) −1 E (x ):= min preexpansion (res) res∈results-at (x ) One last property is invoke-nf , which encodes if the numerical computations need to be interrupted and the results of the normal form need to be invoked. First, we define abstractly when this is necessary, namely on the stable manifold of the origin. That is, the set of all points, which tend to the origin in infinite time. We restrict our attention to the part of the stable manifold whose trajectories do not intersect Σ for positive time. 123 A Verified ODE Solver and the Lorenz Attractor 103 Definition 15 Γ := {x |[0, ∞] ⊆ ex-ivl x ∧ (∀t > 0.φ(x , t)/ ∈ Σ) ∧ (φ (x , t ) −→ (0, 0, 0))} t →∞ When invoke-nf is true, the computations will be interrupted once the reachable sets arrive at the small cube L =[−0.1, 0.1]×[−0.1, 0.1]×[−0.1, 0.1] inside which the normal form estimates are valid. In our computations, solutions are guaranteed to enter the cube L through a rectangle T and the tangent vectors are in the cone that contains DT : ⎛ ⎞ [0.8, 1.7] ⎝ ⎠ T:= ([−0.1, 0.1], [−0.00015, 0.00015], 0.1) × [0.0005, 0.002] That is, sets are very slim in the y-direction, and the expanding direction is closely around the x axis. From Tucker’s analysis ([48, Proposition 3.1]), we devised the following bounds for the sets E , E2 (and corresponding cones in DE , DE ) through which solutions emanating 1 1 2 from T exit the cube L: ⎛ ⎞ ⎝ ⎠ E := ([−0.12, −0.088], [−0.024, 0.024], [−0.012, 0.13]) × [−0.56, 0.56] [−0.6, −0.08] ⎛ ⎞ ⎝ ⎠ E := ([0.088, 0.12], [−0.024, 0.024], [−0.012, 0.13]) × [−0.56, 0.56] [0.08, 0.6] When we interrupt computations close to L, we check that the sets entering L do so within T and continuous computations from E ∪ E . Since we have not verified Tucker’s normal 1 2 form theory, we need to trust the following assumption: Assumption 14 (Normal Form Theory Bounds) T  (E ∪ E ) 1 2 7.2 Checking the Input Data In the previous section, we only defined what the input-data encodes. Now we check if the numerical bounds prescribed by the input-data are actually correct. This involves three steps: First, we need to find a suitable setup to be able to use the algorithm poincare,which computes derivatives and not cones. Second, we set up the check that a single element of the input-data is correct. Third, we check all elements of the input-data , from which we conclude the formal counterparts of Theorems 1 and 2. 7.2.1 Representation of Cones Concerning the checking of cone conditions, first note that C res is an infinite cone, i.e., an unbounded set of vectors. In contrast to that, all of our numerical algorithms are tailored towards bounded enclosures. We therefore perform the computations with the line segment connecting the two tangent vectors with unit length. matrix-segment x y x y e encodes 1 1 2 2 a line segment (parameterized by e) in a matrix (such that it can be used as matrix initial condition DX of poincare , compare Theorem 13). mat-seg-of-deg uses this to define the line segment between the endpoints of unit vectors with given angles u,v to the x axis. A cone can therefore represented with the help of mat-seg-of-deg : 123 104 F. Immler Lemma 2 (Matrix Representation of Cone) ⎧ ⎛ ⎞ ⎫ ⎨ m ⎬ (1,1) − + ⎝ ⎠ C(res) = cone hull m m ∈ mat-seg-of-deg (angle res)(angle res) (2,1) ⎩ ⎭ with ⎛ ⎞ x + e · (x − x ) 00 1 2 1 ⎝ ⎠ matrix-segment x y x y2 e:= y + e · (y − y ) 00 1 1 2 1 2 1 mat-seg-of-deg u v:= matrix-segment (cos (rad u))(sin (rad u))(cos (rad v))(sin (rad v))[0, 1] 7.2.2 Checking a Single Result Element Algorithm 3 Check Result 1: function check-line-c1 (res) 2: X ← N (res) − + 3: DX ← mat-seg-of-deg (angle res)(angle res) 4: RES ← poincare X DX Σ 0 0 5: (P × DP ) ← split-along NRES i i i − − + + 6: RET ← get-results (retx res, rety res)(retx res, rety res) 7: return ∀i. ∃ret ∈ RET . returns-within res X DX ret i i Algorithm 3 outlines how to check that a single result element res ∈ input-data represents correct numerical bounds. It works as follows: X is the initial rectangle, DX the initial 0 0 data for the derivatives, which encodes the associated cone with angles angle res and angle res. Then the ODE solver returns with a union of return images RES, which are split along the boundaries of the individual rectangles making up N. This splitting ensures that each individual element (X , DX ) resulting from the splitting is contained in exactly i i on individual element of N. We write singleton parts of the result of this splitting X , DX . i i In RET , there are all elements of the input-data within which res is specified to return. The final check makes sure that every part X , DX of the splitting returns within one element i i ret of the collection RET . It is defined as follows and precisely formulates that X and DX, which emanate from a result res and hit the result ret, satisfy the prescribed bounds on cones and expansion. returns-within res X DX ret:= X ⊆ N (ret ) ∧ − + check-cone-bounds (angle res)(angle res) XDX ∧ −1 DX≥ E (res) ∧DX≥ E (ret ) check-cone-bounds is checked using affine arithmetic: It assumes that u and u are on the x y line segment encoding a cone according to mat-seg-of-deg , therefore checks that u = 0and ignores the other entries of the argument matrix. It further checks that the segment is on the right side (0 < u ) and that the boundary angles L and U (given in degrees) also represent a 123 A Verified ODE Solver and the Lorenz Attractor 105 cone pointing to the right side. The main purpose is in the last line, the check that the angle of the vector (u , u ) with the horizontal axis is between L and U. x y ⎛ ⎞ ⎛ ⎞ x u v w x x x ⎝ ⎠ ⎝ ⎠ check-cone-bounds LU y u v w := y y y z u v w z z z − 90 < L ∧ L ≤ U ∧ U < 90 ∧ 0 < u ∧ u = 0 ∧ x z u u y y tan(rad L) ≤ ∧ ≤ tan(rad U ) u u x x Correctness of check-line-c1 states that the set N (res) is mapped into the part return-of res of the forward invariant set. Vectors in the cone C(res) are mapped by the derivative DP into the cone field with the prescribed expansion estimates. The theorem states that the derivative exists and is defined when approaching x within Σ ={(x , y, z) | z ≤ 27}. Theorem 15 (Correctness of check-line-c1) check-line-c1 (res) = return True ⇒ ∀x ∈ N (res) − Γ. ∀dx ∈ C(res). returns-to Σ x ∧ P(x ) ∈ N (return-of res) ∧ (∃DP.(P has-derivative DP)(at x within Σ ) (DP(dx )≥ E (res) ·dx ) ∧ (∃ret ∈ return-of res. −1 P(x ) ∈ N (ret ) ∧ DP(dx ) ∈ C(ret ) ∧DP(dx )≥ E (ret ) ·dx )) The theorem follows rather directly from the definition of algorithm 3 and the specifications and definitions of the occurring functions. 7.2.3 Checking All Results We have indeed the theorem that all input-data is correct: Theorem 16 (Global Numerical Results) ∀res ∈ input-data . check-line-c1 res = return True We prove formally that under the Assumption 14, Theorem 16 implies Theorems 1 and 2, which is the main result of this article. It follows from combining the individual instances of Theorem 15 in a suitable way. Theorem 16 is proved by computing check-line-c1 (res) for every res ∈ input-data.The computations are carried out using by evaluating the statement Parallel.forall (λres. check-line-c1 res) input-data with Isabelle/HOL’s evaluation engine eval. Parallel.forall results in parallel processing of the 400 individual elements of input-data . Further parallelism is introduced when enclosures are split during reachability analysis. Split sets can be processed in parallel until they reach the next (intermediate) Poincaré section, where they might be (partially merged) upon resolving the intersection (Sect. 6.5). Figure 11 shows the plot of an enclosure for the Lorenz attractor resulting from the verified computation. The plot hints at the intermediate Poincaré sections that were manually set up 123 106 F. Immler (for some initial rectangles) at about z = 27, z = 30, x =±5, x =±1.5, x =±1, x =±0.75, x =±0.1, and z = 0.1. The black part of Fig. 10 is an enclosure for P(N ) resulting from these computations, and it is as expected and verified contained in N. The timing results of a computation on a machine with 22 cores are given below: – Threads: 22 – Elapsed Time: 6 h 33 min 9 s – CPU Time: 131 h 52 min 40 s – Parallelization Factor: 20.13 – Garbage Collection Time: 42 min 36 s To compare this with Tucker’s C++ program, I compiled Tucker’s program in a Vir- tual Machine running Ubuntu 4.20 and gcc version 3.3.4 on a machine with a 2,6 GHz Intel® Core™ i7 CPU and 16 GB RAM. Tucker’s program finished after a total CPU time of 30h and 24min. The algorithms and data structures are very different, so a direct comparison is hard. But with regard to the total CPU time (131 h) of our algorithm, a factor of less than five compared to a C++ program signifies reasonable performance for a verified algorithm. In earlier developments [20], an enclosure for the Lorenz attractor was computed with neither derivative nor cones. This earlier version verified an enclosure for the Lorenz attractor in about 7000 CPU hours. With the present version algorithms, such a computation (without derivatives and cones) can be performed in about 3 CPU hours. The speedup compared to the earlier work is mostly due to less aggressive splitting of reachable sets, and a smaller number of intermediate Poincaré sections: In the earlier work [20], intermediate Poincaré sections were introduced heuristically on-the-fly, and in the present work only where they are really effective. This is beneficial, because resolving the intersection incurs overapproximations. 8 Conclusion This article presented an overview of the diversity of algorithms, abstract results and tech- niques for the formal verification of a general-purpose ODE solver and its application to the rigorous numerical part of Tucker’s proof. 8.1 Future Work Future work would be to lift the assumption 14 on the normal form theory to formal founda- tions. This involves in particular multivariate formal power series and a number of analytic estimations, proving existence and convergence of specifically constructed formal power series. This involves more than 24 pages of Tucker’s article and is therefore a larger formal- ization effort. It also includes computer assisted parts: smalldiv.cc and coeffs.cc help devising the normal form, but neither their specification nor their implementation is particularly challenging; they essentially only evaluate a large number of fixed arithmetic expressions. Another direction for future work would be to formally conclude from the numerical results that the Lorenz equations support a robust singular hyperbolic attractor. One part would be to prove (similar to Tucker’s expansion.cc) that the expansion estimates given by E are sufficiently large to guarantee the locally eventually onto property. We validated this with a non-verified computation for the expansion estimates prescribed by our input-data.Further Intel Xeon CPU E5-2699 v4 @ 2.20 GHz. 123 A Verified ODE Solver and the Lorenz Attractor 107 + + Fig. 11 Enclosure of φ{(x , [0; τ(x )]) | x ∈ N }, during the computation of P(N ) mathematical foundations are required in order to conclude from the computed forward invariant conefield C and the expansion bounds that there exists a stable foliation, and that this foliation can be used to reduce the two-dimensional dynamics on Σ to a one-dimensional map. This requires in particular the formalization of differentiable manifolds, and theorems like the existence of differentiable invariant sections for fiber contractions [41]. 8.2 Size of Code Table 1 shows some statistics on the size in terms of lines of code of several programs related to this verification. RODES is the rigorous ODE solver used by Tucker, it consists of 3800 lines of C++ code and builds on a library for interval arithmetic (Profil/BIAS) of about twice the size. Similar to the sum of those two is the size of the generated SML code. The verification required more effort, but the largest part is generic: the part specific to the Lorenz attractor makes up less than 10% of the total number of lines of code. 8.3 Trust Base We use the evaluation oracle eval in Isabelle/HOL. This is common practice to speed up rewriting. Isabelle/HOL equations are mapped to function definitions in Standard ML. These are compiled and evaluated in PolyML. We also use the common extension http://polyml.org/. 123 108 F. Immler Table 1 Size of code and formalization Sections Language Lines of code/proof RODES C++ 3800 Profil/BIAS C++ 8852 generated ODE solver SML 13200 Flow, Poincaré map 2 Isabelle/HOL theory 12000–16500 Affine Arithmetic 3 8500 Intersection 4 5000 Refinement/Enclosures 3 5000 Reachability Analysis 6 10000 Lorenz Attractor 7 3000 (HOL-Library.Code_Target_Numeral in Isabelle2017) of the code generator that maps Isabelle/HOL integers to the integer type IntInf.int of PolyML, which can be based on either PolyML’s own implementation of arbitrary precision integers or GMP [8]. This setup means that the trusted code base is increased: The translation of Isabelle/HOL terms to SML code is not verified. One needs to trust PolyML and its compiler, but PolyML is Isabelle’s implementation language and therefore anyhow part of the trusted code base. Reducing the trusted code base is an orthogonal issue: there is ongoing work [16,17]for verified code generation to CakeML [28], a verified implementation of ML. Isabelle/HOL’s eval speeds up evaluation by translating terms to the implementation language of the proof checker (PolyML). In view of this, it is more similar to Coq’s native_compute [2], which evaluates terms after translation to Coq’s implementation language OCaml, than to Coq’s virtual machine [9]. 8.4 Related Work Integrals and Differential Equations in Proof Assistants Spitters and Makarov [33]imple- ment Picard iteration to calculate solutions of ODEs in the interactive theorem prover Coq, but restricted to relatively short existence intervals. Boldo et al. [3] approximate the solution of one particular partial differential equation with a C-program and verify its correctness in Coq. Mahboubi and Sibut-Pinote [32] compute rigorous approximations of integrals with Taylor models. Rigorous Numerics in Proof Assistants Rigorous numerical approximation of arithmetic expressions has been done in Coq [34] for different types of enclosures (Taylor models, inter- vals, centered forms). Muñoz and Lester [37] use rational interval arithmetic in PVS. Rigorous numerics with first order interval approximations has been implemented by Solovyev for the Flyspeck project [44] in HOL Light. This work is remarkable in that it is not relying on code generation but uses only primitive inference rules of HOL Light’s kernel. Computational Geometry Pichardie and Bertot [40] were the first to formalize the ccw system and verify a functional convex hull algorithm in Coq. Meikle and Fleuriot [35] formalized an imperative algorithm and verified it using Hoare logic in Isabelle/HOL. Brun et al. [5]verify an algorithm based on hypermaps to compute the convex hull. 123 A Verified ODE Solver and the Lorenz Attractor 109 Acknowledgements I would like to thank Jeremy Avigad for drawing my attention to this particular appli- cation of rigorous ODE solving. I am very grateful to Tobias Nipkow for encouraging and supporting me to pursue this multi-faceted project. The anonymous reviewers provided helpful and constructive feedback. This research was financially supported by DFG RTG 1480 (PUMA) and DFG Koselleck Grant NI 491/16-1. I would like to thank Johannes Hölzl for supervising my student project on the Picard–Lindelöf theorem. I thankfully acknowledge Christoph Traut’s work for his student project on formalizing the variational equation. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. Back, R.J., Wright, J.: Refinement Calculus: A Systematic Introduction. Springer, New York (2012) 2. Boespflug, M., Dénès, M., Grégoire, B.: Full reduction at full throttle. In: International Conference on Certified Programs and Proofs, pp. 362–377. Springer, New York (2011) 3. Boldo, S., Clément, F., Filliâtre, J.C., Mayero, M., Melquiond, G., Weis, P.: Wave equation numerical resolution: a comprehensive mechanized proof of a C program. Journal of Automated Reasoning 50(4), 423–456 (2013). https://doi.org/10.1007/s10817-012-9255-4 4. Bouissou, O., Chapoutot, A., Djoudi, A.: Enclosing temporal evolution of dynamical systems using numerical methods. In: Brat, G., Rungta, N., Venet, A. (eds.) NASA Formal Methods, LNCS, vol. 7871, pp. 108–123. Springer, New York (2013). https://doi.org/10.1007/978-3-642-38088-4_8 5. Brun, C., Dufourd, J.F., Magaud, N.: Designing and proving correct a convex hull algorithm with hyper- maps in Coq. Computational Geometry 45(8), 436–457 (2012). https://doi.org/10.1016/j.comgeo.2010. 06.006.(Geometric Constraints and Reasoning) 6. de Figueiredo, L., Stolfi, J.: Affine arithmetic: concepts and applications. Numerical Algorithms 37(1–4), 147–158 (2004) 7. Girard, A., Le Guernic, C.: Zonotope/hyperplane intersection for hybrid systems reachability analysis. In: Egerstedt, M., Mishra, B. (eds.) Hybrid Systems: Computation and Control, LNCS, vol. 4981, pp. 215–228. Springer, New York (2008). https://doi.org/10.1007/978-3-540-78929-1_16 8. Granlund, T., et al.: GNU MP 6.0 Multiple Precision Arithmetic Library. Samurai Media Limited, Hong Kong (2015) 9. Grégoire, B., Leroy, X.: A compiled implementation of strong reduction. ACM SIGPLAN Not. 37(9), 235–246 (2002) 10. Haftmann, F., Krauss, A., Kuncar ˇ , O., Nipkow, T.: Data Refinement in Isabelle/HOL, pp. 100–115. Springer, Berlin (2013). https://doi.org/10.1007/978-3-642-39634-2_10 11. Haftmann, F., Wenzel, M.: Constructive type classes in Isabelle. In: International Workshop on Types for Proofs and Programs, pp. 160–174. Springer, New York (2006) 12. Harrison, J.: A HOL theory of Euclidean space. In: Hurd, J., Melham, T. (eds.) 18th International Confer- ence on Theorem Proving in Higher Order Logics (TPHOLs 2005), Lecture Notes in Computer Science, vol. 3603, pp. 114–129 (2005) 13. Hölzl, J.: Proving inequalities over reals with computation in Isabelle/HOL. In: Reis, G.D. Théry L. (eds.) Programming Languages for Mechanized Mathematics Systems (ACM SIGSAM PLMMS’09), pp. 38–45 (2009) 14. Hölzl, J., Immler, F., Huffman, B.: Type classes and filters for mathematical analysis in Isabelle/HOL. In: Blazy, S., Paulin-Mohring, C., Pichardie, D. (eds.) ITP 2013, pp. 279–294. Springer, New York (2013) 15. Huffman, B., Kuncar ˇ , O.: Lifting and transfer: a modular design for quotients in Isabelle/HOL. In: Inter- national Conference on Certified Programs and Proofs, pp. 131–146. Springer, New York (2013) 16. Hupel, L.: Dictionary Construction. Archive of Formal Proofs (2017). http://isa-afp.org/entries/Dict_ Construction.html, Formal proof development 17. Hupel, L., Nipkow, T.: A verified compiler from Isabelle/HOL to CakeML. Springer LNCS (2018). https:// lars.hupel.info/pub/isabelle-cakeml.pdf 18. Immler, F.: Formally verified computation of enclosures of solutions of ordinary differential equations. In: Badger, J.M., Rozier, K.Y. (eds.) NFM 2014, pp. 113–127. Springer, New York (2014) 123 110 F. Immler 19. Immler, F.: A verified algorithm for geometric zonotope/hyperplane intersection. In: CPP 2015, pp. 129– 136. ACM (2015) 20. Immler, F.: A verified enclosure for the Lorenz attractor (rough diamond). In: Urban, C., Zhang, X. (eds.) ITP 2015, pp. 221–226. Springer, Berlin (2015) 21. Immler, F.: Verified reachability analysis of continuous systems. In: Baier, C., Tinelli, C. (eds.) TACAS 2015, pp. 37–51. Springer, Berlin (2015) 22. Immler, F.: Affine Arithmetic. Archive of Formal Proofs (2017). http://isa-afp.org/entries/Affine_ Arithmetic.shtml, Formal proof development 23. Immler, F., Hölzl, J.: Ordinary Differential Equations. Archive of Formal Proofs (2017). http://isa-afp. org/entries/Ordinary_Differential_Equations.shtml, Formal proof development 24. Immler, F., Hölzl, J.: Numerical analysis of ordinary differential equations in Isabelle/HOL. In: Interna- tional Conference on Interactive Theorem Proving, pp. 377–392. Springer, New York (2012) 25. Immler, F., Traut, C.: The flow of ODEs: Formalization of variational equation and Poincaré map. J. Autom. Reason. (2018). https://doi.org/10.1007/s10817-018-9449-5 26. Immler, F., Traut, C.: The flow of ODEs. In: Blanchette, J.C., Merz, S. (eds.) ITP 2016, pp. 184–199. Springer, Cham (2016) 27. Knuth, D.: Axioms and Hulls. Springer, Berlin (1992). Number 606 in Lecture Notes in Computer Science 28. Kumar, R., Myreen, M.O., Norrish, M., Owens, S.: CakeML: a verified implementation of ML. In: ACM SIGPLAN Notices, vol. 49, pp. 179–191. ACM (2014) 29. Lammich, P.: Refinement for monadic programs. Archive of Formal Proofs (2012). http://afp.sf.net/ entries/Refine_Monadic.shtml, Formal proof development 30. Lammich, P.: Automatic data refinement. In: International Conference on Interactive Theorem Proving, pp. 84–99. Springer, Heidelberg (2013) 31. Lorenz, E.N.: Deterministic nonperiodic flow. J. Atmos. Sci. 20(2), 130–141 (1963). https://doi.org/10. 1175/1520-0469(1963)020<0130:DNF>2.0.CO;2 32. Mahboubi, A., Melquiond, G., Sibut-Pinote, T.: Formally Verified Approximations of Definite Integrals, pp. 274–289. Springer, Cham (2016). https://doi.org/10.1007/978-3-319-43144-4_17 33. Makarov, E., Spitters, B.: The Picard algorithm for ordinary differential equations in Coq. In: Blazy, S., Paulin-Mohring, C., Pichardie, D. (eds.) Interactive Theorem Proving, LNCS, vol. 7998, pp. 463–468. Springer, Berlin (2013) 34. Martin-Dorel, E., Melquiond, G.: Proving tight bounds on univariate expressions with elementary func- tions in Coq. J. Autom. Reason. 57(3), 187–217 (2016). https://doi.org/10.1007/s10817-015-9350-4 35. Meikle, L., Fleuriot, J.: Mechanical theorem proving in computational geometry. In: Hong, H., Wang, D. (eds.) Automated Deduction in Geometry, Lecture Notes in Computer Science, pp. 1–18. Springer, Berlin (2006). https://doi.org/10.1007/11615798_1 36. Morales, C., Pacifico, M., Pujals, E.: Singular hyperbolic systems. Proc. Am. Math. Soc. 127(11), 3393– 3401 (1999) 37. Muñoz, C., Lester, D.: Real number calculations and theorem proving. In: Hurd, J., Melham, T. (eds.) Theorem Proving in Higher Order Logics (TPHOLs 2005), LNCS, pp. 195–210 (2005) 38. Nipkow, T.: Order-sorted polymorphism in Isabelle. In: Logical Environments, pp. 164–188. Cambridge University Press (1993) 39. Nipkow, T., Paulson, L.C., Wenzel, M.: Isabelle/HOL: A Proof Assistant for Higher-Order Logic, LNCS. Springer, Heidelberg (2002) 40. Pichardie, D., Bertot, Y.: Formalizing convex hull algorithms. In: Boulton, R., Jackson, P. (eds.) Theorem Proving in Higher Order Logics, Lecture Notes in Computer Science, vol. 2152, pp. 346–361. Springer, Berlin (2001). https://doi.org/10.1007/3-540-44755-5_24 41. Robinson, C.: Dynamical Systems—Stability, Symbolic Dynamics, and Chaos. CRC Press, Boca Raton (1999). https://doi.org/10.1007/978-1-4613-0003-8 42. Rump, S.M., Kashiwagi, M.: Implementation and improvements of affine arithmetic. Nonlinear Theory Appl. IEICE 6(3), 341–359 (2015). https://doi.org/10.1587/nolta.6.341 43. Smale, S.: Mathematical problems for the next century. Math. Intell. 20(2), 7–15 (1998). https://doi.org/ 10.1007/BF03025291 44. Solovyev, A., Hales, T.C.: Formal verification of nonlinear inequalities with Taylor interval approxima- tions. In: NASA Formal Methods Symposium, pp. 383–397. Springer, New York (2013) 45. Sparrow, C.: The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors. No. 41 in Applied Mathematical Sciences. Springer. New York. https://doi.org/10.1007/978-1-4612-5767-7 46. Tucker, W.: My thesis: The lorenz attractor exists. http://www2.math.uu.se/~warwick/main/pre_thesis. html 47. Tucker, W.: The Lorenz attractor exists. Comptes Rendus de l’Académie des Sciences-Series I- Mathematics 328(12), 1197–1202 (1999) 123 A Verified ODE Solver and the Lorenz Attractor 111 48. Tucker, W.: A rigorous ODE solver and Smale’s 14th problem. Found. Comput. Math. 2(1), 53–117 (2002) 49. Viana, M.: What’s new on Lorenz strange attractors? Math. Intell. 22(3), 6–19 (2000)

Journal

Journal of Automated ReasoningSpringer Journals

Published: Jan 22, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off