Numerical mathematics is concerned with the representation of general mathematical objects (e.g. the solution of some differential equation) in terms of a limited collection of simpler objects like the basic arithmetic functions on R. Traditionally, the correctness of such representations is rigorously investigated but when it comes to the implementation, error-prone representation changes to data structures are not checked in a formal manner. To avoid this problem, mathematical models of data structures and algorithms are required to enable rigorous arguments down to the level of actual implementation steps. We have carried out this process for various problems in the area of differential equations. In order to keep proofs about the correctness of implementations manageable, intermediate data structures with high reusability are constructed.
In this article, we use the word model to denote a description of a system (i.e. a collection of objects with
mutual relations) in a formal language. For example, a mathematical model of the solar system S describes
the objects (the sun and planets) in terms of scalar (mass) parameters and vector valued functions of a scalar
(time) variable for positions and velocities. The mutual relations are then coded as a set of ordinary differential
equations (Newton’s law with gravitational forces). From a mathematical point of view, the model (call it Mphy)
appears as a theory which, according to [
In this context, an approximation of the ordinary differential equation by a discrete equation can be considered as a model (call it Mnum) of Mphy. Differences between Mphy and Mnum are due to additional (discretization) parameters, changes in the equations which relate the objects and possibly also in the object representation (the domain of the (time) variable may be restricted to a finite set, for example).
Finally, if a solver for the discrete equations is implemented, the resulting program P is also a model, because it describes Mnum in a formal language: The mathematical objects are represented by appropriate data structures and the mutual relations are coded in the solution algorithm. The resulting cascade of models is depicted in Figure 1.
Since Mphy involves various abstractions and approximations of S, modeling errors can occur which means that predictions generated from Mphy differ from corresponding observations in S. The modeling error between Mnum and Mphy is traditionally denoted approximation error and we would call Mnum a reasonable model, if the approximation error converges to zero under suitable assumptions. Finally, modelling errors between P and Mnum may be due to neglected integer overflow or round-off errors (if mathematical integers and reals are modeled by representations of finite bit length), or due to approximations of transcendental by rational functions. Apart from these technical errors also undetected typos and logical mistakes may be part of this modeling error. 3
In the modeling process described above, the best error control is obtained between Mnum and Mphy when the convergence proof is verified by sufficiently many independent mathematicians or by some trusted proof checker. Besides excluding logical mistakes, the proof also serves as a very careful explanation, how Mnum models Mphy.
For the link between Mphy and S, such an explanation is generally not possible, because the language describing S is not formal enough to support strict arguments. Nevertheless, the chance of logical errors may be reduced if the concepts used to describe S are closely matched by notions available in the mathematical language so that the translation is simplified. Conversely, the bigger the difference in vocabulary, the more (uncontrolled) reasoning is required to explain the relation between S and Mphy which impedes comprehensibility and increases the chance of mistakes (in view of this, one educational objective in our modeling classes is the design of comprehensible models: First, mathematical notions are defined which translate key words of the problem – then the problem is reformulated using these terms).
The same idea applied to the link between Mnum and P calls for some library with data structures and
algorithms adapted to the counterparts in Mnum. The use of such a library may turn the program into a
comprehensible model of Mnum and thus reduce the chance of mistakes. To increase confidence in the correctness
of P , time consuming tests may be carried out. However, they are often incomplete and therefore much less
satisfactory [
In numerical math courses, algorithms are generally presented in some pseudo-code which combines mathematical notation with a small selection of control structures common to most programming languages. Data types like int, real or list are identified with Z; R and Cartesian powers respectively. In this way, restriction to a single language is avoided and notation stays closer to common mathematical usage.
We pick up the latter approach by first describing certain aspects of (a family of) computer languages L in terms of a mathematical model Mcmp (here the objects are data structures and algorithms). Then (see figure 2), we reformulate Mnum using these abstract objects and prove mathematically, how the resulting model Malg corresponds to Mnum. Finally, if Malg-expressions, which are based solely on Mcmp objects, are automatically translated into some L-coherent programming language, even translation errors can be avoided.
Since the modeling process is very flexible, we can include high level algorithms into Mcmp which simplifies the formulation of Malg but reduces mathematical control. Conversely, rigor and effort in the formulation of Malg is increased if the model is based on low level algorithms. Similarly, the model assumptions on the data structures can range from IEEE 754 rules to a mere identification of computer and mathematical reals if round-off errors are considered less critical.
We have worked out an exemplary version of Mcmp which starts with the fundamental object rawData representing a set whose elements are interpreted as data objects in L.
In order to model the binding of raw data to mathematical concepts, we could postulate a mapping from rawData to objects defined in MATh-theories. However, such a global binding map would cause some difficulty: Since the development of numerical algorithms generally comes along with the definition of new data structures for relevant mathematical objects, new axioms on the binding map are required during the development. In contrast to usual mathematical practice, it would be impossible to collect all axioms of Mcmp before drawing conclusions.
To avoid this, we favor a less global approach and introduce representations which bind rawData elements to mathematical meaning in records with a data field for the raw data and an obj field for its mathematical interpretation. Formally, the MATh-definition of this notion would be part of Mcmp.
representation := r : record with dom(r)
f’data’; ’obj’g; r:data : rawData Collections of representations whose mathematical interpretations range in some set X with no rawData object connected to more than one interpretation will be called data types of X. More specifically, dataType(X) := T with T
representation; (x 2 T 7! x:data) : injective; 8x 2 T holds x:obj 2 X Using this notion, the abstract Mcmp-parameter real is introduced with the axiom real : dataType(R). This assumption does not entail that each real number is represented by some raw data which would, however, follow from the axiom objRange(real) = R. Here objRange(T ) is defined as the set of all objects corresponding to the elements of the dataType T .
To describe functions from rawData to rawData which are realizable in L, the notion rawMapping is introduced as a recursively defined set: It contains identity and constant maps as explicit elements and certain rules allow the creation of new raw mappings from given ones (e.g. by composition). An algorithm is then defined as a representation-valued mapping on some dataType with the property that the induced map between the data fields of arguments and images is a rawMapping. In this way, the recursive construction principle is passed on from raw mappings to algorithms.
While algorithms give rise to raw mappings when representations are restricted to their data fields, they may also be considered as functions between mathematical interpretations. This requires, however, that representations of the same object are mapped to a common mathematical interpretation. In this case, we speak of algorithmic representations
algorithmicRepresentation(f ) := f : algorithm with 8x 2 dom(f) holds f(x):obj = f (x:obj) Altogether, there are three levels on which an algorithmic representation can operate. On the level of rawDataobjects as a rawMapping, on the level of representation-objects as an algorithm and on the level of purely mathematical objects as a mathematical function.
In our model, algorithmic representations are used to relate abstract algorithms like plus on real to their mathematical counterparts (e.g. with the axiom plus : algorithmicRepresentation(plus) the relation to the addition function plus on R is established). Moreover, in Malg, certain algorithms are proved to be algorithmic representations of functions in Mnum which serves as careful explanation of the relation between Malg and Mnum.
When modeling compound data types of L like lists or records, one could either follow the approach above or use related types in MATh directly. For example, a simple way to axiomatize record formation in L is to consider a MATh-record r as rawData element if and only if the images of r (i.e. the indexed values) are rawData.
While such a close relation between L- and MATh-records alleviates the use in MATh expressions, it may complicate automatic translations to L-compatible languages. This effect is more obvious in the case of lists because the integer indices are distinguished objects in L and MATh and 1 as the minimal list index in MATh may differ from the choice in other languages. For this reason, we have modeled L-lists without resorting to the MATh counterparts directly. Instead, list has been introduced as an abstract set of representations in Mcmp where each element x of list represents a MATh-list (which is a function on a set f1; : : : ; ng). In order to enable list usage, additional L objects are required like length which maps each x 2 list to a non-negative int and eval which connects x with an algorithmic representation of x.obj. Finally, the creation of a list from a finite number of representations requires the abstract object explicitList. The precise interplay of these objects and other fundamental list operations like joining, composing, repeating, transposing or flattening are described with appropriate axioms. 5
Our experience with algorithms for the approximate solution of ordinary and partial differential equations in
bachelor and master courses is that linear algebra aspects like setting up coordinate vectors, functionals and
matrices turn out to be primary sources for logical errors. This is particularly the case when ad-hoc
representations of finite difference grids, finite element meshes, finite dimensional vector spaces as well as linear, non-linear
and discrete functions are formulated in terms of simple list structures. In comparison, round-off errors are
generally less critical so that details of floating point arithmetic have not been introduced in our model. Instead,
we identified int and real with Z and R (see [
To give an example, we consider an equidistant rectangular grid with distance h > 0 between nodes in both space directions as it may appear in a finite difference discretization. In this case, natural indices are integer pairs (i; j) 2 Z2 to refer to the grid nodes h (i; j) 2 G and the vector space of grid functions G ! R (also denoted RG) plays a central role. A useful basis of its dual space consists of the evaluation functionals at the grid nodes so that the natural indexing is also in terms of pairs (i; j). Finite difference operators are simple linear combinations of these basis functionals which typically involve indices of geometrically close nodes.
Since basic arrays in programming languages are indexed by f1; : : : ; ng or f0; : : : ; n 1g respectively, a bijective index-mapping (resulting, for example, from row- or column-wise node enumeration) has to be employed which renders the relevant expressions more technical and error-prone. This situation is aggrevated if systems of equations with time dependence are considered. In fluid dynamics, for instance, the unknowns may consist of a velocity vector and a pressure value at each node and at each time point so that the relevant function space is ((R2 R)G)T (here T is the temporal grid). To realize elements of such structured spaces, we do not translate them directly to lists but introduce reusable intermediate data structures and algorithms: representations of finite sets with operations like product, intersection or union, representations of functions between finite sets with evaluation and composition, representations of functions from finite sets into general data types and representations of vector space operations on arbitrary finite products of real.
Establishing all these structures and proving their relationship to the mathematical counterparts amounts to the construction of a generalized software library with a very precise documentation and high reusability. The actual formulation of the target algorithms uses this library and the associated mathematical theory alleviates the proof of the theorems which connect the algorithms to the original mathematical problems. Apart from teaching code generation which is particularly amenable for program verification, we are also planning a tool which creates code automatically from logically verified MATh-algorithms in different programming languages. So far, we have tested this translation process manually with Scala implementations of the algorithms.