Two finite-difference schemes that preserve the dissipation of energy in a system of modified wave equations

TW O FINITE-DIFFERENCE SCHEMES THA T PRESER VE THE DISSIP A TION OF EN E R GY IN A SYSTEM OF MODIFIED W A VE EQUA TIONS J. E. MA C ´ IAS-D ´ IAZ AND S. JEREZ-GALIANO Abstract. In this work, w e presen t t w o n umerical methods to appro ximate solutions of systems of dissipative sine-Gordon equations that arise in the study of one-dimensional, semi-infinite arra ys of Josephson junctions coupled through superconducting wires. Also, we presen t schemes for the total ene rgy of such systems in asso ciation with the finite-difference sc hemes used to approximate the s olutions. The prop osed methods ar e conditionally stable techniques that yield consisten t approximations not only in the domains of the solution and the total energy , but also in the appro ximation to the rate of c hange of energy with r espect to time. The metho ds are employ ed in the estimation of the threshold at which nonlinear supratransmiss ion takes place, in the presence of parameters such as in ternal and external damping, generalized mass, and generalized Josephson curr en t. Our results are qualitatively in agreement with the corr esponding problem i n mechanical ch ains of coupled oscill ators, under the pr esence of the same parameters. 1. In troduction The well-kno wn sine-Gordon equation is an impor tant diff erential equation in ph ysics due to its multiple applicatio ns . F or example, it appea rs in the study of long Josephson junctions betw een superco nductors [1], in the in vestigation o f discrete arrays of Josephs o n junctions coupled through s uper conducting wir es [2], in the mo deling of the motion of a damp ed str ing in a non- Hoo kean medium or, more generally , in the description of the motion of rigid p endula attached to a stretched wire [3], and a s a model for rapidly rotating baro c linic fluids [4]. Mo dified versions of the sine-Go rdon eq uation often a ppea r in the form of Klein-Gordon equa tions, in pr oblems suc h as the study of fluxons in Josephso n transmiss ion lines [5], in the statistical mechanics of nonlinear coherent structures such as solita ry wa ves (see [6] pp. 29 8–30 9 ), in the study of Alfv en wa ves in nonuniform media [7], o r as T aylor approximations to pro blems tha t in volv e the sine-Gordon model. The contin uous form of the sine-Gor don eq ua tion p ossesses many interesting ana- lytical prop erties. F or example, the s ine - Gordon equation is an integrable equation, a characteristic which do es not sha re, for instance, with the nonlinea r K lein-Gordon equation. The existence of tr av eling-wa ve solutions, so liton solutions and, more- ov er, nonlinea r intrinsic mo des (ca lled moving br e athers ), ar e also well-kno wn facts asso ciated to this equation [8 ]. In summary , many interesting pro per ties ar e alr eady av ailable for the analy sis of contin uous sine - Gordon systems; ho wev er, the theory 2010 Mathematics Subje ct Classific ation. (P ACS) 45.10.-b; 05.45.- a; 02.30.Hq. Key wor ds and phr ases. finite-difference sche mes; sine-Gordon equation; discrete Josephson- junction arrays; nonlinear supratransmission; stabili t y analysis. 1 2 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO underneath the foundations o f this differential equation is ra ther far from being completely understoo d, whence it follows that the determination of novel results related to con tinuous s ine-Gordon media, as well as the study of new physical ap- plications, are hig hly transited high w ays in the literature of applied mathematics and ph ysics. F r om a numerical per sp ective, the study of the classic al, contin uous sine-Go rdon equation has b een a popular source o f in terest since its inception [9]. F or instance, Josephson tunnel junctions were studied n umerically via the sine-Gordo n equation in [5], while the numerical study of wav e tra nsmission in long Jo sephson structures sub j ect to harmonic driving was b egun b y Olsen and Samuelsen [10]. O n the other hand, the (1 + 1 )-dimensional Langevin equation was solved throug h num erica l sim- ulations in [11]. Later o n, Strauss and V´ azquez designed a numerical tec hnique to approximate radia lly s ymmetric solutions o f conserv ative Klein-Gor don equatio ns [12], one of the most imp or tant features of their n umerical metho d b eing that the discrete energy asso ciated with the differen tial equation is conserved. It must be men tioned tha t, after Stra uss and V´ azquez’s pio ne e r ing work, the developmen t of numerical techniques with prop erties of inv ar iance in so me physical domains bec ame a very fruitful topic of resea rch. Particularly , many finite-difference schemes (implicit and explicit) with pr op erties of inv aria nce in the energy domain were propo sed for classe s of nonlinear , co nserv ativ e Klein-Go rdon equations [12, 13, 1 4, 15]. How ev er, it is worth ment ioning that the nonco nserv ativ e case has bee n left a side, a nd that its study has been recently initiated [16, 17 , 18, 19], motiv a ted by the study of the nonlinear phenomeno n of supratransmissio n of energy in c ertain dissipativ e s ystems, where the sine-Go rdon eq uation clearly dominates the theoretical [2, 20, 21] and the applied [22, 23] scena rios. This article pres ent s tw o conditionally stable, finite-difference schemes with pro p- erties of consistency in the doma ins of the solution and the energy . The numerical techn iques are used in the inv estigation of the pro cess of nonlinea r supr atransmis- sion in discrete arrays o f Josephson junctions sub ject to harmonic p erturbatio ns in one end, for which the metho ds pr e sented her e seem to be natural c hoices in the study of the problem. Section 2 of this pap er introduces the mathematical mo del under study and the energy expressio ns emplo y ed in the analysis of supra transmission. Sec tio n 3 presents t w o n umerical schemes ass o ciated to o ur problem and some pra ctical ob- serv ations on the implementation of the metho ds. The next section e s tablishes the mo s t imp or tant numerical pr op erties o f the computatio nal techniques; her e, we show that the pro p o sed metho ds pr ovide co nsistent approximations in the energy domain — a highly desirable characteristic in the analysis of the pro cess of supra- transmission —, and establish stabilit y prop erties. In Section 5, our techniques are applied to the determination of the supratrans mis s ion thresho ld in the pres ence of several parameter s. Finally , the a rticle closes with a section of co ncluding r emarks and directions of further research. 2. P rel iminaries 2.1. Mathematical mo de l. The ph ysical motiv ation of the present work is the study o f discrete, linear arrays of parallel Josephson j unctions coupled through su- per conducting wires, sub ject to harmonic driv ing on one end. The model co nsiders the effects of int ernal and external damping, r elativistic mass, and a gener alized TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 3 Josephson curre n t. More co ncretely , w e consider a s e quence of N co upled Jo s ephson junctions initially at rest, and study the assoc ia ted initial- v a lue problem (1) ¨ u 1 −  c 2 + β D t  ∆ 2 x u 1 + γ ˙ u 1 + m 2 u 1 + sin u 1 = J + φ ( t ) , ¨ u n −  c 2 + β D t  ∆ 2 x u n + γ ˙ u n + m 2 u n + sin u n = J, 1 < n < N , ¨ u N +  c 2 + β D t  ∆ 2 x u N − 1 + γ ˙ u N + m 2 u N + sin u N = J − I , s. t. ( u n (0) = 0 , 1 ≤ n ≤ N , du n dt (0) = 0 , 1 ≤ n ≤ N . Here, D t is the deriv ative opera tor with respect to time t , (2) ∆ 2 x u n = ( u n +1 − u n ) χ [1 ,N ] ( n + 1 ) + ( u n − 1 − u n ) χ [1 ,N ] ( n − 1) , and χ [1 ,N ] ( n ) is the characteristic function on the set [1 , N ] ev alua ted at n , which is equal to 1 if n b elongs to [1 , N ], and e qual to 0 otherwise. Notationally , c is a no nnegative constant called the c oupling c o efficient , and β and γ are also nonnega tive co nstants called the internal and external damping c o efficients , resp ectively . The cons ta nt m — called the mass of the sy s tem — is a real or pure- imaginary complex num ber tha t has b een included in the system of equations in order to sug g est p o ssible, further applications of our technique to phi- four theories and sup e r symmetry [24]. Meant ime, the function φ is called the input intensity function , and it is a har mo nic disturbanc e irr adiating a t a freq uenc y in the forbidden band-ga p of the s ystem, that takes the concrete form φ ( t ) = A sin(Ω t ). Finally , the constants J and I a re called, resp ectively , the gener alize d Josephson curr ent and the output cu rr ent intensity of the system. It has b een establis hed [25] that the inclusion of tw o conv enient no des, one lo cated at the b eginning and the other at the end of the ar r ay o f the N Josephson junctions in (1), tra nsforms our initial-v alue pr oblem into the equiv a lent initial- v a lue problem with boundary data (3) ¨ u n −  c 2 + β D t  ∆ 2 x u n + γ n ˙ u n + m 2 u n + sin u n = J, 1 ≤ n ≤ N , s. t.    u n (0) = 0 , ˙ u n (0) = 0 , 1 ≤ n ≤ N , c 2 ( u 0 − u 1 ) + β ( ˙ u 0 − ˙ u 1 ) = φ, t ≥ 0 , u N +1 − u N = 0 , t ≥ 0 , where γ n = γ for every n < N , and γ N = γ + 1 / R . The num ber R is ca lled the output r e ad ing r esistanc e of the system, and it is rela ted to the output current int ensity thr ough Ohm’s la w: I = ˙ u N /R . Let us cons ide r the ca s e β = γ = J = 0. Then, the linear disp ersio n rela tio n asso ciated with the linea rized form of the differential e q uations in (3) is g iven by ω 2 ( k ) = m 2 + 1 + 2 c 2 (1 − cos k ). It is clear that 1 − cos k ≥ 0 for every k , so that ω 2 ( k ) ≥ m 2 + 1 . It follows that the forbidden band-gap is pr ovided b y the inequality Ω < √ m 2 + 1. 2.2. Energy expressions . It is worth noticing that the Hamiltonian of the n -site in the conserv a tive version of (3) is given b y (4) H n = 1 2  ˙ u 2 n + c 2 ( u n +1 − u n ) 2 + m 2 u 2 n  + V ( u n ) − J u n , with V ( u n ) b eing the p otential function for the classica l sine-Gordo n equation ev a luated at u n , namely , V ( u n ) = 1 − c os u n . The inclusion of the p otential b etw een 4 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO the co upling o f the first t w o sites in the chain of Josephso n junctions results in the following expres sion for the tota l energy of the system: (5) E = N X n =1 H n + c 2 2 ( u 1 − u 0 ) 2 . As it is customary in the analysis of supra transmission, it is conv enient to co n- sider semi-infinite, discr e te arrays of no de s , in order to determine the critical thresh- old at whic h the sy stem under study starts to propaga te wa v e signa ls in the form of lo caliz e d mo des . T o do that, one may even tually need to cons ide r infinite se- quences ( u n ) ∞ n =1 of rea l functions satis fying the differential equations o f problem (3), for which the sequences of their co rresp onding deriv atives are members of ℓ 2 ( R ). More precisely , o ne needs to consider seq uences ( u n ) ∞ n =1 for which P ∞ n =1 ( ˙ u n ) 2 is conv ergent a t an y time. This type of sequenc e s will b e called squar e-summable . Prop ositio n 1 (Mac ´ ıas-D ´ ıaz and Pur i [17]) . L et ( u n ) N n =1 b e a solution of the system of differ ent ial e quations in (3 ). Then, the instant ane ous r ate of change of the total ener gy of the system with r esp e ct to t ime is given by (6) dE dt = c 2 ( u 0 − u 1 ) ˙ u 0 − β " N X n =1 ( ˙ u n − ˙ u n − 1 ) 2 + ( ˙ u 1 − ˙ u 0 ) ˙ u 0 # − γ N X n =1 ˙ u 2 n . Mor e over, if ( u n ) ∞ n =1 is a squar e-summable solut ion of an infi nite class of c ouple d differ ent ial e quations describ e d by (3), then the instantane ous r ate of change of total ener gy of the system is obtaine d by taking the limit when N tends to infinity in the right-hand side of (6).  Corollary 2 (Geniet a nd Leon [21]) . L et ( u n ) N n =1 b e a solution of the undamp e d version of (3). Then, the total ener gy of the system is given by (7) E ( t ) = − c 2 Z t 0 ˙ u 0 ( s )( u 1 ( s ) − u 0 ( s )) ds.  3. Comput a tional techniques 3.1. A finite-diffe rence sc heme. Consider a class of N differential equatio ns satisfying (3), and let us take a reg ular partition 0 = t 0 < t 1 < · · · < t M = T of a time interv al [0 , T ], with time step equal to ∆ t . F o r each k = 0 , 1 , . . . , M , let us represent by u k n the approximate so lution to our problem on the n th lattice site at time t k , and the actual v alue o f φ at the k th time step by φ k . Let us conv ey now that (8) δ t u k n = u k +1 n − u k − 1 n , δ 2 t u k n = u k +1 n − 2 u k n + u k − 1 n , δ 2 x u k n = u k n +1 − 2 u k n + u k n − 1 , ¯ δ u k n = u k +1 n + u k − 1 n . TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 5 ✻ t ✲ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ n n − 1 n + 1 t k t k +1 t k − 1 ✉ ✉ ✉ ✉ × × × Figure 1. F orward-difference stencil for the approximation to the n th differe ntial equation in (3) at time t k , using scheme (9). The black cir cles repre s ent known approximations to the actual solu- tions at times t k − 1 and t k in s cheme (9), and the cr osses denote the unknown a pproximations at time t k +1 . Then, problem (3) tak es the discrete form (9) δ 2 t u k n (∆ t ) 2 −  c 2 2 ¯ δ + β 2∆ t δ t  δ 2 x u k n + γ 2∆ t δ t u k n + m 2 2 ¯ δ u k n + V ( u k +1 n ) − V ( u k − 1 n ) u k +1 n − u k − 1 n = J,        1 ≤ n ≤ N , 1 ≤ k ≤ M , s. t.    u 0 n = 0 , u 1 n = 0 , 1 ≤ n ≤ N , c 2 ¯ δ ( u k 0 − u k 1 ) + β δ t ( u k 0 − u k 1 ) / ∆ t = 2 φ k , 1 ≤ k ≤ M , u k N +1 − u k N = 0 , 1 ≤ k ≤ M . The forw ard- differ ence s tencil for this method is presented in Fig. 1. In the unbounded situation (that is, when we consider a semi-infinite system of coupled Jo sephson junctions), we choose a large sy stem of N Jose phs on junctions following (3 ), in which the external damping co efficient includes b o th the effect of external damping and a simulation of an a bsorbing b oundar y slowly increasing in magnitude o n the la st N − N 0 oscillator s. Ex plicitly , we let γ b e the sum of the actual v alue of the external damping co efficient and the function (10) γ ′ ( n ) = 0 . 5  1 + tanh  2 n − N 0 + N 6  , where N 0 = 50 and N ≥ 200 for our computations. It is impo rtant to posses discrete schemes to approximate the lo cal ener gy den- sity and the tota l energy of the pr oblem under study . W e asso cia te the fo llowing expression to finite-difference scheme (9), in order to approximate the v alue of the 6 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO ✻ t ✲ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ n n − 1 n + 1 t k t k +1 t k − 1 ✉ ✉ ✉ ✉ ✉ ✉ × × × Figure 2. F orward-difference stencil for the approximation to the n th differen tial equation in (3) a t time t k , using scheme (13). The black cir cles repre s ent known approximations to the actual solu- tions at times t k − 1 and t k in s cheme (13), and the cro sses denote the unknown a pproximations at time t k +1 . Hamiltonian of the n th site at the k th time: (11) H k n = 1 2  u k +1 n − u k n ∆ t  2 + c 2 8 h  u k +1 n +1 − u k +1 n  2 +  u k +1 n − 1 − u k +1 n  2 +  u k n +1 − u k n  2 +  u k n − 1 − u k n  2 i + m 2 2 ( u k +1 n ) 2 + ( u k n ) 2 2 + V ( u k +1 n ) + V ( u k n ) 2 − J u k +1 n + u k n 2 . Moreov er, the total ener gy of the system is estimated b y (12) E k = N X n =1 H k n + c 2 4 "  u k +1 1 − u k +1 0  2 +  u k 1 − u k 0  2 2 # . 3.2. A se cond fin i te-difference scheme. L e t δ 2 u k n = u k +1 n + 2 u k n + u k − 1 n . A second, finite-differ ence scheme to approximate solutions to pr oblem (3) is present ed next. Its forward-difference stencil is depicted in Fig. 2: (13) δ 2 t u k n (∆ t ) 2 −  c 2 4 δ 2 + β 2∆ t δ t  δ 2 x u k n + γ 2∆ t δ t u k n + m 2 2 ¯ δ u k n + V ( u k +1 n ) − V ( u k − 1 n ) u k +1 n − u k − 1 n − J = 0 ,        1 ≤ n ≤ N , 1 ≤ k ≤ M , s. t.    u 0 n = 0 , u 1 n = 0 , 1 ≤ n ≤ N , c 2 ¯ δ ( u k 0 − u k 1 ) + β δ t ( u k 0 − u k 1 ) / ∆ t = 2 φ k , 1 ≤ k ≤ M , u k N +1 − u k N = 0 , 1 ≤ k ≤ M . It is worth mentioning that finite-differe nc e schemes (9) and (13) ar e mo dified versions of the o ne presented in [17], whic h in turn is a spa tially discretized and adapted version of a co mputational metho d to appr oximate radially symmetric TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 7 solutions of initial-v a lue proble ms in volving a general clas s of contin uous, three- dimensional Klein-Gordo n equa tions [16]. F o r a ppr oximations o btained via (13 ), the Hamiltonian of the n th site o f the system will be co mputed using the express ion (14) H k n = 1 2  u k +1 n − u k n ∆ t  2 + c 2 4   u k +1 n +1 + u k n +1 2 − u k +1 n + u k n 2 ! 2 + u k +1 n − 1 + u k n − 1 2 − u k +1 n + u k n 2 ! 2   + m 2 2 ( u k +1 n ) 2 + ( u k n ) 2 2 + V ( u k +1 n ) + V ( u k n ) 2 − J u k +1 n + u k n 2 . With this notation at ha nd, the tota l energ y of the system will b e approximated via (15) E k = N X n =1 H k n + c 2 4 u k +1 1 + u k 1 2 − u k +1 0 + u k 0 2 ! 2 . It is imp ortant to mention that there a re several computational r easons to pr efer an implicit scheme ov er an explicit one, in order to approximate so lutions to (3). First o f all, it has bee n noted that so me explicit schemes used to a pproximate solutions to mo dified versions of our proble m hav e pr oved to b e hig hly unstable [12]; on the o ther hand, the use of an implicit scheme seems inevitable in or der to approximate the term D t ∆ 2 x u n with a consistency of order at least O (∆ t ) 2 . Also, we must mention that the prop osed numerical metho ds ar e consistent with r esp ect to the pr oblem under study , nonlinear , and require an applica tion o f Newton’s method for systems of nonlinear equations, together with an application of Crout’s reduction tec hnique for tridiagonal systems [26]. Finally , it is worthwhile mentioning her e that the stopping criter io n used for Newton’s method in our computations was (16) max {| u k +1 n − u k n | : 1 ≤ n ≤ N } < 1 × 10 − 5 . 4. N u merical prope r ties In this section, w e prove that finite-difference schemes (9) and (13) p o s ses pro p- erties of co nsistency that make them useful metho ds in the s tudy of the pro cess of nonlinear supratra nsmission in chains of coupled J o sephson junctions. With this aim in mind, we pr esent here the most impo r tant results concerning the energy prop erties of our methods. Prop ositio n 3. The discr ete r ate of change of ener gy of a se quenc e ( u k n ) satisfying finite-differ enc e scheme (9) is given by (17) E k − E k − 1 ∆ t = − β ( N X n =1  δ t u k n − δ t u k n − 1 2∆ t  2 +  δ t u k 1 − δ t u k 0 2∆ t  δ t u k 0 2∆ t ) − γ N X n =1  δ t u k n 2∆ t  2 + c 2  ¯ δ u k 0 − ¯ δ u k 1 2  δ t u k 0 2∆ t . 8 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO Pr o of. W e provide a sk etch of the pro of in view that the task is tedious but, after all, a mere algebraic w ork. T erm by ter m, it is straight -forward to chec k that N X n =1 "  u k +1 n − u k n ∆ t  2 −  u k n − u k − 1 n ∆ t  2 # = N X n =1  δ 2 t u k n ∆ t 2  δ t u k n , (18) N X n =1  ( u k +1 n ) 2 + ( u k n ) 2 2 − ( u k n ) 2 + ( u k − 1 n ) 2 2  = N X n =1  u k +1 n + u k − 1 n 2  δ t u k n , (19) N X n =1  [ V ( u k +1 n ) + V ( u k n )] − [ V ( u k n ) − V ( u k − 1 n )]  = N X n =1  V ( u k +1 n ) − V ( u k − 1 n ) u k +1 n − u k − 1 n  δ t u k n . (20) F o r the sake of conv enience, w e let (21) B k = c 2 8 N X n =1 h  u k +1 n +1 − u k +1 n  2 +  u k +1 n − 1 − u k +1 n  2 +  u k n +1 − u k n  2 +  u k n − 1 − u k n  2 i . Using finite-difference scheme (9) a nd the discrete, b oundary co nditions imp osed on the problem, it is straight-forw ard to chec k that (22) B k − B k − 1 = c 2 8 N X n =1  − 2  ¯ δ δ 2 x u k n   δ t u k n  +  ¯ δ u k n +1 − ¯ δ u k n   δ t u k n +1  −  ¯ δ u k n − ¯ δ u k n − 1   δ t u k n  +  ¯ δ u k n +1 − ¯ δ u k n   δ t u k n  −  ¯ δ u k n − ¯ δ u k n − 1   δ t u k n − 1  . Notice that the last four terms under the summation sym b ol form t w o telescoping series, which may b e readily simplified. Next, the difference E k − E k − 1 is computed. After some simplification, w e obtain that (23) E k − E k − 1 ∆ t = β N X n =1  δ t δ 2 x u k n 2∆ t   δ t u k n 2∆ t  − γ N X n =1  δ t u k n 2∆ t  2 + c 2  ¯ δ u k 0 − ¯ δ u k 1 2  δ t u k 0 2∆ t . The result follows now after s implifying the term multiplied by β , follow ed by an application of the discrete version of Green’s Firs t Identit y [17] to the sequence ( u k +1 n − u k − 1 n ) N +1 n =0 .  Prop ositio n 4. The discr ete r ate of change of ener gy of a se quenc e ( u k n ) satisfying finite-differ enc e scheme (13) is given by (24) E k − E k − 1 ∆ t = − β ( N X n =1  δ t u k n − δ t u k n − 1 2∆ t  2 +  δ t u k 1 − δ t u k 0 2∆ t  δ t u k 0 2∆ t ) − γ N X n =1  δ t u k n 2∆ t  2 + c 2  δ 2 u k 0 − δ 2 u k 1 4  δ t u k 0 2∆ t  . Pr o of. Here, we e mploy some of the identities used in the pro of o f Pr op osition 3. Notice first of all tha t if w e define C k by the for mu la (25) C k = c 2 4 N X n =1   u k +1 n +1 + u k n +1 2 − u k +1 n + u k n 2 ! 2 + u k +1 n − 1 + u k n − 1 2 − u k +1 n + u k n 2 ! 2   , TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 9 then (26) C k − C k − 1 = c 2 4 (  δ 2 u k 0 − δ 2 u k 1 4   δ t u k 0 + δ t u k 1  − N X n =1  δ 2 δ 2 x u k n 2  δ t u k n ) . W e mimic here the pro of of P rop osition 3. As a c o nsequence of the fact that the sequence ( u k n ) satisfies (13) and after an applica tion of the discrete version of Green’s First Identit y [17], it is r eadily chec ked that the discrete rate of change of the energy betw een times k a nd k − 1 is given by (27) E k − E k − 1 ∆ t = − β ( N X n =1  δ t u k n − δ t u k n − 1 2∆ t  2 +  δ t u k 1 − δ t u k 0 2∆ t  δ t u k 0 2∆ t ) − γ N X n =1  δ t u k n 2∆ t  2 + c 2 2  δ 2 u k 0 − δ 2 u k 1 4   δ t u k 1 + δ t u k 0 2∆ t  + c 2 2  δ 2 u k 0 − δ 2 u k 1 4   δ t u k 0 − δ t u k 1 2∆ t  , whence the result follo ws after an easy simplification.  The metho ds described b y finite-difference sc hemes (9) a nd (13) a re clearly c o n- sistent of order O (∆ t ) 2 . Moreover, in the following, we g ive some necessar y and sufficient conditions in or der for these metho ds to b e stable. Notice that when m is a r eal n umber, pro blem (3) repr esents a dis crete, no nlinear mo dification o f the classical linear Klein-Gordon equation. Prop ositio n 5. Consider finite-differ enc e scheme (13) wi th V ′ = 0 , J = 0 and m 2 a nonn e gative c onstant. Then, ∆ t ≤ 1 c is a ne c essary c onditio n for metho d (13) t o b e line arly stable. Pr o of. W e a pply the von Neumann s tability analysis to metho d (13) and we obtain the following augmentation matrix: (28) A ( ξ ) = ˆ f ( ξ ) ˆ g ( ξ ) − ˆ h ( ξ ) ˆ g ( ξ ) 1 0 ! , where (29) ˆ f ( ξ ) = 2 − 2 c 2 (∆ t ) 2 sin 2  ξ 2  , ˆ g ( ξ ) = 1 +  c 2 (∆ t ) 2 + 2 β ∆ t  sin 2  ξ 2  + m 2 (∆ t ) 2 + γ ∆ t 2 , and ˆ h ( ξ ) = 1 +  c 2 (∆ t ) 2 − 2 β ∆ t  sin 2  ξ 2  + m 2 (∆ t ) 2 − γ ∆ t 2 . Since a ll norms in a finite-dimens io nal space are equiv alent a nd the spectr a l radius ρ ( A ) of a square matrix A is the g r eatest low er b ound for the set o f matrix-induces norms of A , then a neces sary condition for linear stability is that ρ ( A ) ≤ 1. The eigenv alues of A ( ξ ) are: (30) λ ± ( ξ ) = ˆ f ( ξ ) ± q ˆ f 2 ( ξ ) − 4 ˆ g ( ξ ) ˆ h ( ξ ) 2 ˆ g ( ξ ) . 10 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO 2 2.5 3 3.5 4 4.5 5 5.5 0 1 2 3 4 5 6 7 8 9 10 x 10 5 Amplitude Total Energy Figure 3. Gr a ph of total ene r gy versus driving amplitude for a system (3) s ub ject to a harmonic frequency of 0 . 8, with coupling co efficient equal to 5, and a time perio d of 100 00. Different v alues of external damping were chosen: γ = 0 (solid), 0 . 1 (dashed), 0 . 2 (dash-dotted), and 0 . 3 (dotted). In order to show tha t the sp ectra l radius o f A ( ξ ) is less than or equal to 1, we will prov e that λ ± ( ξ ) ≤ 1 for all ξ ∈ [0 , π ] by c o nsidering tw o ca ses: (1) If λ ± ( ξ ) ∈ C then ˆ f 2 ( ξ ) − 4 ˆ g ( ξ ) ˆ h ( ξ ) ≤ 0 or, equiv alen tly , that 4 ˆ g ( ξ ) ˆ h ( ξ ) ≥ ˆ f 2 ( ξ ), whence it follows that ˆ h ( ξ ) is nonnegative. Moreover, (31) | λ ± ( ξ ) | = s ˆ f 2 ( ξ ) + (4 ˆ g ( ξ ) ˆ h ( ξ ) − ˆ f 2 ( ξ )) 4 ˆ g 2 ( ξ ) = s ˆ h ( ξ ) ˆ g ( ξ ) ≤ 1 . (2) If λ ± ( ξ ) ∈ R a nd ∆ t ≤ 1 c , then we obse rve fir st of all that ˆ f ( ξ ) is non- negative with ˆ f 2 ( ξ ) ≤ 4. This case leads to the chain of inequalities 1 ≥ ˆ h ( ξ ) ≥ 2 − ˆ g ( ξ ) ≥ ˆ f ( ξ ) − ˆ g ( ξ ). Clear ly , 2 ˆ g ( ξ ) ≥ ˆ f ( ξ ), so that (32) | λ ± ( ξ ) | ≤ ˆ f ( ξ ) + q ˆ f 2 ( ξ ) − 4 ˆ g ( ξ ) ˆ h ( ξ ) 2 ˆ g ( ξ ) ≤ ˆ f ( ξ ) + q (2 ˆ g ( ξ ) − ˆ f ( ξ )) 2 2 ˆ g ( ξ ) = 1 .  Prop ositio n 6 . Consider fin ite-differ enc e scheme (13) with V ′ = 0 , J = 0 and m 2 a nonn e gative c onstant. Th is scheme is line arly s table in the infi nite norm if 2 γ < ∆ t < √ 2 c . Pr o of. A sufficient condition for scheme (13) to be linearly stable is that || A ( ξ ) || ≤ 1 for all ξ ∈ [0 , π ], for s o me norm k · k ; in our case, we will consider the infinite norm TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 11 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 12 Frequency Amplitude Figure 4. Diagr am of sma lle st driving amplitude a t which supra- transmission occurs v ersus driving frequency for a system (3) sub- ject to a har monic frequency with co upling c o efficient equal to 5 . Different v alues of externa l damping were chosen: γ = 0 (so lid), 0 . 1 (dashed), 0 . 2 (dash-dotted), and 0 . 3 (dotted). k · k ∞ applied to the augment ation matrix (28) . Using the same notation as in the previous result, observe that the inequality k A ( ξ ) k ∞ ≤ 1 is satisfied in case that (33)    2 ˆ g ( ξ )    +    − ˆ h ( ξ ) ˆ g ( ξ )    ≤ 1 . In tur n, condition (33) holds if 2 ≤ ˆ g ( ξ ) + ˆ h ( ξ ) and 2 ≤ ˆ g ( ξ ) − ˆ h ( ξ ). But the former inequality is alwa ys true, while the latter is satisfied if and only if (34) 1 ≤ 2 β ∆ t s in 2  ξ 2  + γ ∆ t 2 , ∀ ξ ∈ [0 , π ] .  Corollary 7. The finite-differ en c e scheme (9) is line arly stable if and only if c on- dition 2 γ ≤ ∆ t ≤ 1 c holds. Pr o of. W e a pply a gain the v on Neumann stabilit y a nalysis for the sc heme (9). The augmentation ma trix for this scheme is the same as (28) when w e set ˆ f ( ξ ) = 2 for all ξ ∈ [0 , π ] and the pro of of Prop ositions 5 and 6 a re v alid for this cas e .  5. A pplica tion Before we provide an application of the numerical techniques presented in this work, w e must declare that the res ults o f this se ction were obtained by mea ns of (9) 12 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 12 Frequency Amplitude Figure 5. Diagr am of sma lle st driving amplitude a t which supra- transmission occurs v ersus driving frequency for a system (3) sub- ject to a har monic frequency with coupling co efficient equal to 5. Diff erent v alues of the Jose phs o n current were chosen: J = 0 (solid), 0 . 1 (dashed), 0 . 2 (dash-dotted), and 0 . 3 (dotted). and (13). Simulations in b oth the domain o f solutions a nd the energ y domain have shown that the metho ds a g ree excellently . Moreover, it is imp ortant to mention that we have p erfor med several numerical sim ulations (not shown here) with b o th metho ds, for several nonnegative v alues o f β , γ and m 2 , several v alues of the p ositive constant c 2 and b oth V ′ and J equa l to zer o, observing stability of the metho ds when the corresp onding stabilit y conditions are satisfied. The v a lidit y of the implementation o f our methods is c heck ed next against the prediction of the no nlinear supra transmission thr e shold for a semi-infinite, linear array of J o sephson junctions in parallel a nd coupled thro ugh sup erco nducting wires, for which the mathematical mo del is given by (3). Throughout, we employ a coupling co efficient equal to 5, a nd follow the systematic pro cedure us e d in [1 7], namely: (1) F or a fixe d frequency Ω in the forbidden band- g ap, we exa mine the b ehav- ior of solutions of o ur problem for different driving a mplitudes, us ing the finite-difference schemes prop osed in this article. Acco rding to the theor y , there exis ts a cr itical amplitude ab ov e which the qualitative nature o f the solutions drastically changes. (2) This drastic change in the b ehavior of the solutions is check ed next in the energy doma in, by mak ing use of the e ner gy s chemes asso cia ted to each finite-difference scheme. A sudden increase in the total ene r gy injected into the medium by the driving bo undary must be obser ved above the pre dicted supratrans mission threshold. TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 13 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 Frequency Amplitude Figure 6. Diagr am of sma lle st driving amplitude a t which supra- transmission occurs v ersus driving frequency for a system (3) sub- ject to a har monic frequency with co upling c o efficient equal to 5 . Different v alues o f internal da mping were chosen: β = 0 (so lid), 0 . 1 (dashed), 0 . 3 (dash-dotted), and 0 . 6 (dotted). (3) Finally , this pro cedure is a pplied to a re g ular partition { Ω i } l i =0 of the fre- quency interv a l [0 , √ m 2 + 1], obtaining thus a second arr ay { A i } l i =0 con- taining, for each driving freq uency Ω i , the co rresp onding cr itical amplitude A i at which s upratransmiss ion starts. In this wa y , a g raph of o ccur rence of critical amplitude versus driving frequency ma y be constructed. (4) This pro cedure may be rep eated for differen t c hoices of the parameters of the mo del s tudied. Let us fix a driving freq uency Ω = 0 . 8 in the fo rbidden band-gap of system (3), let c = 5, and set all the other scalar parameters of the model equal to zero. Select a fixed p erio d o f time of 10000, during which the system will be p e rturb ed by the harmonic driving φ ( t ) = A s in(Ω t ), wher e the v alue of A will b e chosen inside an int erv al that contains the predicted thre s hold. According to [2], the critical v alue of the con tinuous-limit ca se is pro vided b y the expression (35) A s = 2 c (1 − Ω 2 ) , which yields a critical amplitude equal to 3 . 6 in our case. In these circ ums ta nces, Fig. 3 pre sents the gra ph o f the total energy of the system for v alues of the driving amplitude b etw een 2 and 5 . 5, in the form of a s olid line. A drastic increa se in the total energy is o bserved around an amplitude of 3 . 75, in g o o d agreement with the contin uous-limit case. Let us denote by E 1 = E 1 ( A ) the total energy of the system describ ed in the pre- vious para graph ass o ciated with the amplitude A ∈ [2 , 5 . 5] at time 1 0000, obtained 14 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0 2 4 6 8 10 12 14 Frequency Amplitude P 4 P 2 P 1 P 0 P −1 P −2 P −4 Figure 7. Diagr am of sma lle st driving amplitude a t which supra- transmission occurs v ersus driving frequency for a system (3) sub- ject to a har monic frequency with co upling c o efficient equal to 5 . F o r every ℓ = 0 , ± 1 , ± 2 , ± 4, the plot P ℓ corres p o nds to a system with mass m satisfying p m 2 ℓ + 1 = 1 + ℓ 40 . using finite-difference scheme (9). Simila rly , let E 2 = E 2 ( A ) b e the corr esp onding energy obtained us ing scheme (13). It is imp ortant to ment ion that our computa- tions show that (36) max { | E 1 ( A ) − E 2 ( A ) | : A ∈ [2 , 5 . 5] } < 3 × 10 − 10 . Moreov er, for a mplitudes below the critical threshold, (37) max { | E 1 ( A ) − E 2 ( A ) | : A ∈ [2 , 3 . 5] } < 8 × 10 − 16 . Next we wish to check the effect of external damping on the o ccurrence of the supratrans mission v alue. T o this effect, we fix three different v alues of exter nal damping: γ = 0 . 1, 0 . 2 and 0 . 3; the re s ults are displayed in Fig. 3 , in which the co rresp onding gr aphs are presen ted as dashed, das hed-dotted and dotted lines, resp ectively . It is worth noticing that the pr esence of external damping pro duce s a right shift in the occurr ence of the supratra ns mission v alue and, at the same time, a decr e ase in the total energy of the system. Similar r emarks have been made for the same problem with Diric hlet bo undary data [17]. In a next sa ge, we wish to compute diagrams for the o ccurrence of supra transmis- sion in the pr esence of all pa rameters. So, we co mpute gra phs of dr iving amplitude at which supratransmissio n fir st star ts versus driv ing frequency for system (3). Un- der the same computational setting as b efor e, Fig s. 4 and 5 present such dia grams for differen t v alues of the ex ternal damping co efficient and the generalize d Jo seph- son current, resp ectively . A quick c omparison with similar results o btained using a different finite- differ ence s cheme [25] lea ds us to cer tify the v alidity of our metho d. TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 15 Additionally , we provide diagra ms which evidence the physical implications of in- cluding a rela tivistic mass and internal damping (see Figs . 7 for mas s, and 6 for int ernal damping), the results b e ing in qualitative a greement with the Dir ichlet case, in the follo wing senses: • The phenomenon of harmonic phonon quenching still appea rs in the pres- ence of external and in ternal damping. • The discrepancy reg ion due to phonon quenching is shor tened as the exter- nal damping co efficient is increa sed. • The thr e shold v alue at which supra transmissio n fir st o ccurs, for fixed fre- quencies outside the discr epancy r egion, increases as the da mping co efficient (in ternal or external) increase s . • The diagr am, for any fixed v alue of β ≥ 0, is approximately equal to the corres p o nding diagra m for the undamp ed system, shifted β horizontal units to the right. • A hor izontal shift of √ 1 + m 2 − 1 units in the diagra m of a sine- Gordon system of mass m with resp ect to the corres po nding massless sy stem, is observed for small mas ses and frequencies outside the dis crepancy region. 6. Concl u s io ns In this article, we hav e presented tw o n umerical metho ds to appr oximate solu- tions of a s ystem of differen tial equations , which mo dels discrete arrays consisting of Josephson junctions coupled with super conducting wires, in whic h each site is de- scrib ed by a mo dified, s patially discrete version of the one- dimens ional sine-Gordon equation. B o th methods ar e consistent and conditionally stable. As an imp ortant part of the metho ds intro duce d, the w ork provides co nsistent, discrete schemes for the lo cal e nergies of the links and the to tal energy of the system. These sc hemes a r e such that the discrete rates of change of total energy of the system are consistent approximations of the corresp onding cont inuous rates of change, s o that these techniques a r e ideal metho ds in the study of the pr o cess of nonlinear supratransmissio n. As applications, we obtained diagr ams of dr iving amplitude versus driving fre- quency , in order to establish the cr itical amplitude at whic h supratr ansmission starts. This s tudy was carr ied out in the presence of the parameters under study , that is , internal a nd external damping, mass and generalized Jose phson current. Our qua litative results are essentially similar to thos e o btained for the corre s po nd- ing analysis of the Diric hlet case, corrob ora ting thus the v alidit y of our technique. Of course, several directions of further resea rch still remain op en. F or instance, the design of computatio nal techniques that employ symplectic metho ds in the analysis o f the transmissio n of energy in discr ete ar r ays of Josephson junctions, is a topic that merits close a tten tion. Particularly , one may use numerical metho ds that split the vector field y ′ as the sum of a conserv ative co nt ribution f ( y ) — solved in the in terv a l [ t, t + h/ 2 ] using a symplectic metho d for the sine-Gordon equation [27, 2 8, 29] —, and a linear dissipative par t g ( y ) — solv ed exactly in [ t, t + h ] with an expo nent ial integrator, for instance ). The symplectic dyna mics may b e so lved then for the next ha lf step, obtaining th us a sec o nd-order, s plitting metho d that is linearly implicit, and symplectic in the unperturb ed cas e. 16 J. E. MAC ´ IAS-D ´ IAZ AND S. JERE Z-GALIANO Ac kno wledgments. One of us (JEMD) wishes to expr ess his deep est g ratitude to Dr. F. J . ´ Alv ar ez Ro dr ´ ıguez, dea n of the F aculty of Science of the Univ ersidad Aut´ ono ma de Agua s calientes, M´ exico, and to Dr. F. J. Avelar Gonz´ a le z, de a n of the Office for Resear ch and Graduate Studies o f the same university , for unin terestedly providing the res ources to pro duce this ar ticle. The present work represents a set of partial res ults under re search pro ject PIM07-2 at this university , and is dedicated to P atricia Arias Mu˜ no z on the occasio n of her 38 th academic anniv ersar y . References [1] M. Remoissenet. Waves Called Solitons . Spr inger-V erlag, New Y or k, third edition, 1999. [2] D. Chevrieux, R . Khomeriki, and J. Leon. Theory of a Josephson junction parallel array detect or s ensitiv e to very weak signals. Phys. R ev. B , 73 :214516, 2006. [3] A. Barone, F. Esp osito, C. J. Magee, and A .C. Scott. Theory and applications of the si ne- Gordon equation. Riv. Nuovo Cim. , 1 :227–267, 1971. [4] J.D. Gibb on, I.N. James, and I.M. Mor oz. The sine-Gordon equat ion as a model for a r apidl y rotating baro clinic fluid. Phys. Script. , 20 : 402–408, 1979. [5] P .S. Lomdahl, O.H. Soerensen, and P .L. Christiansen. Soliton excitations in Josephson tunnel junctions. Phys. R ev. B , 25 (9):5737–5748, 1982. [6] V.G. Makhanko v, A.R. Bishop, and D.D. H olm, D. D., editor. Nonline ar Evolution Equa- tions and Dynamic al Syste ms Ne e ds ’94; L os A lamos, NM, USA 11-18 Septemb e r ’94: 10th International Workshop . W orld Scient ific Pub. Co. Inc., Singapore, first edition, 1995. [7] Z.E. Musielak, J.M. F ontenla, and R.L. Mo ore. Klein-Gordon equation and reflection of Alfven w av es in nonun iform media. Phys. Fluids B , 4 : 13–18, 1992. [8] D. Zwi llinger. Handb o ok of Differ ential Equations . Academic Press, Boston, MA, thir d edi- tion, 1997. [9] J.K. P erring and T. H.R. Skyrm e. A m o del unified field equation. Nucl. Phys. , 3 1 :550–555, 1962. [10] O.H. Olsen and M. R. Sam uelsen. Hysteresis in rf-dri v en large-area Josephson junctions. Phys. R ev. B , 3 4 :3510–3512, 1986. [11] F.J. Al exander and S. Habib. Statistical mec hanics of kinks i n 1 + 1 dimensions. Phys. R ev. L ett. , 71 :955–958, 1993. [12] W.A. Strauss and L. V´ azquez. Numeri cal solution of a nonlinear Kl ein-Gordon equation. J. Comput. Phys. , 28 :271–278, 1978. [13] G. Ben-Y u, P .J. Pascual, M.J. Rodr ´ ıguez, and L. V´ azquez. Numerical solution of the sine- Gordon equation. J. Math. Comput. , 18 :1–14, 1986. [14] Z. F ei and L. V´ azquez. Two energy conserving numerical schemes for the si ne-Gordon equa- tion. Appl. Math. Comput. , 45 :17–30, 1991. [15] S. Li and L. V u-Quo c. Finite difference calculus inv ari an t s tructure of a class of algorithms for the nonlinear Klein-Gordon equation. SIAM J. Num. Anal. , 45 : 1839–187 5, 1991. [16] J.E. M ac ´ ıas-D ´ ıaz and A. Puri. A n umerical met ho d for computing r adially symmetric solu- tions of a dissipative nonlinear modified Klein-Gordon equation. Num. Meth. Part. Diff. Eq. , 21 :998–1015, 2005. [17] J. E. M ac ´ ıas-D ´ ıaz and A. Pur i . An energy-based computational metho d in the analysis of the transmission of energy i n a chain of coupled oscillators. J. Comput. Appl. Math. , 214 :393–405, 2008. [18] J.E. Mac ´ ıas-D ´ ıaz. Numerical study of th e transmission of energy in discrete arr a ys of sine- Gordon equations in tw o space di mensions. Phys. R ev. E , 77 : 016602, 2008. [19] J. E. Mac ´ ıas- D ´ ıaz. Bi t propagation in (2+1)-dimensional systems of coupled si ne-Gordon equations. Commun. Nonline ar Sci. Numer. Simul. , 14 :1025–1031, 2009. [20] F. Geniet and J. Leon. Energy transmissi on in the forbidden band gap of a nonlinear ch ain. Phys. R ev. L ett. , 89 :134102, 2002. [21] F. Ge niet and J. Leon . Nonlinear supratransmission. J. Phys.: Condens. Matter , 15 :2933– 2949, 2003. [22] J.E. M ac ´ ıas-D ´ ıaz and A. Puri. An application of nonlinear supratransmission to the propa- gation of bi nary signals in weakly damp ed, mechanica l systems of coupled oscill ators. Phys. L ett. A , 3 66 :447–450, 2007. TW O E NERG Y-PRESER VING F INITE-DIFFERENCE SCHEMES 17 [23] J.E. M ac ´ ıas-D ´ ıaz and A. Pur i. On the propag ation of binary signals in damped mec hanical systems of oscillators. Physic a D , 228 :112–121, 2007. [24] A.E. Kudrya vtsev. Solitonlike solutions for a Higgs scalar field. JETP L ett. , 22 (3):82–83, 1975. [25] J.E. M ac ´ ıas-D ´ ıaz and A. Pur i. On the transmi ssion of binary bits in discrete Josephson- junction arrays. Phys. L et t. A , 2008. Submitted. [26] R.L. Burden and J.D . F air es. Numeric al Analysis . PWS-KENT Publi shing, Boston, fourth edition, 1988. [27] B. Leimkuhler and S. Reich . Simulating Hamiltonian Dynamics . Cambridge Monographs on Applied and Computational Mathe matics. Camb ridge Unive rsity Press, Great Br itain, first edition, 2005. [28] E. Hairer, C. Lubich, and G. W anner. Ge ometric Numeric al Inte gr ation: Structur e- Pr eserving Algorithms for Or dinary Differ ential Equations . Springer-V erlag, Berlin, Ger- man y , second edition, 2006. Springer Series i n Computational Mathematics. [29] R. M cLac hlan. Symplectic integrat ion of Hamiltonian w av e equations. Numer. Math. , 6 6 :465– 492, 1993. Dep ar t amento de Ma tem ´ aticas y F ´ ısica, Universidad Aut ´ onoma de Aguascalientes, A venida Universidad 940, Ciudad Universit aria, Aguascalientes, Ags. 20100, Mex ico E-mail addr ess : jemacias@ correo.uaa.m x (CIMA T) Centro de Investigaci ´ on en Ma tem ´ aticas, A. C., Jalisco S/N, Colonia V a - lenciana, Guanaju ato , Gto. 36240 , Mex ico E-mail addr ess : jerez@cim at.mx