A Constraint Satisfaction Approach to Parametric Differential Equations

1 Introduction

Initial value problems for ordinary differential equations arise naturally in many applications in science and engineering, including chemistry, physics, molecular biology, and mechanics to name only a few. An ordinary differential equation (ODE)  is a system of the form

often denoted in vector notation by D9 BC B4D8B5 BP CUB4D8BN D9B4D8B5B5 or D9 BC BP CUB4D8BN D9B5. In addition, in practice, it is often the case that the parameters and/or the initial values are not known with certainty but are given as intervals. Hence traditional methods do not apply to the resulting parametric ordinary differential equations since they would have to solve infinitely many systems. Interval methods, pioneered by Moore [Moo79], provide an approach to tackle parametric ODEs as they return reliable enclosures of the exact solution at different points in time. Interval methods typically apply a one-step Taylor interval method and make extensive use of automatic differentiation to obtain the Taylor coefficients [Moo79]. Their major problem however is the explosion of the size of the boxes at successive points as they often accumulate errors from point to point and lose accuracy by enclosing the solution by a box (this is called the wrapping effect). Lohner's AWA system [Loh87] was an important step in interval methods Yves Deville INGI-UCL Louvain-La-Neuve, Belgium which features coordinate transformations to tackle the wrapping effect. More recently, Nedialkov and Jackson's IHO method [NJ99] improved on AWA by extending a HermiteObreschkoff's approach (which can be viewed as a generalized Taylor method) to intervals. Note that interval methods inherently accommodate uncertain data. Hence, in this paper, we talk about ODEs to denote both traditional and parametric ODEs.

This research takes a constraint satisfaction approach to ODEs. Its basic idea [DJVH98; JDVH99] is to view the solving of ODEs as the iteration of two processes: (1) a forward process that computes initial enclosures at given times from enclosures at previous times and bounding boxes and (2) a pruning process that reduces the initial enclosures without removing solutions. The real novelty in our approach is the pruning component. Pruning in ODEs however generates significant challenges since ODEs contain unknown functions.

The main contribution of this paper is to show that an effective pruning technique can be derived from a relaxation of the ODE, importing a fundamental principle from constraint satisfaction into the field of differential equations. Three main steps are necessary to derive an effective pruning algorithm. The first step consists in obtaining a relaxation of the ODE by safely approximating its solution using, e.g., generalized Hermite interpolation polynomials. The second step consists in using the mean-value form of this relaxation to solve the relaxation accurately and efficiently. Unfortunately, these two steps, which were skeched in [JDVH99], are not sufficient and the resulting pruning algorithm still suffers from traditional problems of interval methods. The third fundamental step consists in globalizing the pruning by considering several successive relaxations together. This idea of generating a global constraint from a set of more primitive constraints is also at the heart of constraint satisfaction. It makes it possible, in this new context, to address the problem of dependencies (and hence the accumulation of errors) and the wrapping effect simultaneously. 1

Theoretical and experimental results show the benefits of the approach. The theoretical results show that the pruning step can always be implemented to make our approach faster than existing methods. In addition, it shows that our approach should be significantly faster when the function CU is very complex. Experimental results confirm the theory. They show that the approach often produces many orders of magnitude improvements in accuracy over existing methods while not increasing computation times. Alternatively, at similar accuracy, other approaches are significantly slower. Of particular interest is the fact that our approach is not dependent on high-order Taylor coefficients contrary to other methods.

1 Global constraints in ordinary differential equations have also been found useful in [CB99]. The problem and the techniques in [CB99] are however fundamentally different.

The rest of the paper is organized as follows. Section 2 introduces the main definitions and notations. Sections 3, 4, and 5 describe the pruning component. Sections 7 and 8 present the theoretical and experimental analyses.

2 Background and Definitions

Basic Notational Conventions Small letters denote real values, vectors and functions of real values. Capital letters denote real matrices, sets, intervals, vectors and functions of intervals.
CA denotes the set of all closed intervals AI CA whose bounds are floating-point numbers. A vector of intervals BW BE
CA D2 is called a box. If D6 BE CA, then D6 denotes the smallest interval
BE
CA such that D6 BE
. If D6 BE CA D2 , then D6 BP B4D6 BD BN BM BM BM BN D6 D2 B5. In this paper, we often use D6 instead of D6 for simplicity. If BT AI CA D2 , then BEBT denotes the smallest box BW BE
CA D2 such that BT AI BW. We also assume that D8 BC BN BM BM BM BN D8 CZ , D8 CT and D8 are reals, D9 BC BN BM BM BM BN D9 CZ are in CA D2 , and BW BC BN BM BM BM BNBW CZ are in
CA D2 . Finally, we use D8 CZ to denote CWD8 BC BN BM BM BM BN D8 CZ CX, D9 CZ to denote CWD9 BC BN BM BM BM BN D9 CZ CX and BW CZ to denote CWBW BC BN BM BM BM BNBW CZ CX.

We restrict attention to ODEs that have a unique solution for a given initial value. Techniques to verify this hypothesis numerically are well-known [Moo79; DJVH98]. Moreover, in practice, the objective is to produce (an approximation of) the values of the solution of  at different points D8 BC BN D8 BD BN BM BM BM BN D8 D1 . This motivates the following definition of solutions and its generalization to multistep solutions.

Definition 1 (Solution of an ODE) The solution of an ODE  is the function D7 BM CA A2 CA D2 AX CA AX CA D2 such that D7B4D8 BC BN D9 BC B5 satisfies  for an initial condition D7B4D8 BC BN D9 BC B5B4D8 BC B5 BP D9 BC .

Definition 2 (Multistep solution of an ODE) The multistep solution of an ODE  is the partial function D1D7 BM BT AI CA CZB7BD A2 B4CA D2 B5 CZB7BD AX CA AX CA D2

where D7 is the solution of  and is undefined otherwise.

Since multistep solutions are partial functions, we generalize the definition of interval extensions to partial functions.

Definition 3 (Interval Extension of a Partial Function)

The interval function BZ BM
CA D2 AX
CA D1 is an interval extension of the partial function CV BM BX AI CA D2 AX CA D1 if

where CVB4BTB5 BP CUCVB4DCB5 CY DC BE BTCV.

Finally, we generalize the concept of bounding boxes, a fundamental concept in interval methods for ODEs, to multistep methods. Intuitively, a bounding box encloses all solutions of an ODE going through certain boxes at given times.

Definition 4 (Bounding Box) Let  be an ODE system, D1D7 be the multistep solution of , and CUD8 BC BN BM BM BM BN D8 CZ CV AI CC BE
CA. A box BU is a bounding box of D1D7 over CC wrt (D8 CZ ,BW CZ ) if, for all D8 BE CC , D1D7B4D8 CZ BN BW CZ B5B4D8B5 AI BU.

3 The Constraint Satisfaction Approach

The constraint satisfaction approach followed in this paper was first presented in [DJVH98]. It consists of a generic algorithm for ODEs that iterates two processes: (1) a forward process that computes initial enclosures at given times from enclosures at previous times and bounding boxes and (2) a pruning process that reduces the initial enclosures without removing solutions. The forward process also provides numerical proofs of existence and uniqueness of the solution. The intuition of the successive integration steps is illustrated in Figure 1. Forward components are standard in interval methods for ODEs. This paper thus focuses on the pruning process, the main novelty of the approach. Our pruning component is based on relaxations of the ODE as described in the next section. To our knowledge, no other approach uses relaxations of the ODE to derive pruning operators and the only other approach using a pruning component [NJ99; Rih98] was developed independently.

4 Pruning

The pruning component uses safe approximations of the ODE to shrink the boxes produced by the forward process. To understand this idea, it is useful to contrast the constraint satisfaction to nonlinear programming [VHMD97] and to ordinary differential equations. In nonlinear programming, a constraint CRB4DC BD BN BM BM BM BN DC D2 B5 can be used almost directly for pruning the search space (i.e., the Cartesian products of the intervals
CX associated with the variables DC CX ). It suffices to take an interval extension BVB4CG BD BN BM BM BM BNCG D2 B5 of the constraint. Now if BVB4
BC BD BN BM BM BM BN
BC D2 B5 does not hold, it follows, by definition of interval extensions, that no solution of CR lies in
BC BD A2 BM BM BM A2
BC D2 . The interval extension can be seen as a filter that can be used for pruning the search space in many ways. For instance, Numerica uses box(CZ)-consistency on these interval constraints [VHMD97]. Ordinary differential equations raise new challenges. In an ODE BK D8 BM D9 BC BP CUB4D8BN D9B5BN functions u and u' are, of course, unknown. Hence it is not obvious how to obtain a filter to prune boxes.

One of the main contributions of our approach is to show how to derive effective pruning operators for parametric ODEs. The first step consists in rewriting the ODE in terms of its multistep solution ms to obtain

Let us denote this formula BK D8 BM CUD0B4D8 CZ BN D9 CZ B5B4D8B5BM This rewriting may not appear useful since ms is still an unknown function. However it suggests a way to approximate the ODE. Indeed, we show in Section 6 how to obtain interval extensions of ms and BSD1D7 BSD8 by using polynomial interpolations together with their error terms. This simply requires a bounding box for the considered time interval and safe approximations of ms at successive times, both of which are available from the forward process. Once these interval extensions are available, it is possible to obtain an interval formula of the form

which approximates the original ODE. The above formula is still not ready to be used as a filter because D8 is universally quantified. The solution here is simpler and consists of restricting attention to a finite set CC of times to obtain the relation

which produces a computable filter. Indeed, if the relation BYB4D8 CZ BNBW CZ B5B4D8B5 does not hold for a time D8, it follows that no solution of D9 BC BP CUB4D8BN D9B5 can go through boxes BW BC BN BM BM BM BNBW CZ at times D8 BC BN BM BM BM BN D8 CZ . The following definition and proposition capture these concepts more formally.

Definition 5 (Multistep Filters) Let  be an ODE and s its solution. A multistep filter for  is an interval relation BY BM CA CZB7BD A2 B4
CA D2 B5 CZB7BD AX CA AX BUD3D3D0 satisfying

Proposition 1 (Soundness of Multistep Filters) Let  be an ODE and let BY be a multistep filter for . If BYB4D8 CZ BNBW CZ B5B4D8B5 does not hold for some D8, then there exists no solution of  going through BW CZ at times D8 CZ .

How can we use this filter to obtain tighter enclosures of the solution? A simple technique consists of pruning the last box produced by the forward process. Assume that BW CX is a box enclosing the solution at time D8 CX (BC AK CX AK CZ) and that we are interested in pruning the last box BW CZ . A subbox BW AI BW CZ can be pruned away if the condition

does not hold for some evaluation point D8 CT . Let us explain briefly the geometric intuition behind this formula by considering what we call natural filters. Given interval extensions MS and DMS of ms and BSD1D7 BSD8 , it is possible to approximate the ODE D9 BC BP CUB4D8BN D9B5 by the formula

In this formula, the left-hand side of the equation represents the approximation of the slope of D9 while the right-hand represents the slope of the approximation of D9. Since the approximations are conservative, these two sides must intersect on boxes containing a solution. Hence an empty intersection means that the boxes used in the formula do not contain the solution to the ODE system. Figure 2 illustrates the intuition. It is generated from an actual ordinary differential equation, considers only points instead of intervals, and ignores error terms for simplicity. It illustrates how this technique can prune away a value as a potential solution at a given time. In the figure, we consider the solution to the equation that evaluates to D9 BC and D9 BD at D8 BC and D8 BD respectively. Two possible points D9 BE and D9 BC BE are then considered as possible values at D8 BE . The curve marked KO describes an interpolation polynomial going through D9 BC BN D9 BD BN D9 BC BE at times D8 BC BN D8 BD BN D8 BE . To determine if D9 BC BE is the value of the solution at time D8 BE , the idea is to test if the equation is satisfied at time D8 CT . (We will say more about how to choose D8 CT later in this paper). As can be seen easily, the slope of the interpolation polynomial is different from the slope specified by CU at time D8 CT and hence D9 BC BE cannot be the value of the solution at D8 BE . The curve marked OK describes an interpolation polynomial going through D9 BC BN D9 BD BN D9 BE at times D8 BC BN D8 BD BN D8 BE . In this case, the equation is satisfied at time D8 CT , which means that D9 BE cannot be pruned away. The filter proposed earlier generalizes this intuition to boxes. Both the left- and the right-hand sides represent sets of slopes and the filter fails when their intersection is empty. Traditional consistency techniques and algorithms based on this filter can now be applied. For instance, one may be interested in updating the last box produced by the forward process using the operator BW CZ BP BECUD6 BE BW CZ CY BYB4D8 CZ BN CWBW BC BN BM BM BM BNBW CZA0BD BN D6CXB5B4D8 CT B5CVBM The following definition is a novel notion of consistency for ODEs to capture pruning of the last D6 boxes.

Definition 6 (Forward Consistency of Multistep Filters)

Let FL be a multistep filter. FL is forward(D6)-consistent wrt D8 CZ , BW CZ and D8 if

The algorithm used in our computational results enforces forward consistency of the filters we now describe.

5 Filters

Filters rely on interval extensions of the multistep solution and its derivative wrt D8. These extensions are, in general, based on decomposing the (unknown) multistep solution into the sum of a computable approximation D4 and an (unknown) error term CT, i.e.,

There exist standard techniques to build D4BN BSD4 BSD8 and to bound CTBN BSCT BSD8 . Section 6 reviews how they can be derived from generalized Hermite interpolation polynomials. Here we simply assume that they are available and we show how to use them to build filters.

In the previous section, we mention how natural multistep filters can be obtained by simply replacing the multistep solution and its derivative wrt D8 by their interval extensions to obtain

It is not easy however to enforce forward consistency on a natural filter since the variables may occur in complex nonlinear expressions. In addition, in ODEs, the boxes are often small which makes mean-value forms particularly interesting.

5.1 Mean-Value Filters

The mean-value theorem for a function CV states that

A mean-value interval extension BZ of CV in interval
can then be defined as

for some D1 in
, where BWBZ is an interval extension of the derivative of CV. In order to obtain a mean-value filter, consider the ODE

where the multistep solution has been expanded to make the approximations and the error terms explicit. We are thus interested in finding a mean-value interval extension wrt D9 CZ of the function defined as

However, as mentioned earlier, CT is not a computable function and we cannot compute its derivative wrt D9 CZ . Fortunately, it is possible to construct a mean-value filter by considering the error terms as constants and taking a mean-value interval extension of the function CVB4D8 CZ BN D9 CZ BN CTBN CSCTB5B4D8B5 defined as

If BZ is a mean-value interval extension of CV, then the ODE is approximated by an equation of the form

where BX and BWBX are interval extensions of the error term CT and of its derivative wrt D8. In the following, we call the resulting filter an implicit mean value filter.

The implicit mean-value filter for a time D8 produces an interval constraint of the form

where each CG CX and BU are vectors of D2 elements and each BT CX is an D2A2D2 matrix. Assuming that we are interested in pruning the last box BW CZ , we obtain an operator as

The implicit mean-value filter can be turned into explicit form by inverting the matrix BT CZ .

Definition 7 (Explicit Mean-Value Filter) Let  be an ODE D9 BC BP CUB4D8BN D9B5 and rewrite an implicit mean-value filter for  as

where BT CX is an D2A2D2 matrix and BU a vector of D2 elements. An explicit mean-value filter for  is given by

It is simple to make an explicit mean-value filter forward(1)- consistent, since the variables in CG CZ have been isolated.

5.2 Global Filters

The explicit mean-value filter produces significant pruning of the boxes but it still suffers from a dependency problem typical of interval methods. Consider the application of the filter at two successive time steps of the method. We obtain expressions of the form

Note that the second expression considers all tuples of the form CWDA BD BN BM BM BM BN DA CZ CX BE CWCG BD BN BM BM BM BNCG CZ CX to compute CG CZB7BD . But only a subset of these tuples satisfy the first equation for a given DA BC BE CG BC . The key idea to remove this dependency problem is to use a global filter which clusters a sequence of CZ mean-value filters together. Figure 3 gives the intuition for CZ BP BF. A global filter is then obtained by predicting and pruning CZ boxes simultaneously. More precisely, the idea is to cluster the CZ successive filters

This filter can be turned into an explicit form by substituting CG CZ in the second equation, CG CZ and CG CZA0BD in the third equation and so on. The resulting system can be written as a system

where  is a square D2CZ A2 D2CZ matrix and BU is a vector of D2CZ elements. It is important to mention that such a global filter is forward(CZ)-consistent wrt D8 BECZA0BD BNBW BECZA0BD and D8, since the variables CG CZ BM BM BMCG BECZA0BD have been isolated. The dependency problem has thus been reduced to the well-known wrapping effect which has standard solutions in interval methods.

5.3 TheWrapping Effect

Approximating the set of solutions to an equation by a box may induce a significant loss of precision. This problem is known as the wrapping effect and its standard solution is to choose an appropriate coordinate system to represent the solutions compactly. Let us review briefly how to integrate this idea with our proposals. Consider a global filter that we write as CI BD BP BT BC CI BC B7 BU BC . The idea is to introduce a coordinate transformation CH BD BP  BD CI BD to obtain a system  BD CI BD BP  BD BT BC CI BC B7 BD BU BC where  BD BT BC is diagonally dominant, i.e., the non-diagonal elements are very small compared to the diagonal. If  BD BT BC is diagonally dominant, then the wrapping effect induces no or little loss of precision in computing CH CX . Our pruning step can easily accommodate this idea. Iteration CX defines CH CX BP  CX CI CX and it sends CI CX , CH CX and  CX to iteration CX B7 BD of the main algorithm. Iteration CX B7 BD will have to solve the system

A coordinate transformation is applied once again to obtain

so that  CXB7BD BT CXB7BD  A0BD CX is diagonally dominant. Thematrices  CX are computed using QR-factorizations [Loh87].

6 Hermite Filters

We now show how to build interval extensions of D4, BSD4BPBSD8, CT, BSCTBPBSD8 by using Hermite interpolation polynomials [SB80]. Informally, a Hermite interpolation polynomial approximates a continuously differentiable function CU and is specified by imposing that its values and the values of its successive derivatives at various points be equal to the values of CU and of its derivatives at the same points. Note that the number of conditions at the different points may vary.

Definition 8 (Hermite(AR) Interpolation Polynomial) Let

D8 CX BE CABN D9 B4BCB5 CX BP D9 CX BE CA D2 and D9 B4CYB5 CX BP CU B4CYA0BDB5 B4D8 CX BN D9 CX B5BN CX BP BCBN BM BM BM BN CZBN CY BP BCBN BM BM BM BN AR CX A0 BD. Consider the ODE D9 BC BP CUB4D8BN D9B5. Let AR BE  CZB7BD BN AR CX BIBP BCBN CX BP BCBN BM BM BM BN CZBN AR D7 BP  CZ CXBPBC AR CX . The Hermite(AR) interpolation polynomial wrt CU and B4D8 CZ BN D9 CZ B5, is the unique polynomial D5 of degree AK AR D7 A0 BD satisfying

It is easy to take interval extensions of a Hermite interpolation polynomial and of its derivatives. The only remaining issue is to bound the error terms. The following standard theorem (e.g., [SB80], [Atk88]) provides the necessary theoretical basis.

Theorem 1 (Hermite Error Term) Let D4B4D8 CZ BN D9 CZ B5 be the Hermite(AR) interpolation polynomial wrt CU and B4D8 CZ BN D9 CZ B5. Let D9B4D8B5 AH D1D7B4D8 CZ BN D9 CZ B5B4D8B5, CC BP BECUD8 BC BN BM BM BM BN D8 CZ BN D8CV, AR D7 BP  CZ CXBPBC AR CX and DBB4D8B5 BP CZ CXBPBC B4D8 A0 D8 CX B5 AR CX . Then, for CX BP BDBN BM BM BM BN D2,

• AF BL AO CX BE CC BM CT CX B4D8 CZ BN D9 CZ B5B4D8B5 BP BD AR D7 AX CU B4AR D7 A0BDB5 CX B4AO CX BN D9B4AO CX B5B5DBB4D8B5;

How to use this theorem to bound the error terms? It suffices to take interval extensions of the formula given in the theorem and to replace AO CX BN AO BDBNCX BN AO BEBNCX by CC and D9B4AO CX B5BN D9B4AO BDBNCX B5BN D9B4AO BEBNCX B5 by a bounding box for the ODE over CC . In the following, we call Hermite(AR) filters, filters based on Hermite(AR) interpolation and we denote a global Hermite(AR) filter by GHF(AR). The following novel result is fundamental and specifies the speed of convergence of Hermite filters. It shows that the order of natural and mean-value Hermite(AR) filters is the sum of the elements in AR, i.e., the number of conditions imposed on the interpolation.

Proposition 2 (Order) Let BY be a natural or mean-value Hermite(AR) filter for ODE D9 BC BP CUB4D8BN D9B5. Let D8 CT A0 D8 CZ BP A2B4CWB5, D8 CXB7BD A0D8 CX BP A2B4CWB5BN CX BP BCBN BM BM BM BN CZA0BD, AR D7 BP  CZ CXBPBC AR CX andBW CZ BP BECUD9 BE CA D2 CY BYB4D8 CZ BN CWD9 BC BN BM BM BM BN D9 CZA0BD BN D9CXB5B4D8 CT B5CV. Then, under some weak assumptions on CU and BY, DBB4BW CZ B5 BP B4CW AR D7 B7BD B5BM

7 Theoretical Cost Analysis

We now analyse the cost of our method and compare it to Nedialkov's IHO(D4BN D5) method [NJ99], the best interval method we know of. We use the following assumptions. At each step, the forward process uses Moore's Taylor method and the pruning component applies a global Hermite filter together with coordinate transformations (using Lohner's QRfactorization technique). For simplicity of the analysis, we assume that (the natural encoding of) function CU contains only arithmetic operations. We denote by  BD the number of A3BN BP operations in CU , by  BE the number of A6 operations, and by  the sum  BD B7 BE . We also assume that the cost of evaluating BSCU B4D6B5 BPBSD9 is D2 times the cost of evaluating CU B4D6B5 . We define AR D1 as maxB4ARB5, AR D7 BP AR BC B7BM BM BMB7AR CZ BN D4B7D5B7BD BP AR D7 BN D5 BE CUD4BN D4B7BDCV. We also report separately interval arithmetic operations involved in (1) products of a real and an interval matrix (Cost1) and (2) the generation of Jacobians (Cost-2). Table 7 reports the main cost of a step in the IHO(D4BN D5) method (IHO in the table) and our method GHF(AR) (GHF in the table). It also shows the complexity of two particular cases of GHF(AR). The first case (GHF-1) corresponds to a polynomial with only two interpolation points (CZ BP BDBN CYAR BD A0AR BC CY AK BD), while the second case corresponds to a polynomial imposing two conditions on every interpolation points (AR BC BP BM BM BM BP AR CZ BP BE). Note that the methods are of the same order in all cases.

The first main result is that GHF-1 is always cheaper than IHO, which means that our method can always be made to run faster by choosing only two interpolation points. (The next section will show that substantial improvement in accuracy is also obtained in this case). GHF-2 is more expensive than GHF-1 and IHO when CU is simple. However, when CU contains many operations (which is often the case in practical applications), GHF-2 can become substantially faster because Cost-1 in GHF-2 is independent of CU and Cost-2 is substantially smaller in GHF-2 than in GHF-1 and IHO. It also shows the versatility of the approach that can be taylored to the application at hand.

8 Experimental Analysis

We now report experimental results of our C++ implementation on some standard benchmarks [HNW87; Loh87] and two molecular biology problems, which are real-life parametric ODEs given to us by a biologist:

Table 2 compares GHF(AR) and IHO(D4BN D5) methods of the same order. It reports the precision at the last step and execution time (in seconds) of both methods for the same (constant) stepsize. The experimental results follow the same assumptions as in the theoretical analysis. The forward process uses Moore's Taylor method of order D5 B7 BD (same order as the predictor used in IHO(D4BN D5)) and a Taylor series method of order AR D7 to compute a bounding box, except for 2BP and 3BP where we use a series of order 1. The choice of the evaluation time D8 CT involved in GHF(AR) has not been discussed yet. So far we have no theoretical result about the optimal choice of D8 CT . We use a simple binary search algorithm to determine a good value for D8 CT at the beginning of or during the integration. In our experiments, we chose D8 CT between the last two interpolation points, keeping the distance constant throughout the integration. Our results could be further improved by using a variable distance.

The results indicate that our method produces orders of magnitude improvements in accuracy and runs faster than the best known method. The gain in precision is particularly significant for lower orders. The theoretical results are also confirmed by the experiments. When CU contains many operations (e.g. in 3BP), using many interpolation points is particularly effective. For very complex functions, the gain in computation time could become substantial. When CU is simple, using few interpolation points becomes more interesting.

As a consequence, we believe that a constraint satisfaction approach to parametric ordinary differential equations is a very promising avenue that complements well traditional approaches. The main open issue at this point is the choice of an optimal evaluation point. This will be the topic of our research in the coming months.

Acknowledgment

This research is partially supported by the actions de recherche concert´ee ARC/95/00-187 and an NSF NYI award. Special thanks to Philippe Delsarte for interesting discussions and for his detailed comments on this paper.