The Nonlinear Analysis of Perturbation Solution for a Parabolic Differential System

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 timedependent 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.


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 International Journal of Mathematics and Mathematical Sciences 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.

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 : where u and v are kinetic variables characterizing the dynamical system usually propagation and recovery ones , D 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 x 1 , x 2 from stationary state of 2.1 , namely, Considering 2.1 and 2.2 , we have the exact model consisting of nonlinear equations for perturbations here the stationary state u 0 , v 0 is determined analogously 1, 6, 14

2.3
International Journal of Mathematics and Mathematical Sciences 3 Let us suppose that x 1,2 are small enough and apply the normal modes representation as follows: where eigen wavenumber k m πm/l, 0 ≤ r ≤ l.This decomposition is similar to Fourier series expansion, but on one variable k m , while ω m f k m 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 C m 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: where coefficients which provide the temporal dependence are written as

2.7
Here ω m ω m i ω m , while ω m and ω m describe damping i.e., stability properties and oscillation processes, respectively.In 2.5 , the dependence f m on spatiotemporal parameters provides complete information about linear operators involved by system 2.1 , it is expressed as follows: where 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 C m 1 , m 0, 1, 2, . . ., M, 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 International Journal of Mathematics and Mathematical Sciences coincidence with our supposition concerning small perturbations.So, we introduce the more general form of 2.5 as follows: where 10 and A m 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 , 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 ω m on k m reflecting related frequency and dispersion curves.In this work, we suppose that ε i 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 r ∈ 0, l for derivation of TDS from 2.9 .In this work, we consider quadratic nonlinearities involved by ς M 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 cosk n r, where n 0, 1, 2, . . ., M, 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 A n determination: where D n/2 determines the delayed functional dependence for an even component and can be written as follows: where the digital function η n is determined in accordance with A.5 .The sums B n − and B n consist of terms which remain as nonzero ones after integration due to fulfillment of conditions |k i − k j | k n and k i k j k n , respectively, those are determined by means of the expressions: where index boundary is determined accordingly A. n − 2i and decreases by two samples for every subsequent pair of terms.Equation 2.12 form the system at n 1, 2, . . ., M, while for n 0 the following equation is valid, it is derived analogously above ones: Derivation of TDS 2.12 -2.15 allows us an opportunity to consider M 1 independent equations instead of one nonlinear equation 2.9 and thereby increases its solvability for determination of M 1 complex variables A n .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.

Matrix Representations and Modulo Operations Involved by Nonlinear Integration Procedure
Regarding representation of quadratic nonlinearities described by 2.9 , the following matrix is useful: Λ M A i A j , where i, j 1, 2, . . ., M. 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 Λ M form the dependence D n/2 in 2.12 , while relevant pairs of mixed multiplications constitute B n − and B n .Evidently, Λ M can be obtained by direct multiplication of the following matrices: where Q M consists of equal elements in its columns, that is, As it follows from the structure of Λ M , elements with equal difference of indices forming B n − constitute a pair of lateral diagonals in Λ M , and by neglecting symmetric repetitions, one can conclude that the triangle matrix representing equal index differences can be obtained from 3.1 as follows: where the matrix Q l r is constructed so that every last column is just the upper first lateral diagonal in Λ M i.e., in reverse order .In turn, the sign of transformation − − − → represents the transition to index difference within a pair of elements in Q l r .On the other hand, the matrix of index difference can be obtained using direct record of left and right indices in Λ M Q left and Q right , taking into account only terms upper the diagonal of Λ M for every pair of products as follows: and thus Q − Q right −Q left provides index difference for every product in its lateral diagonals.
In turn, consideration of Q Q right Q left for exploration of elements with equal sum of indices makes sense only if Q i,j ≤ M, that is, for every "cross" diagonal in the upper half i.e., subtriangle of Q .As a result of consideration with respect to 3.3 -3.5 , one can conclude that with increase of n in 2.12 by one sample, the quantity of elements products in B n − decreases by one, while that in B n increases by one only if n increases by two samples it can be shown by counting the number of terms in "cross" diagonals in the upper half of Q , 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 and take multiplication of every A i by A j at i < j 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 Λ M ⊗ sh s Q M , where sh s • means the s-step cyclic shift of matrix regarding columns in our consideration ; that is, if for Q M the number of a column is where the shift index s 1, 2, . . ., M, and the modulo operation is defined similarly to 20 as follows: Evidently, if k s / M, then k ⊕ s k s Mod M , while at k s M such operation does not coincide with the conventional one.

International Journal of Mathematics and Mathematical Sciences
Thus, the columns of sh s Q M can be obtained by the adding operation of indices accordingly 3.6 , and therefore we have the following expressions: and so on, in 3.7 we show only indices of A i , similarly to 3.3 , and hence we imply the spatiotemporal dependence 2.10 -2.11 in the represented record.Evidently, sh M Q M Q M .It is worth noting that a single line of cubic elements A 3 i forms by Q M ⊗ Λ M in its main diagonal, while mixed products including square terms such as A i A 2 j also form due to this operation, namely, M M−1 terms after excluding cubic ones.The rest M−1 multiplications result in 2M M − 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 3M M − 1 such terms altogether.
Alternatively, the process of matrix product derivation for cubic nonlinearities can be represented through M successive one-step direct multiplications of Λ M by the square matrix consisting of similar elements, that is, A s , where s 1, 2, . . ., M. As a result, we have one cubic term and 3 M − 1 mixed square terms at every multiplication, namely, M − 1 diagonal terms and 2 M − 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 A i A j A k , where i / j / k 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 M − 1 M − 2 , 1/2 M − 2 M − 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 M 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: where 1 3 corresponds to A 1 cosk 1 r 3 , 1, 2 2 − A 1 cosk 1 r A 2 cosk 2 r 2 , and so on.

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 M m 1 A m cosk m r 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 : where l k , l p determine index relations between equation number n and current index values i, j for every successive sum; the total quantity of the sums for each term in 4.1 that contains l k , l p equals the upper value of k and P, respectively, and l k , l p are determined in accordance with the following expressions:

4.2
In turn, the step function η 1 n is expressed as follows similarly to the above-defined function η n : where n, m 1, 2, . . ., M, and terms describing square nonlinearities are determined analogously above the considered ones 4.4

International Journal of Mathematics and Mathematical Sciences
It is worth noting that in 4.1 the Goldstone mode A 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. 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: 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 where characteristics for nth mode are as follows:

4.7
and dissipative term displays the temporal dependence and simultaneously reflects changes with the mode order as α d e −2ω m t .After some transformations, we obtain from 4.6 the expression for amplitude estimation of the form International Journal of Mathematics and Mathematical Sciences 11 At numerical simulations, we use the dispersion and dissipation dependencies obtained before for linearized analysis of 2.1 , namely, see 6

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 D 1 1, D 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 u 0 −1.0328 6, 14 and temporal scale parameters ε 0.02, ε 3 0.01, l 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.

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 T * similarly to 1 , one should take into consideration that related temporal functions e ω m t 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 ω k .For example, for conservative systems where ω m 0 and using related operations of multiplication and integration as above, we obtain for a term of TDS the following expression: where the variable ω σ includes all frequencies involved by related eigenvalues of the term, that is, and similarly A Σ includes amplitude values related to the term as follows: 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π/T * that for every eigenvalue ω n ∼ q n ω * is valid, where q n is an integer.Such approximation seems to be suitable for conservative systems International Journal of Mathematics and Mathematical Sciences 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: , where ϕ m ω n t − kr, and from 4.8 at ω m 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.
Taking into account the above quasiorthogonal approximation, we have for dissipative system that provides more exact consideration than conservative one because exp 2πi ω n / ω * ∼ 1 and |ω n | > ω * , ω n < 0. Therefore, similarly to 5.4 , one can obtain from the TDS 2.12 the following system: If we take TDS 2.12 at some fixed t 0 e.g., t 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 |A n |.
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 ω m , as well as dispersion dependency ω m , 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 |A n | 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.

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 reactiondiffusion 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 A m cosk m r dr A 0 A n , 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 : |i ± j| n, on one hand, and m/2 n, on the other hand.Let us note that A.4 -A.7 is valid for n 1, 2, . . ., M, and the related expression for n 0 is derived analogously.It differs only with some coefficients.
6 .Evidently, both B n − and B n represent just the terms describing the mode interactions; the sum B n − is constructed by constant range of indices n, and in B n the range of indices is

bFigure 1 :
Figure 1:The dependencies of the amplitude level A m on the mode order reflecting the influence of cubic nonlinearity on the self-organization process evolution, a corresponds to t 5, b is related to t 50, thus the obtained curves represent an asymptotic estimation of amplitude distribution.

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

1 A m cosk m r M m 1 A m cosk m r 2 , A. 2 where 1 A m cosk m r 2 M m 1 A m cosk m r 2 2 M 5 Internationalwhere
m cosk m r dr l 2 A n f n , A.1 where A.1 is valid for n 1, 2, . . ., M. For integration of ς M 2 , let us rewrite it with separation of Goldstone terms as follows: , in turn, considering mixed multiplications which describe cross-modal interactions, one can write M m j cosk i r cosk j r, A.3 and integration similarly to A.1 yields the following expressions for nonlinear terms: indicates that the consideration is implemented for the even coefficient η n 1, if n 2m, 0, otherwise, A.Journal of Mathematics and Mathematical Sciences and the bound for the sum with change of contracting range is expressed as follows: • r means rounding to ∞, and for the second term of A.2 the following is valid: