Regularization of inclusions of differential equations solutions based on the kinematics of a vector field in stability problems

This paper presents the new results of computing the inclusions of sets of solutions of ordinary differential equations corresponding to perturbations acting on the solutions of the system. The regularization of algorithms for the inclusion of solutions is investigated. Inclusions of solutions are used to study the stability over a finite time interval under perturbations of the system parameters. The boundaries of solutions sets using methods that construct symbolic formulas that characterize the behavior of the system are computed. In this case, the Influence of permanent disturbances on the solutions is taken into account. It is proposed to use the parameters of the vector field of the problem in order to compensate for a strong growth of the computable solution boundaries, which is often encountered in many methods. This means the regularization of the problem of estimating inclusions of solution sets.


The practical stability or the stability over a finite time interval
The problem of studying the trajectories of sets of mechanical systems has appeared in the published works of motion stability. The concept of stability of a set of solutions has common features with the well-known definitions of the Lyapunov stability, the practical stability, and the stability over a finite time interval [1]- [6], but it has many differences. The practical or technical stability over a finite time interval means that solutions are uniformly bounded with respect to a set of initial values and a set of disturbing influences. The Lyapunov direct method, which implements the idea of small deviations of solutions of a differential equation over a time interval with small variations of the initial data of this solution, is not sufficient for solving such problems.
The practical or technical stability over a finite time interval means that solutions are uniformly bounded with respect to a set of initial values and a set of disturbing influences. The Lyapunov direct method, which implements the idea of "small" deviations of solutions of a differential equation over a time interval [0, +∞) with "small" variations of the initial data of this solution, is not sufficient for solving these problems. A real phenomenon or a process often operates on a finite time interval and is influenced by arbitrary perturbations (large or small), often for these perturbations only the boundaries of possible values are known. With perturbations of the initial data for the stability analysis, it is necessary to track the proximity of solutions y 1 t; t 0 , y 0 1 and y 2 t; t 0 , y 0 2 with the initial data. Stable solutions should slightly differ from one another for any t > t 0 , if y 0 1 and y 0 2 are sufficiently close. An essential here 2 requirement is, that the solutions are close for all values t > t 0 , since if the consideration is limited to a finite interval of time variation, then the indicated property of preserving proximity will be a direct consequence of the classical theorem on the continuous dependence of solutions of a system of differential equations on the initial data.
In this paper, we consider the problems of motion stability under initial conditions given by their boundaries [6]. Consider the Cauchy problem for equations with initial data given by their boundaries: where y 0,l , y 0,u ∈ R n are the vectors of the upper and the lower bounds of the initial values y 0 , f (t, y) ∈ C R + × R n , R {n , and f (t, 0) = 0 for all t ≥ t 0 .
We assume that the motions of system (1) under initial conditions (2) are defined for all t ∈ R + .
In what follows, Y (t) is a set of trajectories of system (1) with initial data given by their boundaries, i.e.
The study of a set of trajectories is complicated by the fact that, in addition to stable zero solutions of system (1), a set of trajectories can contain singular points, separatrices, and limit cycles. For systems of equations, whose dimension is n > 2, system (1) may have chaotic attractors. All this is the reason for the behavior of the trajectories to be hard to predict near the vicinity of the attractor.
In [6] matrix inequalities to determine stability are used. There, is a reference to the article [7], as an example of the successful application of methods of numerical integration of motion problems with multivalued parameters. The publication [6] also provides a conclusion on the need to continue the research on this issue. This is also important because after the publication of [7], new results were obtained. This paper is aimed at the description the current methods for construction of external boundaries or inclusion boundaries of solutions or two-sided estimates of the sets of solutions of ordinary differential equations based on symbolic formulas in the stability problems.
Let us make a small survey of the stability studies with initial data given by estimates. The number of works is quite large, although there are no researches analyzing various methods for including solutions sets (bilateral boundaries of these sets). Perhaps one of the few exceptions is article [8].
The computation of the boundary of the ODE solutions was carried out by specialists in the use of the Taylor interval polynomials [9]. There are many publications describing the interval Taylor polynomials, but the methods do not solve this problem with sufficient accuracy [10]. The interval methods compute the boundaries of solutions with exponential growth for no reason [10], [9].
The names of their methods usually use the term "Taylor series". For the Taylor series, there is a convergence of estimates for solutions of certain classes of ODE systems. In fact, computations are performed with Taylor polynomials, that is, a finite number of terms of the Taylor series. This means that convergence must be proven. Moreover, all the justifications in these studies ignore this fact. Special terms wrapping, shrink wrapping have appeared, which allow for many explanations and translation options: thermo-packaging, thermo-shrinking material, stretch film, wrapping [10]. The meaning of this term is that there is a strong growth of the boundaries of the sets of solutions, which does not correspond to the actual behavior of the sets of such solutions.

Description of the symbolic algorithm
The problem is being considered, where y, f ∈ R n for each t ∈ T, y 0 ∈ Y 0 ⊂ R n . A set of all possible solutions (trajectories) are: Unfortunately, this set is usually impossible to accurately compute for many reasons. Instead, we compute approximations to this set. In order to be able to check the stability of the system, in any form, this approximation must be an external inclusion.
To eliminate most of the difficulties in the implementation of methods for estimating solution sets, it is proposed to use symbolic formulas. Symbolic formulas are combinations of signs of operations, constants and variable values that have an independent meaning and represent the algorithm for computing the value of an expression. A possible definition of symbolic formulas is to write down expressions to determine the range of values of a quantity through its parameters.
An efficient numerical algorithm based on the symbolic formulas in the class of problems under consideration is determined by the fact that the computation of the range of values of the expression (or their estimates) by this formula is performed with the least time and the number of operations.
An algorithmic scheme for estimating the sets of solutions of the ODE system Input: To choose an effective symbolic formula Sym [0,δ] (t i , Y 0 ) for plotting the range of a set of solutions, it is often necessary to know the geometry of decision plots, and other properties of the range. In this case, the most time-consuming operations will be operations for organizing and storing the symbolic formulas involved in computing estimates of ranges of values. The total number of these operations will be significantly less than the number of operations performed when computing the range of solutions.  Figure 1. External bounds of a solution set as well as internal bounds.
Let Sym [0,δ] (t i , Y 0 ) be a symbolic vector formula that reflects the dependence of the solution y(t) = ξ(t; y 0 ) on variable t of the original ODE system and on the components of the initial data vector y 0 , varying within the given regions (defined by their boundaries).
Then the expression denotes the range of values of the function y(t, y 0 1 , . . . , y 0 n ) in the specified range of the variables. Obviously, under the conditions of continuity and continuous differentiability of the symbolic formula, this set will be connected with a smooth boundary.
To effectively operate with these sets, it is possible to approximate them by the regions defined using polyhedra, parallelepipeds, polygons, balls, ellipses in the dimensional space. Figure 1 schematically depicts a set of solutions, a set inside and around this set and shows the boundaries of external inclusions, as well as internal inclusions.
Let, in the case of a linear system of ODEs the entries a ij (t) of the matrix A and u i t) of the vector u(t) be defined and continuous at t ≥ 0. To represent the symbolic formula for the solution, we will use the Cauchy formula. The fundamental matrix of system (4) is a nondegenerate square matrix function Φ (t) that satisfies for all t ≥ 0 the matrix systems of differential equations For system (7) the initial conditions are given: Then the Cauchy formula holds for the solution of the ODE system (7) with the initial conditions (8).
If the matrix of system (1) is constant, then the fundamental matrix has the following form, and the solution formula can be written down as Although the Cauchy formula is a symbolic formula for solutions, its direct application requires the implementation of some algorithm based on the Cauchy formula. Just writing an exponential matrix function is not practical.
The algorithm for finding a symbolic formula for a solution consists of the following stages. 1) Construct a symbolic formula for the fundamental matrix of solutions to system (7). Note that the matrix function will be a solution to the matrix system of equation (7) if and only if any of its columns is a solution to a linear homogeneous system.
2) Find symbolic formula of the inverse matrix for the fundamental solution matrix.
3) Multiply the symbolically inverse matrix for the fundamental matrix by the vector of the complementary term on the right-hand side. 4) Find the antiderivative of the resulting expression. 5) Completely collect formula (9) 6) Compite the estimate of the range of values of formula (9). The class of algorithms based on the symbolic formulas used, among other things, for problems of practical stability on a finite time interval, is also described in [11].
For system (1) with a nonlinear right-hand side, it is possible to construct a "perturbed" system associated with system (1) with changed initial data, the solutions of these systems satisfy the relation [12], [13] ∂y ∂y 0 (s, t, z(s))g(s, z(s))ds. (10) The difficulties arising in the problems of estimating the sets of solutions of ODE systems have led to the need analyzing this situation and creating the necessary algorithms. A strong increase in the bounds of estimates is observed in most assessment methods. It is possible to eliminate its influence in algorithms based on symbolic decision formulas.
3. Regularization of the symbolic algorithm based on the kinematic parameters of the vector field A.N. Tikhonov in [14] and in many other publications formulated the notion of a regularizing algorithm (regularizing operator, regularizer) for an ill-posed problem as a one-parameter family of operators that in a special way approximates the inverse operator and provides a stable restoration of the desired solution.
For the Hadamard-correct problems, the inverse operator it it is can be taken as a formal regularizing algorithm. Another thing is that such an algorithm may turn out to be nonconstructive (practically unrealizable). The concept of a regularizing algorithm turned out to be highly efficient [15].
For the methods of estimating sets of solutions, the kinematic parameters of the vector field and its derivative affinor of the solution domain as a regularizer are proposed [17]. Having estimated these kinematic parameters, it is possible to determine the reference points on the boundary of the solution sets, control the geometric properties of the solution domain, and approximate the inverse operator.
Let us describe the kinematic parameters of the initial vector field (a mapping that sets each point M of the phase space a vector with the origin at this point), given by the formula The region Ω in which the vector field is specified will be assumed to be filled with some mobile deforming medium, for example, a liquid. Let the vector f (M ) describe the velocity with which a particle of the liquid moves, which is at a given moment at a point M . Each point, in [16], is called a "particle of liquid", describes a certain trajectory over time In this case, the projections of the velocity vector onto the coordinate axes are, as is known, On the other hand, the velocity vector coincides at each point M with the field vector f (M ), and its projections are equal f i (M ) = f i (y 1 (t), y 2 (t), y 3 (t)). As a result dy i dt = f i (y 1 , y 2 , y 3 ), i = 1, 2, 3.
A liquid particle located at a given point M with respect time, after an infinitely small time interval ε = ∆t will a shift by the vector ∆t · f (M ), if we neglect the infinitely small higher order. Since the velocity of movement of a fluid depends on its position, then in the process of movement the velosity will, generally speaking, change. However, in an infinitely small period of time, it manages to change, starting from the value f (M ) only by an infinitely small amount. As a result ∆t · f (M ), it gives a displacement with an error of an infinitesimal higher order (roughly speaking, an infinitesimal domain velocity ∆f is multiplied by Err = ∆t · ∆f.) Let M ′ be some point in a ball of infinitesimal radius ρ centered at M (a certain solution is chosen in the domailn of solutions). All points of the ball to the specified displacement τ · f (M ) along the path of the solution, in particular, the point M, M ′ go to some points L, L ′ .
In other words, infinitesimal vectors outgoing from the center of the neighborhood of solutions of radius ρ transform to vectors outgoing from the center of the displaced (and deformed) neighborhood, subject to the action of an affinor U, i.e. a linear law by which each vector in space is associated with a certain vector, denoted by y = U z, (the matrix of derivatives).
But the action of such an affinor is reduced to a pure deformation generated by the affinor E +∆t·B, and to rotation by means of the affinor E +∆t·C. In this case b ij , c ij − are symmetric In our kinematic interpretation, this means that the volumetric expansion occurs at zero velocity, i.e. any spatial region filled with liquid particles at the initial moment shifts and deforms over time without changing its volume (here this is shown only for infinitely small liquid drops, and the transition to a finite volume would require some more reasoning). The deformation (a change in shape) of the domain of differential equations solutions occurs with increasing time.
In order to find a useful information about a numerical solution, the condition number of the problem is calculated, that is, the coefficient of sensitivity of the solution to the problem perturbation.
Definition 3 (conditionality). Let f be a problem that has the data d and the solution y. Then the problem f is said to be well-conditioned if small changes in the data d lead only to small changes in the solution y, and f is ill-conditioned if small changes in the data lead to large changes in the solution y. The conditionality of a problem is usually characterized by a conditionality number that describes how much the data disturbance in the solution is amplified [17], [18]. One can determine the conditionality number of the problem f.
Thus, we can see that the condition number, like a derivative, measures the rate of a change in solutions relative to changes in the data d (maximum in a certain range of the data D) To regularize the problem of estimating sets of solutions to a system of ODEs, the Alekseev formula (10) can be chosen as the problem f . Using this formula, we estimate the conditionality number and then correct the operator of the problem, changing the boundaries of the set, changing the kinematic characteristics.. Let us consider an example of studying the stability of the motion of a technical system with disturbing influences. Region bounds of oscillations of a technical system with control actions u 1 , u 2 , y 1 , y 3 , i.e. generalized coordinates, y 2 , y 4 , i.e. generalized impulses are computed when solving the ODE system   Figure 3. The computed boundaries of sets of solutions of the system: the projection onto the plane t − y 1 .

Conclusion
A real phenomenon or process often operates on a finite time interval and is influenced by arbitrary perturbations (large or small), often for these perturbations only the boundaries of possible values are known. With perturbations of the initial data for the stability analysis, it is necessary to track the proximity of solutions with the initial data. Stable solutions should differ slightly from one another. The essential here is the requirement that the solutions are close for all values , since if the consideration is limited to a finite interval of time variation, then the indicated property of preserving proximity will be a direct consequence of the classical theorem on the continuous dependence of solutions of a system of differential equations on the initial data.
For applied problems, importance is not only the existence of a number ε for a given number δ that satisfy the definition of Lyapunov stability, but also the verification of the suitability of the estimates in specific conditions of the problem. Therefore, methods for solving stability problems that make it possible to obtain these estimates are the main ones. The symbolic methods considered in this paper effectively solve the problems of motion stability under the initial conditions given by their boundaries.