# An energy-based computational method in the analysis of the transmission of energy in a chain of coupled oscillators

## Abstract.

In this paper we study the phenomenon of nonlinear supratransmission in a semi-infinite discrete chain of coupled oscillators described by modified sine-Gordon equations with constant external and internal damping, and subject to harmonic external driving at the end. We develop a consistent and conditionally stable finite-difference scheme in order to analyze the effect of damping in the amount of energy injected in the chain of oscillators; numerical bifurcation analyses to determine the dependence of the amplitude at which supratransmission first occurs with respect to the frequency of the driving oscillator are carried out in order to show the consequences of damping on harmonic phonon quenching and the delay of appearance of critical amplitude.

###### Key words and phrases:

finite-difference scheme; bifurcation analysis; nonlinear supratransmission###### 2010 Mathematics Subject Classification:

(PACS) 45.10.-b; 05.45.-a; 02.30.Hq; 05.45.-a## 1. Introduction

The study of nonlinear continuous media described by modified Klein-Gordon equations subject to initial conditions is a topic of interest in several branches of the physical sciences [1, 2, 3, 4]. Here the analytical study on features of solutions of Klein-Gordon-like equations as well as the development of new and more powerful computational techniques to approximate them have been the most transited highways in the mathematical research of the area.

The behavior of nonlinear continuous media described by modified Klein-Gordon equations subject to boundary conditions is also an interesting topic of study in mathematical physics. This type of mathematical models have proved to be useful when applied, for instance, to the description of the project of data transmission in optical fibers in nonlinear Kerr regimes [5, 6] or in the study of the property of self-induced transparency of a two-level system submitted to a high-energy incident (resonant) laser pulse [7, 8]. More concretely, the behavior of a continuous medium submitted to a continuous wave radiation constitutes a fundamental problem that has not been studied in-depth, yet it has shown to have potential applications as a model in the study of Josephson transmission lines [9, 10, 11].

Numerical studies on discrete versions of these models [12] followed by experimental results [13] have shown that there exists a bifurcation of wave transmission within a forbidden band gap in certain semi-infinite undamped Klein-Gordon-like chains of coupled oscillators which are periodically forced at the end. This bifurcation is a consequence of the generation of nonlinear modes by the periodic forcing at the end of the chain, and allows energy to flow into the medium via a nonlinear process called *nonlinear supratransmission*, which has been proved numerically not to depend on integrability as long as the model possesses a natural forbidden band gap.

The process of nonlinear supratransmission has been studied numerically by means of computational algorithms already built in standard mathematical packages. Most particularly, the numerical results obtained in [12], for instance, rely on the use of a Runge-Kutta method of high order, which has the advantage of possessing high accuracy and efficiency in the computations, but lacks the consistency in the numerical estimation of the energy flowing into the medium which is desired in the study of supratransmission. With this shortcoming in mind, we present in this paper an alternate numerical formulation to study the process of supratransmission in a nonlinear system of differential equations, that has the advantage of being consistent in the approximation of the solutions to the problem and in the estimation of the continuous energy. More concretely, in the present work we study the process of nonlinear supratransmission in a semi-infinite nonlinear discrete system of coupled oscillators governed by modified sine-Gordon and Klein-Gordon equations with constant internal and external damping. We develop some numerical results for dissipative nonlinear modified Klein-Gordon-like equations, and validate our conclusions against [12]. Our study — rigorous in numerical nature — will soon be succeeded by future applications in forthcoming papers.

In Section II of this paper we introduce the mathematical problem under study and derive the expression of the instantaneous rate of change of the energy transmitted to the system through the boundary. Section III is devoted to introducing the finite-difference scheme; here we present the discrete energy analysis associated with the problem under study and establish in detail the numerical properties of our method. Numerical results are presented in the next section, followed by a brief discussion.

## 2. Analytical results

### 2.1. Mathematical model

In this article we study the effects of the nonnegative constant parameters and , and the real constant on the behavior of solutions to the mixed-value problem with mass ,

(1) |

which describes a system of coupled oscillators with coupling coefficient , and where and clearly play the roles of internal and external damping coefficients, respectively. The number is used to represent the spatial second-difference for every , and the boundary-driving function is assumed continuously differentiable on . In our study, we will let for a chain of coupled sine-Gordon equations, and for zero mass in the case of a Klein-Gordon chain.

It is worth recalling that the differential equation under study has multiple applications in the continuous limit. For instance, a similar equation appears in the study of Josephson junctions between superconductors when dissipative effects are taken into account [1]. A modification of this equation is used in the study of fluxons in Josephson transmission lines [2], and a version in the form of a modified Klein-Gordon equation appears in the statistical mechanics of nonlinear coherent structures such as solitary waves (see [3] pp. 298–309) when no internal damping is present. Besides, our equation clearly describes the motion of a damped string in a non-Hookean medium.

For purposes of this paper, we will consider a system of coupled oscillators initially at rest at the origin of the system of reference, subject to an external harmonic forcing described by , and pure-imaginary or real mass satisfying . Analysis of the undamped linearized system of differential equations in (1) shows that the linear dispersion relation reads in any case. The driving frequency will take on values in the forbidden band gap region , in which case the linear analysis yields the exact solutions , where

Meanwhile, the massless undamped nonlinear case possesses an exact solution in the continuous limit which happens to work well for high values of the coupling coefficient. It has been shown numerically that, for each frequency in the forbidden band gap, the massless medium starts to transmit energy by means of nonlinear mode generation for amplitudes greater than the threshold value .

### 2.2. Energy analysis

In this section we derive the expression of the instantaneous rate of change of total energy in system (1). Here, by a square-summable sequence we understand a real sequence for which is convergent.

###### Lemma 1 (Discrete Green’s First Identity).

If is a square-summable sequence then

###### Proof.

Hölder’s inequality implies that both series in the formula of the lemma converge. Moreover, the sequence defined by for every positive integer converges to zero and the result follows now from the identities

∎

###### Proposition 2.

Let be solutions to (1) such that is square-summable at any fixed time . The instantaneous rate of change of the total energy in the system is given by

###### Proof.

The energy density of the undamped system of coupled equations with a potential energy in the -th oscillator is given by . After including the potential energy from the coupling between the first two oscillators, the total energy of the system at any time becomes

and the fact that tends to as increases implies that the sequence converges to zero pointwisely at any fixed time. Simplifying and rearranging terms of the derivative of the Hamiltonian yields the expression

The result follows now from this identity after computing the derivative of the total energy of the system, using the formula for telescoping series and applying our discrete version of Green’s First Identity. ∎

Observe from the proposition that the expression of the rate of change of energy associated with damped system (1) is in agreement with the undamped formula derived in [12]. Moreover, it is clear that the external damping coefficient contributes decreasing the total amount of energy in the chain system for equal to zero.

## 3. Numerical analysis

### 3.1. Finite-difference scheme

From a practical point of view, we will restrict our study to a system consisting of a finite number of coupled oscillators with constant internal and external damping coefficients, described by the system of ordinary differential equations

for , where includes both the effect of external damping and a simulation of an absorbing boundary slowly increasing in magnitude on the last oscillators. More concretely, we let be equal to zero at all time , and let be the sum of external damping and the function

In practice, we will let , , and .

We proceed now to discretize problem (1) using a finite system of differential equations and a regular partition of the time interval with time step equal to . For each , let us represent the approximate solution to our problem on the -th oscillator at time by . If we convey that , that and that , our problem takes then the discrete form

(2) |

Observe that the proposed numerical method is nonlinear and requires an application of Newton’s method for systems of equations in order to be implemented. Notice also that if at the -th time step is approximated by then the finite-difference scheme becomes linear and an application of Crout’s reduction technique for tridiagonal systems suffices to approximate solutions of (1). In such case, it is readily seen that the vector equation must be satisfied for every , for the three -dimensional tridiagonal matrices and the -dimensional vector

respectively, and constants

Here for every , and is the -th dimensional vector whose -th component is equal to . This latter formulation of our problem will be used for validation purposes only, since we will prefer the nonlinear formulation due to the quadratic order of convergence of Newton’s method and other reasons that will be presented in the next section.

### 3.2. Discrete energy

The total energy in the system at the -th time step will be approximated numerically via

###### Proposition 3.

The following identity holds for every sequence satisfying scheme (2) and every positive index :

###### Proof.

The proof of this result is merely an algebraic task. We need only say that an application of the discrete Green’s First Identity with is indispensable in order to reach the correct expression of the term with coefficient . ∎

### 3.3. Stability analysis

The following result summarizes the main numerical properties of our method.

###### Proposition 4.

###### Proof.

Following the notation in [14], let and for each and . For every and let be the column vector whose components are and . Our problem can be written then in matrix form as

where

Applying Fourier transform to the vector equation we obtain

The matrix multiplying in this equation is the amplification matrix of the problem. It is easy to check that the eigenvalues of for are given by

Suppose for a moment that . If the radical in the expression for the eigenvalues of is a pure real number then . So for every , grows faster than for any constants and . A similar situation happens when the radical is a pure imaginary number, except that in this case represents the usual Euclidean norm in the field of complex numbers. Therefore in order for our numeric method to be stable order it is necessary that , which is what we wished to establish. ∎

## 4. Numerical results

In this section we study the effects of the internal and external damping coefficients on the behavior of solutions to mixed-value problem (1). Particularly, we wish to establish the effect of weak damping on the minimum amplitude value necessary for the phenomenon of supratransmission to take place at a fixed driving frequency. Throughout we validate our code against [12] and against an implementation of our finite-difference scheme using the Runge-Kutta method of order four.

### 4.1. External Damping

For the remainder of this work we consider a semi-infinite coupled chain of oscillators initially at rest in their equilibrium positions, subject to harmonic forcing at the end described by at any time . The functions and are both identically equal to zero and, in order to avoid the shock wave produced by the vanishing initial velocity in the boundary, we let the driving amplitude linearly increase from to in a finite period of time before we start to compute the total energy. In the present section we will let be equal to zero and consider a discrete system of coupled oscillators described by (1) over a typical time period of , with a time step of and a coupling coefficient equal to .

Let us first examine the massless case when no damping is present at all and . In order to verify that our mixed-value problem produces a bifurcation it is necessary to study the qualitative behavior of solutions near the predicted threshold value which, in this case, reads approximately . The first row of Figure 1 shows the function in the solution of a sine-Gordon system for two different values of the amplitude of the driving end, while the second row presents the corresponding graphs in the solution of a Klein-Gordon chain. The graphs evidence the presence of a bifurcation in the behavior of solutions around the amplitude value for both chains.

Naturally, the next step in our investigation will be to determine the behavior of the total energy flow injected by the periodic forcing at the end of the undamped discrete chain of oscillators as a function of the amplitude. The solid lines of Figure 2 represent the graphs of total energy transmitted into the system vs. amplitude for a forcing frequency of and external damping coefficient equal to zero, for a sine-Gordon system in the first column and a Klein-Gordon system in the second. It is worthwhile noticing the abrupt increase in the total energy administered to the system for some amplitude value between and . This feature of the graphs evidences the existence of a bifurcation value around the predicted amplitude , after which the phenomenon of nonlinear supratransmission takes place.

Figure 2 also presents graphs of total energy vs. forcing amplitude for weak constant damping coefficients , and in a sine-Gordon system of oscillators. The graphs show a tendency of the bifurcation value to increase linearly as the external damping coefficient is increased. Another interesting feature of this figure is the decrease of total energy for increasing values of , at least for fixed amplitudes greater than the undamped bifurcation threshold. Needless to point out that similar conclusions are obtained for Klein-Gordon chains of oscillators.

The qualitative effect of in a sine-Gordon system is also of interest in the analysis of solutions of this chain and is numerically carried out in Figure 3 for a chain with external damping equal to and pure-imaginary and real masses, using graphs of total energy administered into the system through the boundary vs. driving amplitude. A horizontal shift in the occurrence of the bifurcation value is readily noticed in these graphs. Indeed, the displacement of the bifurcation amplitude for the system (1) of mass with respect to the bifurcation amplitude of the massless system is a monotone increasing function of . Analogous computational results (not included here) were obtained for a similar Klein-Gordon chain.

### 4.2. Internal damping

Consider again a system of oscillators ruled by mixed-value problem (1) over a time period of , with time step , coupling coefficient equal to and constant external damping equal to zero. In this context, Figure 4 shows graphs of total energy vs. forcing amplitude for internal damping coefficients , , and , for sine-Gordon and Klein-Gordon systems. As in the case of external damping, we observe that the threshold value at which supratransmission starts tends to increase as the value of is increased. Opposite to the case of external damping, though, the minimum value for which supratransmission starts varies with in a nonlinear way.

Next we verify our conclusions on the effect of the mass on the qualitative behavior of solutions of the sine-Gordon system. Figure 5 shows graphs of total energy vs. driving amplitude for a system with no external damping, and driving frequency of , for different pure-imaginary and real masses. As observed before, the graphs evidence a shift on the bifurcation amplitude with respect to the massless bifurcation value which is an increasing function of .

### 4.3. Bifurcation analysis

Both the chain of coupled sine-Gordon equations and the chain of Klein-Gordon equations have numerically proved to undergo nonlinear supratransmission for a driving frequency equal to and different values of the external and internal damping coefficients, thus establishing that the results obtained in this work do not depend on integrability. Naturally we are interested in determining that the process of nonlinear supratransmission happens for any frequency value in the forbidden band gap. With that purpose in mind, we obtained graphs of total energy vs. driving frequency and various amplitude values for undamped discrete systems of coupled oscillators with coupling coefficient , time period of and time step . The -dimensional results for both chains of oscillators are shown in Figure 6 together with a graph of the continuous-limit threshold amplitude vs. driving frequency on the amplitude-frequency plane for comparison purposes.

Several observations may be immediately drawn from Figure 6. First of all, the predicted bifurcation values display an excellent concordance to the corresponding values obtained using the finite difference-schemes associated with the sine-Gordon and the Klein-Gordon chains. Second, the process of supratransmission ceases to appear in the Klein-Gordon case for driving frequencies below . Third, for driving frequencies close to the band gap limit, the bifurcation threshold is not clearly determined from the energy vs. driving amplitude graph of the sine-Gordon chain, as prescribed by [12]. For those frequencies, it is indispensable to increase the time period of approximation at least up to . Fourth, strong numerical proof of the existence of the occurrence of the supratransmission process is at hand and we have now reasons to believe that there exists a bifurcation function for driving amplitude associated with (1).

We proceed then to obtain graphs of amplitude values for which nonlinear supratrasmission starts vs. driving frequency for a massless sine-Gordon chain of coupled oscillators with internal damping coefficient equal to zero and different values of . The numerical results are summarized in Figure 7 together with the plot of the prescribed continuous-limit bifurcation amplitudes . It is worth noticing that the bifurcation threshold increases with for fixed frequencies above , as previously evidenced for a frequency equal to . We must notice also that the discrepancies that appear for frequencies below when , also appear for greater values of , each time bounded in smaller intervals. In fact, we have checked that the discrepancy region — an effect of the phenomenon of harmonic phonon quenching — tends to vanish for higher values of external damping (results not included). In this state of matters, we wish to point out that better numerical approximations to the bifurcation threshold are obtained for larger systems of oscillators at the expense of superior needs in terms of computational time. Likewise, we have confirmed that larger values of the coupling coefficient lead to better approximations to the bifurcation threshold, as prescribed by [12].

Figure 8 shows bifurcation diagrams for a massless sine-Gordon system of oscillators with no external damping and different constant internal damping coefficients. As expected, the bifurcation threshold tends to increase with for a fixed driving frequency. The effects of harmonic phonon quenching are present again in all the bifurcation diagrams, but contrary to the case of external damping, in the case of internal damping the range over which discrepancies occur slightly widens as increases. Also, it is worth observing that the length of the forbidden band gap increases with the parameter apparently in a linear fashion. Moreover, the graphs of bifurcation diagrams for nonzero values of are approximately obtained by shifting horizontally the corresponding graph of the undamped diagram a number of units to the right. More concretely, is a good approximation for .

Finally, we compute bifurcation diagrams for an undamped sine-Gordon system with different real and pure-imaginary masses in order to establish the effect of on the occurrence of the bifurcation threshold. The numerical results are summarized in Figure 9 for some real and pure-imaginary masses, and driving frequencies starting at . Our results lead us to conclude that the bifurcation diagram of the sine-Gordon system of oscillators with mass is approximately equal to the graph of the corresponding massless graph shifted horizontal units, for . In order to test our claim numerically, we obtained graphs of absolute differences between the massless undamped bifurcation diagram, and shifted undamped bifurcation diagrams for several mass values (results not included). We observed that the differences tend to attenuate for high frequency values and that smaller differences are obtained for smaller values of in magnitude.

## 5. Conclusions

In this paper we have developed a numerical method to approximate solutions of a semi-infinite nonlinear chain of coupled oscillators ruled by modified sine-Gordon equations harmonically driven at its end. The proposed finite-difference scheme is consistent order and we provided a necessary condition in order for the method to be stable order . The process of nonlinear supratransmission for a coupled system of oscillators described by sine-Gordon equations was studied numerically under the scope of this numerical technique, and the dependence of supratransmission on damping was analyzed.

Several conclusions can be drawn from our computational experiments on the sine-Gordon system of coupled oscillators (1). First of all, we have shown that the phenomenon of harmonic phonon quenching still appears in the presence of external and internal damping and that the discrepancy region due to phonon quenching is shortened as the external damping coefficient is increased, while it slightly widens as the internal damping coefficient increases. Second, the threshold value at which supratransmission first occurs for fixed frequencies outside the discrepancy region is seen to increase for both external and internal damping as the damping coefficient increases; both conclusions are clearly consequences of the dispersive and dissipative natures of the parameters and , respectively.

Finally, the bifurcation diagram for a value of the parameter is approximately equal to the corresponding diagram for the undamped system shifted horizontal units to the right. Likewise, a horizontal shift of units in the bifurcation diagram of a sine-Gordon system of mass with respect to the corresponding massless system is observed for small masses and frequencies outside the discrepancy region.

#### Acknowledgment

It is our duty to express our deepest gratitude to Dr. Álvarez Rodríguez, dean of the Centro de Ciencias Básicas of the Universidad Autónoma de Aguascalientes, and to Dr. Avelar González, head of the Dirección General de Investigación y Posgrado of the same university, for uninterestedly providing us with computational resources to produce this article. The present work represents a set of partial results under project PIM07-2 at this university.

### References

- Remoissenet, M. Waves Called Solitons. Springer-Verlag, New York, 3rd edition, 1999.
- Lomdahl, P. S.; Soerensen, O. H.; Christiansen, P. L. Soliton excitations in Josephson tunnel junctions. Phys. Rev. B, 25(9):5737–5748, 1982.
- Makhankov, V. G.; Bishop, A. R.; Holm, D. D., editor. Nonlinear Evolution Equations and Dynamical Systems Needs ’94; Los Alamos, NM, USA, 11-18 September ’94: 10th International Workshop. World Scientific Pub. Co. Inc., Singapore, first edition, 1995.
- Barone, A.; Esposito, F.; Magee, C. J.; Scott, A. C. Theory and applications of the sine-Gordon equation. Riv. Nuovo Cim., 1:227–267, 1971.
- Hasegagwa, A.; Tappert, F. Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. I. Anomalous dispersion. Appl. Phys. Lett., 23:142–144, 1973.
- Mollenauer, L. F.; Stolen, R. H.; Gordon, J. P. Experimental observation of picosecond pulse narrowing and solitons in optical fibers. Phys. Rev. Lett., 45:1095–1098, 1980.
- McCall, S. L.; Hahn, E. L. Self-induced transparency. Phys. Rev., 183:457–485, 1969.
- Ablowitz, M. J.; Kaup, D. J.; Newell, A. C. Coherent pulse propagation, a dispersive, irreversible phenomenon. J. Math. Phys., 15:1852–1858, 1974.
- Ustinov, A. V. Solitons in Josephson junctions. Physica D, 123:315–329, 1998.
- Ustinov, A. V.; Mygind, J.; Oboznov, V. A. Phase-locked flux-flow Josephson oscillator. J. Appl. Phys., 72:1203–1205, 1992.
- Ustinov, A. V.; Mygind, J.; Pedersen, N. F.; Oboznov, V. A. Millimeter-wave-induced fluxon pair creation in flux-flow Josephson oscillators. Phys. Rev. B, 46:578–580, 1992.
- Geniet, F.; Leon, J. Energy transmission in the forbidden band gap of a nonlinear chain. Phys. Rev. Lett., 89:134102, 2002.
- Geniet, F.; Leon, J. Nonlinear supratransmission. J. Phys.: Condens. Matter, 15:2933–2949, 2003.
- Thomas, J. W. Numerical Partial Differential Equations. Springer-Verlag, New York, 1995.