On the Phenomenon of Low-Frequency, Large-Amplitude Oscillations in a High-Dimensional Linear Dynamical System Текст научной статьи по специальности «Математика»

Russian Journal of Nonlinear Dynamics, 2024, vol. 20, no. 2, pp. 259-276. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd240203

MATHEMATICAL PROBLEMS OF NONLINEARITY

MSC 2010: 39A21

On the Phenomenon of Low-Frequency, Large-Amplitude Oscillations in a High-Dimensional

Linear Dynamical System

A. I. Gudimenko, A. V. Lihosherstov

This article considers a linear dynamic system that models a chain of coupled harmonic oscillators, under special boundary conditions that ensure a balanced energy flow from one end of the chain to the other. The energy conductivity of the chain is controlled by the parameter a of the system.

In a numerical experiment on this system, with a large number of oscillators and at certain values of a, the phenomenon of low-frequency high-amplitude oscillations was discovered. The primary analysis showed that this phenomenon has much in common with self-oscillations in nonlinear systems. In both cases, periodic motion is created and maintained by an internal energy source that does not have the corresponding periodicity. In addition, the amplitude of the oscillations significantly exceeds the initial state amplitude. However, this phenomenon also has a fundamental difference from self-oscillations in that it is controlled by the oscillation synchronization mechanism in linear systems and not by the exponential instability suppression mechanism in nonlinear systems.

This article provides an explanation of the observed phenomenon on the basis of a complete analytical solution of the system. The solution is constructed in a standard way by reducing the dynamic problem to the problem of eigenvalues and eigenvectors for the system matrix. When solving, we use methods from the theory of orthogonal polynomials. In addition, we discuss two physical interpretations of the system. The connection between these interpretations and the system is established through the Schrödinger variables.

Keywords: linear dynamical system, harmonic chain, high-amplitude oscillations

Received August 21, 2023 Accepted January 23, 2024

Aleksey I. Gudimenko gudimenko@iam.dvo.ru

Institute for Applied Mathematics, Far Eastern Branch, Russian Academy of Sciences ul. Radio 7, Vladivostok, 690041 Russia

Andrey V. Lihosherstov likhosherstov.02@mail.ru

Far Eastern Federal University Universitetskiy pr., Vladivostok, 690091 Russia

1. Introduction

This article is mainly concerned with the analytical explanation of the low-frequency large-amplitude oscillations observed in numerical experiments on a multidimensional linear dynamical system of the form

xi = xi+i - xi-i, l = 1, ..., N, (11) ax0 + x1 = 0, axN + xN+1 =0, a g R, a < 0.

A review of the literature shows that bounded oscillations whose amplitude significantly exceeds the amplitude of the initial state have not been noticed before in linear systems. However, there are a number of aspects of these oscillations that make them interesting to study.

First, because the system (1.1) is autonomous, the large-amplitude oscillations cannot be associated with an external or parametric resonance [1, 2]. Another explanation of this phenomenon is required.

Second, as will be shown below, the system (1.1) is an active oscillatory system, i.e., a system in which oscillations are supported by an internal source of energy [3]. Together with autonomy, this property brings the system closer to the class of self-oscillatory systems, which are additionally characterized by the presence of limit cycles in their phase spaces [4, 5]. A limit cycle is a closed trajectory in the phase space that attracts or repels all other trajectories from some of its neighborhood. It is well known that there are no limit cycles in linear systems. Therefore, the observed phenomenon cannot be explained within the framework of the classical theory of self-oscillations. It should be noted here that the concept of self-oscillations is widely used in the modern science. For example, we refer the reader to [6-11], and for low-frequency self-oscillations, to [12].

Third, it is believed that the behavior of a linear dynamical system is well studied and is completely determined by its Jordan matrix [1]. In other words, the dynamics of the system is entirely determined by the structure of the invariant (stable, unstable, and central) subspaces of the system's phase space [13, 14]. This is, of course, true of the basis of (generalized) eigenvectors of the system. However, the transition to the Jordan basis in multidimensional systems can be very complicated and, as a consequence, the system can exhibit unexpected effects in the original variables, such as the one discussed in this article.

Our interest in the dynamical system (1.1) is also due to its physical interpretations. The transition to them is carried out through the Schrodinger transformation of variables [15-17].

Schrodinger variables. Consider a system of second-order differential equations

Vi - Vi+i +2yi - Vi-i = 0, l = 1,...,L, (1.2)

with boundary conditions

ayo + Vi - Vo = 0 yL+i + a(VL+1 - Vl) = 0 (1-3)

or

aV0 + Vi - Vo = 0, aVL + Vl+i - Vl = (l4)

We call the variables Vi natural variables because they are directly related to the physical variables considered below.

The Schrodinger variables are introduced by the expressions

x2i = yl, x2i+i = Vi+i - Vi. (L5) RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2024, 20(2), 259 276_

In these variables, as a simple check shows, Eq. (1.2) with the boundary conditions (1.3) or (1.4) turns into the system (1.1), in which N = 2L + 1 for the former condition and N = 2L for the latter. Thus, in the Schrödinger variables, two different boundary value problems, (1.2), (1.3) and (1.2), (1.4), correspond to one dynamical system, (1.1), but with N of different parity.

Mechanical interpretation. Boundary value problems (1.2), (1.3) and (1.2), (1.4) describe the oscillations in a finite chain of material points of mass m, connected by springs of stiffness k, with a damper at the left end of the chain and an antidamper at the right end (see Fig. 1).

m

m

A qt (h A ql+1 ql+1

(a)

7 k

%

(b)

m

aVl+i Ql+i AqL qL

-7

■3-

Ql

(c)

(d)

Fig. 1. Chain of coupled oscillators as a dynamical model for the problems (1.2), (1.3) and (1.2), (1.4): (a) internal chain segment, (b) left end of the chain, (c) and (d) variants of the right end of the chain. Here y is damping, ql is the displacement of the Zth oscillator from the equilibrium position, and Aql =

= qi- q-i

Indeed, applying Newton's laws, it is easy to show that the dynamics of the chain segments (a)-(d) of Fig. 1 are described by the equations

mi qi - kAql+1 + kAql = 0, Z = 1,...,L, Yqo - kAqi = 0 mqL+i - YAqL+i = 0 (YqL - k^qL+i = 0).

By normalizing the time variable to y^r and setting a = we get (1.2)-(1.4). For example,

to obtain the first equality of the second line in (1.6), we must equalize two forces acting on the damper piston: the elastic force of the spring, kAq1, and the viscous force of the damper, Yq0. However, to obtain the last equation of this system, in addition to establishing a similar equality of force,

kAqL + mLqL - YqL = 0,

we have to use the first equation of (1.6) at Z = L.

The mechanical interpretation is the simplest and most obvious. However, it is not easy to imagine real mechanisms for the antidampers (c) and (d). In addition, in mechanics, the quantity a is used rarely, if at all. In electrical interpretation, these drawbacks are less pronounced.

Electrical interpretation. The problems (1.2), (1.3) and (1.2), (1.4) describe also the oscillations in the chain of electric oscillatory circuits shown in Fig. 2 by its typical and boundary fragments.

-R

Fig. 2. Chain of oscillatory circuit: (a) a chain segment, (b) the left end of the circuit, (c), (d) variants of the right end of the circuit. In this figure, C is the capacitance, L is the inductance, R is the resistance, and il is the electric current in the Zth circuit segment

Applying Kirchhoff's laws or using the electromechanical analogy [18], we find for the chain fragments (a)-(d) the equations

ÄC^-t1+to = 0>

A dt2 di

1

LC'

(ii_ 1 - 2ii + h+i) = o

L+1

dt

~ j{'L+1 - h)

(1.7)

Normalizing the time to y/LC, we get (1.2)—(1.4), where a = In the theory of electric

circuits, the quantity and its inverse are known as the quality factors of a parallel and

a series RLC circuit, respectively [19].

The principal feature of the boundary circuits (c) and (d) is the presence of a negative resistor. Such an electrical element models an internal energy source (an electric generator), in our case, with a linear current-voltage characteristic [20]. The boundary circuit (b) contains a positive resistor. On this element, energy dissipation occurs.

In the theory of electric circuits, the system (1.7) describes a lossless transmission line with

a load on the left end of the line and a generator on the right end [21-23]. The quantity is called the characteristic impedance of the line. The transmission line is characterized by the that is, by the value of a. If |a| ^ 0, the intensity of the energy transfer in

ratio of R and

the line tends to zero, if |a| increases, this intensity grows.

Physics of the system. It is convenient to unify the physical interpretations given above using the concept of a chain of coupled harmonic oscillators (briefly, a harmonic chain) as a sequence of quantities yl, l = 1, ..., L, whose dynamics obey Eq. (1.2). Thus, using this terminology, we say that the system (1.1) describes oscillations in a harmonic chain in Schrödinger variables.

To better understand the physical background of the system (1.1), we introduce the energy

N 2

function E = J2 '2~- This is just the energy of a harmonic chain expressed in the Schrödinger i=i

o

variables. For example, in the variables of the mechanical interpretation, we have

L

e = £

1=1

™iQi HQi+I ~ Vif

2

2

Let us show that the left and right boundary conditions in (1.1) determine the energy dissipation and pumping, respectively. We have

N

N

e — xixi — xi(xi+i

— x

l-1 )

I _ -1 2 2

xixo + xn Xn+1 — a xi — ax N,

(1.8)

i=1

i=1

where (1.1) is used to obtain the second and fourth equality. We see that the energy dynamics of the system is determined only by the first and last chain oscillators. Since the first term of the resulting expression in (1.8) is always negative and the second is positive, we conclude that, for the harmonic chain defined by the system (1.1), the energy is dissipated at the left end and generated at the right end of the chain. In other words, there is an energy flow in the chain in the direction from right to left.

The specific behavior of the energy of the system is determined by the solution of (1.1). For N = 3, this solution is easy to find explicitly by calculating the matrix exponent:

eAt —

x— r)Atry0 ty - f ty ty ty 1 t 123 0 1 x0 x0 x2, x3 ]T ,

1 a a2 ' a4 — a2 —a —a2

—a —a2 —a3 e(a 1-a)t + a a4 + 1 a3 (cos t + A sin t)

a2 a3 a4 —a2 —a3 1 — a2

(1.9)

Then

E —

eatx

o

eAtx0)

eAT t eAt x0 ™0

e e x , x

(1.10)

2 2 2 From (1.9) and (1.10), it follows that, in general, the energy of the system undergoes quasi-periodic oscillations, decreasing at —1 < a < 0 and increasing without bound when a < -1.

This behavior of energy is consistent with that of the phase volume Q = dx1 A ■■■ A dxN of the system (1.1). It is well known that the evolution of the phase volume in a linear system is determined by the trace of the system matrix. In more detail, for the Lie derivative of Q along the vector field of the system, we have

Lvfi — d(u J fi) — d

N

J2(xi+1 — xi-1)di j (dx1 A ••• A dxN)

i=1 1

= d\a 1x1 dx2 A ■■■ A dxN — (—1)N axN dxl A - ■ ■ A dxN __ J = (a 1 — a) Ü,

where J is the contraction of forms [24]. We see that at —1 < a < 0 the phase volume decreases, which corresponds to the prevalence of energy dissipation over energy pumping. At a = —1, the balance of these processes is observed. At a < —1, the phase volume increases, and therefore, the energy generation exceeds the dissipation.

Thus, as a mechanical system, the system under consideration is a nonconservative system, in contrast, for example, to the systems considered in the classical works of Gantmacher and

Krein [25, 26]. However, in the limit a —0, our system becomes conservative, since in this case, the boundary conditions in (1.1) are replaced by x1 = xN+1 = 0, and, as in (1.8), we find E = —x2x1 + xNxN+1 = 0. These boundary conditions correspond to harmonic chains with a free left end and a free (N is even) or fixed (N is odd) right end.

Purpose and plan of the article. The main purpose of this study is to provide an analytical explanation of the large-amplitude low-frequency oscillations observed in the system (1.1) at a close to —1 and with a special choice of initial conditions. To explain this effect, it is not sufficient to know the individual behavior of the system's eigenmodes. This explanation can only be obtained by analyzing the join behavior of the system's eigenmodes. This is what the main part of our work is devoted to.

The plan of the article is as follows. In Section 2, we specify the problem statement, in Section 3, we analytically solve the system (1.1), in Section 4, we study the behavior of phase curves of the system over the whole range of a < 0, in Section 5, we explain the phenomenon of generation of high-amplitude oscillations, and in Section 6, we draw conclusions.

2. Problem statement

This study considers a linear dynamical system defined by the differential-difference equation

xi = xi+i — xt-i> l = 1, ■ ■ ■,

and boundary conditions

ax0 + x1 = 0, axN + xN+1 = 0, a E R, a < 0■

(2.1)

(2.2)

In matrix form, the system is

where

x1

xN -1

xN

ry> - A rp

«aj --y

(a-1

A =

1 0

-1 0 1

V

-1 0 1 0 -1 -a

Physical interpretations of the system (2.1)-(2.2) are described in the Introduction. Summing up, the system describes oscillations in a harmonic chain under special boundary conditions that ensure energy flow from the right side of the chain to the left side. The parameter a plays a double role. On the one hand, it coordinates the energy loss at the left end of the chain and its generation at the right end. The energy balance is achieved at a = —1. On the other hand, it acts as a relative boundary conductivity, i. e., the ratio of the chain's boundary conductivity to its specific conductivity. For example, the minimum conductivity is realized at a = —0. Energy in such chains is conserved, i.e., there is no absorption or generation of energy at the chain boundaries.

We are interested in the dynamics of the system for different values of a and different initial conditions. We pay special attention to the phenomenon of low-frequency large-amplitude oscillations in the system. This phenomenon is observed when a is close to or equal to —1 and

N

x

2

x

the initial data is rapidly changing in the spatial variable l, i.e., when the spatial frequency of the change in the initial data is sufficiently close to or coincides with the maximum frequency of their change.

The characteristic behavior of oscillations is illustrated by the plots of the numerical solution of the system (2.1)-(2.2) shown in Fig. 3.

1

0.5 0

-0.5 -1

u

n

u

n

x¡ +99

lN

u

100 50 0

-50 -100

n

N = 99,a = -0.25, = cos(ttZ) N = 99, a = -1, xf = cos ($) N = 99, a = -1, œ? = cos(ttZ)

(a)

(b)

(c)

300

H

200 100 0

10 15

n

AT =100, a =-0.997, x? = cos($) N = 100, a = -0.997, z? = cos(7rI) TV = 100, a = -l, a;? = cos(7rZ)

(d)

(e)

(f)

0 2 4 6

-2000 -4000 -6000 -8000

30000

H

20000 10000

2 4 6

N = 99, a = -1.003, x= cos(ttZ) N = 100, a = -1.003, xf = cos(ttZ)

(g)

(h)

Fig. 3. Plots of x15 (t) at different parity of N, values of a < 0, and initial conditions. The amplitude of the initial state is chosen to be the same and equal to 1

We see that for a ^ —1 close to or equal to —1 and initial data changing with the maximum spatial frequency, xl = cos(nZ) = (—1)1, low-frequency high-amplitude oscillations are observed. The amplitude of these oscillations is N times greater than the amplitude of the initial data (see Figs. 3c and 3e). Moreover, if N is even, the shift of oscillations relative to zero can be made arbitrarily large by varying a (see Fig. 3e). When a = —1 and even N are chosen, a linear growth of oscillations is observed, as shown in Fig. 3f. When a < —1, the amplitude of oscillations grows exponentially (see Figs. 3g and 3h). For other values of a < 0 and initial data, the amplitude of oscillations is comparable with the amplitude of spatial change of initial data (Figs. 3a, 3b, and 3d).

In Fig. 3, the function xl(t) is plotted for l = 15. The behavior of oscillations for other values of l is similar (see Fig. 4).

20 40 60 80

20 40 60 80

-200

200J

N = 99, a = -1, xJ> = cos(trl) N = 99, a = -1, xf = cos(ttZ)

Fig. 4. Plots of xt(t) as a function of l for different values of t. The values of t are taken from the time interval T equal to the oscillation period of the function plotted in Fig. 3c

The purpose of this article is to provide an analytical explanation of the oscillation dynamics presented in Figs. 3 and 4. To achieve this, we

1) find in Section 3 an analytical solution to the initial value problem for the system (2.1)-(2.2) with the initial conditions

xi(0)= xl, l = 1,...,N, (2.3)

2) study the behavior of the phase trajectories of the system for the entire range of values of a < 0 and initial data (Section 4),

3) give an analysis of the found solution in the vicinity of a = —1 (Section 5).

3. Solving the dynamical system

We solve the initial boundary value problem (2.1)-(2.3) by reducing it to the eigenvalue and eigenvector problem for the system matrix A [1]. In view of the special (tridiagonal) form of A, the latter problem can be formulated as a boundary value problem for the recurrence relation

= Vi+i — Vi-1, l = 1, ...,N — 1, (3.1)

with boundary conditions

aVo + Vi = 0 aVN + Vn+1 = (3.2)

The recurrence relations of the form (3.1) are studied in the theory of orthogonal polynomials [27-29]. To calculate the eigenvalues and eigenvectors of A, we use the technique adopted in this theory.

We begin with the eigenvalues. Setting

n(2k — N)

uk = 2sin<

™ 2N k rk' (3.3)

r0 = —a, rk = ei(pk, k = 1, ..., N — 1,

we have

Proposition 1. The eigenvalues of A are

\ -1-1

A0 = r0 — rn = a — a, 0 0 0 ' (3.4)

Ak = rk — rk = lwk, k = 1,...,N — 1

Thus, there exists only one eigenvalue, A0, that depends on a. This eigenvalue is real and the other eigenvalues are imaginary. If a = —1 or N is odd, all the eigenvalues are distinct. Otherwise, that is, if a = —1 and N is even, there exists a double eigenvalue A0 = An/2 = 0.

Proof.

Following the approach adopted in the theory of orthogonal polynomials, to find the eigenvalues, we must solve the initial value problem for (3.1) with the initial conditions

a^0 + Vi = 0, Vi = 1 (3.5)

and then impose on the solution the second of conditions (3.2).

Indeed, let Al be the principal leading 1th order submatrix of A and I be the identity matrix of an appropriate order. Expanding the determinant of Al — AI along the last row, we see that the sequence of polynomials in A

Pi = 1, Pi+i = (—1)1 det(Ai — AI), l = 1, ...,N, (3.6)

is a unique solution of the initial problem (3.1), (3.5). On the other hand, from the second of conditions (3.2) and Eq. (3.1) for l = N — 1, it follows that

(a + X)Vn + Vn -i =

When applied to (3.6), this means that pN+i = (—1)N det(AN — AI) = 0, i. e., A is an eigenvalue of the matrix A.

We solve the initial value problem (3.1), (3.5) by the characteristic equation method [30]. For the recurrence relation (3.1), the characteristic equation is

A = r — r-i. (3.7)

This equation has two roots. If the roots are distinct and r is one of them, according to Vieta's formulas, the other is —rki. In this case, the solution of the problem (3.1), (3.5) is represented in the form

Pi = cirl + C2 ( 1)1 rkl. (3.8)

Substituting (3.8) into (3.5), we find ci, c2 and obtain

(rki — a) rl + (a + r)(—1)rkl 0

» = --V + r-) (3-9)

For multiple roots, we look for a solution in the form pl = (ci + c2l)rl. Proceeding similarly to the case of distinct roots, we obtain

Pl={a + r)l~rr\ r2 + 1 = 0. (3.10)

l ar

Now substituting (3.9) and (3.10) into the second of conditions (3.2), we arrive at the system of equations for r:

(r + a)(ar - 1) [r2N - (-1)N] =0, r2 + 1=0, (3.11)

(a ± i)2 = 0, r2 + 1 = 0.

The second equation in this system is unsolvable for a g R. Therefore, the eigenvalues of the matrix A are determined by the first equation. Solving (3.11) and substituting the found roots into (3.7), we get (3.4). □

Remark 1. It is worth noting that according to Proposition 1, the characteristic polynomial of A is represented in the form

Pn+1 = (a - a0W, where the polynomial pN has exactly imaginary roots.

Now, we describe the eigenvector of A. For each k = 0, ..., N - 1, we denote by Yk an eigenvector (perhaps generalized) corresponding to the eigenvalue \k, and by Yk 1, ..., Yk n, the entries of Yk.

Proposition 2.

1. If a = —1 or N is odd, the vectors y0, ■ ■ ■ , Yn- 1 form a complete set of eigenvectors, and their entries can be taken as

Yk,i = (r-1 - a) rk + (rk + a)(-1)1 r-1, l = 1, ..., N. (3.12)

In particular,

Yo,i = - {a~l + a) (-a)l, Yn/2,i = 2. (3.13)

If a = -1 and N is even, to the zero eigenvalue there corresponds one eigenvector Yn/2 and one generalized eigenvector y0■ In this case, the vectors y0, ■■■, Yn- 1 form a complete set of eigenvectors and generalized eigenvectors■ If y1, ■■■, Yn- 1 are of the form (3.12), then

1 ( 1)i

7o,z = l~2+-^L, 1 = 1, ...,N. (3.14)

2. The eigenvectors are orthogonal in the following sense: If a = -1 or N is odd, then

N

^(-1)lYk,iYk',i = Pkh,k', k,kk = 0,...,N - 1, (3.15)

\=\

where

Po = (1 + a2) [(-1)wa2N - 1, pk = 2N(r-1 - a) (rk + a), k = 0.

If a = —1 and N is even, equality (3.15) holds also, with the exception that

N N

e(-1)' i = N2, Y. (-1)1 Yo, iYn/2 , i = 2N. (3.16)

i=i i=i

Proof.

1. If a = —1 or N is odd, all the eigenvalues are distinct, and therefore, the corresponding eigenvectors form a complete set. The expression (3.12) follows from (3.9) and the fact that the sequence pl(Ak), l = 1, ..., N, is a solution of the problem (3.1)-(3.2) for A = Xk. We, however, give a direct proof of (3.12).

We seek the sequence Yk i, ..., Yk N as a solution of the problem (3.1)-(3.2). Since the characteristic equation of relation (3.1) has no multiple roots (see the proof of Proposition 1), we have

Yk,i = Ck, irl + (—1) Ck , 2r-1, (3.17)

where ck 1, ck 2 are arbitrary complex constants. Substituting (3.17) into (3.2) we find

ckii(a + r) + ck,2 (a — r-1) = 0, Ck ,irN(a + r) + ck ,2(—1)Nr-N (a — r-1) = 0.

Solving this system using (3.11), we find ck 1, ck 2 and obtain (3.12).

If a = —1 and N is even, it follows from Proposition 1 that A0 = AN/2 = 0, i.e., the zero eigenvalue is a multiple eigenvalue. Let YN/2 be an eigenvector. A direct check shows that y0 defined by the expression (3.14) is a generalized eigenvector:

Yo,i+i — Yo,i-i = yn/2,^ l = 1,...,N, Yo , o = Yo i Yo, N+i = Yo , N. Note that we can find (3.14) by the formula

Yo

lim

YjV/2 " Yo

AN/2 -A0 a>-1

Indeed, then

AYN/2 — AYo .. AN/2YN/2 — Ao Yo

= llI111 V-\- = llm1 -\-\- = llm1 70 = 7W/2>

a^-i an/2 — Ao a^-i an/2 — Ao a^-i 1

i.e., Yo is a generalized eigenvector attached to the eigenvector yn/2.

2. To obtain equalities (3.15) and (3.16), one only needs to calculate the sums they contain.

□

Now, we establish the main result of this section.

Theorem 1. If a = —1 or N is odd, the solution of the problem (2.1)-(2.3) is

N-i 1

xi = X] —

k=0 pk

If a = —1 and N is even, then N

N

E(—1)^ x° Ykl

.l'=1

(3.18)

— y (-if

2N v '

x=

l' 0

2N

1) x

i'=i

N

Yn/2j' (Yo,/ + tlN/2,i) + ( Yo,/' ~ yYw/2,/' ) In/2,1

+

Ni

k=l Pk

k=N/2

N

E(—1)" x0' Ykii

i'=i

YkieXkt. (3.19)

a

Proof.

For a = —1 or odd N, we have

N-1 N-1 _ N-1

77 ^ht-ÏT _ ~ I \ ^ ^ ^

•'•/ = a-oeAoi7o,/ + I] akeXktlkj + J] H^lki = ^e^So,/ + J] afceAfei7fc,/+ k=1 k=1 k=1 N—1 N—1 N—1 N—1

+ J] = aoeAoi7o,/ + J] afceAfei7fc,/ + J] âW-fceAfci7W = J] VV '"/,./•

k=1 k=1 k=1 k=0

Here, the second equality is obtained by applying the identity 7fc/ = 7which is easily verified using (3.3) and (3.12). In the third equality, we have changed the summation index in the second sum. In the last equality, we have relabeled the indeterminate coefficients. So, we have

N-1

xi bkeXk*Yk, 1. (3.20)

k=0

Note that this expression also follows from the diagonalizability of matrix A.

The case when a = —1 and N is even is similarly considered. As a result, we have

N—1

xl = b0(Y0,l + tYN/2,l) + bN/2,l Yn/2,1 + ^ bkeh * Yk,l• (3-21)

k=1 k=N/2

Applying the orthogonality relations (3.15) and (3.16) (taken at t = 0) to equalities (3.20) and (3.21), we obtain (3.18) and (3.19), respectively. □

4. Invariant subspaces

In a linear dynamical system, the structure of phase space (central, slow, stable, and unstable subspaces) is determined by the eigenvalues of the system matrix [13, 14]). For our system, as a direct consequence of Proposition 1, we have

Theorem 2. If a = —1, the phase space decomposes into a direct sum of the central subspace spanned by the eigenvectors Y1, ■■■, Yn- 1 and stable (—1 < a < 0) or unstable (a < —1) subspace spanned by y0 ■ In this case, if N is even, the central subspace includes the slow subspace spanned

by yn/2 ■

If a = —1, the entire phase space is central and contains the slow subspace spanned by Y0 for odd N, and y0, Yn/2 for even N■

The following theorem describes the central subspaces of the system and indicates the subspaces on which the dynamics is bounded.

Theorem 3. If a = —1, the central subspace is given by the equation

N

E(—1)1 xiY0l = 0. (4.1)

1=1

For odd N, there are no slow subspaces■ For even N, the slow subspace is generated by the vector Yn/2 ■ The dynamics on the central subspace is bounded■

If a = —1 and N is odd, the phase space is fibered into the central subspaces

N

E(—1)1 XYo ,i = C c e R, (4-2)

i=i

along the slow subspace spanned by the eigenvector y0 . In this case, the dynamics is bounded on the entire phase space.

If a = —1 and N is even, the fibering (4.2) also takes place, but along the generalized eigenvector y0 . The subspaces of bounded dynamics are given by the equations

NN

E(—1)1 xiYo,i = E(—1)1 xiYn/2,i = 0, C e R. (4.3)

i=i i=i

Proof.

All the assertions are a consequence of the orthogonality relations (3.15), (3.16). If a = —1, the central subspace is generated by the vectors y1 , ..., YN- 1 and therefore, according to (3.15), is characterized by the orthogonality of these vectors to the eigenvector y0. Applying (3.15), we get (4.1). The boundedness of the dynamics follows, for example, from Theorem 1.

When a = —1 and N is odd, the entire phase space is central, and therefore, any vector is decomposed into y0, ..., Yn-r Applying (3.15) to this decomposition, we get (4.2). The case when a = —1 and N is even is similarly considered. Conditions (4.3) follow from the explicit form of the solution (3.19). □

Theorems 2 and 3 give a general representation of the dynamics of the system at different values of a < 0 and initial data. The dynamics is bounded in the following cases:

(i) —1 < a < 0, or a = —1 and N is odd,

N

(ii) a < —1 and ^ (—1)1 xly0 l = 0,

i=1 i 0 i

N

— m In

i=1

(iii) a = —1, N is even, and ^ (—1)'%iYn/2 i = 0

In all other cases, the dynamics is unbounded. Exponential growth of oscillation occurs at a < — 1, and linear growth is observed at a = —1 and even N.

After substituting the expressions (3.13) for y0 l and YN/21, the sums from (ii) and (iii) are transformed into a single sum

N

a1 = 0. (4.4)

i=i

It follows from (3.18) and (3.19) that (4.4) is a condition for suppressing both the exponential growth (decay) and the linear growth of oscillations. Since (4.4) expresses the orthogonality of y0 and yn/2 to the vectors y1, ..., YN-1, the vector x = (x1, ..., xN) can be considered as a linear combination of the vectors Re(Yk) and Im(Yfc), k = 1, ..., N — 1. Using (3.12) and (3.3), we

Re(Yfc ,i) = [1 + (—1)il cos[0fc(l — 1)] — [1 — (—1)il acos[0fcl],

r S1 r S1 (4.5)

Im(Yfc>i)= 1 — (—1)i sin[0fc(l — 1)] — 1 + (—1)i asin[0fcl], l = 1,...,N. The suppression of the linear and exponential growth of oscillations is illustrated by Fig. 5.

■ n

AT = 100, a = —l

(a)

0 -1

t_ IN

AT = 99, a = -1.003

(b)

Fig. 5. Plots of xig(t) at x0 — Re(Y49 ^). These plots correspond to Figs. 3f and 3g

5. Analysis of the solution near a = — 1

To analytically study the dynamics of the system in the vicinity of a = —1, it is more convenient to use a different form of solution than that in Theorem 1.

Theorem 4. The solution of the problem (2.1)-(2.3) can be represented as follows: for a = —1 or odd N,

ii —2

(0-i-0)t l + a

N

(—a)1 al x° +

(—1)N a2N — 1

i'=1 2N

H--Y^ e""k

+ 2N ^

k=1

0U4k ^ e^fc + a e—i^k- a

e»vfc + —--(_i ye-»vk

N

¿2xï,e—*+-, (5.1)

i'=1

and for a — —1 and even N, N

¿D-d

xl —

i' ^ o

N

1) x

i'

i'=1

l + l' + 2t — N — 1 +

(—1)i

+

+

2N

2N

£

k=1

k=3N/2

eiwk *

eii^fc +

e^k — 1

e-^k + 1

(—1)l e—ii^k

N

J>0 e—ii. (5.2)

i'=1

Proof.

To obtain (5.1), we transform the sum over k in (3.18) as

N1

N1

e k Pk Yk,iYk,i' = e k

k=1

k=1

("D V ( 'i + '>V I +

Tk — a

1

+rk

.i' 'k

Tu - a

rk + a

Tk + (—1) r

i —i

k

N —1

2N-1

k=1

i t-at u a

rk — a

+ Z^ e r2W-fc | I-—r2N_k + (-1) r2W_fc ] =

k=N +1

'2N—k + a 2N-1

2N — 1 /

E e^—f r— U +

k=1

,tf 1 I J , rjL+a ! i\lm-l r—l- a

k '

k

(—1)ir

k

2

1

Here, the first equality has been obtained by substituting (3.12), in the second equality we have split the sum into two and changed the summation index. In the last equality, we applied the easily verifiable identity

-1

r2N-k = —rk ,

k = 1,

N - 1.

The expression (5.2) was obtained from (3.19) similarly.

Corollary 1. If a = —1 and N is even, the solution has the form

(5.3)

□

_ (a 1-a)t_

xi = e

1 + a

-2

N

(—1)N a2N — 1

1

(—'a)' E alx° +

i'=i

1

2N

(—1)' +

1—a

1 + a

N

E(—1)" x?+

l' = 1

2N

+ 2N ^

k=1 k=3N/2

eil^k +

e^fc + a e-l^h — a

(—1)1 e-il^k

N

Ex?' e-il'^. (5.4)

i'=i

Proof.

This follows from (5.1) after substitutions $3N/2 = n and w3N/2 = 0. □

The analytical solution was tested for consistency with the numerical solution over a wide range of values of the time coordinate, parameter a, and initial conditions. The results always coincided within the predetermined accuracy.

Further, we focus on the analysis of the steady-state bounded oscillations in the system. According to (5.1) and (5.2), the amplitude of these oscillations is determined by the joint behavior of the quantities

fi^k + a

DkJ(a) = + _ + (-l)'e-^ e i(pk — a

N

xk ^E

(5.5)

l=i

From the first equality in (5.5), we see that the main contribution to the magnitude of Dk l(a) comes from a and k close to —1 and ^j-, respectively. Variation of I has no significant effect on Dkl. If N is odd, the magnitude of Dkl(a) reaches its maximum at k = ™ ± ^ and, as simple calculations show, the estimate is valid:

lA-,/(«)l L.-3JV+I

'h — 2 =c2

2 /J_

|l + a| + V^2

7T \N

If N is even, we have

2

\Dk,l(a)\\

k=21±ztl

\1 + a\ 2N

+ O

1

+ O

n

if a is sufficiently far from —1,

if a = 1.

if a is sufficiently far from —1,

if a = 1.

(5.6)

(5.7)

The second quantity in (5.5), Xk, is the kth component of the shifted discrete Fourier transform for the sequence of 2N numbers x1, ..., xN, 0, ..., 0:

Xi, —

2N

N

l=i

This transform characterizes the spatial frequency spectrum of the initial state of the system.

1

The behavior of the quantities Dk l(a) and \Xk\ as functions of k at large N and a close to —1 is illustrated by Figs. 6a-6c.

\Dk,i\

n

0.5

0.5 1 1.5 2n

a = -1, 1 = 25

(a)

I* J

n

0.75

I

0.50 0.25 0

0.5 1 1.5 2n

xï=cos{<p3N/2l)

(b)

I* J

n

0.75

I

0.50 0.25 0

0.5 1 1.5 2^

a% = coB(<l>If/2l)

(c)

Fig. 6. The plots of

and ^ as functions of k at N = 99

Note that, the values of k close to correspond to the low-frequency oscillations of the system with eigenfrequencies close to the minimum eigenfrequency. Indeed, according to (5.3) and (3.4), the identity wk = w2N-k is valid, which implies that the frequency wN/2 = 0 corresponds to k = At". For odd N, for example, the minimum eigenfrequency satisfies the estimate

W(N±1)/2 = 2sin 0(N±1)/2 = ±2sin

2N=±N+°

1

?P

The shift of the steady-state bounded oscillations relative to zero, S, can be estimated from the first summands in (5.1) and (5.2), and the second summand in (5.4). We have

N

B-1)

-1)' X0

S=

1=1

N

ivB-1« '=1

1

if a = — 1 and N is odd,

if a = — 1 and N is even,

(5.8)

1

N

1

1 + a

+ O(1)

N

E(-1)'x0 if a = -1 and N i

is even.

'=1

The estimates (5.6)-(5.8) agree in order of magnitude with the numerical results. For example, for N = 99, a = -0.25, and x0 = (—1)', the amplitude of oscillations is estimated as œ 1.3, whereas at a = -1, we have œ 63 (see Fig. 3a,c). In the last case, the shift of oscillation is -N. For N = 100, a = -0.997, and x0 = (-1)', the shift is estimated as œ 333 (see Fig. 3e). To obtain more accurate estimates for the amplitudes of oscillations, one needs to consider other terms of the sum over k in (5.1).

t

6. Conclusion

In this article, we have investigated a linear dynamical system that describes oscillations in a chain of harmonic oscillators under special boundary conditions that ensure a coordinated flow of energy from the right end of the chain to the left. The character of the flow is determined by

the parameter a < 0 and initial conditions xO, l — 1, ..., N. Based on numerical and analytical studies, we draw the following conclusions about the dynamics of the system.

There are bounded and unbounded oscillations in the system. Bounded oscillations occur when at least one of the conditions is met:

(i) —1 < a < 0, or a — —1 and N is odd,

The unbounded ones take place when the corresponding conditions are

(iii) a < —1 and (ii) is false,

(iv) a = —1 and N is even, and (ii) is false.

If the value of a is sufficiently far from —1, the amplitude of bounded oscillations is comparable to the initial one. If a is close to —1, the low-frequency high-amplitude oscillations arise. Their amplitude can be N times greater than the initial amplitude and is proportional to the maximum frequency in the spatial Fourier expansion of the initial state. For even N, the shift of the oscillations relative to zero can be made arbitrarily large by varying a. If condition (iii) or (iv) is valid, the amplitude of oscillations grows exponentially or linearly in time, respectively. Condition (ii) suppresses this growth as well as the exponential decay of oscillations. In this sense, the system becomes conservative, and the problem arises of representing it in Hamiltonian form.

As to the spectrum of the system matrix, there is only one real eigenvalue that can take nonzero values; all other eigenvalues are imaginary. If a = —1 or N is odd, all the eigenvalues are distinct. If a = —1 and N is even, there exists a double, zero, eigenvalue. A numerical experiment, the result of which we do not present in this article, shows that the spectrum is sensitive to a small mismatch of boundary conditions. For example, adding a small negative number to the rightmost entry of the matrix results in a spectrum in which all eigenvalues have nonzero negative real parts. In this case, the system becomes purely dissipative. A detailed study of this issue can be the subject of a subsequent study.

Conflict of interest

The authors declare that they have no conflict of interest.

References

[1] Müller, P. C. and Schiehlen, W. O., Linear Vibrations: A Theoretical Treatment of Multi-Degree-of-Freedom Vibrating Systems, Monographs and Textbooks on Mechanics of Solids and Fluids: Mechanics. Dynamical Systems, vol. 7, The Hague: Martinus Nijhoff Publ., 1985.

[2] Gatti, P. L., Applied Structural and Mechanical Vibrations: Theory and Methods, 2nd ed., London: CRC, 2014.

[3] Pikovsky, A., Rosenblum, M., and Kurths, J., Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge Nonlinear Sci. Ser., vol. 12, New York: Cambridge Univ. Press, 2001.

[4] Andronov, A. A., Vitt, A. A., and Khaikin, S.E., Theory of Oscillators, Oxford: Pergamon, 1966.

[5] Jenkins, A., Self-Oscillation, Phys. Rep., 2013, vol. 525, no. 2, pp. 167-222.

N

(ii) ¿ alx0 = 0.

i=i

[6] Monkewitz, P. A., Bechert, D. W., Barsikow, B., and Lehmann, B., Self-Excited Oscillations and Mixing in a Heated Round Jet, J. Fluid Mech, 1990, vol. 213, no. 1, pp. 611-639.

[7] Zheng, Z. and Qian, Y., Dominant Phase-Advanced Driving Analysis of Self-Sustained Oscillations in Biological Networks, Chin. Phys. B, 2018, vol. 27, no. 1, 018901, 10 pp.

[8] Li, M. and Bao, J., Effect of Self-Oscillation on Escape Dynamics of Classical and Quantum Open Systems, Entropy, 2020, vol. 22, no. 8, Paper No. 839, 10 pp.

[9] Arnaut, L.G. and Ibanez, S., Self-Sustained Oscillations and Global Climate Changes, Sci. Rep., 2020, vol. 10, no. 1, 11200, 19 pp.

[10] Zhang, H. and Detournay, E., An Alternative Formulation for Modeling Self-Excited Oscillations of Rotary Drilling Systems, J. Sound Vibration, 2020, vol. 474, 115241, 22 pp.

[11] Zhai, C. and Sun, W., Analytical Approximation of a Self-Oscillatory Reaction System Using the Laplace-Borel Transform, Chaos Solitons Fractals, 2021, vol. 142, Paper No. 110508, 12 pp.

[12] Balmashnov, A.A., Kalashnikov, A.V., Kalashnikov, V.V., Stepina, S. P., and Umnov, A.M., Self-Excitation of Low-Frequency Oscillations in the Plasma Ring Formed by an ECR Discharge in a Narrow Coaxial Cavity, Plasma Phys. Rep., 2018, vol. 44, no. 6, pp. 626-629.

[13] Chicone, C., Ordinary Differential Equations with Applications, 2nd ed., Texts Appl. Math., vol. 34, New York: Springer, 2006.

[14] Teschl, G., Ordinary Differential Equations and Dynamical Systems, Grad. Stud. Math., vol. 140, Providence: AMS, 2011.

[15] Schrödinger, E., Zur Dynamik elastisch gekoppelter Punktsysteme, Ann. Phys., 1914, vol. 44, pp. 916-934.

[16] Takizawa, É. I. and Kobayasi, K., Heat Flow in a System of Coupled Harmonic Oscillators, Chinese J. Phys., 1963, vol. 1, no. 2, pp. 59-73.

[17] Kobayasi, K. and Takizawa, É. I., Effect of a Light Isotopic Impurity on the Energy Flow in a System of One-Dimensional Coupled Harmonic Oscillators, Chinese J. Phys., 1964, vol. 2, no. 2, pp. 68-79.

[18] Olson, H.F., Dynamical Analogies, New York: D. Van Nostrand, 1943.

[19] Nilsson, J. W. and Riedel, S.A., Electric Circuits, 12th ed., New York: Pearson, 2023.

[20] Groszkowski, J., Frequency of Self-Oscillation, Oxford: Pergamon, 2013.

[21] Pozar, D.M., Microwave Engineering, New York: Wiley, 2011.

[22] Ulaby, F. T. and Ravaioli, U., Fundamentals of Applied Electromagnetics, New York: Pearson, 2023.

[23] Irwin, J. D. and Nelms, R. M., Basic Engineering Circuit Analysis, 12th ed., New York: Wiley, 2020.

[24] Lee, J.M., Manifolds and Differential Geometry, Grad. Stud. in Math., vol. 107, Providence, R.I.: AMS, 2009.

[25] Gantmacher, F.R., The Theory of Matrices: Vol. 1, New York: Chelsea, 1998.

[26] Gantmacher, F. P. and Krein, M.G., Oscillation Matrices and Kernels and Small Vibrations of Mechanical Systems, Providence, R.I.: AMS, 2002.

[27] Atkinson, F. V., Discrete and Continuous Boundary Problems, Math. Sci. Eng., vol. 8, New York: Acad. Press, 1964.

[28] Guseinov, G.Sh., On an Inverse Problem for Two Spectra of Finite Jacobi Matrices, Appl. Math. Comput., 2012, vol. 218, no. 14, pp. 7573-7589.

[29] Akhiezer, N.I., The Classical Moment Problem and Some Related Questions in Analysis, Edinburgh: Oliver & Boyd, 1965.

[30] Elaydi, S., An Introduction to Difference Equations, New York: Springer, 2005.