Introduction

Mathematical models of many applied processes in areas of complex structure are initial and boundary value problems for differential equations with a deviating argument. Such equations have found wide application due to the fact that they well describe many phenomena with consequences in technological processes, radio technical and electrical devices, economic and ecological systems. In many models, delay is introduced as a characteristic of poorly studied processes, which are not taken into account at this stage of model construction. All this became the reason why systems of differentialfunctional equations are a relevant object for research [1][2][3][4][5].

An important task is to study the stability of solutions of differential-difference equations. In connection with numerous practical applications, considerable attention is paid to obtaining stability conditions for linear differential equations with a delay [6][7][8].

Of particular interest are studies that allow using the methods of the theory of ordinary differential equations for the analysis of differential functional equations.

Schemes of approximation of differential-difference equations by special systems of ordinary differential equations, in particular, make it possible to use approaches to the study of ordinary dynamic systems and to solve a number of problems from the theory of stability of linear systems with a delay.

In this work, approximation schemes for initial problems for systems of linear differential-difference equations with many delays are investigated by the sequence of systems of ordinary differential equations, and algorithms for approximate finding of non-asymptotic roots of quasi-polynomials are constructed. Using these algorithms, a method of modeling the stability of solutions of linear systems with a delay was developed, and constructive computational algorithms were proposed for constructing coefficient regions of stability of linear systems with many delays [8][9][10].

Approximation schemes

Consider the initial problem for a system of differential equations with a delay 𝑑𝑥(𝑡) 𝑑𝑡 = 𝐹(𝑡, 𝑥(𝑡), 𝑥(𝑡 − 𝜏)), 𝑡 ≥ 𝑡 0 ,

𝑥(𝑡) = 𝜑(𝑡), 𝑡 ∈ [𝑡 0 − 𝜏, 𝑇 0 ],(1)

where 𝐹(𝑡, 𝑢, 𝑣) is a continuous function that satisfies the Lipshitz condition for 𝑢 + 𝑎𝑣, 𝜑(𝑡) is a given continuous initial function 𝑥 ∈ 𝑅 𝑛 , 𝜏 > 0.

An approximate replacement of the initial problem ( 1)-( 2) with the Cauchy problem for the system of ordinary differential equations is proposed in works [11][12][13]. It is based on the study of a sequence of connected 𝑚 delay elements and their replacement by a sequence of aperiodic links that carry out a shift by the amount ) , 𝑗 = 0, 𝑚. ̅̅̅̅̅̅

Theorem 1 [9]. If the solution of the problem (1)-( 2 ) → 0 when 𝑚 → 0, which means that the solution of the Cauchy problem (4)-( 5) approximates the solution of the initial problem for a system of differential equations with a delay (1)-(2).

Approximation nonasymptotic roots of quasipolynomials

Consider a linear system of differential-difference equations

𝑑𝑥 𝑑𝑡 = ∑ 𝐴 𝑖 𝑘 𝑖=0 𝑥(𝑡 − 𝜏 𝑖 ),(7)

where 𝐴 𝑖 , 𝑖 = 0, 𝑘 ̅̅̅̅̅ , , i=1,k i AB fixet 𝑛 × 𝑛 matrices, 𝑥 ∈ 𝑅 𝑛 , 0 = 𝜏 0 < 𝜏 1 < ⋯ < 𝜏 𝑘 = 𝜏. Qualitative properties of system (7) depend on the location of the zeros of its quasi-polynomial

Φ(𝜆) = 𝑑𝑒𝑡(𝜆𝐸 − ∑ 𝐴 𝑖 𝑘 𝑖=0 𝑒 −𝜆𝜏 𝑖 ) = 0. (8)

The zeros of the quasi-polynomial ( 8) are divided into two groups: asymptotic and non-asymptotic. Asymptotic formulas can be obtained for asymptotic roots. Non-asymptotic roots located near the origin of coordinates, determine the main properties of solutions of system (7) and can be found only with the help of approximate algorithms.

The study of the approximation schemes of linear systems with a delay showed that the nonasymptotic roots of their quasi-polynomials can be effectively approximated by the roots of the characteristic polynomials of the corresponding approximating systems of ordinary differential equations.

According to scheme (4)-( 5), we correspond to the system with delay (7) by the system of ordinary differential equations

𝑑𝑧 0 (𝑡) 𝑑𝑡 = ∑ 𝐴 𝑖 𝑘 𝑖=0 𝑧 𝑙 𝑖 (𝑡), 𝑙 𝑖 = [ 𝜏 𝑖 𝑚 𝜏 ],(9)𝑑𝑧 𝑖 (𝑡) 𝑑𝑡 = 𝜇(𝑧 𝑖−1 (𝑡) − 𝑧 𝑖 (𝑡)), 𝑖 = 1, 𝑚 ̅̅̅̅̅̅ , 𝜇 = 𝑚 𝜏 .

Lemma 1 [9,10]. Equality holds for the characteristic equation of the approximating system ( 9)

Ψ 𝑚 (𝜆) = 𝑑𝑒𝑡 (𝜆𝐸 − ∑ 𝐴 𝑖 𝑘 𝑖=0 ( 𝜇 𝜇+𝜆 ) 𝑙 𝑖 ) (𝜇 + 𝜆) 𝑚𝑛 = 0(10)

and sequence of functions

𝐻 𝑚 (𝜆) = Ψ 𝑚 (𝜆) (𝜇 + 𝜆) 𝑚𝑛 (11)

converges as 𝑚 → ∞ to the quasipolynomial (8).

Remark 2. [9,10]. The zeros of the functions Ψ 𝑚 (𝑥) and 𝐻 𝑚 (𝜆) coincide, so the zeros of the characteristic polynomial (10) approximate the non-asymptotic roots of the quasi-polynomial (8).

To increase the accuracy of the approximation of non-asymptotic roots, a scheme of increased accuracy is proposed

𝑑𝑧 0 (𝑡) 𝑑𝑡 = ∑ 𝐴 𝑖 𝑘 𝑖=0 𝑧 𝑙 𝑖 (𝑡) , (12) 𝑑𝑧 𝑗 (𝑡) 𝑑𝑡 = 𝑧 𝑚−𝑗 (𝑡), 𝑑𝑧 𝑚+𝑗 (𝑡) 𝑑𝑡 = 2𝜇 2 (𝑧 𝑗−1 (𝑡) − 𝑧 𝑗 (𝑡)) − 2𝜇𝑧 𝑚+𝑗 (𝑡), 𝑗 = 1, 𝑚 ̅̅̅̅̅̅ .

Lemma 2 For the characteristic equation of the approximating system of increased accuracy (2), equality holds

Ψ 2𝑚+1 (𝜆) = 𝑑𝑒𝑡 [(𝐴 0 − 𝜆𝐸) (1 + 𝜆𝜏 𝑚 (1 + 𝜆𝜏 2𝑚 )) 𝑚 + ∑ 𝐴 𝑖 𝑘 𝑖=1 (1 + 𝜆𝜏 𝑚 (1 + 𝜆𝜏 2𝑚 )) 𝑚−𝑙 𝑖 ] × × (1 + 𝜆𝜏 𝑚 (1 + 𝜆𝜏 2𝑚

))

𝑚𝑛 (13) and sequence of functions

𝐻 2𝑚+1 (𝜆) = Ψ 2𝑚+1 (𝜆) (1 + 𝜆𝜏 𝑚 (1 + 𝜆𝜏 2𝑚 )) 𝑚𝑛 (14)

converges as 𝑚 → ∞ to the quasipolynomial (8).

Example 1. Consider a scalar differential equation with a delay 𝑑𝑥(𝑡) 𝑑𝑡 = 𝑎𝑥(𝑡) + 𝑏𝑥(𝑡 − 𝜏), (15) where 𝑎, 𝑏, 𝜏 ∈ 𝑅, 𝜏 > 0.

The quasipolynomial for equation (15) has the form Φ(𝜆) = 𝑎 − 𝜆 + 𝑏𝑒 −𝜆𝜏 = 0.

(16) The approximate system of ordinary differential equations corresponding to scheme (9) for equation (15) 𝑑𝑧 0 𝑑𝑡 = 𝑎𝑧 0 + 𝑏𝑧 𝑚 ,

(17) 𝑑𝑧 𝑖 𝑑𝑡 = 𝜇(𝑧 𝑖−1 − 𝑧 𝑖 ), 𝑖 = 1, 𝑚 ̅̅̅̅̅̅ .

For the characteristic equation of system (17), we obtain equality

Ψ 𝑚 (𝜆) = (𝑎 − 𝜆) (1 + 𝜆𝜏 𝑚 ) 𝑚 + 𝑏 = 0.

After making the substitution 𝑆 = 1 + 𝜆𝜏 𝑚 , we obtain the equation.

𝑆 𝑚+1 − (1 + 𝜏 𝑚 𝑎) 𝑆 𝑚 − 𝜏 𝑚 𝑏 = 0, (18)

which is convenient for calculating its roots.

The high-accuracy approximating system of ordinary differential equations corresponding to scheme (12) for equation (15) has the form

𝑑𝑧 0 𝑑𝑡 = 𝑎𝑧 0 + 𝑏𝑧 𝑚 ,(19)𝑑𝑧 𝑗 𝑑𝑡 = 𝑧 𝑚+𝑗 , 𝑑𝑧 𝑚+𝑗 𝑑𝑡 = 2𝜇 2 (𝑧 𝑗−1 − 𝑧 𝑗 ) − 2𝜇𝑧 𝑚+𝑗 , 𝑗 = 1, 𝑚 ̅̅̅̅̅̅ .

For the characteristic equation of the system (19), the representation is valid

Ψ 2𝑚+1 (𝜆) = (𝑎 − 𝜆) (1 + 𝜆𝜏 𝑚 (1 + 𝜆𝜏 2𝑚 )) 𝑚 + 𝑏 = 0. (20)By replacing 𝜆 = 2𝑚 𝜏 (𝑆 −1 2

) in equation (20), the equation can be obtained in the standard form 𝛼 0 𝑆 2𝑚+1 + 𝛼 1 𝑆 2𝑚 + 𝛼 2 𝑆 2𝑚−1 + ⋯ + 𝛼 2𝑚 𝑆 + 𝛼 2𝑚+1 = 0, where

𝛼 2𝑖 = − 2𝑚 𝜏 𝐶 𝑚 𝑖 2 𝑚−2𝑖 , 𝑖 = 0, 𝑚 ̅̅̅̅̅̅ , 𝛼 2𝑖+1 = (𝑎 + 𝑚 𝜏 ) 𝐶 𝑚 𝑖 2 2𝑚−2𝑖 , 𝑖 = 0, 𝑚 − 1 ̅̅̅̅̅̅̅̅̅̅̅ , 𝛼 2𝑚+1 = (𝑎 + 𝑚 𝜏 ) 2 −𝑚 + 𝑏.

Consider the model equation with a delay 𝑑𝑥(𝑡) 𝑑𝑡 = 2𝑥(𝑡) + 𝑥(𝑡 − 1),

the characteristic quasipolynomial of which 𝜆 = 2 + 𝑒 −𝜆 .

(23) The real root of the quasi-polynomial (23) with the largest real part, found by dividing the segment in half, is equal to λ=0.120028.

We find the approximation of this root of the quasi-polynomial (23) using the roots of the characteristic polynomials of the approximating systems (18) and (21). The results of numerical experiments for different 𝑚 are shown in Table 1, where 𝜆 𝑖 𝑠𝑡. is obtained according to the standard scheme (17), and 𝜆 𝑖 𝑖𝑎. is obtained according to the scheme of increased accuracy (19). 1, we see that the scheme of increased approximation accuracy is significantly more efficient than the classical approximation scheme.

Table 1

Stability analysis of linear systems with a delay

The theory of stability for systems of differential-difference equations ( 7) is currently sufficiently well researched.

Theorem 2 [6,[14][15][16]. In order for the zero solution of system (7) to be exponentially stable, it is necessary and sufficient that all slopes of its quasipolynomial (8) lie in the left half-plane 𝑅𝑒 𝜆 < 0.

Thus, the location of the zeros of quasi-polynomials characterizes the stability of linear stationary systems with a delay. In practice, this check is possible only in the simplest cases.

We will use Lemma 1 and Lemma 2 to analyze the location of the roots of quasi-polynomials using the characteristic equations of the corresponding approximating systems of ordinary differential equations.

Theorem 3 [10]. If the zero solution of the system with delay (1) is exponentially stable (not stable), then there is 𝑚 0 > 0 such that for all 𝑚 ≥ 𝑚 0 , the zero solution of the approximating system (3) is also exponentially stable (not stable).

If for all 𝑚 ≥ 𝑚 0 the zero solution of the approximation system (3) is exponentially stable (not stable) then the zero solution of the system with a delay (1) is exponentially stable (not stable).

Remark 3. From Theorem 3, we have that there exists a number 𝑚 0 such that 𝑚 ≥ 𝑚 0 the asymptotic stability of the zero solution of the system with a delay ( 7) is equivalent to the asymptotic stability of the zero solution of the approximating system of linear differential equations (9).

Consider equalityΦ(𝜆) = 𝐻 𝑚 (𝜆) + 𝑅 𝑚 (𝜆),(25)

where 𝐻 𝑚 (𝜆)the function defined by equality (11), which at 𝑚 → ∞ converges to a quasi-polynomial Φ(𝜆), а 𝑅 𝑚 (𝜆) = Φ(𝜆) − 𝐻 𝑚 (𝜆).

The number 𝑚 0 should be chosen so that the inequality holds for 𝑚 ≥ 𝑚 0 min

𝜆 |Φ(𝜆)| > max 𝜆 |𝑅 𝑚 (𝜆)|.

Theorem 3 can be used to study the stability of solutions of linear differential-difference equations with a delay. Algorithms of such research are given in works [10,17].

Estimation of the effect of delay on the stability of solutions of linear equations with delay

In engineering practice, systems with a delay are replaced by systems without a delay, if it is small. Let us consider the method of finding the upper limit of the delay, in which the stability of systems with a delay is equivalent to the stability of the corresponding systems without a delay.

We will consider the mathematical justification of the possibility of replacing differential-difference equations with delay by ordinary differential equations is considered, and consider on finding upper bounds of delay, for which the stability regime of systems with delay is analogous to the stability regime of the corresponding systems without delay, is carried out.

Consider the linear system with a delay (7) and its corresponding system without a delay

𝑑𝑥(𝑡) 𝑑𝑡 = ∑ 𝐴 𝑖 𝑘 𝑖=0 𝑥(𝑡). (26)

Theorem 4 [18,19]. If the zero solution of the system (26) is asymptotically stable, then there exists a constant Δ > 0 such that for 0 < 𝜏 < Δ the zero solution of the system with a delay is also asymptotically stable.

Consider the algorithm for applying the approximation scheme of differential-difference equations (7) by the system of ordinary differential equations (9) to find the upper limit of the delay in system (7), in which it is asymptotically stable [20]:

1. We match the system with a delay (7) to the approximating system (9). 2. We reduce the characteristic polynomial of the approximating system (10) to a form that is convenient for finding its roots.

3. We set the initial limit delay Δ 0 and a step ℎ (example, ℎ = 1). 4. We find the approximate root of the quasi-polynomial 𝛼 with the largest real part 𝛼 at 𝜏 = Δ 0 . 5. If 𝛼 < 0, then for the given Δ 0 the zero solution of system with delay ( 7) is asymptotically stable. We put Δ 𝑖+1 = Δ 𝑖 + ℎ, 𝑖 = 0,1, …, and we go to point 4 until this property is fulfilled.

6. If received 𝛼 > 0, then for this upper bound Δ 𝑖+1 the zero solution of the delay system ( 7) is asymptotically unstable.

7. We obtained the interval [Δ 𝑖 , Δ 𝑖+1 ] of the values of the upper limit of the delay, at the ends of which the asymptotic stability changes. Using the method of division in half, we narrow this interval to the required accuracy and select the left end of the obtained interval for the required value. which is exponentially stable. Using the standard approximation scheme (17) and the increased accuracy approximation scheme (19) for the differential-difference equation (27), we will approximate equation (27) by quasipolynomial (16) with characteristic polynomials of the corresponding approximating systems of ordinary differential equations.

Table 2 shows the approximate values of the upper limits of the delay ∆, which are found by the algorithm proposed above for different values of the dimensions of the approximating systems m according to the standard approximation scheme and the approximation scheme of increased accuracy. The exact value of the quantity ∆ for which the exponential stability of equation ( 15) is preserved, found by the D-partition method [17] is equal to

𝜏 * = cos −1 (− 𝑎 𝑏 ) 𝑏 sin (cos −1 (− 𝑎 𝑏 )) .

In the case of equation ( 27), we have 𝜏 * = 2.4181. It follows from Table 2 that the application of approximation schemes of differential-difference equations makes it possible to find an approximate value of the upper limit of the delay for which exponential stability is preserved. At the same time, the scheme of increased accuracy provides better accuracy of approximation.

Conclusions

The aims to study the schemes of approximation of initial value problems for systems of differential equations with many delays by a sequence of systems of ordinary differential equations and their application to the study of the stability of systems of linear differential-difference equations with many delays and finding the upper bound of delay for which the stability of the system with delay is preserved.

The study of these problems is reduced to checking the conditions for the negativity of the real parts of all zeros of the corresponding quasipolynomials. Since a direct verification of this condition in practice is possible only in the simplest cases, we analyze the roots of characteristic polynomials of the corresponding approximating systems of ordinary differential equations to solve it.

The classical approximation scheme results in algorithms for finding the asymptotic roots of quasipolynomials that are convenient for computer implementation but require a high dimensionality of the corresponding approximating system. To improve the approximation accuracy of the nonasymptotic roots of quasipolynomials, a scheme for increased approximation accuracy is proposed and compared on test model examples.

The implementation of the proposed algorithms for studying the stability of solutions of linear differential equations is demonstrated on model test examples. The analysis of numerical experiments carried out for the test model examples confirms the theoretical results presented in the paper.