On the simulation of the energy transmission in the forbidden band-gap of a spatially discrete double sine-Gordon system

ON THE SIMULA TION OF THE ENERGY TRANSMISSION IN THE F ORBIDDEN BAND-GAP OF A SP A TIALL Y DISCRETE DOUBLE SINE-GORDON SYSTEM J. E. MA C ´ IAS-D ´ IAZ Abstract. In this work, w e presen t a n umerical method to consisten tly ap- prox imate solutions of a spatially discrete, double sine-Gordon c hain which considers the presence o f external damping. In addition to t he finite-difference sc heme employ ed to approximate t he solution of the difference-different ial equations of the mo del under inv estigation, our method provides positivity- preserving sc hemes to appro ximate the lo cal and t he total energy of the system, in suc h wa y tha t the discr ete rate of c hange of the total ene rgy with resp ect to time provides a consistent appro ximation of the corresp onding contin uous rate of change . Simulation s are per f ormed, first of all , to assess the v alidity of the computational tec hnique against kno wn qualitative solutions of coupled sine-Gordon and coupled double sine-Gordon chains. Secondly , the method is used in the inv estigation of the phenomenon of nonlinear transmission of energy i n double si ne-Gordon systems; the qualitative effects of the damping coefficient on the occurrence of the nonlinear pro cess of supratransmiss i on are briefly determined in this work, too. 1. Introduction The well-kno wn sine-Go rdon equa tion is a partial differential equation that ap- pea rs in ma ny applications , either in its or iginal for m o r as a slight mo dification of the classical version. F or instance, a damp ed s ine-Gordon equation app ears in the study of long Josephso n junctions b et ween sup erconductors when dissipative effects a re taken into a ccoun t [1 ]. A s imilar par tial differential equation with dif- ferent no nlinea r term app ears in the study of flux o ns in Josephson transmission lines [2]. Meanwhile, a mo dified Klein-Gor don equation app ears in the sta tistical mechanics of nonlinear co heren t structures —such as s olitary wa ves—, in the form of a L a ngevin equation (see [3], pp. 298–3 09). The spatially dis c rete version o f the sine-Gordon equation also has ma n y imp o r- tant applications. F o r instance, a coupled system of discrete s ine-Gordon equa tions may describ e a chain of harmonic oscillato rs coupled thro ugh Hookean springs [4] or a s ystem of Jose phs on junctions a ttac hed through sup erconducting wir es [5]. In the former case, a system initially at rest, with void initia l velo cities and sinusoidal Dirichlet b oundary conditio n at one end, is employ ed to s tudy the phenomeno n of supratransmiss io n of e nergy [4, 6], which is a nonlinear pro cess characterized by a sudden increa se in the amplitude of w av e signals ge ne r ated by the p erturbe d bo undary , for driving a mplitudes above a cr itical v alue c alled the supr atr ansm ission thr eshold . This phenomenon is als o present in the inv estigation of the tra nsmission 2010 Mathematics Subje ct Classific ation. (P ACS) 02.60.Lj; 02.70.Bf; 45.10.-b. Key wor ds and phr ases. double sine-Gordon ch ain; difference-differen tial equation ; finite- difference sc heme; energy sc heme; nonlinear supratransmission. 1 2 J. E. MAC ´ IAS-D ´ IAZ of ener gy in chains of Josephs o n junctions, except that, in this case, a harmo nic Neumann bo undary co nditio n needs to b e imp osed up on the physical pro blem for the sa k e of meaning fulness [5]. W e must rema rk, in this p oin t, that several other bo unded, nonlinear regimes present the phenomenon of supratr ansmission when they are harmo nically p er- turb ed at one end, such is the case o f Klein-Gordon arr a ys [4], F ermi-Pasta-Ula m systems [7], Bragg media in the no nlinear Kerr regime [8], and even in spatially con- tin uous, bounded media describ ed by unda mped sine-Gordon equatio ns [9]. Mo re- ov er, the presence of the phenomenon of nonlinear supr a transmission is also found in discrete, double s ine - Gordon chains [6]. How ever, the sp ecialized mathematical literature unfortunately lacks studies to approximate the o ccurr ence of the pro cess in these systems. Nevertheless, in this w or k, we study the pro cess of nonlinear supratr ansmission in dissipative, double sine-Gordo n chains, emplo ying a dissipation-pr eserving finite- difference scheme. Of course, there exist many ana lytical [10, 11] and numerical [12, 13, 14] techniques to appr oximate solutions of Klein-Gordo n-lik e equations. The c o mputational metho d e mplo yed in this work distinguishes from many other techn iques av aila ble in the literature in that it co ns isten tly approximates not only the solution of the physical mo del under study , but also the loca l e ne r gy densit y , the total energ y , a nd the rate of change o f the ener gy with r espect to time. Moreover, the nonnega tiv e character of bo th the ener gy density and the total energy —a c har- acteristic which is not preser v ed by other numerical techniques [1 5]—, is pres erv ed by o ur metho d. As we shall see later o n, these qualities are, by definition, highly desirable characteristics of any computational technique employ ed in the study of the phenomeno n o f supratransmiss io n. In Section 2, we present the system of or dinary differen tial equa tio ns that moti- v a tes our study , together with the lo cal ener gy functions ass ociated, and the total energy of the sy s tem. A prop osition which summarizes the expre s sion of the de- riv ativ e of the total ene r gy with r espect to time is provided in this stage, a s well as a br ief description of the pr ocess of no nlinear supratransmiss ion, particularly , in double sine-Gordo n sy stems. Section 3 intro duces the numerical metho d employ ed to approximate the solutions o f the mo del under inv estigation, the lo cal energ y distribution, and the total e ner gy of the system. A subsection o n the numerical prop erties o f the metho d summarizes the pr operties of c o nsistency established for conv enience in the app endices, and ano ther s ubsection presents some remar ks on the co mputational implementation of our technique. Section 4 presents simulations of sine-Go rdon and double sine- Gordon systems, obtained by means of our metho d. The for mer regime is employ ed only for v alidation purp oses, while the s im ulations on the latter (which follow the same metho dology pro p osed in [4]) ar e aimed at establishing the existence of the pro cess o f nonlinear supratra nsmission in this sys- tem. Finally , this work clo s es with a section of co ncluding remark s and further directions o f re s earch. 2. Preliminaries 2.1. Ph ysical mo del. T hr oughout this work, we let Z N = { 1 , 2 , . . . , N − 1 } , for every p ositiv e integer N ; ob vious ly , we will a ssume that N > 1 for the sake o f non-triviality . Moreover, w e le t Z N = Z N ∪ { 0 , N } , that is, Z N = { 0 , 1 , 2 , . . . , N } . ENERGY SIMULA TION IN DOUBLE SINE -GORDON 3 Let N be a p ositive int ege r a nd, for every n = 0 , 1 , . . . , N , le t u n be a rea l function o f time t ≥ 0. Moreover, let c b e a p ositive real num b er, and let γ be a nonnegative n umber. Thro ughout, we consider a mechanical chain of no nlinear os- cillators o beying the sys tem of o rdinary differential equations with initial-b oundary conditions (1) d 2 u n dt 2 − δ (2) x u n + γ du n dt + V ′ ( u n ) = 0 , ∀ n ∈ Z N ,          u n (0) = 0 , ∀ n ∈ Z N , du n dt (0) = 0 , ∀ n ∈ Z N , u 0 ( t ) = φ ( t ) , ∀ t ≥ 0 , u N ( t ) − u N − 1 ( t ) = 0 , ∀ t ≥ 0 . In o ther words, we consider a spatially discr ete, b ounded system initially at r est, with zero initial veloc ities , per turbed at the left end by a function φ which we will assume to be co n tinu ous, and with dis crete Neumann bo undary condition on the rig h t end. The constant γ is immediately identified a s the external damping co efficien t, while c is clearly the coupling co efficient b et ween no des. Here, the spatial, se c o nd-difference o p erator (2) δ (2) x u n = c 2 ( u n +1 − 2 u n + u n − 1 ) , ∀ n ∈ Z N , has b een employed for co n venience. F o r the sa k e o f concreteness, we will co nsider a dr iving function o f the for m (3) φ ( t ) = A sin(Ω t ) , where the dr iving amplitude A and the driving frequenc y Ω are p ositive rea l n um- ber s. Mor eo ver, we co nsider a p otential function o f the form (4) V ( u ) = 1 2 − 1 6 [2 cos u + cos(2 u )] , whence the double sine-Go rdon law V ′ ( u ) = 1 3 [sin u + sin(2 u )] readily results. In- deed, let c = 1 ∆ x . If ∆ x is r elativ ely small (or, equiv alently , c is relatively lar ge), then the system of ordinar y differen tial equations of (1) approximates the spatially contin uous, partial differ en tial equatio n (5) ∂ 2 v ∂ t 2 − ∂ 2 v ∂ x 2 + γ ∂ v ∂ t + V ′ ( v ) = 0 , x ∈ [0 , L ] , where L = N ∆ x , a nd v is a function of space x and time t . This equation is clearly ident ified with the classica l, double sine-Gordo n equation with c o nstan t external damping. Of course, different p oten tials may giv e rise to other imp ortant models in mathe- matical ph ysics . F or instance , V ( u ) = 1 − cos( u ) is the potential for the sine-Gordon regime, while V ( u ) = 1 2! u 2 − 1 4! u 4 + 1 6! u 6 corres p onds to the p oten tial of a nonlin- ear Klein-Gor do n equation. In fact, it is imp ortant to warn the reader that the dissipation-pres erving nu merica l technique presented in this work is v a lid no t only for the double sine-Gor don p otential, but a lso for a n y differentiable function V defined o n the rea l line. 4 J. E. MAC ´ IAS-D ´ IAZ ✻ t ✲ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ n n − 1 n + 1 t k t k +1 t k − 1 ✉ ✉ ✉ ✉ × × × Figure 1. F o rw ard- difference stencil fo r the approximation to the partial differential equa tion (1) a t time t k , using the finite- difference scheme (18). The black circles r epresent known approx- imations to the actual solutions at times t k − 1 and t k , and the crosses deno te the unk nown approximations at time t k +1 . 2.2. Energy of the system. Let n ∈ Z N . F or physical rea sons, it is imp ortant to notice that the lo cal energy of the n th node in the undamped system governed by (1) is provided by the expressio n (6) H n = 1 2 "  du n dt  2 + ( δ x u n ) 2 # + V ( u n ) , ∀ n ∈ Z N , where the spatial, fir st-order differenc e op erator δ x is defined through (7) δ x u n = c ( u n +1 − u n ) , ∀ n ∈ Z N . In these ter ms, the total energy E of the sy stem (1) is o btained by adding the lo cal energies H n , for n ∈ Z N , and the p oten tial due to the coupling in the bo undaries o f the chain. In other words, (8) E = X n ∈ Z N H n + 1 2 ( δ x u 0 ) 2 . The following prop osition is easy to establish. Prop osition 1. The ra te of change of ener gy with r esp e ct to time of a system governe d by (1) is given by (9) dE dt = − c ( δ x u 0 ) du 0 dt − γ X n ∈ Z N  du n dt  2 .  As a coro llary , the system (1 ) co ns erv es the total energ y if no exter nal damping is pr e sen t a nd, either φ is a consta n t function or a void Neumann c o ndition is impo sed o n the left end of the chain. ENERGY SIMULA TION IN DOUBLE SINE -GORDON 5 0 50 100 150 200 −5 0 5 t u 60 (t) 0 100 200 0 5 10 t H 60 (t) 0 50 100 150 200 −5 0 5 t u 60 (t) 0 100 200 0 5 10 t H 60 (t) Figure 2. Gra phs of approximate solution u 60 versus time, of the 60th no de of a sy stem (1) of length N = 2 00, for c = 4, γ = 0 and a po ten tial V ( u ) = 1 − cos u . The system was pe r turbed by means of (3) with Ω = 0 . 9, and t wo different amplitudes w ere used: A = 1 . 7 7 (top graph) a nd A = 1 . 7 8 (right column). The insets depict the corres p onding temp oral e volution of the lo cal energy of the 60 th no de. 2.3. Nonlinear supratransmis sion. As observed in the litera tur e (see [4, 6 , 7]), the double sine-Gordon sys tem (1), a s well as the nonlinea r Klein-Gor don and the sine-Gordon chains, and the class ical β -F er mi-P asta- Ulam systems, presents the phenomenon of s upr atransmission of energy , which is a nonlinear pr oces s charac- terized by a sudden incr ease in the a mplitude of wa ve s ig nals pro pagated into a nonlinear medium by a dr iving source which ir radiates at a fre q uency in the for - bidden ba nd-gap. The mec hanism of this transmission of energy is through the generation of lo calized, nonlinear mo des at the driving b oundary , in the fo r m of moving breather s or soliton so lutions [1 6]. More concretely , co nsider a nonlinear system of any of the types men tioned in the previous par agraph, whic h is p erturbe d at o ne end by a ha rmonic function of the form (3), with Ω a fixed v alue in the for bidden ba nd-gap o f the system. Relatively small driving amplitudes A result in the propaga tion of practically no energy into the sy stem; how ever, as the v alue of A is increa sed, the exis tence of a critical v alue A s , ab ov e which the s y stem b egins to absorb g reat amounts o f energ y from the bo undary , is immediately noticed. The v alue A s int ro duced in the previo us para graph, is called the supr atr ansmis- sion thr eshold , and its ex istence ha s b een analy tica lly pr o ved for discrete [4] a nd contin uous [9] sine-Gordon chains, as well as for systems of anha rmonic o s cillators [7]. Howev er, as it was men tioned b efore, the study of the double sine-Gor don regime has been left practically unexplor ed. F o r our particular study , a simple a nalysis of the undamp ed, linearized system of differential equations in (1) shows that the linear dispe rsion relatio n is given by (10) ω 2 ( k ) = 1 + 2 c 2 (1 − cos k ) . In our simulations, the driving frequency Ω will take on v alues in the forbidden band-gap reg ion Ω < 1. 6 J. E. MAC ´ IAS-D ´ IAZ 0 50 100 150 200 −3 −2 −1 0 1 2 3 t u 60 (t) 0 100 200 0 1 2 t H 60 (t) 0 50 100 150 200 −3 −2 −1 0 1 2 3 t u 60 (t) 0 100 200 0 1 2 t H 60 (t) Figure 3. Gra phs of approximate solution u 60 versus time, of the 60th no de of a sy stem (1) of length N = 2 00, for c = 4, γ = 0 and a po ten tial given b y (4). The sys tem w as p erturbed by means o f (3) with Ω = 0 . 9 , and tw o different amplitudes were used: A = 1 . 03 (top graph) a nd A = 1 . 0 4 (right column). The insets depict the corres p onding temp oral e volution of the lo cal energy of the 60 th no de. 3. Numerical method 3.1. Finite-difference sc heme. Le t N b e a p ositive in teger , and let T be a p osi- tive real n umber. In order to approximate the solutions o f the system (1) at time T , we fix a re gular partition o f the int erv al [0 , T ] of the form 0 = t 0 < t 1 < . . . < t M = T , of no r m ∆ t = T / M . Additionally , we let u k n be the numerical approximation of the actual v alue of u n at time t k , for k = 0 , 1 , . . . , M . More over, in or der to simplify o ur no tation, we define the temp oral differ ences δ t u k n = u k +1 n − u k n ∆ t , (11) δ (1) t u k n = u k +1 n − u k − 1 n 2∆ t , (12) δ (2) t u k n = u k +1 n − 2 u k n + u k − 1 n (∆ t ) 2 , (13) for every n ∈ Z N and k ∈ Z M . F urther more, for s uc h v alues of n and k , we employ the temp oral average op erator (14) µ (1) t u k n = 1 2  u k +1 n + u k − 1 n  , and the discrete deriv ativ e o f V with r espect to u and the time av era g e of V at u , resp ectiv ely: δ (1) u V ( u k n ) = V ( u k +1 n ) − V ( u k − 1 n ) u k +1 n − u k − 1 n , (15) µ t V ( u k n ) = V ( u k +1 n ) + V ( u k n ) 2 . (16) Finally , le t (17) φ k = φ ( t k ) . ENERGY SIMULA TION IN DOUBLE SINE -GORDON 7 0 50 100 150 200 0 100 200 −4 −2 0 2 4 n t u n (t) 0 50 100 150 200 0 100 200 −4 −2 0 2 4 n t u n (t) Figure 4. Graphs of appr o ximate solution versus no de site n and time t , of a sy s tem (1) of length N = 200 , for c = 4, γ = 0 and a po ten tial given b y (4). The system was p erturb ed by means of (3), with Ω = 0 . 9, and tw o different amplitude v a lues: A = 1 . 03 (left) a nd A = 1 . 0 4 (rig h t). With these con ven tions at hand, the numerical method to approximate s o lutions of (1) is summarized a s follows: (18)  δ (2) t − µ (1) t δ (2) x + γ δ (1) t + δ (1) u V  ( u k n ) = 0 , ∀ n ∈ Z N ,        u 0 n = 0 , ∀ n ∈ Z N , u 1 n = 0 , ∀ n ∈ Z N , u k 0 = φ k , ∀ k ∈ Z M , u k N − u k N − 1 = 0 , ∀ k ∈ Z M . F o r conv enience, the forward-difference stencil of this metho d has b een depicted in Fig. 1. 3.2. Energy sc heme. With the s ame notation a s in the previous par agraph, the lo cal energ y o f the system (1) at the n th no de a nd at the k th time step will be approximated by means o f the disc r ete fo r m ula (19) H k n = 1 2    δ t u k n  2 + n X j = n − 1 k +1 X l = k  δ x u l j  2 4   + µ t V ( u k n ) , where n ∈ Z N and k ∈ Z M . Meanwhile, the total ener gy of the system at time t k is calcula ted through the expressio n (20) E k = X n ∈ Z N H k n + 1 2 k +1 X l = k  δ x u l 0  2 4 . Before closing this stage o f our inv estiga tion, it is impor tan t to p oin t out that the lo cal e nergy function H n in (6) is nonnega tiv e for the case of the double sine- Gordon regime; in additio n, its discr ete co un terpart, namely , Eq. (19), is likewise nonnegative. It follows that the total ener g y of the system (1) a s g iv en by (8), as well as the discre te total energy (2 0) a re b oth nonnega tive a t any time. 8 J. E. MAC ´ IAS-D ´ IAZ 0 50 100 150 200 0 100 200 0 1 2 3 n t H n (t) 0 50 100 150 200 0 100 200 0 1 2 3 n t H n (t) Figure 5. Graphs of approximate lo cal energ y v ersus no de site n and time t , of a system (1) of length N = 2 00, for c = 4, γ = 0 and a p otent ial given b y (4). T he system was p erturb ed by means of (3 ), with Ω = 0 . 9 , and tw o different amplitude v alues: A = 1 . 03 (top g raph) and A = 1 . 04 (rig h t column). 3.3. Numerical prop erties. As men tioned previously , the n umerica l metho d pre- scrib ed by the expr essions (18), (19) a nd (20) preserves the p ositivit y character of the lo cal and the tota l energy of the system (1). Moreov er, the finite-difference schemes presented in (18) provide consis ten t solutions of (1) of order the second o r- der in time (see App endix A for a brief discussion of the consis tency of the metho d). The fact that the lo cal ener gy estimate (19) is a consistent approximation of the co n tin uous lo cal e ne r gy (6), and that the discrete total energ y (20), in turn, is a consistent estimation of the cor respo nding co n tinuous expressio n (8), is eviden t. The following r e sult shows that this c onsistency is also preser v ed o n the grounds of the rate of change of energy with resp ect to time. Prop osition 2 . Consider the finite-differ enc e scheme (18), with lo c al ener gy given by (19), and total ener gy (20). Then, the discr ete r ate of change of ener gy of the metho d at time t k − 1 is given by (21) δ t E k − 1 = − c  µ (1) t δ x u k 0   δ (1) t u k 0  − γ X n ∈ Z N  δ (1) t u k n  2 Pr o of. See Appendix B.  3.4. Computational remarks. Clear ly , the finite-difference scheme (18) is non- linear and implicit when V is not a co nstan t function, as it is the cas e of the double sine-Gordon chain. Thus, in order to appr o ximate the solution of the system (1) at time t k +1 when the approximations at times t k and t k − 1 are at hand, we employ Newton’s metho d for nonlinear systems of equations. Once, aga in, let us a dopt the co n ven tions of Se c tion 3.1. F or every k ∈ Z M , let u k = ( u k 0 , u k 1 , . . . , u k N ), and let f n be the left-hand s ide of the n th differ ence equation o f (18), that is, let (22) f n ( u k ) = h δ (2) t − µ (1) t δ (2) x + γ δ (1) t + δ (1) u V i ( u k n ) . ENERGY SIMULA TION IN DOUBLE SINE -GORDON 9 0.8 0.9 1 1.1 1.2 0 2000 4000 6000 8000 10000 12000 Amplitude Total Energy Figure 6. Graph of a ppro ximate total energ y ov er the time p e- rio d [0 , 200] o f the undamp ed system (1) v ersus driving amplitude, sub ject to harmonic driv ing of the fo r m (3) and a p oten tial (4). The para meters c = 4, N = 20 0 and Ω = 0 . 9 w ere employed in the simulations. for every n ∈ Z N . Additionally , let f 0 ( u k ) = u k 0 − φ k , (23) f N ( u k ) = u k N − u k N − 1 . (24) Moreov er, let f = ( f 0 , f 1 , . . . , f N ). Using a re cursiv e pro cess, as sume that the vectors u k and u k − 1 hav e b een prev iously co mputed. Then (25) u k +1 = u k − y , where y is the ( N + 1)-dimensio nal vector which s atisfies the ma trix equation (26) J ( u k ) y = − f ( u k ) . Evidently , the ( N + 1 ) × ( N + 1) matr ix J is the Jaco bian matrix of f , which is given by (27) J ( u k ) =          1 0 0 0 · · · 0 0 0 a d 1 a 0 · · · 0 0 0 0 a d 2 a · · · 0 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 · · · a d N − 1 a 0 0 0 0 · · · 0 − 1 1          , 10 J. E. MAC ´ IAS-D ´ IAZ 0.2 0.4 0.6 0.8 1 0 2 4 0 1 2 3 4 x 10 5 Ω A E 200 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 Ω A −1 −0.5 0 0.5 1 1.5 2 2.5 x 10 5 Figure 7. Graph of approximate total ener gy ov er the time perio d [0 , 200] of the unda mp ed system (1) v ers us driving amplitude and driving frequency (left), sub ject to ha rmonic driving of the fo r m (3) and a p otential (4). The par ameters c = 4 and N = 20 0 were employ ed in the simulations. The left graph is the chec kb oard plot of the top one. where a = − c 2 2 , (28) d n = 1 (∆ t ) 2 + c 2 + γ 2∆ t (29) +  u k +1 n − u k − 1 n  V ′ ( u k +1 n ) + V ( u k − 1 n ) − V ( u k +1 n )  u k +1 n − u k − 1 n  2 , for every n ∈ Z N . The tridiago nal system (26) is solved then employing Crout’s reduction technique with pivoting [17]. Of c o urse, for our simulations, Newton’s metho d requires o f a sto pping criterion in or der to a ppro ximate the vector u k +1 in Eq. (25). Particularly , in this work, this criterio n is g iv en by the condition k y k 2 < ǫ , where the toler ance par ameter ǫ is equal to 1 × 10 − 4 , a nd k · k 2 is the classical Euclidea n nor m in R N +1 . 4. Simula tions Throughout this section, we co nsider a system g o verned by (1 ), where the driving function assumes the sinusoidal form (3). In o rder to avoid the cr e a tion of sho ck wa ves pro duced by the s udden mov ement of the driving b oundary at the initial time, we linearly increase the driving amplitude from zero to its actua l v a lue A during a finite p e riod of time T 0 . P ar ticularly , in the simulations p erformed in this study , we fix T 0 = 50. 4.1. Sine-Gordon c hain. As a means to verify the v alidit y of our metho d, we consider, first of all, a discrete c hain of har monic oscilla tors coupled through iden- tical springs, in which case, the gov erning equations a re given by (1), with V ( u ) = 1 − cos u . Moreover, a ssume tha t the s ystem under study is undamp ed, let c = 4, N = 200 a nd fix a driving frequency of 0 . 9. Acco rding to [4], the supratrans mis sion threshold o f the system occ ur s ar o und the critical v alue A s = 1 . 78. ENERGY SIMULA TION IN DOUBLE SINE -GORDON 11 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 Ω A s Figure 8 . Gr aph o f a ppro ximate driving a mplitude A s ab o ve which supratrans mis s ion o ccurs in the undamp ed sys tem (1) ver- sus driving frequency Ω, sub ject to har monic driving of the form (3) and a p otential (4). The par ameters c = 4 and N = 20 0 were employ ed in the simulations. F r om a computational pe rspective, w e let ∆ t = 0 . 05, and compute appr o xi- mations to the actual solution of the initial- boundary -v alue pr obme (1) and the corres p onding lo cal energy , over a time interv a l of length T = 20 0. Under these circumstances, Fig. 2 presents the temp oral evolution of the solution and the lo cal energy of the 60 th no de of the system fo r t wo v a lues of the driving amplitude: A = 1 . 77 (top gra ph) a nd A = 1 . 78 (b elo w graph). The results show a drastic change in the qualitative b ehavior of s olutions a r ound the prop osed critical ampli- tude A s . The s e results a re clearly in agr eemen t with [4 ], and they a re considered as evidence in fa vor of b oth the v alidity o f our metho d a nd the existence o f supra- transmission in the s ine-Gordon chain. 4.2. Double sine-Gordon c hain. As mentioned befo re, the study of the pheno m- enon of no nlinear supratr ansmission of energ y in the double sine- Go rdon chain is a topic of interest that has b een left aside. In this section, howev er, we pro ceed to compute bifurca tion dia grams similar those co ns tructed to predict the pro cess o f supratrans mission in disc rete sine-Gordon and Klein-Gor don sys tems (see [4 , 6 ]). So, as in the pr evious s ta ge of our inv estigation, we consider a n undamp ed sy stem gov erned by (1), with parameters c = 4 , N = 200 , T = 200, Ω = 0 . 9, and p otential given by (4). Computationa lly , let ∆ t = 0 . 05. With these co ns iderations, Fig. 3 presents the time-dep e nden t gr aphs o f the solution and the lo cal energy of the 6 0th no de o f the system, for tw o different v a lues of the driving amplitude, namely , A = 1 . 03 (top graph) and A = 1 . 04 (bo ttom gr aph). As in the cas e of the discrete sine-Gordo n system, we obser ve a drastic qualitative difference in the b ehavior of the solution and the lo cal ener gy of 12 J. E. MAC ´ IAS-D ´ IAZ 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 4 Amplitude Total Energy γ = 0 γ = 0.01 γ = 0.02 γ = 0.03 Figure 9. Graph of a ppro ximate total energ y ov er the time p e- rio d [0 , 20 0] of the system (1) versus driving a mplitude, sub ject to harmonic driving of the form (3) and a p o ten tial (4). The par ame- ters c = 4, N = 200 and Ω = 0 . 9 w ere employed in the sim ulations, with four differ en t v a lues of the damping co efficient: γ = 0 (solid), γ = 0 . 01 (dashed), γ = 0 . 02 (das h- dotted), γ = 0 . 03 (dotted). the 6 0th no de, aro und the critical v alue A s = 1 . 04. Indeed, this obser v ation is in agreement with the av a ilable litera ture [6 ]. Fig. 4 presents the solution of the s ystem studied in the previous paragraph, with resp ect to no de site n and time t , where t ∈ [0 , 20 0]. The graphs of the solutions clearly change drastically for the t wo driv ing a mplitudes considere d: for A = 1 . 03 , the b oundary obviously do es not propag ate wa ve signals in to the s ystem; on the contrary , the gra ph corr esponding to A = 1 . 04 shows transmission of energy into the medium. This obser v ation is verified in Fig. 5, which pr esen ts the c o rresp onding graphs o f lo cal energ y for the tw o amplitudes co nsidered. Evidently , the qualitative observ ations done in the doma in of the so lutions car ry ov er to the domain of the lo cal energy of the system. F o r the next step in our discussio n, we define the total energy of the system (1) ov er the time interv al [0 , T ] as (30) E T = Z T 0 E ( t ) dt, where E ( t ) is the total energy of the system at time t , given by expressio n (8). Clearly , E T is consis tently approximated by means of the formula (31) E ′ T = M − 1 X k =1 E k ∆ t, where ea ch E k is g iv en by (20). ENERGY SIMULA TION IN DOUBLE SINE -GORDON 13 With this notation, Fig. 6 presents the total energy over the tempo ral in terv al [0 , 200], of a sy s tem (1) with the same parameter s as a bov e, when the driv ing ampli- tude ta kes o n v a lues in the interv al [0 . 8 , 1 . 3]. The graph ev iden tly sho ws the drastic change in the behavior of the to tal energ y of the system b efore a nd after the ampli- tude v alue 1 . 0 4. With this strong ev idence of the existence of supratr a nsmission in the double sine-Gordon chain, the critical v a lue A s = 1 . 04 is immediately identified as the nonlinear supr a transmission threshold for Ω = 0 . 9. Obviously , these r e s ults are in per fect agreement with [6]. W e hav e p erformed similar simulations for several v alues of Ω in the interv a l [0 . 2 , 1], and v alues of A in [0 , 4], a nd we ha ve obtained q ualitativ ely identical r esults. Indeed, Fig . 7 summariz e s our findings , in the form of a graph of total energy over the p erio d of time [0 , 20 0 ], v ersus driving frequency and driving amplitude. Clear ly , for every such driving frequency , there exists a sma llest driving amplitude A s ab o ve which the system b egins to absor b energ y from the b oundary . F ro m here, a graph of A s versus driving freq uency is obtained and presen ted as Fig. 8. It is worth noticing that the results ar e in qualitative ag reemen t with those obtained for sine- Gordon chains [4, 6]. Finally , it must be mentioned that the numerical metho d employed in this work is a lso useful in order to establish the effects of damping in a dis crete double sine- Gordon chain governed by (1). Indeed, co nsider a system consisting of N = 200 no des coupled thro ugh (1), with c = 4, harmonically p erturb ed by the driving function (3 ) with Ω = 0 . 8, ov er an int erv al of time [0 , 20 0], wher e the p o ten tial function V is given by (4). Fig. 9 presents the effect of the driving a mplitude A on the total energ y of the system, for v alues of A in the interv al [1 , 2], a nd three different v alues of the damping co efficient, na mely , γ = 0, γ = 0 . 01, γ = 0 . 02 a nd γ = 0 . 03. The r esults show the exp ected de c rease in the total energy of the system as γ is increased a nd, mor eo ver, they show that the supra transmission threshold is slightly delayed with the presence of damping. 5. Conclusions In this work, w e have employed a numerical metho d in the study o f the o ccur- rence of the pro cess of nonlinear supratransmis s ion in a discrete chain of oscilla - tors coupled with identical spr ings. The method prop osed is consistent of or de r O ((∆ t ) 2 ), and it is ass ociated to a dis crete scheme to a ppro ximate the lo cal energy of the chain, a s well as a scheme for the total ener gy of the s ystem. Both energ y schemes consistently approximate their contin uous counterparts, a nd the metho d has the prop ert y that the discrete r ate of change of energ y a lso approximates the corres p onding c on tinuous r ate o f change. The metho d was qua lita tiv ely tested aga inst known approximations to the o c- currence of the phenomeno n of no nlinear supratra ns mission in discrete sine-Gordon and double sine-Gordo n chains. The simulations obtained with our metho d ar e in- deed in excellent a greement with the results av ailable in the literature. Moreover, the metho d was employed in the construction o f a bifurcation diagram o f smallest driving amplitude at which supratra nsmission starts in the undamp e d system, ver- sus driving frequency . The graph is actua lly in qua litativ e agreement with those found in the literature for dis crete sine-Gor don and K le in-Gordon chains, which are systems with the sa me forbidden ba nd-gap region. Moreov er, when damping is present, our simulations s ho w that the pr ocess o f supratra nsmission is s till pres en t 14 J. E. MAC ´ IAS-D ´ IAZ in the system under in vestigation, and that the app earance of the critical amplitude v a lue is delay ed as the damping co efficien t inc r eases. Of course , many aven ues of r esearch still remain o p en. Thu s, from a pr actical po in t of view, it is impo rtan t to provide applications of the res ults presented in this work. More concretely , following [9], it is interesting to prop ose a pplications of the pro cess o f nonlinear supratrans mis s ion to the des ign o f amplifiers of weak s ignals, or to the fabr ication of detectors of ultra weak pulses, as it has b een done for the Klein-Gordo n equation [5, 9]. A cknow le dgments. The author would like to acknowledge the enlightening com- men ts of the anonymous reviewers, w hich led to improv e the ov erall quality of the final version of this ma n uscript. Also, he would like to thank Dr. F. J. ´ Alv a rez Ro dr ´ ıguez, dea n of the F aculty of Sciences at the Universidad Aut´ onoma de Aguas- calientes, and Dr. F. J. Av elar Gonz´ alez, Director o f the Office for Graduate Studies and Research of the same universit y , for uninterestedly providing the computational resource s to pr oduce this a rticle. This work presents the final re sults of pro ject PIM08-1 at this universit y . Appendix A. Consistency study A brie f cons istency analysis o f the finite- difference schemes (18) r ev eals that the nu merica l method prop osed in this work is c onsisten t o f or der O ((∆ t ) 2 ). In fact, observe that δ (1) t u k n ≈ du n dt ( t k ) + (∆ t ) 2 12 d 3 u n dt 3 ( t k ) , (32) δ (2) t u k n ≈ d 2 u n dt 2 ( t k ) + (∆ t ) 2 12 d 4 u n dt 4 ( t k ) , (33) for every n ∈ Z N and e v ery k ∈ Z M . Moreover, (34) µ (1) t δ (2) x u k n ≈ δ (2) x u k n + (∆ t ) 2 2 δ (2) x d 2 u n dt 2 . Appendix B. Energy consistency F o r the sake of s implification, we introduce the following notation, for e very n ∈ Z N and k ∈ Z M : µ x u k n = 1 2  u k n +1 + u k n  , (35) ι k = 1 4 k +1 X l = k  δ x u l 0  2 2 , (36) h k n = 1 2 n X j = n − 1 k +1 X l = k  δ x u l j  2 4 . (37) Clearly , the term ι k is identified with the indep endent term (the term which is not prescrib ed by the summation over all n ∈ Z N ) to the rig h t-hand side of Eq . (20). Pr o of of Pr op osition 2. It is eas y to c heck that the following identities are v alid for every n ∈ Z N and k ∈ Z M : (38) 1 2  δ t u k n  2 − 1 2  δ t u k − 1 n  2 =  δ (2) t u k n   δ (1) t u k n  ∆ t, ENERGY SIMULA TION IN DOUBLE SINE -GORDON 15 (39) µ t V ( u k n ) − µ t V ( u k − 1 n ) =  δ (1) u V ( u k n )   δ (1) t u k n  ∆ t. It is a tedious a lgebraic task (though rela tiv ely easy ) to verify that the following equalities ho ld, for every n ∈ Z N and k ∈ Z M : (40) δ t h k − 1 n = − 1 2  δ (1) t u k n   µ (1) t δ (2) x u k n  + c 2  δ (1) t u k n +1   µ (1) t δ x u k n  − c 2  δ (1) t u k n − 1   µ (1) t δ x u k n − 1  . As a co nsequence, (41) X n ∈ Z N δ t H k − 1 n = X n ∈ Z N nh δ (2) t − µ (1) t δ (2) x + δ (1) u V i ( u k n ) ·  δ (1) t u k n o − c  µ x δ (1) t u k 0   µ (1) t δ x u k 0  = − γ X n ∈ Z N  δ (1) t u k n  2 − c  µ x δ (1) t u k 0   µ (1) t δ x u k 0  . Moreov er, (42) δ t ι k − 1 = 1 2  δ (1) t δ x u k 0   µ (1) t δ x u k 0  . The conclusio n of Prop osition 2 is now rea c hed by computing δ t E k − 1 and simpli- fying.  References [1] M. Remoissenet, W av es Called Solitons, 3rd Edition, Springer- V er lag, New Y ork, 1999. [2] P . S. Lomdahl, O. H. Soerensen, P . L. Christiansen, Soliton excitations in Josephson tunnel junctions, Ph ys. Rev. B 25 (9) (1982) 5737–574 8. [3] V. G. Makhank ov, A. R. Bishop, D. D. H olm ( Eds.), Nonli near Evolution Equation s and Dy- namical Systems Needs ’94; Los Al amos, NM, USA 11-18 September ’94: 10th In ternational W orkshop, 1st Edition, W or ld Scientific Pub. Co. Inc., Singap ore, 1995. [4] F. Geniet , J. Leon, Energy transmission in the forbi dden band gap of a nonlinear c hain, Ph ys. Rev. Lett. 89 (2002) 134102. [5] D. Chevriaux, R. Khomeriki, J. Leon, Theory of a Josephson junction parallel arra y detector sensitive to v ery weak signals, Phys. Rev. B 73 (2006) 2145 16. [6] F. Geniet , J. Leon, Nonlinear supratransmission, J. Phys.: Condens. Matter 15 (2003) 2933– 2949. [7] R. Khomeriki, S. Lepri, S. Ruffo, Nonlinear supr atransmi ssion and bistability in the F ermi - Pa sta-Ulam mo del, Phys. Rev. E 70 (2004) 066626. [8] S. A. Leon, J., Gap soliton formation by nonlinear supratransmiss ion in Bragg media, Phys. Lett. A 327 (2004) 474–480. [9] R. Khomeriki, L. J., Bistability in sine-Gordon: The ideal switch, Ph ys. Rev. E 71 (2005) 056620. [10] A. R a vi Kant h, K. Aruna, Differential transform metho d for solvi ng the linear and nonlinear Klein–Gordon equat ion, Computer Ph ysics Communications 180 (5) (2009) 708–711. [11] Z. W ang, Discrete tanh method for nonlinear difference-differen tial equations, Computer Ph ysics Communicat ions 180 (7) (2009) 1104–1 108. [12] M. Dehghan, A. Ghesmati, Numerical simulation of t wo-dimensional sine-Gordon soli tons via a l ocal weak meshless tec hnique based on the radial point interpolation method (RPIM), Computer Ph ysics Communicat ions 181 (4) (2010) 772–786 . [13] A. Aydın, B. Karas¨ oze n, Symplectic and multi-symplectic methods for coupled nonlinear Sc hr¨ o dinger equations with p eriodi c solutions, Computer Physics Communications 177 (7) (2007) 566– 583. 16 J. E. MAC ´ IAS-D ´ IAZ [14] N. Cassol-Seewald, M. Copetti, G. Krein, Numerical appro ximation of the Ginzburg–Landau equation with memory effect s in the dynamics of phase transitions, Computer Ph ysics Com- mun ications 179 (5) (2008 ) 297–309. [15] J. E. M ac ´ ıas-D ´ ıaz, A. Puri, A numerical method for computing radially symmetric solutions of a dissipative nonlinear modi fied Klein-Gordon equation, Numer. Meth. Part. Diff. Eq. 21 (2005) 998– 1015. [16] A. L. F abian, R. Kohl, A. Biswas, Perturbation of top ological soli tons due to s ine-Gordon equation and i ts t yp e, Communications in Nonlinear Science and Numerical Simulation 14 (4) (2009) 1227 –1244. [17] R. L. Burden, J. D. F aires, Numerical Analysis, 4th E di tion, PWS-KENT Publishing Com- pan y , Boston, MA, 1989. Dep ar t am ento de Ma tem ´ aticas y F ´ ısica, Un iversid ad A ut ´ onoma de Aguascalientes, A venida Universidad 940, Ciudad Universit aria, Aguascalientes 2013 1, M exico E-mail addr ess : jemacias@co rreo.uaa.mx