The Nonlinear Analysis of Perturbation Solution for a Parabolic Differential System

Abstract

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 [1–6], 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 [8–10]) 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, 11–13]). 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) [8–10]. 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]: 𝜕𝑢𝜕𝑡=𝑢−𝑢3³−𝑣+𝐷¹𝜕𝑣Δ𝑢,𝜕𝑡=𝜀(𝑢−𝛾𝑣+𝛽)+𝐷²Δ𝑣,(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 𝑥₁,𝑥₂ from stationary state of (2.1), namely, 𝑢=𝑢₀+𝑥₁,𝑣=𝑣₀+𝑥₂.(2.2) Considering (2.1) and (2.2), we have the exact model consisting of nonlinear equations for perturbations (here the stationary state (𝑢₀,𝑣₀) is determined analogously [1, 6, 14]) 𝜕𝑥₁𝜕𝑡=−𝑥3¹³−𝑢⁰𝜕𝑥₂𝑥₁₂+1−𝑢₀₂𝑥₁−𝑥₂+𝐷₁Δ𝑥₁,𝜕𝑡=𝜀𝑥¹−𝛾𝑥₂+𝐷₂Δ𝑥₂.(2.3)

Let us suppose that 𝑥_1,2 are small enough and apply the normal modes representation as follows: ⎛⎜⎜⎝𝑥₁𝑥₂⎞⎟⎟⎠=^𝑀_𝑚₌₀𝐶^(𝑚)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𝑘_𝑚𝑟+𝑢₀𝜎_𝑀²+13𝜎^𝑀³=0,(2.5) where 𝜎_𝑀=^𝑀_𝑚₌₀𝐶^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.8) where 𝜇_𝑚=𝐷₁𝑘_𝑚2−1+𝑢₀₂. In turn, (𝜎_𝑀)² and (𝜎_𝑀)³ 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𝑘_𝑚𝑟+𝜀₂𝑢₀𝜍_𝑀²+𝜀³3𝜍^𝑀³=0,(2.9) where 𝜍_𝑀=^𝑀_𝑚₌₀𝐴𝐴_𝑚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 𝜀₂=𝜀₃=1, 𝐴_𝑚≡𝐶^1(𝑚), while at 𝜀₂=𝜀₃=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, 17–19]. 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 (𝜍_𝑀)² 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.12) where 𝐷^(𝑛/2) determines the delayed functional dependence for an even component and can be written as follows: 𝐷^(𝑛/2)=𝜂(𝑛)2𝐴^𝑛/2²,(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:𝐴₀𝑓₀+𝜀²𝑢₀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, 𝑄_𝑀=⎛⎜⎜⎜⎜⎜⎜⎝𝐴¹𝐴₂⋯𝐴_𝑀𝐴₁𝐴₂⋯𝐴_𝑀𝐴⋯⋯⋯⋯₁𝐴₂⋯𝐴_𝑀⎞⎟⎟⎟⎟⎟⎟⎠.(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⋯𝐴₁𝐴₃𝐴₁𝐴₂0𝐴₂𝐴_𝑀⋯𝐴₂𝐴₄𝐴₂𝐴₃⋯0⋯0⋯⋯𝐴⋯_𝑀−2⋯𝐴_𝑀𝐴_𝑀−2𝐴_𝑀−100⋯0𝐴_𝑀−1𝐴_𝑀⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠₍₋₎−−→⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝(𝑀⋯000−1)((𝑀𝑀⋯00−−22))⋯⋯⋯⋯⋯(((⋯2220)))((((⋯1111))))⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(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=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝⋯1000⋯1200⋯⋯⋯⋯⋯𝑀⋯120−2𝑀𝑀⋯12−−21⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(3.4)_right=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝⋯2000⋯3300⋯⋯⋯⋯⋯𝑀𝑀𝑀⋯0−−−111𝑀𝑀𝑀𝑀⋯⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(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: 𝑘⊕𝑠=𝑘𝑘++𝑠𝑠,−𝑀,iiff𝑘𝑘++𝑠𝑠≤>𝑀𝑀.,(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: 𝑠ℎ𝑠ℎ⁽¹⁾𝑄^ind_𝑀(2)=⎛⎜⎜⎜⎜⎜⎜⎝(((⋯222)))(((⋯333)))⋯⋯⋯⋯(((𝑀𝑀𝑀⋯−−−111)))(((𝑀𝑀𝑀⋯)))(((⋯111)))⎞⎟⎟⎟⎟⎟⎟⎠,𝑄^ind_𝑀=⎛⎜⎜⎜⎜⎜⎜⎝(((⋯333)))(((⋯444)))⋯⋯⋯⋯(((𝑀𝑀𝑀⋯)))(((⋯111)))(((⋯222)))⎞⎟⎟⎟⎟⎟⎟⎠,(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: 𝑄_(2−3)com=⎛⎜⎜⎜⎜⎜⎜⎝1³2,1₂1,2²1,3²⋯1,𝑀²𝑀⋯,1₂⋯2³2⋯,3²⋯⋯2,⋯𝑀²𝑀,2²𝑀,3²⋯𝑀³⎞⎟⎟⎟⎟⎟⎟⎠,(3.8) where (1³) corresponds to (𝐴₁cos𝑘₁𝑟)³, (1,2²)−𝐴₁cos𝑘₁𝑟(𝐴₂cos𝑘₂𝑟)², 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𝑘_𝑚𝑟)³ 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𝐴⁰𝐴_𝑛+𝐵^(𝑛)+𝜀₃⎧⎪⎨⎪⎩14𝐴^𝑛3+𝜂¹(𝑛)3𝐴^𝑛/3³+21⎡⎢⎢⎢⎣𝐴^𝑛𝑀_{𝑖=1𝑖≠𝑛}𝐴_𝑖2+12^𝑀_{𝑖=1𝑖≠𝑙𝑘}𝐴_𝑖2𝐴_𝑙𝑘+^𝑀_𝑖=1,_{𝑗=2𝑖<𝑗<𝑙𝑝}𝐴_𝑖𝐴_𝑗𝐴_𝑙𝑝+𝐴₀⎤⎥⎥⎥⎦𝐴⁰𝐴_𝑛+𝐵^(𝑛)⎫⎪⎬⎪⎭=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𝑖−±2𝑛𝑖,,iiff𝑘𝑘==12,,_𝑝=⎧⎪⎨⎪⎩𝑛𝑛(𝑖±±+𝑗((𝑖𝑖)+−−𝑗𝑗𝑛)),,,iiifff𝑝𝑝𝑝===123,,.(4.2) In turn, the step function 𝜂₁(𝑛) is expressed as follows (similarly to the above-defined function 𝜂(𝑛)): 𝜂₁(𝑛)=10,,iofth𝑛e=rw3is𝑚e,,(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 𝐴₀ 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. [1–3, 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: 𝑓_𝑛+𝜀³4𝐴^𝑛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𝛼^𝑑𝜀₃𝐴_𝑛2sin2𝜔_𝑛𝜇_𝑚+𝜔_𝑚+𝜀̃_𝑛𝛾_𝑛,𝑡=−4𝛼^𝑑𝜔_𝑛𝜀₃𝐴_𝑛21−𝜀̃_𝑛,(4.6) where characteristics for 𝑛th mode are as follows: 𝛾_𝑛=𝜀𝜔̃_𝑛𝑛+𝛾𝜀+𝐷₂𝑘_𝑛2,=𝜀𝜔^2𝑛+𝛾_n²,(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𝛼_𝑑𝜀₃⁴𝜇_𝑛+𝜔_𝑛+𝜀̃_𝑛𝛾_𝑛²+𝜔^2𝑛𝜀̃_𝑛−1².(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.(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, 𝐷₂=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 𝑢₀=−1.0328 [6, 14] and temporal scale parameters 𝜀=0.02,𝜀₃=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.

(a)

(b)

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: 𝑥₁≅∑^𝑀𝑚=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 ^𝑇∗₀𝐴_𝑛𝑑𝑡=𝐴^𝑛𝜔_𝑛exp2𝜋𝜔𝜔^𝑛_∗−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𝜂𝜔(𝑛)_𝑛/2𝐴^𝑛/2²+^𝑀−𝑛_𝑖=0𝐴_𝑖𝐴_𝑛+𝑖𝜔_𝑖+𝜔_𝑛+𝑖+^𝑠_𝑖=0𝐴_𝑖𝐴_𝑛−𝑖𝜔_𝑖+𝜔_𝑛−𝑖=0.(5.5) If we take TDS (2.12) at some fixed 𝑡₀ (e.g., 𝑡₀=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.

Appendix

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𝑘_𝑛𝑟^𝑀_𝑚₌₀𝐴_𝑚𝑓_𝑚cos𝑘_𝑚𝑟𝑑𝑟=𝑙2𝐴^𝑛𝑓_𝑛,(A.1) where (A.1) is valid for 𝑛=1,2,…,𝑀.

For integration of (𝜍_𝑀)², let us rewrite it with separation of Goldstone terms as follows: 𝜍_𝑀²=𝐴₀₂+2𝐴^0𝑀_𝑚=1𝐴_𝑚cos𝑘_𝑚𝑟+^𝑀_𝑚=1𝐴_𝑚cos𝑘_𝑚𝑟²,(A.2) where, in turn, considering mixed multiplications which describe cross-modal interactions, one can write ^𝑀_𝑚=1𝐴_𝑚cos𝑘_𝑚𝑟²=^𝑀_𝑚₌₁𝐴^𝑚cos𝑘_𝑚𝑟²+2^𝑀_𝑖_{,𝑗=1𝑖≠𝑗}𝐴_𝑖𝐴_𝑗cos𝑘_𝑖𝑟cos𝑘_𝑗𝑟,(A.3) and integration similarly to (A.1) yields the following expressions for nonlinear terms: ^𝑙0cos𝑘_𝑛𝑟^𝑀_𝑚₌₁𝐴_𝑚cos𝑘_𝑚𝑟²𝑑𝑟=𝑙⎡⎢⎢⎢⎣𝜂(𝑛)𝐴^𝑛/22²+12⎛⎜⎜⎜⎝^𝑀−𝑛_{𝑖=1𝑛<𝑀}𝐴_𝑖𝐴_𝑖+𝑛+^𝑠_𝑗=1𝑛>2𝐴_𝑗𝐴_𝑛−𝑗⎞⎟⎟⎟⎠⎤⎥⎥⎥⎦,(A.4) where 𝜂(𝑛) indicates that the consideration is implemented for the even coefficient 𝜂(𝑛)=10,,iofth𝑛e=rw2is𝑚e,,(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𝐴₀^𝑙0cos𝑘_𝑛𝑟^𝑀_𝑚₌₁𝐴_𝑚cos𝑘_𝑚𝑟𝑑𝑟=𝐴₀𝐴_𝑛,(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.