**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–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.

**Citation:** Morozov, A.Y.; Zhuravlev, A.A.; Reviznikov, D.L. Sparse Grid Adaptive Interpolation in Problems of Modeling Dynamic Systems with Interval Parameters. *Mathematics* **2021**, *9*, 298. https://doi.org/ 10.3390/math9040298

Academic Editor: Mikhail Posypkin Received: 18 December 2020 Accepted: 28 January 2021 Published: 3 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Papers [11–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 (*p* + 1) *<sup>d</sup>* nodes, where *p*—a polynomial degree for each dimension and *d*—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 (*p* + 1) *d* . 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–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.
