CC BY
8Russian Journal of Nonlinear Dynamics, 2020, vol. 16, no. 3, pp. 479-490. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd200306
MATHEMATICAL PROBLEMS OF NONLINEARITY
MSC 2010: 65G40, 65L05, 65L07, 65L11
Modeling of Dynamic Systems with Interval Parameters in the Presence of Singularities
In solving applied and research problems, there often arise situations where certain parameters are not exactly known, but there is information about their ranges. For such problems, it is necessary to obtain an interval estimate of the solution based on interval values of parameters. In practice, the dynamic systems where bifurcations and chaos occur are of interest. But the existing interval methods are not always able to cope with such problems. The main idea of the adaptive interpolation algorithm is to build an adaptive hierarchical grid based on a kd-tree where each cell of adaptive hierarchical grid contains an interpolation grid. The adaptive grid should be built above the set formed by interval initial conditions and interval parameters. An adaptive rebuilding of the partition is performed for each time instant, depending on the solution. The result of the algorithm at each step is a piecewise polynomial function that interpolates the dependence of the problem solution on the parameter values with a given precision. Constant grid compaction will occur at the corresponding points if there are unstable states or dynamic chaos in the system; therefore, the minimum cell size is set. The appearance of such cells during the operation of the algorithm is a sign of the presence of unstable states or chaos in a dynamic system. The effectiveness of the proposed approach is demonstrated in representative examples.
Keywords: interval ODE, bifurcations, instability, dynamic chaos, adaptive interpolation algorithm, kd-tree
Received March 10, 2020 Accepted May 19, 2020
Alexander Yu. Morozov morozov@infway.ru
Federal Research Center Computer Science and Control of the Russian Academy of Sciences
ul. Vavilova 44, bld. 2, Moscow, 119333 Russia
Moscow Aviation Institute (National Research University)
Volokolamskoe sh. 4, Moscow, 125993 Russia
Dmitry L. Reviznikov reviznikov@mai.ru
Moscow Aviation Institute (National Research University) Volokolamskoe sh. 4, Moscow, 125993 Russia
Federal Research Center Computer Science and Control of the Russian Academy of Sciences ul. Vavilova 44, bld. 2, Moscow, 119333 Russia
A. Yu. Morozov, D. L. Reviznikov
1. Introduction
In solving applied and research problems, there often arise situations where certain parameters are not exactly known, but there is information about their ranges. For such problems, it is necessary to obtain an interval estimate of the solution based on interval values of parameters. There has been a large amount of research devoted to the relevant methods [4, 5, 13, 20]. There are methods based on the representation of a set of solutions through geometric primitives: parallelepipeds and ellipses [3, 8]; methods based on symbolic calculations [11, 18], as well as stochastic (Monte-Carlo) methods.
In practice, problems where bifurcations and chaos take place [6, 10] are of interest. But the existing interval methods are not always able to cope with them successfully.
In this paper, a modification of the adaptive interpolation algorithm is considered, which allows us to model dynamic systems with interval parameters in the presence of singularities. Also, during the operation of the algorithm, these singularities are identified. The main idea of the algorithm is to build a piecewise polynomial function which interpolates the dependence of the solution on the values of parameters specified by intervals. Bifurcation is usually understood as a qualitative change in the behavior of a system with smooth change of parameters [1]. In the vicinity of the bifurcation point a feature arises that does not allow the solution to be expanded as a power series of the parameters. A distinctive feature of chaos is a significant dependence on initial conditions. This means that initially close trajectories will diverge dramatically after some time. In both of these cases, an attempt to construct a function that interpolates the dependence of the solution on parameters or initial conditions will fail. The proposed modification of the algorithm is based on this idea.
A complete description of the adaptive interpolation algorithm for modeling dynamic systems with interval parameters is given in [14-17]. A theoretical basis of the algorithm is given in [14]. In [15] we show that the global error of the algorithm is directly proportional to the kd-tree height. The effectiveness of the algorithm is demonstrated by problems of chemical kinetics. In [16], various aspects of parallelization and implementation of the algorithm using CUDA technology are described. The present work is devoted to the application of the adaptive interpolation algorithm to dynamical systems with bifurcations and deterministic chaos.
The main idea of the algorithm is to build an adaptive hierarchical partitioning (Adaptive Mesh Refinement, AMR) based on the kd-tree over the uncertainty area. Each cell of kd-tree contains an interpolation grid. Depending on the features of the solution, a rebuilding of the partition is performed for each time instant. This is how the adaptability manifests itself. The result of the algorithm at each time instant is a piecewise polynomial function that interpolates the dependence of the solution on the point values of the parameters with a given precision.
An additional parameter, the minimum size of the partition cell, is introduced into the adaptive interpolation algorithm. If cells of this size appear during the operation of the algorithm, this is a sign of the presence of unstable states in the dynamic system. If the number of such cells is constantly increasing, it is possible to judge the presence of dynamic chaos in the system.
2. Adaptive interpolation algorithm
The process of modeling a dynamic system containing interval parameters can be represented as a sequence of mappings:
Vk+i = Fk (Vk"),
yo e yo = (y1,y20,...,ym,vm+l,---,Vo)T;
k = 0,...,N - 1;
where Fk: Rra ^ Rra, y0 E IRra are initial interval conditions. We note that the vector function Fk satisfies the conditions that ensure the existence and uniqueness of the solution for all y0 E y0. Hereinafter, the intervals are marked in bold, and their set is denoted by IR = {x = = [*£) x E R}.
For a discrete dynamic system, Fk describes the law of transition from one state to another.
A dynamic system can also be specified using a system of ordinary differential equations (ODE). The interval formulation of the Cauchy problem in this case takes the form
V = f (u),
u(o) = uo e uo = (u0, u20,...,um ,umm+1 ,...,u )
t E [0,îjV].
Here Fk can be interpreted as a solver of the Cauchy problem, which "moves" the solution from the kth time layer to the (k + 1)th time layer.
For further description of the algorithm, we will assume that the dynamic system is autonomous and the intervals are contained in the initial conditions only. Otherwise, additional dummy equations are added to the system and the intervals, together with the independent variable, are transferred to the initial conditions.
The purpose of the algorithm under consideration is to build for each k a piecewise polynomial function that interpolates the dependence of yk on specific values of interval parameters. This function will allow us to calculate a solution quickly for any starting point and, as a result, determine the interval estimate y&.
Let us suppose that at the kth step we know the solution yk(y0). We build an interpolation grid over the space formed by initial interval conditions. Let us suppose that it is regular and uniform:
yr-tm = (yl + + ... + y^K^vr1, • • • ,yo)T,
yk1i2-im = yk (y0li2 -im),
where p is a degree of the interpolation polynomial.
We calculate for each yk i2-im the value of yk+T"^ = Fk (yk i2-im) = = Fk(Fk-T(... F0(yt01 i2-im)...)). According to the obtained table function
G2+1 = {(»r-<m,»K3rim)|ii = M, i=T^}.
we build the interpolation polynomial, Pk+ T(y), in the Lagrange form [19]:
P0+i(y) = E... it yk+i"imli1i2"'im(y),
il=0 im =0
where
j=0 l=0,ii=j 11 - j If a posteriori interpolation error
error = maxOTey0 ||yk+i(y0) - p°+i(y0)|| _RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2020, 16(3), 479 490
is greater than a given value, e, then G0 is split into two grids G\ and G\, so that their estimate of interpolation error is smaller than error. All the same actions could be performed for Gk and Gk, and they could be split as well, if necessary. If a yk+1(y0) is continuously p + 1 times differentiable by y0, then it can be shown that the splitting process is finite. In general, from a theoretical point of view, the interpolation error is estimated using higher derivatives with proportionality coefficients, which characterize the interpolation grid and its size [14]. With each splitting, these coefficients decrease and, as a result, the interpolation error decreases. In practice, the minimum cell size is determined, upon reaching which the fragmentation does not occur any further, and the corresponding area of space is marked as "discontinuity area". If cells of this size appear during the operation of the algorithm, it is a sign of the presence of unstable states in a dynamic system. And if the number of such cells is constantly increasing, then we can judge the presence of dynamic chaos in the system. When building an interval estimation of the solution, linear interpolation is used in the "discontinuity regions".
As a result, a kd-tree and the corresponding piecewise polynomial function interpolating the solution with a given accuracy will be obtained at the (k + 1)th step. The process of building a kd-tree is illustrated in Fig. 1. It is not necessary to build a kd-tree from scratch at each step; instead, the tree obtained at the previous step is used, and depending on the estimate of the interpolation error, it is rebuilt.
The process of splitting vertices always occurs at the previous step (dotted lines), because when new vertices are created, interpolation of the values associated with their nodes is performed, which must be carried out at the moment when the error is still acceptable. If for a vertex and all its descendants the interpolation error becomes acceptable, then the descendants are removed and the vertex itself becomes a leaf node.
k = 0
k = 1
k = 2
k = 3
2/o
Gh
Vk
Gl
Gl Gf
Gj GS
@
/ \
G' 1S| ifi^ „ 0,1
y>
' s-i 3 I / 41 1(15 I I r< 6 >
g| g27 g|
Gl gs
' 1 < Ci ^
g| g30 g32
g3 g31
g3 gs
I r® 1 int-OK '/"111 1/^12*
'J V "î; >^3;
Fig. 1. Illustration of the algorithm operation.
Additionally, the maximum absolute value for phase variables is determined. If all solutions that are in the same interpolation grid increase without bound, then the vertex containing this grid is removed from the kd-tree, and the corresponding area of space is marked as "the area of infinity".
In practice, the estimation of the interpolation error is performed only for some points from the region that correspond to each particular vertex of the tree. We can distinguish two approaches to their choice. The first is to add randomly a test set of points to each vertex when it is created (Fig. 2a). The second approach is based on lowering the interpolation order, when some nodes and the values at those nodes are obtained with reduced order interpolation using the values at the remaining nodes (Fig. 2b). Although the error value in the second case is overestimated, it requires a smaller amount of computational resources, so we turn our attention to this alternative.
O O
o o
o o
-9-
I I I
----Q----<j>----6----O
A
I l I
-o-
(a) (b)
Fig. 2. Selection of the points at which the interpolation error is estimated.
Figure 3 shows the geometric interpretation of interpolation error in the Euclidean norm for a single vertex. All lines are obtained by interpolation over shaded points and correspond to the lines in Fig. 2b.
Fig. 3. Geometric interpretation of an interpolation error for one vertex. RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2020, 16(3), 479-490
-w
3. Results
We consider a one-dimensional dynamic system described by the equation
/ \ 3 5
x — Axx I x x ,
where A is an external control parameter. For this system, it is possible to find all points of bifurcations analytically and construct a bifurcation diagram (Fig. 4). There is a severe loss of stability. The dashed lines indicate unstable states.
-1 -0.5 0 0.5 1
A
Fig. 4. Bifurcation diagram.
Using the adaptive interpolation algorithm, we solve the following Cauchy problem:
(x' — Ax + x3 — x5; \x(0) — xo e [—2, 2],
where A e [—1,1].
Figure 5 shows the dependence of the solution on the interval parameters of the problem. The upper figures show the dependence of the solution on real values x0 and A, the lower figures show the resulting splits of the uncertainty region during the operation of the algorithm. The cells with the minimum size (hereinafter referred to as "bad" cells) are highlighted in black. The resulting black lines correspond to the dotted lines in the bifurcation diagram and indicate unstable states.
Thus, it is possible to identify the presence of unstable states in a dynamic system. A Duffing oscillator is the simplest model of an oscillator with reactive nonlinearity, which is given by the equation
x'' + ax' + w0(1+ x2)x — F (t), (3.1)
where a is a dissipation rate, w0 is a natural frequency, and F(t) is an external force.
According to Ref. [2], this equation can be obtained by considering oscillations of a mathematical pendulum at small deflection angles, oscillations of a weight suspended on a spring with a nonlinear restoring force located on a horizontal surface or when describing the motion of a particle in two potential wells, for example. Reference [9] shows that a Duffing oscillator is a system with an excitation and two stable equilibrium states that are separated by a single separatrix. At certain parameters, chaotic solutions arise near the separatrix.
t= 1
t = 5
x(xo,\)
t= 10
£ = 40
Fig. 5. The dependence of the solution on the initial conditions and parameters.
Substituting certain values a, and F(t) in (3.1), we write the Cauchy problem with two interval initial conditions [12]:
x = y,
y' = x - 0.25y - x3 + 0.3 cos(t), x(0) G [-2, 2], y(0) G [-2, 2], t G [0, 8].
(3.2)
Figure 6 shows a variety of solutions at different time points. At the initial time, the solution is a square. During the process of integration, the structure of the set changes greatly; it twists and decreases in size. Figure 7 shows the time dependence of the interval solution. Kinks in the graphs are clearly visible here. This is because the point corresponding to the upper or lower bound of the solution region at a certain time instant jumps from one place to another and its trajectory is discontinuous. Figure 8 shows the resulting grids. Here there is a pronounced compaction of cells, which increases in size.
If the number of minimum size cells is constantly increasing during the computational process, then we can assume that there is chaos in the dynamic system.
Let us consider a model of convection currents in the atmosphere, which was investigated by Edward Lorenz in 1961 where a strange attractor arises (Lorenz attractor) [7]. The model is described by the ODE system:
' X = a (-x + y ), y' = rx - y - xz, z' = -bz + xy,
x(0) = xo G [1,1.1], y(0) = yo G [0,0.1], z(0) = 0, t G [0, 29],
(3.3)
H 0
Fig. 6. A set of solutions to the system (3.2) at different time instants.
Fig. 7. The dependence of the upper and lower estimates for the system solution (3.2) on time.
where b = 8/3, a =10 and r = 28. Figure 9 shows the projections of the solution region (3.3) on the coordinate planes at different points in time. The black dots indicate all nodes of the interpolation grids; the grey lines indicate portions of the solution region that correspond to cells of normal size.
Figure 10 shows the resulting grid obtained in the computational process, where the "bad" cells are highlighted in black.
£ = 0.0
£ = 0.5
£ = 1.0
Xo Xq XQ
Fig. 8. Space splits resulting in the process of solving the system (3.2).
"Bad" cells form "split lines", and their number exponentially increases over time. Solutions that are on different sides to these lines drastically diverge from each other in the phase space. This indicates a substantial dependence on the initial conditions, and it is the main feature inherent in chaotic dynamics.
4. Conclusion
In this paper, the modification of the adaptive interpolation algorithm is considered. The modification allows us to model dynamic systems with interval parameters in the presence of singularities. During the operation of the algorithm, these singularities are identified. The basic idea is that it is impossible to build a polynomial function to interpolate the dependence of a solution on the parameters at the points where bifurcations occur. A modification of the adaptive interpolation algorithm consists in limiting the cell size of the adaptive grid. The emergence of the minimum size cells is a sign that there are unstable states in the system.
i = 25.4
£ = 26.8
i = 28.1
-30 -15 0 15 30 y
-30 -15 0 15 30 y
-30 -15 0 15 30 y
Fig. 9. The set of solutions to the system (3.3) at different points in time.
¿ = 24.0 ¿ = 25.4 ¿ = 26.1
Xq XQ XQ
Fig. 10. Space splits resulting from the solution of the system (3.3).
If consolidations of such cells constantly increase in size and do not move in space, this indicates a significant dependence of the solution on the parameters, which, in turn, is a characteristic feature of dynamic chaos. The operation of the modified adaptive interpolation algorithm is illustrated by some classical examples.
References
[1] Arnold, V. I., Afrajmovich, V. S., Il'yashenko, Yu. S., and Shil'nikov, L.P., Bifurcation Theory and Catastrophe Theory, in Dynamical Systems V. I. Arnold (Ed.), Encyclopaedia Math. Sci., vol. 5, Berlin: Springer, 1999.
[2] Astakhov, V. V., Koblyansky, S.A., and Shabunin, A. V., Duffing Oscillator, Saratov: Saratov. Gos. Univ., 2007 (Russian).
[3] Chernousko, F.L., State Estimation for Dynamic Systems, Boca Raton, Fla.: CRC, 1993.
[4] Dobronets, B.S., Interval Mathematics, Krasnoyarsk: Krasnoyarsk. Gos. Univ., 2004 (Russian).
[5] Eijgenraam, P., The Solution of Initial Value Problems Using Interval Arithmetic: Formulation and Analysis of an Algorithm, Mathematical Centre Tracts, No. 144, Amsterdam: Stichting Mathematisch Centrum, 1981.
[6] Kubyshkin, E. P. and Moriakova, A. R., Features of Bifurcations of Periodic Solutions of the Ikeda Equation, Russian J. Nonlinear Dyn., 2018, vol. 14, no. 3, pp. 301-324.
[7] Kronover, R. M., Fractals and Chaos in Dynamical Systems: Fundamentals of Theory, Moscow: Postmarket, 2000 (Russian).
[8] Lohner, R. J., Enclosing the Solutions of Ordinary Initial- and Boundary-Value Problems, in Comput-erarithmetic. Scientific Computation and Programming Languages, E. Kaucher, U. Kulisch, Ch. Ullrich (Eds.), Stuttgart: Teubner, 1987, pp. 255-286.
[9] Loskutov, A. Yu., Dynamical Chaos: Systems of Classical Mechanics, Physics-Uspekhi, 2007, vol.50, no. 9, pp. 939-964; see also: Uspekhi Fiz. Nauk, 2007, vol. 177, no. 9, pp. 989-1015.
[10] Maglevanny, 1.1., Smolar, V. A., and Karyakina, T.I., Weak Signals Amplification through Controlled Bifurcations in Quasi-Two-Dimensional Electron Gas, Russian J. Nonlinear Dyn., 2018, vol. 14, no. 4, pp. 453-472.
[11] Makino, K. and Berz, M., Verified Computations Using Taylor Models and Their Applications, in Numerical Software Verification (NSV 2017), A. Abate, S.Boldo (Eds.), Lect. Notes Comput. Sci., vol.10381, Cham: Springer, 2017, pp. 3-13.
[12] Makino, K. and Berz, M., Rigorous Reachability Analysis and Domain Decomposition of Taylor Models, in Numerical Software Verification (NSV 2017), A. Abate, S.Boldo (Eds.), Lect. Notes Comput. Sci., vol.10381, Cham: Springer, 2017, pp. 90-97.
[13] Moore, R.E., Interval Analysis, Englewood Cliffs, N.J.: Prentice-Hall, 1966.
[14] Morozov, A.Yu. and Reviznikov, D.L., Adaptive Interpolation Algorithm Based on a KD-Tree for Numerical Integration of Systems of Ordinary Differential Equations with Interval Initial Conditions, Differ. Equ., 2018, vol.54, no. 7, pp. 945-956; see also: Differ. Uravn., 2018, vol.54, no. 7, pp. 963-974.
[15] Morozov, A.Yu., Reviznikov, D.L., and Gidaspov, V. Yu., Adaptive Interpolation Algorithm Based on a KD-Tree for the Problems of Chemical Kinetics with Interval Parameters, Math. Models Comput. Simul., 2019, vol.11, no. 4, pp. 622-633.
[16] Morozov, A.Yu. and Reviznikov, D.L., Modelling of Dynamic Systems with Interval Parameters on Graphic Processors, Programmnaya Ingeneria, 2019, vol. 10, no. 2, pp. 69-76.
[17] Morozov, A.Yu. and Reviznikov, D.L., Methods of Computer Simulation of Dynamic Systems with Interval Parameters, Moscow: MAI, 2019 (Russian).
[18] Rogalev, A. N., Guaranteed Methods of Solving of Ordinary Differential Equations on the Basis of Symbolical Formulas Transformation, Vychislitel'nye Tekhnologii, 2003, vol.8, no. 5, pp. 102-116 (Russian).
[19] Sauer, Th. and Xu, Y., On Multivariate Lagrange Interpolation, Math. Comp., 1995, vol. 64, no. 211, pp.1147-1170.
[20] Shary, S. P., Finite Dimensional Interval Analysis, Novosibirsk: XYZ, 2019 (Russian).
