About this Journal Submit a Manuscript Table of Contents
International Journal of Mathematics and Mathematical Sciences
Volume 2012 (2012), Article ID 605687, 17 pages
Research Article

The Nonlinear Analysis of Perturbation Solution for a Parabolic Differential System

United Institute of Informatics Problems of the National Academy of Sciences of Belarus, Surganov Street 6, 220012, Minsk, Belarus

Received 9 April 2012; Revised 19 June 2012; Accepted 30 July 2012

Academic Editor: Yuri Latushkin

Copyright © 2012 Victor F. Dailyudenko. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


By investigation of perturbation solution for nonlinear reaction-diffusion system, we derive related differential model for perturbations that involves weak nonlinearities up to third order. For a first time, this model is shown to result in derivation of the system for amplitude distribution by means of nonlinear integration on orthogonal basis in spatial region. The obtained time-dependent system (TDS) contains all possible functional relations between the modes of wave train under consideration along with delayed relations, and after numerical simulation it provides some conclusions concerning the natural frequency of the investigated self-organization process in active medium. The related matrix and modulo operations which substantiate the derivation of the TDS are also considered.

1. Introduction

Though the linearized analysis of differential systems is used widely for exploration of stability problems and has been developed for a long time [16], additional investigations of nonlinearities occurring in corresponding analytical models provide an opportunity to obtain important information that is impossible to find restricting with the use of linearizations only. So, in the present work, in difference from [6], we also take into consideration nonlinear terms of the second order (and even higher ones) and show that such exploration allows us to consider the amplitude distribution of the modes constituting the wave train and, moreover, to obtain some estimations for temporal changes of such distribution.

As a result of numerous experimental and theoretical investigations, it was shown that self-organization patterns described by reaction-diffusion systems of differential equations [3, 7, 8] (or proper propagating-recovery model [810]) are just long-living stationary structures with continuous radiation that can be represented through the wave trains propagating in active media (e.g., spiral waves [8, 9, 1113]). Such periodic wave train can have any period (or frequency) within a certain range, and to each period corresponds a particular wavelength (or phase velocity) [810]. So, representation of the wave train through normal modes [1] is suitable both from experimental and theoretical viewpoints (yet allowing us an opportunity to link approaches of linear analysis with considerations of different levels of nonlinearities), and we also use expansion through the orthogonal basis in spatial range along with unique spatiotemporal functional dependence at derivation of propagation solution with supposition of weak nonlinearity.

For a first time, we develop the approach of nonlinearities integration taking into consideration all possible cases of mode coincidence that provides us with complete set of functional relations between the modes of wave train and, moreover, reveals delayed dependencies on related modes with lower order. As a result, we obtain exact expressions of time-dependent system (TDS) describing amplitude distributions over investigated frequency region (this region can also be represented by means of related wavelength range).

The importance of Goldstone mode is also shown at description of amplitudes pertaining to every separate mode. The numerical implementations show the validity of derived mathematical models and allow us to deduce some conclusions concerning temporal evolution of the autowave process and its natural frequency. The obtained results proved to be useful for development of the control methods with respect to investigated self-organizing structures.

2. Derivation of Time-Dependent System for Square Nonlinearities

Let us consider the nonlinear system of reaction-diffusion equations of FitzHugh-Nagumo (FHN) type that can be written as follows [3, 7, 8, 10, 11]: 𝜕𝑢𝑢𝜕𝑡=𝑢33𝑣+𝐷1Δ𝑢,𝜕𝑣𝜕𝑡=𝜀(𝑢𝛾𝑣+𝛽)+𝐷2Δ𝑣,(2.1) where 𝑢 and 𝑣 are kinetic variables characterizing the dynamical system (usually propagation and recovery ones), 𝐷1,2-diffusion coefficients for related variables, Δ-Laplacian, 𝜀is a small positive parameter that defines relation of temporal scales for kinetic variables, 𝛽 and 𝛾 determine excitability of the medium. Introduce the perturbations 𝑥1,𝑥2 from stationary state of (2.1), namely, 𝑢=𝑢0+𝑥1,𝑣=𝑣0+𝑥2.(2.2) Considering (2.1) and (2.2), we have the exact model consisting of nonlinear equations for perturbations (here the stationary state (𝑢0,𝑣0) is determined analogously [1, 6, 14]) 𝜕𝑥1𝑥𝜕𝑡=133𝑢0𝑥12+1𝑢02𝑥1𝑥2+𝐷1Δ𝑥1,𝜕𝑥2𝑥𝜕𝑡=𝜀1𝛾𝑥2+𝐷2Δ𝑥2.(2.3)

Let us suppose that 𝑥1,2 are small enough and apply the normal modes representation as follows: 𝑥1𝑥2=𝑀𝑚=0𝐶(𝑚)1,2𝑒𝜔𝑚𝑡cos𝑘𝑚𝑟,(2.4) where eigen wavenumber 𝑘𝑚=𝜋𝑚/𝑙,0𝑟𝑙. This decomposition is similar to Fourier series expansion, but on one variable 𝑘𝑚, while 𝜔𝑚=𝑓(𝑘𝑚) determines a dispersion curve and instability dependence, in general essentially nonlinear.

At linearized analysis of (2.1) (see [1, 6, 14]), it is impossible to express 𝐶(𝑚)1,2 from related equations since it disappears at dispersion curve determination. Here, we show that nonlinear consideration allows us to overcome this shortage. From (2.3)-(2.4), we obtain the following equation describing the amplitude distribution: 𝑀𝑚=0𝐶1(𝑚)𝑓𝑚cos𝑘𝑚𝑟+𝑢0𝜎𝑀2+13𝜎𝑀3=0,(2.5) where 𝜎𝑀=𝑀𝑚=0𝐶1(𝑚)cos𝑘𝑚𝑟,(2.6) coefficients which provide the temporal dependence are written as 𝐶1(𝑚)=𝐶1(𝑚)𝑒𝜔𝑚𝑡.(2.7) Here 𝜔𝑚=𝜔𝑚+𝑖𝜔𝑚, while 𝜔𝑚 and 𝜔𝑚 describe damping (i.e., stability properties) and oscillation processes, respectively. In (2.5), the dependence 𝑓𝑚 on spatiotemporal parameters provides complete information about linear operators involved by system (2.1), it is expressed as follows: 𝑓𝑚=𝜇𝑚+𝜔𝑚+𝜀𝜔𝑚+𝛾𝜀+𝐷2𝑘𝑚2,(2.8) where 𝜇𝑚=𝐷1𝑘𝑚21+𝑢02. In turn, (𝜎𝑀)2 and (𝜎𝑀)3 express nonlinearities on 𝐶1(𝑚) and provide all functional relations that arise at cross-modal interactions, such nonlinear interactions are studied in many applications [3, 4, 15]. It is worth noting that (2.5) itself is not sufficient for determination of 𝐶1(𝑚),𝑚=0,1,2,,𝑀, as well as for dispersion curves derivation, but we show below that by using some integral transforms we can obtain a system with complete solvability. Such system will be also numerically proved in Section 4.

Equation (2.5) also does not provide enough flexibility concerning the order of nonlinearities and does not allow us to use the “weak” level of nonlinearity that is important for retaining the constant (i.e., time-independent) level of amplitudes [5, 7] and is in good coincidence with our supposition concerning small perturbations. So, we introduce the more general form of (2.5) as follows: 𝑀𝑚=0𝐴𝑚𝑓𝑚cos𝑘𝑚𝑟+𝜀2𝑢0𝜍𝑀2+𝜀33𝜍𝑀3=0,(2.9) where 𝜍𝑀=𝑀𝑚=0𝐴𝑚cos𝑘𝑚𝐴𝑟,(2.10)𝑚=𝐴𝑚𝑒𝜔𝑚𝑡,(2.11) and 𝐴𝑚 is also a constant amplitude which pertains to the mth mode but it is suitable to different levels of nonlinearities. Evidently, in (2.5) we have that at 𝜀2=𝜀3=1, 𝐴𝑚𝐶1(𝑚), while at 𝜀2=𝜀3=0, the completely linear case on levels of deviation from stationary state takes place. This case was considered in [6, 14] and, in spite of its linearization on perturbation levels, it yet allows us the nonlinear dependence of 𝜔𝑚 on 𝑘𝑚 reflecting related frequency and dispersion curves. In this work, we suppose that 𝜀𝑖 can be in general arbitrary but small for weak nonlinearities that is used in some examples at numerical simulations.

In this paper, we confirm that admission of weak nonlinearities is significant since it allows us to apply variety of developed linear analytical methods (which were also proved experimentally for centuries) to investigation of nonlinear phenomena and derivation of new functional relations. Really, linear analysis methods can be considered as background for investigation of nonlinear phenomena, including application of matrix theory, integration and differentiation methods, and also spectral approaches (e.g., Fourier integral and series expansion). In particular, estimation of power spectrum is used widely in related modeling [16], for example, for analysis of spiral waves [7, 13] thereby reflecting some important features of investigated processes which are essential for resonant control [13, 1719]. So, we also use expansion on orthogonal basis for description of functional spatiotemporal dependencies characterizing the propagating kinetic variable.

Fourier integral transform is also widely used for reduction of complexity of partial differential equations [5, 7], but it is suitable mostly for a single pulse propagation (e.g., soliton-like solution) often followed by extrapolation of solution in spatial region [5]. But self-organization processes described by a system of reaction-diffusion equations are characterized by continuous periodic radiation represented by related wave trains, so expansion in a series is more useful than representation by integral transform. Another reason in favor of series representation is that we consider the process in restricted spatial region, in such a case that the eigenvalue spectrum is of discrete structure as it was shown in [1].

Similarly to methods couched in monograph [1], we use integration over all range of 𝑟[0,𝑙] for derivation of TDS from (2.9). In this work, we consider quadratic nonlinearities involved by (𝜍𝑀)2 followed by some generalizations to cubic nonlinearities through related matrix operations. The derivation of TDS is based on successive multiplication of (2.9) by cos𝑘𝑛𝑟, where 𝑛=0,1,2,,𝑀, followed by integration over the above-mentioned range.

At such operations, we use the filtering properties which result from orthogonality of spatial functions (see the appendix for details) and obtain from expression (2.9) the following equations for 𝐴𝑛 determination: 𝐴𝑛𝑓𝑛+𝜀2𝑢0𝐷(𝑛/2)+𝐵(𝑛)+𝐵+(𝑛)=0,(2.12) where 𝐷(𝑛/2) determines the delayed functional dependence for an even component and can be written as follows: 𝐷(𝑛/2)=𝜂(𝑛)2𝐴𝑛/22,(2.13) where the digital function 𝜂(𝑛) is determined in accordance with (A.5). The sums 𝐵(𝑛) and 𝐵+(𝑛) consist of terms which remain as nonzero ones after integration due to fulfillment of conditions |𝑘𝑖𝑘𝑗|=𝑘𝑛 and 𝑘𝑖+𝑘𝑗=𝑘𝑛, respectively, those are determined by means of the expressions: 𝐵(𝑛)=𝑀𝑛𝑖=0𝐴𝑖𝐴𝑛+𝑖,𝐵+(𝑛)=𝑠𝑖=0𝐴𝑖𝐴𝑛𝑖,(2.14) where index boundary is determined accordingly (A.6).

Evidently, both 𝐵(𝑛) and 𝐵+(𝑛) represent just the terms describing the mode interactions; the sum 𝐵(𝑛) is constructed by constant range of indices n, and in 𝐵+(𝑛) the range of indices is 𝑛2𝑖 and decreases by two samples for every subsequent pair of terms.

Equation (2.12) form the system at 𝑛=1,2,,𝑀, while for 𝑛=0 the following equation is valid, it is derived analogously above ones:𝐴0𝑓0+𝜀2𝑢02𝐴02+𝐵(0)2=0.(2.15) Derivation of TDS (2.12)–(2.15) allows us an opportunity to consider (𝑀+1) independent equations instead of one nonlinear equation (2.9) and thereby increases its solvability for determination of (𝑀+1) complex variables 𝐴𝑛. The representation of related functional dependencies in (2.12)–(2.15) can be clarified through matrix consideration along with analysis of corresponding modulo operations, we implement this in the subsequent section. This investigation is also useful for derivation of systems involving nonlinearities of higher order, while (2.12)–(2.15) express nonlinearities of the second order.

3. Matrix Representations and Modulo Operations Involved by Nonlinear Integration Procedure

Regarding representation of quadratic nonlinearities described by (2.9), the following matrix is useful: Λ𝑀𝐴=(𝑖𝐴𝑗), where 𝑖,𝑗=1,2,,𝑀. Here, we do not include the spatial dependence explicitly, but it can be easily represented with related indices accordingly (2.10), and one can conclude that then diagonal elements of Λ𝑀 form the dependence 𝐷(𝑛/2) in (2.12), while relevant pairs of mixed multiplications constitute 𝐵(𝑛) and 𝐵+(𝑛). Evidently, Λ𝑀 can be obtained by direct multiplication of the following matrices: Λ𝑀=𝑄𝑀𝑄𝑀𝑇,(3.1) where 𝑄𝑀 consists of equal elements in its columns, that is, 𝑄𝑀=𝐴1𝐴2𝐴𝑀𝐴1𝐴2𝐴𝑀𝐴1𝐴2𝐴𝑀.(3.2) As it follows from the structure of Λ𝑀, elements with equal difference of indices forming 𝐵(𝑛) constitute a pair of lateral diagonals in Λ𝑀, and by neglecting symmetric repetitions, one can conclude that the triangle matrix representing equal index differences can be obtained from (3.1) as follows: 𝑄𝑙(𝑟)=𝐴1𝐴𝑀𝐴1𝐴𝑀1𝐴1𝐴3𝐴1𝐴20𝐴2𝐴𝑀𝐴2𝐴4𝐴2𝐴3𝐴00𝑀2𝐴𝑀𝐴𝑀2𝐴𝑀1𝐴000𝑀1𝐴𝑀(),(𝑀1)(𝑀2)(2)(1)0(𝑀2)(2)(1)00(2)(1)000(1)(3.3) where the matrix 𝑄𝑙(𝑟) is constructed so that every last column is just the upper first lateral diagonal in Λ𝑀 (i.e., in reverse order). In turn, the sign of transformation () represents the transition to index difference within a pair of elements in 𝑄𝑙(𝑟). On the other hand, the matrix of index difference can be obtained using direct record of left and right indices in Λ𝑀 (𝑄left and 𝑄right, taking into account only terms upper the diagonal of Λ𝑀) for every pair of products as follows: 𝑄left=𝑄1111022200𝑀2𝑀2000𝑀1,(3.4)right=23𝑀1𝑀03𝑀1𝑀00𝑀1𝑀000𝑀,(3.5) and thus 𝑄=𝑄right𝑄left provides index difference for every product in its lateral diagonals. In turn, consideration of 𝑄+=𝑄right+𝑄left (for exploration of elements with equal sum of indices) makes sense only if (𝑄+)𝑖,𝑗𝑀, that is, for every “cross” diagonal in the upper half (i.e., subtriangle) of 𝑄+. As a result of consideration with respect to (3.3)–(3.5), one can conclude that with increase of 𝑛 in (2.12) by one sample, the quantity of elements (products) in 𝐵(𝑛) decreases by one, while that in 𝐵+(𝑛) increases by one only if 𝑛 increases by two samples (it can be shown by counting the number of terms in “cross” diagonals in the upper half of 𝑄+, with those diagonals involving equal sums of indices).

On the other hand, for investigation of mixed products properties (especially for higher-order nonlinearities) the following scheme is useful. We consider the string {𝐴𝑚}𝑀𝑚=1 and take multiplication of every 𝐴𝑖 by 𝐴𝑗 at 𝑖<𝑗 (it is said to be string by string scheme (SBS)), then all resulting terms can be obtained by consideration of corresponding rows in (3.4) and (3.5), with each pair of those containing equal number of nonzero terms.

For investigation of cubic nonlinearities, we cannot restrict our analysis to one-step matrix multiplication as for quadratic case since the more complicated tensor structure is formed; it can be represented through M square matrices which result from direct multiplication of Λ𝑀𝑠(𝑠)(𝑄𝑀), where 𝑠(𝑠)() means the s-step cyclic shift of matrix (regarding columns in our consideration); that is, if for 𝑄𝑀 the number of a column is 1𝑘𝑀, then for 𝑠(𝑠)(𝑄𝑀) it is 𝑘(𝑠)=𝑘𝑠, where the shift index 𝑠=1,2,,𝑀, and the modulo operation is defined similarly to [20] as follows: 𝑘𝑠=𝑘+𝑠𝑀,if𝑘+𝑠>𝑀,𝑘+𝑠,if𝑘+𝑠𝑀.(3.6) Evidently, if 𝑘+𝑠𝑀, then 𝑘𝑠=𝑘+𝑠(Mod𝑀), while at 𝑘+𝑠=𝑀 such operation does not coincide with the conventional one.

Thus, the columns of 𝑠(𝑠)(𝑄𝑀) can be obtained by the adding operation of indices accordingly (3.6), and therefore we have the following expressions: 𝑠(1)𝑄ind𝑀=(,(2)(3)(𝑀1)(𝑀)(1)2)(3)(𝑀1)(𝑀)(1)(2)(3)(𝑀1)(𝑀)(1)𝑠(2)𝑄ind𝑀=,(3)(4)(𝑀)(1)(2)(3)(4)(𝑀)(1)(2)(3)(4)(𝑀)(1)(2)(3.7) and so on, in (3.7) we show only indices of 𝐴𝑖, similarly to (3.3), and hence we imply the spatiotemporal dependence (2.10)-(2.11) in the represented record. Evidently, 𝑠(𝑀)(𝑄𝑀)=𝑄𝑀. It is worth noting that a single line of cubic elements 𝐴𝑖3 forms by 𝑄𝑀Λ𝑀 in its main diagonal, while mixed products including square terms such as 𝐴𝑖𝐴𝑗2 also form due to this operation, namely, 𝑀(𝑀1) terms after excluding cubic ones. The rest (𝑀1) multiplications result in 2𝑀(𝑀1) mixed square terms, those constitute two lateral diagonals and the main one as it can be shown with (3.1) and (3.7); therefore, we obtain 3𝑀(𝑀1) such terms altogether.

Alternatively, the process of matrix product derivation for cubic nonlinearities can be represented through 𝑀 successive one-step direct multiplications of Λ𝑀 by the square matrix consisting of similar elements, that is, 𝐴𝑠, where 𝑠=1,2,,𝑀. As a result, we have one cubic term and 3(𝑀1) mixed square terms at every multiplication, namely, (𝑀1) diagonal terms and 2(𝑀1) in related row and column forming two crossing lines, threefold recurrence also takes place in this scheme for the previously considered mixed square terms. At such consideration, the resulting the previous matrix preserves its symmetry that allows us to make additional conclusions which are necessary for derivation of correct expressions describing nonlinear integration analogously (2.12) derivation. It should be noted that new mixed terms such as 𝐴𝑖𝐴𝑗𝐴𝑘, where 𝑖𝑗𝑘 appear between the main diagonal and the row which involve square terms, and the common quantity of new mixed terms obtained at every multiplication decreases in accordance with the series (1/2)(𝑀1)(𝑀2), (1/2)(𝑀2)(𝑀3),,1. It can be shown using analytical approaches similarly to (3.3)–(3.5) and (3.7). As for terms involving square and cubic powers, their common quantity formed after 𝑀 multiplication can be represented through the following matrix that contains original terms only (i.e., without recurrence); it is obtained similarly to the previously mentioned SBS scheme and is of the following form: 𝑄(23)com=131,221,321,𝑀22,12232,322,𝑀2𝑀,12𝑀,22𝑀,32𝑀3,(3.8) where (13) corresponds to (𝐴1cos𝑘1𝑟)3, (1,22𝐴)1cos𝑘1𝐴𝑟(2cos𝑘2𝑟)2, and so on.

4. Generalization to Nonlinearities of the Third Order along with Numerical Simulation

Comparison of (3.8) with results of above considerations provides us with opportunity to derive the correct form of (𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟)3 for subsequent nonlinear integration, and thus we obtain from (2.9) the following expression considering cubic nonlinearities (along with those of lower order too) for amplitude distribution (analogously TDS (2.12)): 𝐴𝑛𝑓𝑛+𝜀2𝑢02𝐴0𝐴𝑛+𝐵(𝑛)+𝜀314𝐴𝑛3+𝜂1(𝑛)3𝐴𝑛/33+12𝐴𝑛𝑀𝑖=1𝑖𝑛𝐴𝑖2+12𝑀𝑖=1𝑖𝑙𝑘𝐴𝑖2𝐴𝑙𝑘+𝑀𝑖=1,𝑗=2𝑖<𝑗<𝑙𝑝𝐴𝑖𝐴𝑗𝐴𝑙𝑝+𝐴0𝐴0𝐴𝑛+𝐵(𝑛)=0,(4.1) where 𝑙𝑘,𝑙𝑝 determine index relations between equation number 𝑛 and current index values 𝑖,𝑗 for every successive sum; the total quantity of the sums for each term in (4.1) that contains 𝑙𝑘,𝑙𝑝 equals the upper value of k and P, respectively, and 𝑙𝑘,𝑙𝑝 are determined in accordance with the following expressions: 𝑙𝑘=𝑙𝑛2𝑖,if𝑘=1,2𝑖±𝑛,if𝑘=2,𝑝=𝑛±(𝑖+𝑗),if𝑝=1,𝑛±(𝑖𝑗),if𝑝=2,(𝑖+𝑗)𝑛,if𝑝=3.(4.2) In turn, the step function 𝜂1(𝑛) is expressed as follows (similarly to the above-defined function 𝜂(𝑛)): 𝜂1(𝑛)=1,if𝑛=3𝑚,0,otherwise,(4.3) where 𝑛,𝑚=1,2,,𝑀, and terms describing square nonlinearities are determined analogously above the considered ones 𝐵(𝑛)=𝐷(𝑛/2)+𝑀𝑛𝑖=1𝐴𝑖𝐴𝑛+𝑖+𝑠𝑖=1𝐴𝑖𝐴𝑛𝑖.(4.4) It is worth noting that in (4.1) the Goldstone mode 𝐴0 makes additional functional dependence and provides an opportunity to take into account the form of square nonlinearities in the third-order one that is expressed in the last term of (4.1). Thus, this mode takes particular influence on the process that is in good coincidence with results obtained in [21], where it was shown that the Goldstone mode is really of great significance in nonlinear models of wave propagation.

Evidently, the system of (4.1) provides additional delayed functional dependencies and cross-modal interactions, as well as it includes estimation of the total energy of the process. Investigation of the cubic nonlinearity is of great importance for self-organization phenomena modeling because the presence of cubic nonlinearity (or that of higher order) in a related differential equation is the necessary condition of a self-organization phenomenon origin. This fact was proved in [1] by means of comprehensive analysis of related phase curves and singular points and has many confirmations among mathematical models describing self-organization phenomena of various physical nature (see e.g., nonlinear Schrodinger equation, Ginsburg-Landau model, etc. [13, 7]).

In the present work, we reveal influence of cubic nonlinearity itself on the investigated process, for this aim we reduce the general TDS (4.1) to the system where cross-modal relations and delayed dependencies are neglected and therefore obtain the following reduced model: 𝑓𝑛+𝜀34𝐴𝑛2=0.(4.5) Considering (2.8) and (2.11) and dividing (4.5) into real and imagine parts, one can write the expressions for related trigonometric functions as follows: cos2𝜔𝑛𝑡=4𝛼𝑑𝜀3𝐴𝑛2𝜇𝑚+𝜔𝑚+̃𝜀𝑛𝛾𝑛,sin2𝜔𝑛𝑡=4𝛼𝑑𝜔𝑛𝜀3𝐴𝑛21̃𝜀𝑛,(4.6) where characteristics for 𝑛th mode are as follows: 𝛾𝑛=𝜔𝑛+𝛾𝜀+𝐷2𝑘𝑛2,̃𝜀𝑛=𝜀𝜔2𝑛+𝛾n2,(4.7) and dissipative term displays the temporal dependence and simultaneously reflects changes with the mode order as 𝛼𝑑=𝑒2𝜔𝑚𝑡. After some transformations, we obtain from (4.6) the expression for amplitude estimation of the form 𝐴𝑛=2𝛼𝑑𝜀34𝜇𝑛+𝜔𝑛+̃𝜀𝑛𝛾𝑛2+𝜔2𝑛̃𝜀𝑛12.(4.8) At numerical simulations, we use the dispersion and dissipation dependencies obtained before for linearized analysis of (2.1), namely, (see [6]) 𝜔𝑛=12𝜇4𝜀𝑛𝛾𝜀2,𝜔𝑚𝜇=𝑚+𝛾𝜀2.(4.9) So, based on fundamental results and considerations obtained by leading authors in the field of spiral waves modeling [8, 10, 11], we also suggest for numerical considerations that 𝐷1=1, 𝐷2=0. The parameters of excitation are chosen to be as 𝛾=0.5 and 𝛽=0.7, while propagating kinetic variable at stationary state is determined as 𝑢0=1.0328 [6, 14] and temporal scale parameters 𝜀=0.02,𝜀3=0.01,𝑙=100.

Results of numerical calculations of dependence (4.8) are shown in Figure 1; they show that the cubic nonlinearity described previously provides the growth of amplitude; this growth attains its peak in the region of modes with maximal order which simultaneously possess temporally oscillating property, as well as lend the system with maximal temporal stability. This allows us to conclude that the natural frequency of investigated autowave process is located within the modes of the largest wavenumber where frequency (and phase velocity) is minimal at high stability (both convective and temporal [6]) that confirms the conjecture deduced in [6, 14] from linearized analysis. This is also in good coincidence with results obtained after numerical simulations of spiral wave evolution [11, 19], where the natural frequency was shown to be of very small order of magnitude that is characteristic of biological active media. At the same time, it is worth noting that experimental detection of such low frequency can cause some problems since such signals are of low propagating ability.

Figure 1: The dependencies of the amplitude level 𝐴𝑚 on the mode order reflecting the influence of cubic nonlinearity on the self-organization process evolution, (a) corresponds to 𝑡=5, (b) is related to 𝑡=50, thus the obtained curves represent an asymptotic estimation of amplitude distribution.

5. Evaluation of Obtained Results Regarding General Solvability of the Model

In this section, we also consider the conservative approach that is a paradigm of classic physics of wave processes based on Maxwell equations [5] and explore the problem of solvability of the above-obtained TDS. In a case of direct integration of TDS (4.1) or (2.12) on some large enough temporal period 𝑇 (similarly to [1]), one should take into consideration that related temporal functions 𝑒𝜔𝑚𝑡 are essentially nonorthogonal. In other words, the order properties which are intrinsic due to orthogonality are lost here, in difference from spatial functions, due to nonlinearity of dispersion dependence 𝜔(𝑘). For example, for conservative systems where 𝜔𝑚=0 and using related operations of multiplication and integration as above, we obtain for a term of TDS the following expression: 𝑇/2𝑇/2𝐴Σ𝑒𝑖𝜔𝜎𝑡𝐴𝑑𝑡=Σ𝜔𝜎𝜋sin𝜔𝜎𝜔,(5.1) where the variable 𝜔𝜎 includes all frequencies involved by related eigenvalues of the term, that is, 𝜔𝜎=𝜔𝑛1+𝜔𝑛2++𝜔𝑛𝜎,(5.2) and similarly 𝐴Σ includes amplitude values related to the term as follows: 𝐴Σ=𝐴𝑛1𝐴𝑛2𝐴𝑛Σ.(5.3) Evidently, operation of exact integration (5.1) leads to transformation of TDS into a system of transcendental equations which is cumbersome both for analytical and numerical analyses; the similar result is obtained also for systems with dissipation. So, in this paper we use some approximation concerning orthogonality in temporal region, namely, we suppose that there exists such small enough 𝜔=2𝜋/𝑇 that for every eigenvalue 𝜔𝑛𝑞𝑛𝜔 is valid, where 𝑞𝑛 is an integer. Such approximation seems to be suitable for conservative systems only, but below we show its applicability also for dissipative ones, but the order of such quasiorthogonality is quite different in comparison with regular orthogonality of spatial functions. For conservative systems consideration, the normal mode analysis is reduced to representation of perturbations as follows: 𝑥1𝑀𝑚=0𝐵(𝑚)𝑒𝑖𝜑𝑚𝑡, where 𝜑𝑚=𝜔𝑛𝑡𝑘𝑟, and from (4.8) at 𝜔𝑚=0 we have dependence shown in Figure 2, it confirms that the amplitude level increases with the mode order, but the accuracy of conservative approach is restricted with “long-wave” consideration only [5]. On the other hand, obtained dependence complies with the conjecture concerning natural frequency couched previously.

Figure 2: The dependence of the amplitude distribution for conservative system confirms reliable growth for initial modes.

Taking into account the above quasiorthogonal approximation, we have for dissipative system that provides more exact consideration than conservative one 𝑇0𝐴𝑛𝐴𝑑𝑡=𝑛𝜔𝑛𝜔exp2𝜋𝑛𝜔𝐴1𝑛𝜔𝑛,(5.4) because exp(2𝜋𝑖(𝜔𝑛)/(𝜔))1 and |𝜔𝑛|>𝜔,𝜔𝑛<0. Therefore, similarly to (5.4), one can obtain from the TDS (2.12) the following system: 𝑓𝑛𝐴𝑛𝜔𝑛+𝜀2𝑢0𝜂(𝑛)2𝜔𝑛/2𝐴𝑛/22+𝑀𝑛𝑖=0𝐴𝑖𝐴𝑛+𝑖𝜔𝑖+𝜔𝑛+𝑖+𝑠𝑖=0𝐴𝑖𝐴𝑛𝑖𝜔𝑖+𝜔𝑛𝑖=0.(5.5) If we take TDS (2.12) at some fixed 𝑡0 (e.g., 𝑡0=0, similarly to determination of coefficients in [22]), then the obtained equation along with (5.5) forms the system with sufficient solvability for derivation of |𝐴𝑛|.

For estimation of obtained results, it is worth noting that in the present work we considered the behavior of the nonlinear system near its stationary (fixed) point that is similar to approaches developed in [1, 3]. But, in difference from [1], where the linearized operator for nonlinear Brusselator was analysed followed by bifurcation and instability analysis (preferably near the stationary state), the aim of our work is to derive the nonlinear system (taking into account weak nonlinearities) for estimation of amplitude distribution of the wave train. The derived TDS also includes information about instabilities (expressed through 𝜔𝑚), as well as dispersion dependency (𝜔𝑚), but detailed analysis of those is beyond the frames of our work (that leads to transcendental equations), while such analysis for linearized case of reaction-diffusion equations was implemented in our previous papers (see [6, 14]). The considered FHN system contains the essential information about physical processes in the active media under investigation and allows the derivation of models of various physical structures, such as spiral waves and three-dimensional scrolls as it was shown in [7, 8, 11].

As to limitation that arises from exponential representation, it is worth noting that it cannot really catch all features of physical behavior, but only near the stationary state. Nevertheless, in accordance with conclusions proved in the monograph [1] (where such representation is also used), it should be noted that the system with self-organization tends to stationary state (or nearly that) with elapsing the time, so the method developed allows essential properties of the system. Moreover, the method developed in this paper includes features of spatiotemporal nonlinear interactions as one can conclude from (2.9).

In some works, the decomposition in spatial region is also used through the multiple time series [23]. Then, several authors showed that such approach (exposing also nonadiabatic phenomena) results in locking effect which arises from the interaction of the large-scale envelope of the kinetic variable with the small scale underlying the spatial periodic solution (considered also in [24, 25]). At the same time, the main difference of our method is that the derived TDS does not include a spatial dependence in an explicit form due to the nonlinear integration in spatial region which gives the TDS for amplitudes |𝐴𝑛| determination and simultaneously retains the temporal dependence. So, the previously mentioned locking effect cannot appear in our model since the lack of explicit spatial dependence and nonlinear interaction between different modes exposes through including different constant levels of amplitudes (pertained to every separate mode) and related temporal dependencies. Due to nonlinearities involved, the amplitude distribution describes the self-organization structure under investigation.

The main aim of the present paper is investigation of self-organization phenomena in nonlinear systems (namely, the model that can describe spiral waves). In accordance with the theorem proved in [1], such phenomena can take place only in the systems which are described by equations with nonlinearity not less than third order, so FHN model fits for such investigation. Again, the method developed in this paper is designated for the systems describing spiral waves distribution, so more simple models might fail for modeling such self-organization phenomena (as a rule, those model a wave front of plane form [25]), and thus those were not considered, though application of the method developed is also possible for them.

The system with four spatiotemporal parameters is possible for three-dimensional consideration. But we considered radial wave propagation only since the spiral wave (far enough from the center where it is usually measured in biological systems) can be regarded as usual pacemaker (followed by a target wave) and angular changes can be neglected. Nevertheless, if we use also angular considerations, the sense of the method does not change.

6. Conclusions

Thus, the TDS obtained after integration of nonlinearities is shown to be sufficient for derivation of amplitude distribution of investigated wave train and, moreover, provides the opportunity to increase the accuracy of dispersion relation. The related matrix and modulo operations are also considered; those confirm validity of TDS derivation method.

The derived model allows us to find the mode with maximal amplitude that is important for resonance control of the related self-organization process [13, 17]. Again, exploration of this mode explains why the natural frequency of investigated synergetic process is of small level (preferably in biological systems that are modeled with reaction-diffusion equations [3]) and determines the measured frequency at which the wave train propagates. So, the derived mathematical models in this paper can be applicable to exploration of nonlinear processes of different physical nature.


Filtering Properties at Nonlinear Integration on Orthogonal Basis

For linear terms, this operation is similar to that used in Fourier series consideration and yet widely applied in digital signal processing 𝑙0cos𝑘𝑛𝑟𝑀𝑚=0𝐴𝑚𝑓𝑚cos𝑘𝑚𝑟𝑙𝑑𝑟=2𝐴𝑛𝑓𝑛,(A.1) where (A.1) is valid for 𝑛=1,2,,𝑀.

For integration of (𝜍𝑀)2, let us rewrite it with separation of Goldstone terms as follows: 𝜍𝑀2=𝐴02𝐴+20𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟+𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟2,(A.2) where, in turn, considering mixed multiplications which describe cross-modal interactions, one can write 𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟2=𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟2+2𝑀𝑖,𝑗=1𝑖𝑗𝐴𝑖𝐴𝑗cos𝑘𝑖𝑟cos𝑘𝑗𝑟,(A.3) and integration similarly to (A.1) yields the following expressions for nonlinear terms: 𝑙0cos𝑘𝑛𝑟𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟2𝐴𝑑𝑟=𝑙𝜂(𝑛)𝑛/222+12𝑀𝑛𝑖=1𝑛<𝑀𝐴𝑖𝐴𝑖+𝑛+𝑠𝑗=1𝑛>2𝐴𝑗𝐴𝑛𝑗,(A.4) where 𝜂(𝑛) indicates that the consideration is implemented for the even coefficient 𝜂(𝑛)=1,if𝑛=2𝑚,0,otherwise,(A.5) and the bound for the sum with change of contracting range is expressed as follows: 𝑛𝑠=2𝑟1,(A.6) where []𝑟 means rounding to +, and for the second term of (A.2) the following is valid: 2𝐴0𝑙0cos𝑘𝑛𝑟𝑀𝑚=1𝐴𝑚cos𝑘𝑚𝑟𝐴𝑑𝑟=0𝐴𝑛,(A.7) while integration of the first term of (A.2) provides zero. So, (A.4) and (A.7) contain all nonzero terms which remain after nonlinear second-order considerations with subsequent integration, those terms appear only if the following is valid concerning (A.3): |𝑖±𝑗|=𝑛, on one hand, and 𝑚/2=𝑛, on the other hand. Let us note that (A.4)–(A.7) is valid for 𝑛=1,2,,𝑀, and the related expression for 𝑛=0 is derived analogously. It differs only with some coefficients.


  1. G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems, Wiley-Interscience, New York, NY, USA, 1977. View at Zentralblatt MATH
  2. R. A. Barrio and C. Varea, “Non-Linear Systems,” Physica A, vol. 372, no. 2, pp. 210–223, 2006. View at Publisher · View at Google Scholar
  3. M. C. Cross and P. C. Hohenberg, “Pattern formation outside of equilibrium,” Reviews of Modern Physics, vol. 65, no. 3, pp. 851–1107, 1993. View at Publisher · View at Google Scholar
  4. M. Hirota, T. Tatsuno, and Z. Yoshida, “Resonance between continuous spectra: secular behavior of Alfvén waves in a flowing plasma,” Physics of Plasmas, vol. 12, no. 1, p. 012107, 11, 2005. View at Publisher · View at Google Scholar
  5. R. K. Dodd, J. C. Eilbeck, J. D. Gibbon, and H. C. Morris, Solitons and Nonlinear Wave E quations, Academic Press, London, UK, 1982.
  6. V. F. Dailyudenko, “Instabilities on normal modes of perturbation solution for a nonlinear single-diffusive system,” International Journal of Nonlinear Science, vol. 11, no. 2, pp. 143–152, 2011. View at Zentralblatt MATH
  7. H. Mori and Y. Kuramoto, Dissipative Structures and Chaos, Springer-Verlag, New York, NY, USA, 1998. View at Publisher · View at Google Scholar
  8. A. T. Winfree, “Alternative stable rotors in an excitable medium,” Physica D, vol. 49, no. 1-2, pp. 125–140, 1991.
  9. M. Bar and L. Brusch, “Breakup of spiral waves caused by radial dynamics: Eckhause and finite wavenumber instabilities,” New Journal of Physics, vol. 6, pp. 1–22, 2004. View at Publisher · View at Google Scholar
  10. J. J. Tyson and J. P. Keener, “Spiral waves in a model of myocardium,” Physica D, vol. 29, no. 1-2, pp. 215–222, 1987.
  11. H. Sakaguchi and T. Fujimoto, “Forced entrainment and elimination of spiral waves for the FitzHugh-Nagumo equation,” Progress of Theoretical Physics, vol. 108, no. 2, pp. 241–252, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. B. Sandstede and A. Scheel, “Absolute versus convective instability of spiral waves,” Physical Review E. Statistical, Nonlinear, and Soft Matter Physics, vol. 62, no. 6, part A, pp. 7708–7714, 2000. View at Publisher · View at Google Scholar
  13. H. Zhang, B. Hu, and G. Hu, “Suppression of spiral waves and spatiotemporal chaos by generating target waves in excitable media,” Physical Review E, vol. 68, no. 2, part 2, p. 026134, 2003.
  14. V. F. Dailyudenko, “Analytical models of instabilities on normal modes of perturbation solution for a parabolic single-diffusive system,” Advances and Applications in Mathematical Sciences, vol. 8, no. 2, pp. 141–166, 2011. View at Zentralblatt MATH
  15. H. Ji, J. Qiu, K. Zhu, and A. Badel, “Two-mode vibration control of a beam using nonlinear synchronized switching damping based on maximization of converted energy,” Journal of Sound and Vibration, vol. 329, no. 14, pp. 2751–2767, 2010. View at Publisher · View at Google Scholar
  16. H. G. Schuster, Deterministic Chaos, VCH Verlagsgesellschaft mbH, Weinheim, Germany, 2nd edition, 1988.
  17. A. P. Munuzuri, M. Gomez-Gesteira, V. Perez-Munuzuri, V. V. Krinsky, and V. Peres-Villar, “Parametric resonance of a vortex in an active medium,” Physical Review E, vol. 50, no. 5, pp. 4258–4261, 1994. View at Publisher · View at Google Scholar
  18. S. A. Vysotsky, R. V. Cheremin, and A. Loskutov, “Suppression of spatio-temporal chaos in simple models of re-entrant fibrillations,” Journal of Physics, vol. 23, pp. 202–209, 2005.
  19. K. Martinez, A. L. Lin, R. Kharrazian, X. Sailer, and H. L. Swinney, “Resonance in periodically inhibited reaction-diffusion systems,” Physica D, vol. 168-169, pp. 1–9, 2002. View at Publisher · View at Google Scholar
  20. V. F. Dailyudenko, “The integrated and local estimations of instability for a class of autonomous delay systems,” Chaos, Solitons & Fractals, vol. 30, no. 3, pp. 759–768, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. E. Simbawa, P. C. Matthews, and S. M. Cox, “Nikolaevskiy equation with dispersion,” Physical Review E, vol. 81, no. 3, pp. 036220.1–03622011, 2010.
  22. A. N. Tikhonov and A. A. Samarsky, Equations of Mathematical Physics, Nauka, Moscow, Russia, 1977.
  23. D. Bensimon, B. I. Shraiman, and V. Croquette, “Nonadiabatic effect in convection,” Physical Review A, vol. 38, no. 10, pp. 5461–5464, 1988. View at Publisher · View at Google Scholar
  24. M. G. Clerc, C. Falcon, and E. Tirapegui, “Additive noise induced front propagation,” Physical Review Letters, vol. 94, no. 14, p. 148302, 2005. View at Publisher · View at Google Scholar
  25. M. G. Clerc and C. Falcon, “Localized patterns and hole solutions in one-dimensional extended systems,” Physica A, vol. 356, p. 48, 2005. View at Publisher · View at Google Scholar