Sparse Grid Adaptive Interpolation in Problems of Modeling Dynamic Systems with Interval Parameters

: The paper is concerned with the issues of modeling dynamic systems with interval parameters. In previous works, the authors proposed an adaptive interpolation algorithm for solving interval problems; the essence of the algorithm is the dynamic construction of a piecewise polynomial function that interpolates the solution of the problem with a given accuracy. The main problem of applying the algorithm is related to the curse of dimension, i.e., exponential complexity relative to the number of interval uncertainties in parameters. The main objective of this work is to apply the previously proposed adaptive interpolation algorithm to dynamic systems with a large number of interval parameters. In order to reduce the computational complexity of the algorithm, the authors propose using adaptive sparse grids. This article introduces a novelty approach of applying sparse grids to problems with interval uncertainties. The efﬁciency of the proposed approach has been demonstrated on representative interval problems of nonlinear dynamics and computational materials science.


Introduction
Problems related to inaccurately specified data arise in many modern fields of science and technology. When applied to non-stationary processes, they are often formulated as dynamic systems with interval parameters. The result of solving such problems is an interval estimate of the set of possible system states depending on the uncertainties in the parameters. Basic methods of interval analysis are presented in books [1][2][3][4][5]. There are known methods based on the representation of a set of solutions through geometric primitives: parallelepipeds and ellipses [6,7], methods based on symbolic computation [8,9], as well as stochastic methods [10], such as Monte-Carlo methods. Methods based on classical interval arithmetic are subject to the so-called wrap effect [1], which manifests itself in an unlimited increase in the width of the obtained interval estimates of solutions. This effect arises due to the replacement of the exact form of the set of solutions by a simpler form, and for iterative methods, the divergence of intervals' boundaries is often exponential. Existing methods that are not subject to this effect, or weakly susceptible to it, often have exponential complexity with respect to the number of interval parameters. It concerns symbolic methods operating in series, Monte-Carlo methods, and the adaptive interpolation algorithm [11]. Therefore, there is a need for efficient approaches to reduce the computational complexity of methods that are not affected by the wrapping effect.
While solving a considered class of problems, the main idea is to construct an explicit dependence of the solution to the corresponding non-interval problem on the point values of the interval parameters. If such dependence is available, finding an interval estimate would be reduced to solving a certain number of constrained optimization problems for explicitly given functions. 3

of 17
The research methodology is based on methods of mathematical modeling, computational mathematics, and differential calculus. The statement of the problem is formulated in the form of the Cauchy problem for a system of ordinary differential equations with interval parameters. The method is tested on a representative set of problems.
The following sections give a description of the adaptive interpolation algorithm on sparse grids, present the results of testing the algorithm on a number of model problems of nonlinear dynamics, and solve an important problem of computational materials science, namely the determination of an interval stress tensor based on molecular dynamics modeling.

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: , y 2 (t), ..., y n (t)), 1 Hereinafter, the underline denotes the lower bound of the interval, and the overlinethe 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 , ..., f n ) T meets all conditions ensuring the uniqueness and existence of a solution for all y i (t 0 ) ∈ y 0 i , y 0 i , 1 ≤ i ≤ m. The goal is, for each moment of time t k , to construct a piecewise polynomial vector function P k y 0 1 , y 0 2 , ..., y 0 m , where y i (t 0 ) ∈ y 0 i , y 0 i , 1 ≤ i ≤ m, which interpolates the dependence of the solution on the interval parameters with controlled accuracy. If the function P k 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 1 , y 0 2 , ..., y 0 m is known at the moment of time t k , where y 0 i ∈ y 0 i , y 0 i , 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 t k+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+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: The interpolation is performed on a piecewise linear hierarchical basis using the hat function: Define a set of grids G l on a unit interval [0, 1], where l is the level that determines the grid width h l = 2 −l . The grid points x l, i are given as: Families of basis functions ϕ l, i (x) are generated based on the obtained sets of points, using the stretching and transfer of the hat Equation (2): 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: Define a set of grids l G on a unit interval [0, 1] , where l is the level that determines the grid width 2 l l h   . The grid points , li x are given as: Families of basis functions , () li x  are generated based on the obtained sets of points, using the stretching and transfer of the hat Equation (2): (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: 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 ( 1) k  :  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): In this case, the interpolation polynomial given by Equation (4) takes the following form: Next, consider the multidimensional interpolation of a smooth function f (x 1 , x 2 , ..., x d ) using d-dimensional unit cube Ω = [0, 1] d , provided that f | ∂Ω = 0. A multidimensional basis is constructed by the direct product of hierarchical one-dimensional bases: In this case, the interpolation polynomial given by Equation (4) takes the following form: Next, consider the multidimensional interpolation of a smooth function   The interpolation polynomial is written as follows: The interpolation polynomial is written as follows: 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, . By analogy, for multidimensional interpolation, it follows that if l j = 0, then i j = 0, 1 in Equation (6) and ∆ j = 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. The interpolation polynomial is written as follows: In the case when the interpolated function has a nonzero value at the boundary, the one-dimensional basis is supplemented by two additional functions: Figure 4). Two values are added to the polynomial given by Equation (5): . By analogy, for multidimen- (6) and 0 j  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.  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 Consider some examples. Figure 6 shows several functions 2  and the resulting adaptive grid, Figure 7 shows grids for functions 3  . 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 2d descendants, which correspond to the basis functions of the next level. This process continues recursively until the values a l,i 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 6 shows several functions R 2 → R and the resulting adaptive grid, Figure 7 shows grids for functions R 3 → R . predetermined value, then each vertex creates 2d descendants, which correspond to the basis functions of the next level. This process continues recursively until the values , a li 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 6 shows several functions 2  and the resulting adaptive grid, Figure 7 shows grids for functions 3  .  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 ( Figure 6c and Figure 7c), but only in subsets of lower dimension that correspond to these subgroups ( Figure 6a and Figure 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 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:  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 0 i ∈ y 0 i , y 0 i , 1 ≤ i ≤ m is transformed with the help of displacement and stretching into a mdimensional 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: The vector norm a k 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+1 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 ∂ Σα i y k (y 0 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.

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: where C k 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, n check points are randomly generated: For the initial conditions obtained, with the help of a numerical integration method, solutions are constructed at the final moment of time t N . The relative posterior global estimate of the error is written as follows: Let us integrate several ODE interval systems using the described approach. The calculation is performed for two values of ε = 10 −3 and ε = 10 −5 (ε 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: 25].
(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.   Table 1 shows a comparison of computational costs and error estimates for different approaches. When set to a low precision ( ), adaptive sparse grids work a little faster than the classical adaptive interpolation algorithm, and twice as fast as conventional sparse grids. However, for , the classical algorithm wins due to the ap-  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. Next, the Volterra-Lotka model with interval initial conditions and one interval coefficient is considered. The Cauchy problem in the case has the form: 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 ε. 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  . Figure 10. The interval solution of system given by Equation (10) at different times. Figure 10. The interval solution of 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: x(0), y(0), z(0), δ 1 , δ 2 , γ 1 , γ 2 ∈ [1.0, 1 (11) Figure 12 shows the dependences of the interval estimates of solutions on time.
is the distance between two bodies,

v 
is the initial velocity of bodies,   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  For a reasonable time, the solution was obtained using only adaptive sparse grids. For ε = 10 −3 , the obtained result was I = 81, 566.1 and error = 1.2 × 10 −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: where  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 I = 133830.9 and error = 2.6 × 10 −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 10 7 . 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

Computation of an Interval Stress Tensor for Materials with a Covalent Chemical Bond
Let us consider an applied problem of computational materials science, within the framework of which the stresses arising during the deformation of an ideal crystal are calculated [12]. Different angles are possible between the orientation of the crystal lattice and the direction of stretching with a fixed stretch value. The stress tensor thus becomes interval. This problem is solved using molecular dynamics simulation. The motion of atoms is described by the classical equations of dynamics: where r i is the radius vector of the atom with the number i, v i is its velocity, m i is its mass, and F i is the force acting on the atom, in this case F i = −∇ i E, where E is the total energy of the system, and ∇ i is the gradient along the position of the atom with the number i.
In this problem, materials with a covalent interatomic bond are considered. The total energy of interaction between atoms of such materials is well described using the Tersoff potential [25]: where E is total system energy, V ij is the contribution to the interaction energy of atoms with numbers i and j, r ij is the distance between atoms with numbers i and j, f C (r) is a cut-off function, f R (r) and f A (r) are the repulsion and attraction potentials, respectively, and R, D, A, B, n, m, λ 1 , λ 2 , λ 3 , β, γ, c, d and cos(θ 0 ) are potential parameters that are selected in order to reproduce the properties of the simulated material. Methods of parametric identification of the Tersoff potential parameters are considered in papers [26,27]. The initial positions of atoms and their number are determined by the structure of the crystal lattice and the restriction on the minimum size of the simulated space is specified by the structure of the potential.
Consider crystalline silicon as a typical material. The simulated sample is represented by eight unit cells of a diamond crystal lattice making up a cube of 2 × 2 × 2 in size, with periodic boundary conditions; each unit cell contains eight unique atoms (Figure 14), so a total of 64 atoms are involved in the simulation. Initial speeds are considered to be zero. The initial conditions for a dynamical system can be represented as follows: To take into account deformation, 4 additional variables are introduced: one of them reflects the elongation value, and three more are responsible for the angle between the orientation of the lattice and the direction of stretching. In this case, the elongation is set to be fixed, and the variables responsible for the rotation are taken as interval. Rotation is To take into account deformation, 4 additional variables are introduced: one of them reflects the elongation value, and three more are responsible for the angle between the orientation of the lattice and the direction of stretching. In this case, the elongation is set to be fixed, and the variables responsible for the rotation are taken as interval. Rotation is generated evenly using quaternions [28].
The final system looks like this: Integration of the resulting ordinary differential system (13) was carried out using the Verlet velocity method with a constant integration step of 10 −15 s. As a result, the interval stress tensor was obtained (values are given in Pascals): Note that the possibilities of the proposed approach are not limited to a specific type of interatomic interaction potential in a material. The method can be applied to the simulation of the stress-strain state of materials with various types of chemical bonds, including the modeling of composite materials.

Discussion
In the previous sections, the proposed approach was tested on representative interval problems of nonlinear dynamics and computational materials science. It is found that, thanks to sparse grids, it is possible to integrate ODE systems with a large number of interval uncertainties in a reasonable time. To estimate the computational costs, a criterion was used that is equal to the equivalent number of sampling points from the original uncertainty region. Table 1 shows estimates of the computational costs for the ODE system given by Equation (9) describing a conservative oscillator with two interval initial conditions for two values ε. For ε = 10 −3 , the approach proposed in the paper works 1.5 to 5 times faster than the classical adaptive interpolation algorithm. Table 2 shows the computational costs when integrating the ODE system given by Equation (10), describing the Lotka-Volterra model with two interval initial conditions and one interval parameter. Here, for ε = 10 −3 , the use of adaptive sparse grids gives an acceleration of 1.9 − 2.8 times compared to the classical algorithm, and for ε = 10 −5 , the acceleration is from 1.25 time to 15 time.
For ODE systems given by Equations (11) and (12), it was possible to obtain a solution in a reasonable time only using the approach proposed in the paper since the number of interval uncertainties is quite large. To solve system given by Equation (11), the equivalent number of sampling points was about 80 thousand, and in the case of using the classical algorithm with a degree of polynomial 4, the value would be about 10 million. Thus, using sparse grids in this problem gives an acceleration of at least 125 times.
The tables show that adaptive sparse grids work faster than regular sparse grids, and even faster than full grids. This fact is in line with the theoretical foundations. The classical adaptive interpolation algorithm in example (9) with two interval uncertainties showed itself slightly better only when ε = 10 −5 and p = 4. This is primarily due to the chosen value of p. It is known that the greater the degree of the interpolation polynomial, the faster the error decreases with increasing mesh density. Therefore, it seems promising to use sparse grids on a nonlinear basis. We should also note the possibilities of applying adaptive grids to ODE systems not only with interval uncertainties but also with stochastic uncertainties, including applied nonlinear systems.

Conclusions
The adaptive interpolation algorithm allows simulating dynamic systems with interval parameters. In the course of the algorithm operation, a regular grid is constructed in the parameter space. The number of grid nodes depends exponentially on the number of interval parameters, which limits the scope of the algorithm. A typical situation is when the degree of influence of different parameters and their combinations on the solution can differ significantly. This can be used in adaptive interpolation. The paper presents an adaptive interpolation algorithm on sparse grids, which allows for reducing the exponential complexity when solving multidimensional problems in parameter space. The efficiency of the proposed approach has been demonstrated on representative interval problems of nonlinear dynamics and computational materials science. It is shown that, for most variants, adaptive sparse grids are more efficient than the classical adaptive interpolation algorithm in terms of computational costs. With the suggested method, it was possible to solve problems with up to 10 interval parameters in a reasonable amount of time. At the same time, the classical algorithm of adaptive interpolation failed to do so. Taking into account that an increase in the degree of the interpolation polynomial in the classical adaptive interpolation algorithm leads to higher accuracy and lower computational costs, we can outline the use of sparse grids with a nonlinear basis as a direction for further research.