**2. Algorithm for Adaptive Interpolation Using Sparse Grids**

Dynamic systems with uncertainties in parameters arise in many practical areas. Traditionally, interval problems for dynamic systems are formulated in the form of the Cauchy problem for a system of ordinary differential equations (ODE) with interval initial conditions or parameters. It is necessary to obtain an interval estimate of the solution based on interval values of the parameters.

Consider the Cauchy problem with *m* interval initial conditions:

$$\begin{cases} \frac{dy\_i(t)}{dt} = f\_i(y\_1(t), y\_2(t), \dots, y\_n(t)), \; 1 \le i \le n\_\prime\\ y\_i(t\_0) \in \left[ \underline{y}\_i^0, \overline{\underline{y}\_i^0} \right], \; 1 \le i \le m\_\prime\\ y\_i(t\_0) = \underline{y}\_i^0, m+1 \le i \le n\_\prime\\ t \in [t\_0, t\_N]. \end{cases} \tag{1}$$

Hereinafter, the underline denotes the lower bound of the interval, and the overline the upper bound of the interval.

If the ODE system is not autonomous or contains interval parameters, then fictitious equations are added to the system so that it would take the form of system Equation (1). A vector function **f** = (*f*1, *f*2, ..., *fn*) *<sup>T</sup>* meets all conditions ensuring the uniqueness and existence of a solution for all *yi*(*t*0) ∈ 4 *y*0 *<sup>i</sup>* , *<sup>y</sup>*<sup>0</sup> *i* 5 , 1 ≤ *i* ≤ *m*.

The goal is, for each moment of time *tk*, to construct a piecewise polynomial vector function **P**k *y*0 <sup>1</sup>, *<sup>y</sup>*<sup>0</sup> <sup>2</sup>, ..., *<sup>y</sup>*<sup>0</sup> *m* , where *yi*(*t*0) ∈ 4 *y*0 *<sup>i</sup>* , *<sup>y</sup>*<sup>0</sup> *i* 5 , 1 ≤ *i* ≤ *m*, which interpolates the dependence of the solution on the interval parameters with controlled accuracy. If the function **P**<sup>k</sup> is available, finding the interval estimate of the solution (finding the left and right boundaries of the intervals) should be reduced to solving constrained optimization problems for an explicitly given function.

Suppose that the solution to **y**k *y*0 <sup>1</sup>, *<sup>y</sup>*<sup>0</sup> <sup>2</sup>, ..., *<sup>y</sup>*<sup>0</sup> *m* is known at the moment of time *tk*, where *y*<sup>0</sup> *i* ∈ 4 *y*0 *<sup>i</sup>* , *<sup>y</sup>*<sup>0</sup> *i* 5 , 1 ≤ *i* ≤ *m*. An adaptive sparse grid is constructed for the set formed by the interval initial conditions. Each grid point has a corresponding solution to the noninterval system (1) at pointwise values of interval parameters that correspond to the position of a node. To obtain an interval solution at the moment of time *tk*+1, the transfer of all non-interval solutions contained in the grid nodes to the time layer (*k* + 1) is performed, followed by the adaptation of the grid and the construction of an interpolation polynomial **P**k<sup>+</sup>1.

A short description of sparse grid interpolation according to the works [20,21] is given below.

Consider the interpolation of a smooth function *f*(*x*) of one variable on the unit interval [0, 1]. For the sake of simplicity, it is assumed that the function is equal to zero at the boundary points: *f*(0) = *f*(1) = 0.

The interpolation is performed on a piecewise linear hierarchical basis using the hat function:

$$\mathfrak{g}(\mathfrak{x}) = \begin{cases} 1 - |\mathfrak{x}| \, \mathfrak{x} \in [-1, 1] \\ 0 \text{, otherwise} \end{cases} \text{.} \tag{2}$$

Define a set of grids *Gl* on a unit interval [0, 1], where *l* is the level that determines the grid width *hl* = 2−*<sup>l</sup>* . The grid points *xl*, *<sup>i</sup>* are given as:

$$\mathbf{x}\_{l,i} = i \cdot h\_{l'} \ 0 \le i \le 2^I.$$

Families of basis functions *ϕl*, *<sup>i</sup>*(*x*) are generated based on the obtained sets of points, using the stretching and transfer of the hat Equation (2):

$$\varphi\_{l,i}(\mathbf{x}) = \varphi\left(\frac{\mathbf{x} - \mathbf{i} \cdot \mathbf{h}\_l}{\mathbf{h}\_l}\right). \tag{3}$$

.

A nodal basis is formed for each given *l* of Equation (3). Here, the common piecewise linear interpolation (Figure 1) is applied, and the corresponding polynomial is written as follows:

$$P(\mathbf{x}) = \sum\_{i=1}^{2^l - 1} a\_{l,i} q\_{l,i}(\mathbf{x}), \ a\_{l,i} = f(\mathbf{x}\_{l,i}). \tag{4}$$

**Figure 1.** Interpolation on a nodal basis.

Let us make the transition to a hierarchical basis (Figure 2). The basis functions given by Equation (3) is expressed with even level indices *k* in terms of the basis level functions (*k* − 1):

$$\varphi\_{k,2i}(\mathbf{x}) = \varphi\_{k-1,i}(\mathbf{x}) - \frac{1}{2}(\varphi\_{k,2i-1}(\mathbf{x}) + \varphi\_{k,2i+1}(\mathbf{x})), \ 1 \le i \le 2^{k-1} - 1.$$

In this case, the interpolation polynomial given by Equation (4) takes the following form:

$$P(\mathbf{x}) = \sum\_{k=1}^{l} \sum\_{\substack{i=1 \\ i = 1, \\ i \text{ odd}}}^{2^k - 1} a\_{k,i} \varphi\_{k,i}(\mathbf{x}), \; a\_{k,i} = f(\mathbf{x}\_{k,i}) - \frac{1}{2} (f(\mathbf{x}\_{k,i-1}) + f(\mathbf{x}\_{k,i+1})) \tag{5}$$

Next, consider the multidimensional interpolation of a smooth function *f*(*x*1, *x*2, ..., *xd*) using *d*—dimensional unit cube Ω = [0, 1] *d* , provided that *f* |*∂*<sup>Ω</sup> = 0. A multidimensional basis is constructed by the direct product of hierarchical one-dimensional bases:

$$\varphi\_{\mathbf{l},\mathbf{i}}(\mathbf{x}) = \prod\_{j=1}^{d} \varphi\_{l\_j,\mathbf{i}\_j}(\boldsymbol{\pi}\_j) , \ 1 \le \boldsymbol{i}\_j \le 2^{l\_j} - 1 \; \ i\_j \text{ odd} ,$$

where **l** = (*l*1, *l*2, ..., *ld*) are the levels of the corresponding one-dimensional grids, **i** = (*i*1, *i*2, ..., *id*) is the basis function multi-index, **x** = (*x*1, *x*2, ..., *xd*). If *d* ∑ *j*=1 *lj* ≤ *n* + *d* − 1, there is a sparse grid of *n* level; if max *j*=1,*d lj* ≤ *n*, there exists a complete grid (Figure 3). The number of nodes in a sparse grid is estimated as *O p*(log2 *p*) *d*−1 , and the interpolation error is estimated as *O h*2 *<sup>n</sup>*(log2 *p*) *d*−1 ; for a full grid the respective number of nodes is *O pd* , and the error is *O h*2 *n* , where *<sup>p</sup>* <sup>=</sup> <sup>2</sup>*<sup>n</sup>* <sup>−</sup> 1 is the number of nodes in each dimension [20].

**Figure 2.** Interpolation on a hierarchical basis.

**Figure 3.** Sparse grid of the third level: black dots—sparse grid; black and grey dots—full grid: (**a**) sets of points corresponding to basis functions of the same level; (**b**) combining all points into one grid.

The interpolation polynomial is written as follows:

$$P(\mathbf{x}) = \sum\_{\mathbf{l}, \mathbf{i}} a\_{\mathbf{l}, \mathbf{i}} \varphi\_{\mathbf{l}, \mathbf{i}}(\mathbf{x}) , \left| \mathbf{l} \right|\_{1} \le n + d - 1 , \ 1 \le i\_{\mathbf{j}} \le 2^{l\_{\mathbf{j}}} - 1 , \ i\_{\mathbf{j}} \text{odd} \tag{6}$$

where

$$a\_{\mathbf{l},\mathbf{i}} = \sum\_{\Delta\_{\mathbf{l}},...,\Delta\_{\mathbf{l}}} \left(-\frac{1}{2}\right)^{\sum\_{j=1}^{d}|\Delta\_{\mathbf{j}}|} f\left(\mathbf{x}\_{l\_1j\_1+\Delta\_{\mathbf{l}}\prime},\mathbf{x}\_{l\_2j\_2+\Delta\_{\mathbf{l}}\prime},...,\mathbf{x}\_{l\_dj\_d+\Delta\_{\mathbf{l}}}\right)\_{\prime}, -1 \le \Delta\_{\mathbf{j}} \le 1, \ 1 \le j \le d. \tag{7}$$

In the case when the interpolated function has a nonzero value at the boundary, the one-dimensional basis is supplemented by two additional functions: *ϕ*0, 0(*x*) and *ϕ*0, 1(*x*) (Figure 4). Two values are added to the polynomial given by Equation (5): *a*0,0*ϕ*0,0(*x*) and *a*0,1*ϕ*0,1(*x*), where *a*0,0 = *f*(0), *a*0,1 = *f*(1). By analogy, for multidimensional interpolation, it follows that if *lj* = 0, then *ij* = 0, 1 in Equation (6) and Δ*<sup>j</sup>* = 0 in Equation (7). Allowance for boundary values in the multidimensional case can be considered as the construction of sparse grids for all faces of lower dimensions. Figure 5 shows a sparse grid, which takes into account the boundary values.

**Figure 4.** Additional basis functions taking into account boundary values.

**Figure 5.** Sparse grid of the level *n* = 4, which takes into account boundary values.

In addition, there are adaptive sparse grids for which a general tree can be used to perform structuring. Each vertex of the tree corresponds to a certain basis function *ϕ***l**,**i**. If the value of the corresponding coefficient *a***l**,**i**/max(*f*(*x***l**,**i**), 1) > *ε*, where *ε* is some predetermined value, then each vertex creates 2*d* descendants, which correspond to the basis functions of the next level. This process continues recursively until the values *a***l**,**<sup>i</sup>** at all leaf vertices become less than *ε*. With this approach, it is important to make sure that there is no duplication of vertices.

Consider some examples. Figure <sup>6</sup> shows several functions <sup>R</sup><sup>2</sup> <sup>→</sup> <sup>R</sup> and the resulting adaptive grid, Figure <sup>7</sup> shows grids for functions <sup>R</sup><sup>3</sup> <sup>→</sup> <sup>R</sup>.

**Figure 6.** Examples showing interpolation using functions of two variables. (**a**) Linear combination of univariate functions. (**b**) Linear combination of univariate functions and a function of two variables *x* and *y* with a small coefficient. (**c**) Linear combination of univariate and a function of two variables *x* and *y*.

**Figure 7.** Examples of grids for functions of three variables. (**a**) Linear combination of univariate functions. (**b**) Linear combination of univariate functions and a function of two variables *y* and *z*. (**c**) Nonlinear function of three variables.

It can be seen from the figures that if the initial dependence is a linear combination of functions determined by certain subgroups of variables, then the adaptive sparse grid will become more dense not in the entire set (Figures 6c and 7c), but only in subsets of lower dimension that correspond to these subgroups (Figures 6a and 7a,b). The subsets for grid construction are determined by those subgroups of parameters, for which the mixed derivatives are nonzero, and the grid density directly depends on the values of these derivatives (Figure 6b,c).

To build a solution for the system given by Equation (1), the uncertainty area *y*<sup>0</sup> *i* ∈ 4 *y*0 *<sup>i</sup>* , *<sup>y</sup>*<sup>0</sup> *i* 5 , 1 ≤ *i* ≤ *m* is transformed with the help of displacement and stretching into a *m*dimensional unit cube. Taking into account that solving the problem requires interpolating *n* functions at once (*n* is the number of phase variables of the system), Equations (6) and (7) will take the following form:

$$\mathbf{P}^{\mathbf{k}}\left(\mathbf{y}^{0}\right) = \sum\_{\mathbf{l},\mathbf{i}} \mathbf{a}\_{\mathbf{l},\mathbf{i}}^{\mathbf{k}} \boldsymbol{\varphi}\_{\mathbf{l},\mathbf{i}}\left(\mathbf{y}^{0}\right).$$

where

$$\mathbf{a}\_{\mathbf{l},\mathbf{i}}^{\mathbf{k}} = \sum\_{\Delta\_{l},\ldots,\Delta\_{\mathbf{m}}} \left(-\frac{1}{2}\right)^{\sum\_{j=1}^{m}|\Delta\_{j}|} \mathbf{y}^{\mathbf{k}} \Big(y\_{1,l\_{1}i\_{1}+\Delta\_{1}}^{\mathbf{0}} y\_{2,l\_{2}i\_{2}+\Delta\_{2}}^{\mathbf{0}}, \ldots, y\_{m,l\_{m}i\_{m}+\Delta\_{m}}^{\mathbf{0}}\Big), -1 \le \Delta\_{\mathbf{j}} \le 1, \, 1 \le j \le m \tag{8}$$

The vector norm **a**<sup>k</sup> **l,i** (for example, the maximum one) can be used as a criterion for adapting the grid.

Construct an interpolation polynomial **P**k+1 **y**0 . All the solutions that participated in the calculation of the coefficients given by Equation (8) are transferred to the (*k* + 1)-th time layer using some numerical integration method, after which a new set of **a**k+<sup>1</sup> **l,i** coefficients is calculated and the adaptation of the grid is performed. When compacting the grid, the addition of new basis functions occurs at the *k*-th time layer and the solutions involved in computing the corresponding weight coefficients are transferred to the next layer.

The efficiency of the considered approach will be noticeable when many mixed derivatives of the solution with respect to the parameters *<sup>∂</sup>*Σ*α<sup>i</sup>* **<sup>y</sup>**<sup>k</sup>(*y*<sup>0</sup> <sup>1</sup>, *<sup>y</sup>*<sup>0</sup> <sup>2</sup>,..., *<sup>y</sup>*<sup>0</sup> *m*) (*∂y*<sup>0</sup> 1) *<sup>α</sup>*<sup>1</sup> (*∂y*<sup>0</sup> 2) *<sup>α</sup>*<sup>2</sup> ...(*∂y*<sup>0</sup> *m*) *<sup>α</sup><sup>m</sup>* , max 1≤*i*≤*m α<sup>i</sup>* ≤ 2, *y*0 *i* ∈ 4 *y*0 *<sup>i</sup>* , *<sup>y</sup>*<sup>0</sup> *i* 5 , 1 ≤ *i* ≤ *m* are negligible or equal to zero. Particularly, this takes place, if the solution to the ODE system can be represented as a linear combination of functions determined by a certain subset of interval parameters.

Thus, the scope of application of the proposed approach is rather wide and includes various dynamic systems. In the next section, it is demonstrated how the method is applied to some representative problems.

#### **3. Approbation of the Algorithm for Nonlinear Dynamics Problems**

To characterize computational costs, a criterion is determined, which is equal to the average number of integrated non-interval ODE systems at a time step in the computational process:

$$I = \frac{1}{N} \sum\_{k=1}^{N} C\_{k\prime}$$

where *Ck* is the number of nodes at the *k* step. A similar criterion exists for the classical adaptive interpolation algorithm [11]. The *I* value is equivalent to the number of sampling points from the original region of uncertainty.

To estimate the posterior interpolation error at the initial moment of time, *ncheck* points are randomly generated:

$$y\_i^j(t\_0) = rand \left[ \underline{y\_{i'}^0} \, \overline{y\_i^0} \right], 1 \le i \le m, \, 1 \le j \le n\_{check}.$$

For the initial conditions obtained, with the help of a numerical integration method, solutions are constructed at the final moment of time *tN*. The relative posterior global estimate of the error is written as follows:

$$error = \max\_{\substack{1 \le j \le n\_{check\lambda'} \\ 1 \le i \le n}} \frac{\left| P\_i^N \left( y\_1^j(t\_0), y\_2^j(t\_0), \dots, y\_m^j(t\_0) \right) - y\_i^j(t\_N) \right|}{\max \left( \left| y\_i^j(t\_N) \right|, 1 \right)}.$$

Let us integrate several ODE interval systems using the described approach. The calculation is performed for two values of *ε* = 10−<sup>3</sup> and *ε* = 10−<sup>5</sup> (*ε* imposes a restriction on the values of the weight coefficients of the basic functions when constructing an adaptive sparse grid). First, let us take into account an ordinary differential system with two interval initial conditions, which describes a conservative oscillator:

$$\begin{cases} \; \mathbf{x}\prime = \mathbf{y}\prime, \mathbf{y}\prime = -\sin(\mathbf{x})\prime\\ \; \mathbf{x}(0) = \mathbf{x}\_0 \in [-1, 1], \; \mathbf{y}(0) = \mathbf{y}\_0 \in [0, 1], \; t \in [0, 25]. \end{cases} \tag{9}$$

Figure 8 shows a set of solutions for the system given by Equation (9) at different moments of time; it twists into a spiral structure during the integration. Figure 9 shows the grid resulting from applying the algorithm. The points in these two figures correspond to each other.

**Figure 8.** The interval solution of system given by Equation (9) at different moments of time.

**Figure 9.** The grids obtained in the process of solving system given by Equation (9).

Table 1 shows a comparison of computational costs and error estimates for different approaches. When set to a low precision (*ε* = 10−3), adaptive sparse grids work a little faster than the classical adaptive interpolation algorithm, and twice as fast as conventional sparse grids. However, for *ε* = 10−5, the classical algorithm wins due to the application of an interpolation polynomial of a high degree. The levels of the grids were adjusted to obtain approximately the same error as in other approaches.


**Table 1.** Comparison of approaches for system given by Equation (9).

Next, the Volterra-Lotka model with interval initial conditions and one interval coefficient is considered. The Cauchy problem in the case has the form:

$$\begin{cases} \text{ x'} = 4x - \frac{5}{4}xy - ax^2, \, yt = -2y + \frac{1}{2}xy - \frac{1}{20}y^2, \\\ x(0) = x\_0 \in [4, 5], \, y(0) = y\_0 \in [2.8, 3.2], \, t \in [0, 25], \end{cases} \tag{10}$$

where *α* ∈ [−0.05, 0.05].

This model describes predator–prey interactions. A feature of the system is the fact that at *α* < 0 there is an unstable focus and the amplitude of fluctuations in the population of species grows, and at *α* > 0 the focus is stable and the state of the system tends to be stationary over time.

Figure 10 shows the set of solutions for the system at different points in time. The following picture is clearly observed here: some part of the set converges to a point, which corresponds to a stable focus, and another part of the set increases in its size, which corresponds to an unstable focus. Figure 11 shows the resulting grid. Due to the fact that uncertainty is present in the parameters, the set of solutions on the phase plane is only a projection of the three-dimensional set onto the two-dimensional phase space. The additional dimension corresponds to the interval parameter *α*.

Table 2 shows a comparison of the different approaches. Similar to the previous task, adaptive sparse grids are effective with lower accuracy *ε*.


**Table 2.** Comparison of approaches on system given by Equation (10).

**Figure 10.** The interval solution of system given by Equation (10) at different times.

**Figure 11.** The grid obtained in the process of solving the system given by Equation (10).

Consider an ordinary differential system presenting the expanded Volterra-Lotka model with three interval initial conditions and seven interval parameters:

$$\begin{cases} \mathbf{x}' = \mathbf{x}(\delta\_1 - y - \varepsilon \mathbf{x}), & \mathbf{x}(0), y(0), z(0), \delta\_1, \delta\_2, \gamma\_{1\prime}, \gamma\_2 \in [1.0, 1.01], \\ \mathbf{y}t = -\gamma\_1 y(\delta\_2 - \mathbf{x} + z) - \eta y^2, & \varepsilon, \eta \in [-0.0005, 0.0005], \\ zt = -\gamma\_2 z(\mathbf{a} - \mathbf{y}), & \mathbf{a} \in [0.9, 0.91]. \end{cases} \tag{11}$$

Figure 12 shows the dependences of the interval estimates of solutions on time.

**Figure 12.** Time dependence of the upper and lower bounds for the solution of system given by Equation (11): (**a**) *x*(*t*); (**b**) *y*(*t*); (**c**) *z*(*t*).

For a reasonable time, the solution was obtained using only adaptive sparse grids. For *<sup>ε</sup>* <sup>=</sup> <sup>10</sup>−3, the obtained result was *<sup>I</sup>* <sup>=</sup> 81, 566.1 and *error* <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup>−2.

Consider a model describing the motion of interacting bodies. The problem can be formulated as a dynamic system with interval initial velocities. The system of ordinary differential equations in dimensionless variables is as follows:

$$\begin{cases} \left(v\_{i}^{\rm x}\right)' = \sum\_{\substack{j=1, j\neq i}}^{4} m\_{j} (\mathbf{x}\_{j} - \mathbf{x}\_{i}) \tau\_{i,j}^{-3}, \left(v\_{i}^{\rm y}\right)' = \sum\_{\substack{j=1, j\neq i}}^{4} m\_{j} (y\_{j} - y\_{i}) \tau\_{i,j}^{-3}, \left(v\_{i}^{\rm z}\right)' = \sum\_{j=1, j\neq i}^{4} m\_{j} (z\_{j} - z\_{i}) \tau\_{i,j}^{-3},\\ \mathbf{x}\_{i}' = v\_{i}^{\rm x}, y\_{.f}' = v\_{i}^{\rm y}, z\_{i}' = v\_{i}^{\rm z}, 2 \le i \le 4,\\ \mathbf{x}\_{1}(0) = y\_{1}(0) = z\_{1}(0) = v\_{1}^{\rm x}(0) = v\_{1}^{\rm y}(0) = v\_{1}^{\rm y}(0) = 0,\\ \mathbf{x}\_{2,3}(0) = \pm 1, y\_{2,3}(0) = z\_{2,3}(0) = 0, v\_{2,3}(0) = \left(0 \quad \pm v \quad 0\right)^{\rm T} + \Delta v\_{2,3}^{\rm T}\\ y\_{4}(0) = 1, \mathbf{x}\_{4}(0) = z\_{4}(0) = 0, v\_{4}(0) = \left(0 \quad 0 \quad v\right)^{\rm T} + \Delta v\_{4}^{\rm x}\\ t \in [00, 0.02] \end{cases} \tag{12}$$

where *ri*,*<sup>j</sup>* = ! *xj* − *xi* <sup>2</sup> + *yj* − *yi* <sup>2</sup> + *zj* − *zi* <sup>2</sup> is the distance between two bodies, *v* = 316.23 is the initial velocity of bodies, *m*<sup>1</sup> = 105, *m*2, 3, 4 = 10−<sup>5</sup> are body masses, Δ*v*2, 3, 4 = ([−2, 2], [−2, 2], [−2, 2]) are the interval uncertainties in body velocities.

Figure 13 shows graphs for the dependence of the interval estimates of the 2nd body coordinates and velocities on time. Similar to the previous problem, the solution was calculated only using adaptive sparse grids. For *ε* = 10−3, the obtained result was *<sup>I</sup>* <sup>=</sup> 133830.9 and *error* <sup>=</sup> 2.6 <sup>×</sup> <sup>10</sup>−2.

This system is demonstrative because the uncertainty in the speed of a particular body mainly affects the position and speed of that particular body and has little effect on other bodies. In this regard, the solution of the system will have a specific form, as most of the mixed derivatives will be close to zero.

Note that the classical adaptive interpolation algorithm for systems given by Equations (6) and (7) constructs sets of solutions with fewer integrations of the corresponding non-interval ordinary differential systems since it uses nonlinear interpolation. However, when the number of interval parameters increases (systems given by Equations (8) and (9)), the use of adaptive sparse grids becomes more efficient. When increasing the dimension of the problem, it is practically impossible to increase the degree of the interpolation polynomial in the adaptive interpolation algorithm to obtain higher accuracy due to the exponential growth of the number of nodes in the grid. Therefore, for high dimensional-problems, it is suitable to use methods that have lower accuracy, but at the same time allow reasonable computational costs; in particular, adaptive sparse grids.

The examples above demonstrate that by using sparse grids it is possible to simulate dynamic systems with ten interval uncertainties in a reasonable time. When solving system given by Equation (11), the equivalent number of sampling points was about 80 thousand, and in the case of using classical interpolation with the degree of polynomial equal to 4, the value would be of order 107. A lower estimate of the computational cost can be obtained. It

follows from Equation (8) that the number of solved non-interval ODE systems cannot be less than 3*m*. The upper estimate of the computational costs, in the general case, essentially depends on the features of the ODE system being solved, in particular on the values of the mixed derivatives of the solution with respect to the point values of the interval parameters. For comparison, the classical adaptive interpolation algorithm requires at least (*p* + 1) *m* points, and the method proposed in Reference [15] requires (*m*+*p*)! *<sup>m</sup>*!*p*! points, where *p* is the degree of the interpolation polynomial.

**Figure 13.** Time dependencies of upper and lower estimates of the 2nd body coordinates and velocities: (**a**) *x*2(*t*); (**b**) *y*2(*t*); (**c**) *z*2(*t*); (**d**) *v<sup>x</sup>* <sup>2</sup> (*t*); (**e**) *v y* <sup>2</sup>(*t*); (**f**) *<sup>v</sup><sup>z</sup>* <sup>2</sup>(*t*).
