1. 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.
Papers [
11,
12,
13,
14] describe the adaptive interpolation algorithm in detail. 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 theoretical basis of the algorithm is given in References [
11,
12]. The algorithm has a number of essential features: it is not subject to the wrapping effect [
11]; efficiently parallelizes on GPUs; able to simulate rigid systems [
13]; determines the presence of bifurcations and chaos in the system [
14]. It has been tested on applied problems of chemical kinetics [
13], gas dynamics and celestial mechanics, and complex dynamics with bifurcations and chaos [
14]. Nevertheless, there is some drawback. The algorithm uses multidimensional interpolation on a regular grid, which requires
nodes, where
—a polynomial degree for each dimension and
—the number of dimensions. With a large number of interval parameters, the application of the algorithm becomes difficult. However, a typical situation is when the degree of influence of different parameters and their combinations on a solution can differ significantly; therefore, it naturally follows to use approaches that take into account these features and, as a consequence, reduce the computational complexity.
Within the framework of modeling methods for dynamical systems with interval parameters, it is worth noting the work [
15], which describes a method based on the polynomial approximation of the solution, which requires points in the sample less than
. This method has been successfully applied to the problem of modeling a rotating system with both random and interval variables [
16]. This class of problems is significant from an applied point of view.
The curse of dimension, that is, the exponential growth of the number of calculations, is a critical problem. Typically, this situation arises when studying multidimensional functions presented in the form of a black box. The general tactic for reducing computational complexity is to determine and take into account the features of the function under consideration.
Sparse grids [
17] are numerical methods for representing, integrating, or interpolating multidimensional functions based on a hierarchical basis [
18,
19] and reducing the curse of dimension. This approach was first presented by the Russian mathematician Smolyak in 1963. Classic sparse grids result from computational cost optimization for approximating functions with bounded mixed derivatives [
20]. This fact is important since it imposes certain restrictions on the solution’s dependence on interval parameters. Interpolation using sparse grids requires significantly fewer nodes than standard full grid interpolation.
There are many works devoted to sparse grids [
21,
22,
23,
24]. Reference [
21] gives an initial introduction to sparse grids and the technique of combining them. It provides a program code in the Python programming language. In Reference [
22], some parallelization issues are considered; Reference [
23] provides an overview of the foundations and applications of sparse grids, with particular attention to the solution of partial differential equations.
The behavior of the solution to the ordinary differential equations (ODE) system can differ significantly depending on the parameters and initial conditions. Adaptive grids can drastically reduce computational costs by condensing nodes in regions with strong solution dependence on parameters and rarefaction in areas with weak dependence. Besides such adaptation, additional properties of the solution can be taken into account using sparse grids. This approach is effective when the interpolated function has a weak dependence on subsets of variables. For example, if the solution to an ODE system can be represented as a linear combination of functions from certain subsets of parameters and initial conditions, then it is sufficient to consider only the corresponding subsets and construct a grid only from them. Sparse grids are especially effective in multidimensional problems and can significantly reduce computational costs.
The main problem is the high computational costs when solving problems with uncertainties. The main goal of this work is to apply the previously proposed adaptive interpolation algorithm to the case of dynamical systems with a large number of interval parameters. The novelty lies in the application of sparse grids to problems with interval uncertainties, including problems of molecular dynamics.
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.
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
interval initial conditions:
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 meets all conditions ensuring the uniqueness and existence of a solution for all , .
The goal is, for each moment of time , to construct a piecewise polynomial vector function , where , , which interpolates the dependence of the solution on the interval parameters with controlled accuracy. If the function 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 is known at the moment of time , where , . 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 , the transfer of all non-interval solutions contained in the grid nodes to the time layer is performed, followed by the adaptation of the grid and the construction of an interpolation polynomial .
A short description of sparse grid interpolation according to the works [
20,
21] is given below.
Consider the interpolation of a smooth function of one variable on the unit interval . 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
on a unit interval
, where
is the level that determines the grid width
. The grid points
are given as:
Families of basis functions
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
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
in terms of the basis level functions
:
In this case, the interpolation polynomial given by Equation (4) takes the following form:
Next, consider the multidimensional interpolation of a smooth function
using
—dimensional unit cube
, provided that
. A multidimensional basis is constructed by the direct product of hierarchical one-dimensional bases:
where
are the levels of the corresponding one-dimensional grids,
is the basis function multi-index,
. If
, there is a sparse grid of
n level; if
, there exists a complete grid (
Figure 3). The number of nodes in a sparse grid is estimated as
, and the interpolation error is estimated as
; for a full grid the respective number of nodes is
, and the error is
, where
is the number of nodes in each dimension [
20].
The interpolation polynomial is written as follows:
where
In the case when the interpolated function has a nonzero value at the boundary, the one-dimensional basis is supplemented by two additional functions:
and
(
Figure 4). Two values are added to the polynomial given by Equation (5):
and
, where
. By analogy, for multidimensional interpolation, it follows that if
, then
in Equation (6) and
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 . If the value of the corresponding coefficient , where is some predetermined value, then each vertex creates descendants, which correspond to the basis functions of the next level. This process continues recursively until the values 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
and the resulting adaptive grid,
Figure 7 shows grids for functions
.
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
-dimensional unit cube. Taking into account that solving the problem requires interpolating
functions at once (
is the number of phase variables of the system), Equations (6) and (7) will take the following form:
where
The vector norm (for example, the maximum one) can be used as a criterion for adapting the grid.
Construct an interpolation polynomial . All the solutions that participated in the calculation of the coefficients given by Equation (8) are transferred to the -th time layer using some numerical integration method, after which a new set of coefficients is calculated and the adaptation of the grid is performed. When compacting the grid, the addition of new basis functions occurs at the -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 , , , 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:
where
is the number of nodes at the
step. A similar criterion exists for the classical adaptive interpolation algorithm [
11]. The
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,
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
. 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
and
(
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:
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 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
.
This model describes predator–prey interactions. A feature of the system is the fact that at there is an unstable focus and the amplitude of fluctuations in the population of species grows, and at 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
.
Consider an ordinary differential system presenting the expanded Volterra-Lotka model with three interval initial conditions and seven interval parameters:
Figure 12 shows the dependences of the interval estimates of solutions on time.
For a reasonable time, the solution was obtained using only adaptive sparse grids. For , the obtained result was and .
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
is the distance between two bodies,
is the initial velocity of bodies,
,
are body masses,
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
, the obtained result was
and
.
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
thousand, and in the case of using classical interpolation with the degree of polynomial equal to 4, the value would be of order
. 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
. 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
points, and the method proposed in Reference [
15] requires
points, where
is the degree of the interpolation polynomial.
4. 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
is the radius vector of the atom with the number
,
is its velocity,
is its mass, and
is the force acting on the atom, in this case
, where
is the total energy of the system, and
is the gradient along the position of the atom with the number
.
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
is total system energy,
is the contribution to the interaction energy of atoms with numbers
and
,
is the distance between atoms with numbers
and
,
is a cut-off function,
and
are the repulsion and attraction potentials, respectively, and
,
,
,
,
,
,
,
,
,
,
,
,
and
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
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:
where
is the linear size of a unit cell of a silicon crystal.
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:
where
is the mass of atoms,
is the relative elongation of the sample,
,
, and
are stretching direction parameters,
,
,
,
,
,
,
,
,
,
,
,
, and
are the parameters of the potential.
Integration of the resulting ordinary differential system (13) was carried out using the Verlet velocity method with a constant integration step of
. As a result, the interval stress tensor was obtained (values are given in Pascals):
For , the obtained result was and .
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.
5. 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
, the approach proposed in the paper works
to
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
, the use of adaptive sparse grids gives an acceleration of
times compared to the classical algorithm, and for
, the acceleration is from
time to
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 thousand, and in the case of using the classical algorithm with a degree of polynomial , the value would be about million. Thus, using sparse grids in this problem gives an acceleration of at least 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 and . 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.