NONLINEAR NON-AUTONOMOUS BOUSSINESQ EQUATIONS

. We study solitary wave solutions for a nonlinear and non-autonomous Boussinesq system with initial conditions. Since the variable coefficients introduce distortions and modulations of the solution amplitudes, we implement a multiple-scale approach combining various modes in order to capture the coupling between the nonlinear evolution and the effect of the variable coefficient. The differential system is mapped into a solvable system of nonlinear and non-autonomous ODE which is integrable by recursion procedures. We show that even in the limiting autonomous case, the multiple-scale approach gives a new possibly integrable dispersionless coupled envelope system, which deserves further study. We validate our theoretical results with numerical simulations, and we study their stability.


Introduction
There is a high and sustained interest for the scientists and engineers to study and understand the nonlinear waves and soliton propagation under variable conditions which occur in real coastal and oceanographic applications [9,13,18,34].The mathematical modeling of such systems, even in the two-dimensional case, is inherently difficult because it combines the complexity of solving nonlinear and non-autonomous equations.To our knowledge there are no exact results providing integrability and constructing solitons or other nonlinear waves for such equations with coefficients depending on space.In consequence, there is a large amount of approximate theories developed by using various techniques and making various hypotheses, such as linearization, slowly varying waves, as well as numerical, or semi-numerical, methods that have generated a large amount of data.Over the last decade, a large body of literature has evolved attempting to determine the most appropriate analytical or numerical approach to understand the dynamics of nonlinear waves governed by non-autonomous equations.Among such approaches we mention, for example, the use of conformal-mapping spectral method in the study of run-up of waves over vertical walls and breakers [20], the (Saint-Venant) nonlinear shallow water model for waves over significantly varying seabeds [14], Benjamin-Bona-Mahony model and smoothed particle hydrodynamics method for understanding tsunami generation [13,31], pseudo-spectral method for the Babenko equation and Petviashvili iterations for the study of waves over arbitrary depth [8], the Dirichlet-to-Neumann operator for the dissipative Boussinesq problem [15], or are derived from a completely integrable Boussinesq system by a limiting procedure.In section 5 we present numerical results for the non-autonomous nonlinear Boussinesq system under consideration.In subsection 5.1 we describe the numerical algorithms.In section 5.2 we present relevant examples of numerical solutions, for various combinations of parameters and initial conditions, and we analyze how the predictions of the theoretical results presented in Section 4 can apply or justify these numerical solutions.In subsection 5.3 we discuss the stability of the solitary waves obtained numerically.

Boussinesq non-autonomous nonlinear system
We consider a non-autonomous and nonlinear Boussinesq-type of differential system in the form for the solutions z(x, t), u(x, t) where (x, t) ∈ (−L, L)×[0, ∞) and the space domain can be arbitrary extended L to∞.Subscripts x, t represent differentiation.The two parameters α, β ∈ [0, 1] control the nonlinearity and dispersion, respectively, and the variable coefficient q(x) is a time-independent conveniently smooth and L s (R) bounded function.System (2.1) represents the (1+1) conservative (evolutionary) version of the surface-variable Boussinesq system for surface waves [46].While the β dispersion parameter is not qualitatively relevant for this system because it can be absorbed in a scaling transformation, handling the other two parameters α, q(x) can help analyzing asymptotic limiting situations of this Boussinesq system: the nonlinear autonomous limit, and the linear non-autonomous limit.
With this system we associate regular initial Cauchy conditions, z(x, 0) = z 0 (x), u(x, 0) = u 0 (x), z 0 , u 0 ∈ H s (R), where H s (R) is the Sobolev space W s,2 (R) for some s > 1.In fact, in this paper we will use for the initial conditions only one-soliton solutions of the autonomous limit of (2.1) which obey the requested constraints being rapidly decreasing functions in L p (R) of sech types, p ≥ 1, in (3.8).

2.1.
Well posed problem for the nonlinear and non-autonomous Boussinesq system.The autonomous (q = 1) version of the Boussinesq system (2.1) together with the initial conditions (2.2) was shown to be a linearly well posed problem [1], and actually locally nonlinearly well posed for the case when β > 0 [2] if the initial conditions functions z 0 , u 0 belong to a Sobolev space H s with s > 1.
For some generalized Boussinesq systems, which include our case, if the system has Hamiltonian form, the problem becomes even globally well posed in the physically relevant realm of small-amplitude, long-wavelength disturbances.Consequently, these types of autonomous Boussinesq systems represent a good set of models for the propagation of long-crested waves in the small amplitude, long waves with Stokes number of order O(1) regime (Boussinesq regime) with satisfactory mathematical theories, at least as regards the pure initial-value problems (2.2).Our system (2.1) falls into the so-called C-1 category of well posed problems from [1,2] because the conditions for the coefficients of the equations in these works a ≤ 0, c ≤ 0, d ≥ 0, b ≥ 0 are fulfilled in our case, namely This can be proved by a simple substitution in the first equation in (2.1) which turns the term βh 3 For the nonlinear non-autonomous Boussinesq case, in [2,30] it is shown that for the initial-value problem (2.2) under the restriction C-1 (meaning for physically relevant initial disturbances) the problem is globally well posed in time.On the other hand, numerical simulations for this case [2] indicate that the equations do feature singularity formation in finite time for large initial data, just as happens for KdV-type unidirectional models in the same long-wave regime.Moreover, they found that non-homogeneous boundary conditions imposed at finite spatial positions often intrude just as they do for unidirectional KdV models [1].
To apply these results to our non-autonomous case, we need to generalize the boundness criteria to a weighted Sobolev space, in order to include the variable coefficient q(x).The norm of a function f ∈ C k (R) in a weighted Sobolev space W p,k w is given by the sum of the Lebesgue integrals over R of the p th power of the absolute values of products between a weight function w(x) and partial derivatives f, f x , . . ., f xx...x up to order k, denoted ∥f ∥ L p w .The weight function must be a strictly positive locally integrable function on R [40].By Hölder's inequality, if a C k (R) function f is Sobolev weighted, and the weight function is L p (R), then f is also in W p,k w [40,7].It results that the well posed criteria in [1,2,30] can be applied to the non-autonomous Boussinesq system if the variable coefficient fulfills the condition to be a weight functions.It results that for smooth and locally integrable coefficient functions q(x) with values close to 1 our Boussinesq system represents a well posed problem, at least for certain finite time interval following the initial moment.Since numerical simulations used for the validations of the analytical results (see section 5) are invariably performed on bounded domains, if we use homogeneous boundary conditions placed relatively far away from the support of the initial condition functions, we conclude that the system (2.1)-(2.2) represents at least a locally well posed problem.

Asymptotic approach
3.1.Autonomous nonlinear limit.If we take q(x) → 1, Equations (2.1) become autonomous, and we obtain the traditional Boussinesq nonlinear system In this limit (2.1) reduce to the Broer-Kaup (BK) nonlinear system which is integrable, and has three independent Hamiltonian structures [24,46].Integrability for the BK system was proved using the inverse scattering transformation (IST), but there are also integrability proofs using the Bäcklund transformation [21,23], or the Darboux transformation [28].The flat-bottom Boussinesq system is also a member of the Ablowitz-Kaup-Newell-Segur hierarchy (AKNS), [24], having exact rational solutions relevant to the occurrence of rogue waves, [9], a multi-soliton solution that is expressed in a closed implicit form, [46,13,18], as well as exact solutions obtained by the Painlevé method [5].Our first step is to obtain one-or multi-soliton solutions for the autonomous system (3.1), in order to build the perturbative solutions for the non-autonomous system (2.1).By applying a Bäcklund transformation to (2.1) When q = 1, equation (2.1) become (3.1) which are a BK IST-integrable system [24,28] as a consequence of the compatibility condition (zero-curvature equation) for the associated linear spectral problem.In order to demonstrate the integrability of this BK system (3.2) we consider the linear spectral problem for the vector Ψ(X, t) where From the compatibility condition Ψ X,t = Ψ t,X we obtain the zero-curvature equation Pt − NX + [ N , P ] = 0, (3.6) which generates the BK system (3.2) and proves its IST integrability.Technically, to construct the explicit solutions we can follow [46] and perform a Miura transform [46] q = e udx , r = which maps (3.2) into the first member of AKNS hierarchy system.For example, a right-moving one-soliton solution of the autonomous version (3.1) has the form [46] where the last two expressions are the traveling velocity v of this one-soliton, and its half-width L. By applying other types of Darboux transforms we can obtain other multi-soliton solutions.This soliton solution has the property that its group velocity increases with the amplitude of the soliton, and decreases with the increasing of the half-width.The stronger the coefficient of nonlinearity α, the narrower becomes the soliton.This one-soliton has a single peak when its amplitude is less than 2/α and double-peak when the wave amplitude is larger than 2/α, thus having some remarkable features.
3.2.Non-autonomous linear limit.In the linear limit α → 0 of the nonautonomous case, (2.1) become the differential system where we denote by z lin , u lin the solutions of this linear approximation, to discern them from the solutions z, u of the full nonlinear system.While there is no general analytic solution for the system (3.9) for an arbitrary coefficient function q(x) we can always write the solutions in terms of a Fourier integral with respect to time, and a Fourier series with respect to the bounded variable x (which approaches a Fourier integral in the limit L → ∞).This linear problem for the non-autonomous Boussinesq system was solved in literature [17,44,43] for various type of variable coefficients q(x).When this coefficient function is periodic, the linearized solutions for (3.9) can be expressed as series in terms of Floquet solutions in the Fourier time representation.For all these cases, the generic solution for the system (3.9) can be written in the Fourier time representation in the form where B n , U n are the amplitudes to be determined by recursion relations and boundary conditions, and the constants ω, µ ∈ (0, ∞) parameterize the space of the linear solutions with respect to the Fourier component of frequency, and an arbitrary damping coefficient, respectively.Actually, if the coefficient function q is periodic, the parameter µ is related to the Floquet (Lyapunov) exponent [44,43].
4. Multiple-scales method for the Boussinesq system 4.1.Amplitude modulation in the autonomous case: a dispersionless system.In this section we analyze the effect of amplitude modulation for the system (3.1) in the autonomous case, by using multiple-scales method [45].In the next section we will extend this procedure for the non-autonomous case.If we consider u(x, t), z(x, t) to be small, then we can neglect the nonlinear terms, and from the linear ones, when q = 1, we obtain the dispersion relation and accordingly, u(x, t), z(x, t) will be small monochromatic waves.Starting with such small monochromatic solutions one can ask what will be the effect of nonlinearity.That includes the developing of harmonics and variation (modulation) of amplitudes on long spatial scales.Even slow variations of the coefficient function q(x) can induce such modulations, but we will consider this case in the next section.
The autonomous solutions of (3.1) expanded in harmonics will have the form where z n , u n are amplitudes, and we consider the sum to extend over the whole integer set in order to insure reality of the solutions z a , u a .In order to include the modulation on the spatial long scale, we introduce the following stretched variables where the smallness parameter δ and the scaling factor υ are arbitrary parameters at this point.They will be determined from the balancing harmonics, and of the powers of δ.If we consider the amplitudes depending on these stretched variables, then the amplitudes are modulated according to Equating to zero the coefficients of every order n exponential independently, we obtain an infinite set of equations for u n , For n = 0 and the order O(δ) we have For n = 2 and order O(δ 2 ) we obtain For n = 1 we reproduce the dispersion relation in order O(δ).In order O(δ 2 ) we re-obtain the derivative of the dispersion relation with respect to k (υ = dω/dk the group velocity), and in O(δ 3 ) we obtain the following nonlinear coupled system By replacing in these equations the expressions for Z 0 , U 0 , Z 2 , U 2 from above, we obtain the equations describing the envelope nonlinear waves over flat bottom From (4.3)-(4.4)we see that the autonomous multiple-scale approximation (3.1) of the non-autonomous Boussinesq system (2.1) is dispersionless.The signature of this effect is visible even in the linear approximation, which case gives a twisted coupling of equations containing Z n and U n , without dispersion.Accordingly, we may expect that the weak modulation of nonlinear traveling solutions will drive them toward instability and breaking, even in the autonomous limit (of course this is not mandatory; there are dispersionless systems supporting stable solitons).However, since the system (3.1) can be transformed into the integrable Broer-Kaup (BK) system, we expect that the dispersionless complex envelope system (4.3)-(4.4) to be completely integrable, and consequently deserves further study.

4.2.
Multi-scale analysis of the non-autonomous nonlinear system.In this section we introduce a decomposition of the z, u solutions in a series of amplitudes Z nj , U nj which obey a recursive system of linearized equations.Here we extend the calculations presented in the previous section 4.1 from one sequence of phases, to a set of coupled phases.To obtain the hierarchy of solutions for the nonlinear non-autonomous system (2.1), we present below a generalization of the multiplescale expansion method [45,29,22].Such a multiple-scale approach should be able to circumscribe the coupling between nonlinearity and non-autonomous coefficient function contributions.The procedure is to build a linear combination of traveling modes with different phases.This superposition should be realized at least through a three-waves mixing (like where one phase arises from the linearized non-autonomous solutions (section 3.2) solution and another must be introduced to build the multiple-scales (section 3.1).However, since the contribution of the variable coefficient q(x) is time-independent (no frequency parameter associated to this interaction), it becomes difficult to balance the variable coefficient contribution using a classical multiple-scale procedure.Consequently, we have to use an N -waves mixing procedure [36,22].We proceed with a generalized Zakharov-Kuznetsov multiple-scale expansion procedure [45].Following (3.10) we introduce a sequence of coupling phases given by ψ j (x, t) = (j − iµ)x + ωt, for arbitrary positive parameters µ, ω.We use a perturbation of the general form of the linear non-autonomous solutions from (3.10) for some given arbitrary amplitudes B j , U j .Based on these linear solutions we build the solution of the system (2.1) in the form where we introduced the new mixing amplitudes Z −n,j = Z * n,j and U −n,j = U * n,j , and we are using the same smallness parameter 0 < δ < 1 as in the previous section 4.1, and a similar re-scaling of coordinates into the stretched variables (4.2) by x, t → χ = δ(x + υt), θ = −δ 2 t.
In (4.5), (4.6) we have the exponent γ n = |n| if n ̸ = 0, and γ 0 = 2.The stretched variables facilitate the coupling between the scales of linear solutions and the scales of modulation because of the nonlinearity: ∂ t → δυ∂ χ − δ 2 ∂ θ + inω, ∂ x → δ∂ χ , while the label n provides the mixing of scales ψ j → nψ j weighted by the unknown functions Z n,j , U n,j .
We plug the double series in the first equation in (4.5), (4.6) into (2.1) and this system maps into a power series with respect to δ, for the unknown amplitudes U n,j , B n,j , with all other parameters known (q, α, β, U j , B j ).According to the procedure in [45], we approach the limit δ → 0 and we cancel the terms at the minimal (for every n) power of δ.The system (2.1) is reduced, as δ → 0, to explicit expressions for the corresponding Z n,j , U n,j , except for n = 1 at the order O(δ).The summation over j can be extended to the highest level of accuracy, as needed for the solution.
For the order n = 0 the non-zero relevant terms are obtained for O(δ 2 ) and result in j U 0,j U j = 0. (4.7) Equation (4.7) for the amplitudes U 0,j is homogeneous, and the indeterminacy of its solutions can be physically related to the non-conservation of longitudinal momentum because of the interaction with the variable coefficient.The term with n = 1 is identical zero in O(δ).In O(δ 2 ) the equations of interest for this term is where q ′ and Z n,j,χ are χ-derivatives of q and Z n,j , respectively.The term for n = 2 has the first non-zero term in O(δ 2 ) and reads Equations (4.8), (4.9) provide the coupling between wave shape amplitudes Z n,j and velocity amplitudes U n,j , which already occurs at the lowest order.At O(δ 3 ), the equation for n = 1 is The triple summation in the last term in (4.10) manifests the nonlinear coupling between various modes and scales of the amplitudes Z n,j and U n,j , and their coupling with the modes generated by the variable coefficient q.Equation (4.10) is the first one in the hierarchy to contain nonlinear quadratic terms.This equation has the same structure as the vector AKNS system [29].This result is expected, since a similar model, the Fokas-Lenells non-autonomous system was also proved to be integrable towards the vector AKNS system [47].
We use the same procedure for the second equation from (2.1).For the order n = 0 the first non-trivial relation is obtained at O(δ 3 ) For n = 1, the terms of O(δ) cancel, since they satisfy the linearized equation.At O(δ 2 ), we obtain for n = 1 and for n = 2, Although this hierarchy of differential equations, obtained for various n, j, is still non-autonomous, the resulting system of equations (4.7)-(4.13)and further for higher n, can be solved by a recursion procedure, once a maximum value of j is chosen for practical computation of the summation.This solvability feature is possible because at each order n the solutions can be obtained from the previous step of order n − 1, by either simple algebraic relationships, like in the case of system formed by (4.7), (4.9), (4.11), 4.13), or integrating quadrature, like in the case of systems formed by (4.8), (4.10), (4.12), etc.. To illustrate the consistency of this recursion multiple-scale procedure we can choose for example j = 0.By choosing n = 0 in (4.7) we obtain U 0,0 = 0. From (4.8) obtained at n = 1, i.e.O(δ 2 ) we can integrate the 2 × 2 ODE system given by (4.8) and (4.12), F 1 (q)Z 1,0 + qU 1,0,χ = 0, ψ 0 qU 1,0,χ + Z 1,0,χ = 0 by one quadrature, to find U 1,0 , Z 1,0 as functions of q.For the term n = 2 from (4.9) which is order O(δ 2 ), we have to solve an algebraic system given by (4.9) and (4.13), to find U 2,0 , Z 2,0 , and so on.All symbols F k (q) are known functions of q(x) = q(χ/δ + υθ/δ 2 ).The system becomes more complicated when j ̸ = 0 yet |j| ≤ j max ̸ = 0.It is easy to show by direct calculations that for any j max there is always a j 0 < j max such that the components with 0 < j ≤ j 0 of the solutions are arbitrary, and the terms j > j 0 can be obtained from these ones by quadrature.In the recursion system (4.8)-(4.13)we have often the situations when in the same equation there are either only the time derivatives, or only the spatial derivatives, but not mixed derivatives, as in the classical integrable systems.This is not an issue, because in the re-scaled variable χ we have both the space x and time t dependence.These equations can be integrated by an iteration procedure, where the non-homogeneous terms at one step are computed using the solutions from the previous step.In this way the equations reduce to simple (yet tedious) quadratures, so they can be integrated exactly.
The absence of a dispersion relation (or equivalently, zero dispersion) in (4.8)-(4.13) is only apparent.Equatio.(2.1) was mapped into the system (4.7)-(4.13)and (4.7)-(4.13)are AKNS-type integrable (a generalization of the NLS classic system) and consequently one can obtain soliton solutions for them.This means that a real valued dispersion relation is implicitly present in our equations.The difficulty arising from the apparent absence of a dispersion relation is that it may be difficult to verify the integrability, or to find soliton solutions to these equations using the Hirota bilinear formalism.In our case, the dispersion must be real valued, and it is different from the one for the BK system which has imaginary dispersion, as in the case of non-conservative diffusive systems.In a forthcoming paper we will present a new bilinear form for the Kaup-Boussinesq system (found by one of the present authors ASC), which is in fact a Bäcklund bilinear singular form for the Boussinesq system, which allows a faster calculation of the solutions, including multi-solitons.

Numerical solutions
5.1.Numerical algorithm.To obtain numerical solutions for the system (2.1) we define where v = qz and v = qu.We express (2.1) in the form with a ≤ x ≤ b, 0 ≤ t ≤ T .Let ∆x = (b − a)/M and ∆t = T /N .We construct a grid (x i , t n ), with x i = i∆x, i = 0, 1, 2, . . ., M and τ n = n∆τ, n = 0, 1, 2, . . ., N . Let . We solve the above nonlinear advection dispersion system (5.1)numerically using the following Predictor-Corrector (McCormack) scheme [16]. 3) The third derivative v xxx appearing in the dispersion term G n i = G(U n i ) is approximated by using the following second order difference formulas.

Analysis of numerical results.
To interpret the results from the multiplescale analysis presented in section 4, we solved numerically the system (2.1)-(2.2) for several values of the parameters.We choose a region of width L = 150 and impose homogeneous boundary conditions at the ends of this region, namely z(0, t) = z(150, t) = u(0, t) = u(150, t) = 0. We choose these type of boundary conditions because we investigate the evolution of an initial one-soliton solution (3.8) under the perturbation caused by the variable coefficients, and this initial condition is a highly localized function.We also investigate numerically only the first 10 seconds of the solution evolution, such that the generated solitary wave does not reach the boundaries of the interval.In order to validate this hypothesis, we substitute the homogeneous boundary conditions with periodic boundary conditions.For each set of parameters, the numerical solutions did not change with these new boundary conditions.
Concerning the initial conditions, we study only one-soliton solution obtained from the autonomous Boussinesq system, z sol (x), u sol (x) from (3.8), centered at the middle of the space interval, with amplitudes controlled by the nonlinearity parameter α and chosen in the range 0.5 − 0.8 and for two values of the soliton half-width L sol = 30 and 5.We choose a relatively large dispersion parameter β = 0.5-0.7 in all simulations.We compare the evolution of the initial one-soliton function governed by the non-autonomous system (2.1)-(2.2) with its original uniform evolution in shape and velocity if its dynamics would be governed by autonomous system, with q = 1.The variable coefficient is q = 1 + ϵ sin(2x) with amplitude ϵ = 0.1.The initial soliton breaks into smaller multi-soliton solutions, but appears stable within the time frame.When the envelope is modulated by the secondary solitons the height of the wave slightly increases, which shows rudiments of area conservation, even in the non-autonomous case.Once the secondary multi-solitons lag the main one, the height of solution returns to its initial value.   1 and 2. The variable coefficient q = 1 + ϵ sin(2x) has increased amplitude ϵ = 0.2.The larger variation of the coefficient q induces larger amplitude secondary multi-solitons, and larger variations for the maximum value of the solution.The soliton maintains stability after 10 seconds.
In all numerical simulations we use the same periodic signal form for the variable coefficient q(x) = 1 + ϵ sin(2x) with amplitude in the range ϵ = 0.1 to 0.3, but the same wavelength λ = 4π.
From the results of the numerical simulations presented in Figures 1-8 we can understand the perturbations induced by the variable coefficients q(x) of the Boussinesq system on the solitary wave solutions.The perturbation induced in the solitary waves by the periodic variable coefficients are controlled by the relative ration between their space scales, namely between the half-width of the initial soliton  L sol = 5 − 20 and the wavelength of the perturbation coefficient q(x).In general, for small variations of the variable coefficient around 1, for relatively small initial soliton amplitudes A sol = 0.010 − 0.016 representing nonlinearity coefficient α = 0.1 − 0.2, and for relative small values for the dispersion coefficient β = 0.5, the initial soliton propagates with uniform group velocity and generates small amplitude secondary solitons in its trailing region.Also the solitary wave amplitude has some oscillations during its propagation, probably because of a residual effect of the property of area conservation law for the Boussinesq nonlinear autonomous equations.The occurrence of the secondary multi-solitons becomes more intense when the amplitude of oscillation of the variable coefficient increases towards ϵ = 0.3.For larger soliton amplitude A sol > 0.015 (i.e.α > 0.55, larger dispersion coefficient β > 0.6 and for larger variation of the perturbation coefficient ϵ > 0.2, we Figure 6.Numerical solution z(x, t) for Equre (2.1) for α = 0.65, β = 0.7 at five moments of time.The initial condition is the one-soliton solution with A sol = 0.016, L sol = 30.The variable coefficient is q = 1 + ϵ sin(2x) with amplitude ϵ = 0.2.The larger values for the coefficient of nonlinearity and dispersion induce higher frequency and amplitude perturbations in the solution envelope, while the envelope keeps traveling with the same group velocity and the same mean shape.
Figure 7. Numerical solution z(x, t) starting from a larger amplitude one-soliton with A sol = 0.055, L sol = 30 with higher nonlinearity and dispersion parameters α = 0.8, β = 0.7.We also have a larger amplitude perturbation q(x) = 1 + ϵ sin(2x) with ϵ = 0.3.The envelope of the emerging solution develops a very dense modulation by oscillations of higher amplitude.We assume a part of these larger perturbations in the radiation tail is generated by numerical instability.There are no secondary multi-solitons for this configuration, while the envelope of the solution maintains the same mean shape and group velocity after 10 seconds.obtain stronger perturbations of the solitary wave, as expected.In these situations, the secondary multi-solitons have larger phase velocity and occur even in the front Figure 8. Numerical solution z(x, t) for α = 0.65, β = 0.7, with the same function for the coefficients q(x) = 1 + 0.3 sin(2x), this time the initial condition being a narrower soliton solution with L sol = 5, and the same large amplitude A sol = 0.55.This solution becomes unstable very fast.The initial soliton decays quickly in amplitude, it becomes slightly wider, and generates high frequency oscillations in the tail.
of the original solitary wave, and high frequency dispersive oscillations grow in the radiation tail.When α, β and ϵ exceed some critical values, the solitary wave becomes unstable, breaks into high frequency radiation waves, and quickly decreases its amplitude.
In Figures 1-2 we present the time evolution of the numerical solutions z(x, t) and u(x, t), respectively for coefficient oscillations with amplitude ϵ = 0.1.We notice that the envelope of the initial Boussinesq soliton (black curve) travels uniformly, and its trailing slope is modulated by the generation of secondary small amplitude multi-soliton solutions with half-width close to the wavelength as the periodic coefficient q(x).During the evolution t > 0 the periodic modulation decouples from the solitary wave and degenerates into a radiation tail lagging the solitary wave.The modulation effect of secondary solitons is more pronounced in the u(x, t) solution.Comparing Figures 1-2 with Figures 3-4 we observe that the strength of the modulation is proportional to the amplitude ϵ of the variable coefficient: for ϵ = 0.1 the perturbation effect is weaker than in the case of larger coefficient ϵ = 0.2.The soliton velocity, however is not affected by the amplitude of the variable coefficient.The time evolution of the amplitude of the soliton shows the reminiscence of the area conservation law for the integrable autonomous Boussinesq case.When the perturbation affects the soliton envelope, the soliton amplitude has a slight increase to compensate for the loss of area.Once the perturbation is decoupled from the soliton, its amplitude returns to the initial value.
The theoretical results presented in sections 4.1-4.2,concerning the non-dispersive character of the non-autonomous nonlinear system are validated by the numerical results.For relatively small perturbation compared to the contribution of the nonlinear terms of the system, |q − 1| ≪ α, Figures 1-4, the envelope z(x, t) of initial soliton is quasi-stable in time, experiencing small amplitude oscillations around its initial shape, while the secondary multi-soliton generated by the perturbation travel slower and decouple from the solitary wave.The evolution of the solitary wave velocity u(x, t) is affected stroger by the perturbation, and develops in time high amplitude instabilities.
In Figures 5-8 we present the effect of increasing the intensity of the perturbation on solution z(x, t) for several initial solitons (3.8) of different amplitudes and using different dispersion coefficient values.In Figure 5 we present the numerical solution z(x, t) obtained for α = 0.5, β = 0.5, ϵ = 0.2.The perturbation induced by the variable coefficient has a weak effect on the propagating solitary wave in this case, preserving the constant velocity and shape, within small amplitude oscillations caused by the emergence of secondary small amplitude solitons in the trailing part.In Figure 6 we present a case with α = 0.65, β = 0.7, ϵ = 0.2.We notice an increase in the non-autonomous perturbation effect through a larger amplitude modulation of the solitary wave envelope.At the same time, the generation of secondary solitons seem to be overwhelmed by occurrence of high frequency linear waves in the trailing part, while the front of the solitary wave begins to feel a weakly effect of perturbation.Nevertheless, we can state that the solitary wave solution is highly modulated but still stable.In Figure 7 we present a situation with α = 0.8, β = 0.7, ϵ = 0.3.The effect of increasing the amplitude of the perturbative coefficient is not compensated by the increase in the soliton amplitude (stronger value for the perturbation is not compensated by higher order of nonlinearity, even balanced by higher value for the dispersion coefficient).The solitary wave experience large amplitude and high frequency perturbation modulating the whole envelope and traveling together with the solitary wave.For t > 10s this solution becomes unstable.Finally, in Figure 8 we present the extreme case of α = 0.65, β = 0.7, ϵ = 0.3.We notice that the initial soliton becomes quickly unstable, reduces its amplitude and becomes modulated by very strong singular waves.
The phenomenon of coupling of the horizontal space scales described in (4.8)-(4.13) is visible in Figure 8.When we choose the half-width L sol = 5 of the initial soliton smaller than the wavelength 4π of the variable coefficient q(x), the coupling of scales described in section 4 by (4.8)-(4.13)generates in the solution perturbative harmonics of frequencies higher than the fundamental wavelength of the periodic coefficient, exactly as shown in the high frequency oscillations present in the tail of the numerical solution.
These observations are confirmed in literature.In [11] the authors derive a two-dimensional Boussinesq-type system and calculate numerical solutions.They conclude that for large perturbation in the Boussinesq equation, representing in their case large bottom variations, nonlinear effects dominate dispersive ones when the amplitude of bottom variations tends to the shoaling limit.In [33] the authors analyze the propagation of long solitary wave pulse over periodic piece-wise constant topography, in the framework of weakly nonlinear-dispersive theory.They notice a similar behavior of the solitary waves with what we present in Figures 1-6.When the obstacle has width comparable or slightly larger than the pulse width the leading transmitted wave keeps its solitary shape on the average.Similar results are presented on the long-time asymptotic effect of varying bottom over shallow water waves from studies on non-autonomous Euler equations with coefficients depending of the bottom topography [3].Also the Green-Naghdi equations for nonlinear dispersive gravity waves can be mapped into generalized Boussinesq equations in non-autonomous form, with coefficient depending on the boundary conditions [32].Different non-autonomous variations of generalized-Boussinesq equations have been used to study soliton propagation in shallow fluid flow over topography [3,32].Their numerical solutions predict the breaking of the original soliton into a train of upstream propagating smaller solitary waves.However, the balance of dispersion to nonlinearity is maintained through a remarkably large range, a fact which tends to further justify the use of the Boussinesq and KdV approximations in the homogenization limit [11].

Stability of perturbed soliton solution.
In the previous section we discussed how an initial-value problem for the system (2.1)-(2.2) is always locally well posed, for initial conditions provided by smooth L p (R functions.The numerical study above uses for initial conditions the one-soliton solutions (3.8) of the autonomous Boussinesq system which obey these conditions.However, solitary-wave solutions for the non-autonomous system (2.1)-(2.2) are nonlinearly stable only for specific range of their phase speeds [1,2].The one-soliton initial data are stable bidirectional solitons in an autonomous Boussinesq system [46] and generates solitary waves which evolve into global solutions.
We study the stability of these global solutions for periodic coefficients q(x).This perturbation leads to the appearance of multi-solitons for t > 0. In Figures 1 and  3 representing z(x, t) we notice the emergence of multi-soliton solution, where the birth of secondary solitons of a much smaller amplitude is visible in the trail of the original initial condition.This effect become stronger in the case of the component u(x, t) Figures 2, 4, where the secondary solitons emerge also in the front of the solitary wave.
When the initial soliton travel under the perturbed non-autonomous equations with periodic coefficients, a part of the solution is re-directed in opposite direction while the rest of the solution continues traveling at the same velocity.We notice that the solitary wave evolution can be classified into four type of perturbations: (1) propagation with weak distortions, (2) fission of the initial soliton in smaller multisolitons, (3) fission of secondary solitons and peaking of the original soliton, and (4) complete break-up of the soliton in very high frequency oscillations, beginning with the tail.Instabilities of types 2 and 4 were also observed in [19] where the analogous wave patterns generate dispersion chains.
These type of instabilities were also identified in [10] when the authors model the propagation of a soliton wave over a bore in the bottom.Namely, in [10] the height of the initial soliton tend to grow in time.If the width of the obstacle is much narrower, or much larger than the width of the obstacle, the initial soliton is partially reflected back, forming the wave reflection, and the rest passes the obstacle and continues to propagate forward.This is the case of quasi-stable solutions when there is a weak interaction between the soliton and the bottom deformation.For bores with width closer to that of the soliton, the transmitted wave (belongs to mode 4) splits into a large number of sub-harmonics.The secondary solitons break up in dispersive wave chains because of the strong interaction between the nonlinearity and the variable boundaries.
As one can notice in Figures 6-8 the increase in the amplitude of the incoming soliton increases the reflected waves which are waves of radiation.They are highly unstable and decay rapidly with time , which can also be seen from these figures.The waves of radiation in the soliton trail reduce the energy of the initial soliton.This is illustrated in Figure 8.The dispersive radiation waves of small amplitude move to the left because their phase velocity becomes significant for the short waves where [10].
In general is difficult to find analytic solutions when the equation is governed by the nonlinear non-autonomous system of equations.Our numerical results show that smaller amplitude periodic perturbation coefficient q(x) is too weak to affect the wave of large wavelength like a solitary wave L sol = 20 compare to the perturbation q = 1 + ϵ sin(2x) with wavelength λ = 4π.When the initial soliton is narrower, like in Figure 8 perturbation induced by the same form of coefficient scatters both the reflected tail waves and the transmitted front waves.
When both the amplitude of the incoming soliton α and the coefficient of dispersion β are increased, in order to maintain the nonlinear balance and hence the solitary wave stability, the effect of transmitted waves with higher phase velocities is stronger, see Figures 6-7.When this Boussinesq system is used to model variable bathymetry, the wave dispersion scattered by the bottom not only affects significantly the wave deformation, on the primary wave height together with the reflected and the transmitted trail waves, but also generates the occurrence of local vortical flow pattern in the proximity of the bottom deformations [4].Similar behavior of nonlinear waves over variable bed, in the case of weakly nonlinear weakly dispersive Boussinesq-type systems, were obtained when the equations and boundary conditions are formulated in curvilinear coordinates [17].

Conclusions
In this article we investigate solitary wave solutions for a nonlinear non-autonomous Boussinesq system.We study the integrability of this nonlinear system for the corresponding autonomous limit.In addition, we investigate the integrability of the non-autonomous nonlinear problem with the Zakharov-Kuznetsov multiple-scale theory for amplitude modulation of Boussinesq system and we obtain a dispersionless envelope system which is likely to be integrable.We use an extension of the Zakharov-Kuznetsov multiple-scale procedure by combining multiple phases.We generate a hierarchy of differential equations for the solution wave mixing.The resulting differential system of order three, represented by relatively complicated equations, can be always solved by iterations.We solve numerically the nonlinear non-autonomous system and present several examples of solutions in order to validate our theoretical results.These results are also of importance for the field of nonlinear fluid mechanics, because the non-autonomous Boussinesq system is related to models for the dynamics of nonlinear waves over variable bottom in the Boussinesq approximation, which represent a very important field of applications.We study the stability of these numerical solutions, and compare our results with the literature.The presented theoretical approach provides methodical value for the general field of the theory of nonlinear non-autonomous systems, while also underlying the connection between the variable bottom Boussinesq, and the B-K and AKNS integrable systems.

Figure 1 .
Figure 1.Numerical solution z(x, t) (red and blue) for (2.1) for α = 0.5, β = 0.5 at three moments of time, labeled in the frames.The initial condition (black) is a one-soliton solution of the autonomous Boussinesq equation (3.8) with A sol = 0.01, L sol = 30.The variable coefficient is q = 1 + ϵ sin(2x) with amplitude ϵ = 0.1.The initial soliton breaks into smaller multi-soliton solutions, but appears stable within the time frame.When the envelope is modulated by the secondary solitons the height of the wave slightly increases, which shows rudiments of area conservation, even in the non-autonomous case.Once the secondary multi-solitons lag the main one, the height of solution returns to its initial value.

Figure 2 .
Figure 2. Numerical solutions u(x, t) for (2.1) in the same conditions as Figure 1 with initial condition (black) given by the u sol (x) autonomous soliton, second (3.8).For this component of the solution, the amplitude slightly increases in time because of onset of instability.

Figure 3 .
Figure 3. Numerical solutions z(x, t) for (2.1).The initial condition (black) is the same one-soliton solution in (3.8) with A sol = 0.01, L sol = 30 and the equation parameters are the same α = 0.5, β = 0.5 as in Figures1 and 2. The variable coefficient q = 1 + ϵ sin(2x) has increased amplitude ϵ = 0.2.The larger variation of the coefficient q induces larger amplitude secondary multi-solitons, and larger variations for the maximum value of the solution.The soliton maintains stability after 10 seconds.

Figure 4 .
Figure 4. Numerical solutions u(x, t) for (2.1) in the same conditions as Figure 3 with corresponding u sol initial condition (black).The shape of the initial condition is highly perturbed and increases in time.It becomes unstable after 10 seconds, and probably approaches a blow-out singularity, showing that the u component is more sensitive to the effect of variable coefficient.

Figure 5 .
Figure 5.The same numerical solution z(x, t) as in Figure 3, obtained in the same conditions, plotted at five moments of time to emphasize the oscillations of the maximum value of the envelope.