Third International Conference on Applications of Mathematics to Nonlinear Sciences, Electronic Journal of Differential Equations, Conference 27 (2024), pp. 27–47. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.conf.27.l1 NONLINEAR NON-AUTONOMOUS BOUSSINESQ EQUATIONS ANDREI LUDU, HARIHAR KHANAL, ADRIAN STEFAN CARSTEA Abstract. We study solitary wave solutions for a nonlinear and non-autonomous Boussinesq system with initial conditions. Since the variable coefficients intro- duce distortions and modulations of the solution amplitudes, we implement a multiple-scale approach combining various modes in order to capture the cou- pling 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 fur- ther study. We validate our theoretical results with numerical simulations, and we study their stability. 1. 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 con- ditions 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 2010 Mathematics Subject Classification. 35Q51, 35Q53, 35G50, 34E13, 93C70. Key words and phrases. Boussinesq; non-autonomou; nonlinear; multiple-scale; soliton. ©2024 This work is licensed under a CC BY 4.0 license. Published August 20, 2024. 27 28 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 non-hydrostatic σ−model for waves over rapidly varying topography [13]. Non- autonomous equations related to fluid dynamics models are analyzed for the Burg- ers equation [6], and for the MKdV equation [12]. Other approaches for nonlinear waves over a variable bottom include studies of integrability using the four-wave for- malism [36, 22], and special types of variable re-scaling of coordinates [23]. Another traditional approach is to deal with variable boundary conditions using the non- local Dirichlet-Neumann operator [41], and the studies using this operator while performing a conformal transform [17]. A very important model equation for such problems is provided by various types of Boussinesq systems [17, 46]. Such models are used in engineering in the study of the importance of specific water waves in the coastal and harbor design. For these purposes, the mathematical properties of the wave solutions, such as integrability and stability of solitary waves, are crucial. The goal of this article is to build an integrable procedure for the nonlinear non- autonomous Boussinesq system of equations. To accomplish this goal we combine the results from two limiting theories: on the one hand the linear approximation for non-autonomous equations [6, 8, 12, 14, 23, 35, 41], and on the other hand, the theory of integrable Boussinesq and Broer-Kaup (BK) systems with the corre- sponding soliton solutions, where the nonlinearity is considered, but the system is autonomous [21, 22, 24, 26, 28, 29, 36, 45, 47]. The solitary wave solutions to the classical Boussinesq equation (with and without a restoring force) have also been obtained in [27, 42]. This article is organized as follows. In section 2 we present the nonlinear non- autonomous Boussinesq system and the associated initial conditions. In subsection 2.1 we review and elaborate on the extent to which the non-autonomous Boussinesq system under study is a well posed problem. Section 3 presents some qualitative asymptotic analysis of this Boussinesq system. In section 3.1 we study the au- tonomous nonlinear limit, and introduce exact one-soliton solutions which will be used further as initial conditions for the non-autonomous system. In subsection 3.2 we study the other asymptotic limit, the non-autonomous linear case. Section 4 is part of the main core of results of this work. In subsection 4.1 we apply the procedure of amplitude modulation for the autonomous case of the Boussinesq sys- tem, which results in a non-dispersive type of differential system. We describe the property of integrability of the nonlinear system, beginning with the inverse scatter- ing theory integrability problem for the autonomous limiting case (the traditional Boussinesq system). In addition, in this section we investigate the integrability of the nonlinear problem with the Zakharov-Kuznetsov multiple-scale theory for wave amplitude modulation of Boussinesq system (which can be mapped into the inte- grable Broer-Kaup system) and we obtain a dispersionless envelope system which is likely to be integrable. The fully nonlinear and non-autonomous case is studied by using the Zakharov-Kuznetsov multiple-scale procedure in subsection 4.2. To ob- tain explicit analytic solutions we use a generalization of the multiple-scale method of integration by using N -waves mixing procedure where we combine the phases of the autonomous solutions. These calculations result in a recursive hierarchy of differential equations. We noticed that even the autonomous case can produce modulations of solution amplitudes, so this limiting case is instructive. Using this multiple-scale generalized approach we obtain a hierarchy of dispersionless systems of equations for the amplitudes of the waves. This systems are integrable since they EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 29 are derived from a completely integrable Boussinesq system by a limiting proce- dure. In section 5 we present numerical results for the non-autonomous nonlinear Boussinesq system under consideration. In subsection 5.1 we describe the numeri- cal 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. 2. Boussinesq non-autonomous nonlinear system We consider a non-autonomous and nonlinear Boussinesq-type of differential system in the form qzt + (qu+ αzu)x + β 3 (qu)xxx = 0, qut + zx + αuux = 0, (2.1) 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 Ls(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) = z0(x), u(x, 0) = u0(x), z0, u0 ∈ Hs(R), (2.2) where Hs(R) is the Sobolev spaceW 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 Lp(R) of sech types, p ≥ 1, in (3.8). 2.1. Well posed problem for the nonlinear and non-autonomous Boussi- nesq 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 z0, u0 belong to a Sobolev space Hs 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 satisfac- tory 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 a = c = d = 0, b = βh3/3 > 0. 30 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 This can be proved by a simple substitution in the first equation in (2.1) which turns the term βh3(u)xxx into −βh3(u)xxt. 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 physi- cally 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 hap- pens 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 ∈ Ck(R) in a weighted Sobolev space W p,k w is given by the sum of the Lebesgue integrals over R of the pth power of the absolute values of products between a weight function w(x) and partial derivatives f, fx, . . . , fxx...x up to order k, denoted ∥f∥Lp w . The weight function must be a strictly positive locally integrable function on R [40]. By Hölder’s inequality, if a Ck(R) function f is Sobolev weighted, and the weight function is Lp(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. 3. 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 zt + ux + α(zu)x + β 3 uxxx = 0, ut + zx + αuux = 0, (3.1) In this limit (2.1) reduce to the Broer-Kaup (BK) nonlinear system which is in- tegrable, 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 EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 31 system (2.1). By applying a Bäcklund transformation to (2.1) x→ X = x 2 √ 3 β , v(X, t) = −αu, w(X, t) = 1 + αz − 2α √ β√ 3 ux, When q = 1, equation (2.1) become (3.1) which are a BK IST-integrable system [24, 28] vt = 1 2 (2w − vX + v2)X , wt = ( vw + 1 2 wX ) X , (3.2) 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) ΨX = P̂Ψ, Ψt = N̂Ψ, Ψ = (Ψ1,Ψ2) T , (3.3) where P̂ = ( −λ+ v 2 1 −w λ− v 2 ) (3.4) N̂ = ( −λ2 − 1 4 (vX − v2) λ+ v 2 −λw − 1 2 (wX + vw) λ2 + 1 4 (vX − v2) ) (3.5) From the compatibility condition ΨX,t = Ψt,X we obtain the zero-curvature equa- tion P̂t − N̂X + [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 = − ( 1 + z − 1 2 ux ) e− ∫ udx, (3.7) 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] zsol = α(4 + α2) [ 2 + (2 + α2) cosh ( α √ 3(4+α2)(2x−t(2+α2)) 4 √ β )] [ 2 + α2 + 2 cosh (α√3(4+α2)(2x−t(2+α2)) 4 √ β )]2 = 4α 1 + cosh α √ 3(x−t)√ β +O(α2) usol = α(4 + α2) 2 + α2 + 2 cosh (α√3(4+α2)(2x−t(2+α2)) 4 √ β ) (3.8) with usol, zsol ≤ α, vsol = 1 + α2 2 , Lsol = 2 √ β α √ 3(4 + α2) , 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/α 32 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 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 non- autonomous case, (2.1) become the differential system qzlint + (qulin)x + β 3 (qulin)xxx = 0, qulint + zlinx = 0, (3.9) where we denote by zlin, ulin 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 zlin = eiωt−µx ∞∑ n=−∞ Bneinx, ulin = eiωt−µx ∞∑ n=−∞ Uneinx, (3.10) where Bn,Un are the amplitudes to be determined by recursion relations and bound- ary 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 sys- tem. 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 ω(k) = ±k √ βk2 3 − 1, (4.1) 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 nonlin- earity. 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 za = ∞∑ n=−∞ zne in[kx+ω(k)t], ua = ∞∑ n=−∞ une in[kx+ω(k)t], EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 33 where zn, un are amplitudes, and we consider the sum to extend over the whole integer set in order to insure reality of the solutions za, ua. In order to include the modulation on the spatial long scale, we introduce the following stretched variables χ = δ(x+ υt), θ = δ2t, (4.2) 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 u0 = δ2U0(χ, θ), u2 = δ2U2(χ, θ), un = δnUn(χ, θ), U−n = U∗ n, z0 = δ2Z0(χ, θ), z2 = δ2Z2(χ, θ), zn = δnZn(χ, θ), Z−n = Z∗ n. Equating to zero the coefficients of every order n exponential independently, we obtain an infinite set of equations for un, (∂θ + inω)zn + (∂χ + ink)un + α(∂χ + ink) ∑ i∈Z ziun−i + β 3 (∂χ + ink)3un = 0, (∂θ + inω)un + (∂χ + ink)zn + α 2 (∂χ + ink) ∑ i∈Z uiun−i = 0. For n = 0 and the order O(δ) we have Z0 = 2αυ(Z1U ∗ 1 + Z∗ 1U1)− α|Z1|2 2(1− υ2) ≡ A1(Z1U ∗ 1 + Z∗ 1U1) +A2|Z1|2, U0 := 2α(Z1U ∗ 1 + Z∗ 1U1)− αυ|Z1|2 2(υ2 − 1) ≡ A3(Z1U ∗ 1 + Z∗ 1U1) +A4|Z1|2. For n = 2 and order O(δ2) we obtain U2 = 3(2αk2Z1U1 − ωαkU2 1 ) 2(4βh3k4 − 3hk2 + 3ω2) ≡ B1Z1U1 +B2U 2 1 , Z2 = (3αk2h− 4αβh3k4)U2 1 − 6αkωZ1U1 2(4βh3k4 − 3hk2 + 3ω2) ≡ B3Z1U1 +B4Z 2 1 . 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 U1,θ + iαk(U0U1 + U∗ 1U2) = 0, Z1,θ + ikβU1,χχ + iαk(Z0U1 + Z1U0 + Z2U ∗ 1 + Z∗ 1U2) = 0. By replacing in these equations the expressions for Z0, U0, Z2, U2 from above, we obtain the equations describing the envelope nonlinear waves over flat bottom U1,θ + iαk[B2|U1|2U1 + (A3 +B1)Z1|U1|2 +A4U1|Z1|2 +A3Z ∗ 1U 2 1 ] = 0, (4.3) Z1,θ + ikβU1,χχ + iαk((A1 +B2)(Z ∗ 1U 2 1 + Z1|U1|2) + (A2 +A3 +B1)|Z1|2U1 +A3Z 2 1U ∗ 1 +A4|Z1|2Z1 +B4|U1|2U1) = 0. (4.4) 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 34 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 this effect is visible even in the linear approximation, which case gives a twisted coupling of equations containing Zn and Un, 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 Znj , Unj 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 multiple- scale 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 k1 + k2 + k3) 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 Bj ,Uj . Based on these linear solutions we build the solution of the system (2.1) in the form z = ∞∑ n,j=−∞ δγnZn,j(χ, θ)Bjeinψj(χ,θ), (4.5) u = ∞∑ n,j=−∞ δγnUn,j(χ, θ)Ujeinψj(χ,θ), (4.6) 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), θ = −δ2t. 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 Zn,j , Un,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 EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 35 Un,j , Bn,j , with all other parameters known (q, α, β,Uj ,Bj). According to the pro- cedure 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 ex- pressions for the corresponding Zn,j , Un,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 U0,jUj = 0. (4.7) Equation (4.7) for the amplitudes U0,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 q ∑ j U1,j,χUj − 2i 3 βυq′ ∑ j (j − iµ)Z1,jBj + ∑ j [ υq − 2i 3 βq(j − iµ)− 2i 3 βq′ω − 1 3 βυq′′ ] Z1,j,χBj = 0, (4.8) where q′ and Zn,j,χ are χ-derivatives of q and Zn,j , respectively. The term for n = 2 has the first non-zero term in O(δ2) and reads∑ j [2iq(j − iµ) + q′]U2,jUj + ∑ j [4 3 βq(j − iµ)2 + 2iωq + 8 3 βωq′(j − iµ)− 2i 3 βωq′′ ] Z2,jBj + 2iα ∑ j (j − iµ)e−ijχ [ U1,jZ1,−1UjB−1e −iχ + U1,jZ1,0UjB0 + U1,jZ1,1UjB1e iχ ] = 0 (4.9) Equations (4.8), (4.9) provide the coupling between wave shape amplitudes Zn,j and velocity amplitudes Un,j , which already occurs at the lowest order. At O(δ3), the equation for n = 1 is β 3 ∑ j [2iq′(j − iµ) + q′′]Z1,j,θBj − β 3 ∑ j (2υq′ + q)Z1,j,χχBj + iα ∑ j,k,l (j − iµ)Zk,lU1−k,jUjBleik(l−j)χ = q ∑ j Z1,jBj , (4.10) The triple summation in the last term in (4.10) manifests the nonlinear coupling between various modes and scales of the amplitudes Zn,j and Un,j , and their cou- pling 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]. 36 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 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) qυ ∑ j U0,j,χUj + ∑ j Z0,j,χBj + α 2 ∑ j,k,l [k(l − j)Uk,lUk,j + Uk,l,χU ∗ k,j + Uk,lU ∗ k,j,χ]UlUjeik(l−j)χ = 0. (4.11) For n = 1, the terms of O(δ) cancel, since they satisfy the linearized equation. At O(δ2), we obtain for n = 1 υq ∑ j U1,j,χUj + ∑ j Z1,j,χBj = 0, (4.12) and for n = 2, 2iωq ∑ j U2,jUj + 2i ∑ j (j − iµ)Z2,jBj + iα ∑ j,k,l (j − iµ)UklU2−k,jUlUjeik(l−j)χ = 0. (4.13) 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 sys- tem 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 U0,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), F1(q)Z1,0 + qU1,0,χ = 0, ψ0qU1,0,χ + Z1,0,χ = 0 by one quadrature, to find U1,0, Z1,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), F2(q)U2,0 + F3(q)Z2,0 + 2αµU1,0Z1,0 = 0, 2iωU2,0 + 2µZ2,0 = 0, to find U2,0, Z2,0, and so on. All symbols Fk(q) are known functions of q(x) = q(χ/δ + υθ/δ2). The system becomes more complicated when j ̸= 0 yet |j| ≤ jmax ̸= 0. It is easy to show by direct calculations that for any jmax there is always a j0 < jmax such that the components with 0 < j ≤ j0 of the solutions are arbitrary, and the terms j > j0 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 EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 37 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. 5. Numerical solutions 5.1. Numerical algorithm. To obtain numerical solutions for the system (2.1) we define U = [ v v̄ ] , F(U) =  v̄ + α q2 vv̄ 1 q v + α 2q2 v̄2  , G(U) = [ β 3 vxxx 0 ] , where v = qz and v̄ = qu. We express (2.1) in the form ∂U ∂t + ∂F(U) ∂x = G(U), (5.1) with a ≤ x ≤ b, 0 ≤ t ≤ T . Let ∆x = (b − a)/M and ∆t = T/N . We construct a grid (xi, tn), with xi = i∆x, i = 0, 1, 2, . . . ,M and τn = n∆τ, n = 0, 1, 2, . . . , N . Let vni = v(xi, tn), v̄ n i = v̄(xi, tn), U n i = [ vni v̄ni ] , Fni = F(Un i ) and Gn i = G(Un i ). We solve the above nonlinear advection dispersion system (5.1) numerically using the following Predictor-Corrector (McCormack) scheme [16]. U∗ i = Un i − ∆t ∆x ( Fni+1 − Fni ) +∆tGn i , (5.2) Un+1 i = 1 2 (Un i +U∗ i )− ∆t 2∆x ( F∗ i − F∗ i−1 ) + ∆t 2 G∗ i (5.3) The third derivative vxxx appearing in the dispersion term Gn i = G(Un i ) is approx- imated by using the following second order difference formulas. vxxx(xi, tn) = −5vni + 18vni+1 − 24vni+2 + 14vni+3 − 3vni+4 2∆x3 +O(∆x2), i = 0, 1; vxxx(xi, tn) = vni+2 − 2vni+1 + 2vni−1 − vni−2 2∆x3 +O(∆x2), i = 2, 3, . . . ,M − 2; vxxx(xi, tn) = 5vni − 18vni−1 + 24vni−2 − 14vni−3 + 3vni−4 2∆x3 +O(∆x2), i =M − 1,M . 38 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 5.2. Analysis of numerical results. To interpret the results from the multiple- scale 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, zsol(x), usol(x) from (3.8), centered at the middle of the space interval, with amplitudes controlled by the nonlinearity param- eter α and chosen in the range 0.5−0.8 and for two values of the soliton half-width Lsol = 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 gov- erned 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. 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 au- tonomous Boussinesq equation (3.8) with Asol = 0.01, Lsol = 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 mod- ulated 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. EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 39 Figure 2. Numerical solutions u(x, t) for (2.1) in the same condi- tions as Figure 1 with initial condition (black) given by the usol(x) autonomous soliton, second (3.8). For this component of the solu- tion, the amplitude slightly increases in time because of onset of instability. Figure 3. Numerical solutions z(x, t) for (2.1). The initial condition (black) is the same one-soliton solution in (3.8) with Asol = 0.01, Lsol = 30 and the equation parameters are the same α = 0.5, β = 0.5 as in Figures 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 Boussi- nesq system on the solitary wave solutions. The perturbation induced in the soli- tary 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 40 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 Figure 4. Numerical solutions u(x, t) for (2.1) in the same condi- tions as Figure 3 with corresponding usol initial condition (black). The shape of the initial condition is highly perturbed and increases in time. It becomes unstable after 10 seconds, and probably ap- proaches a blow-out singularity, showing that the u component is more sensitive to the effect of variable coefficient. 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. Lsol = 5 − 20 and the wavelength of the perturbation coefficient q(x). In gen- eral, for small variations of the variable coefficient around 1, for relatively small initial soliton amplitudes Asol = 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 ampli- tude 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 equa- tions. 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 Asol > 0.015 (i.e. α > 0.55, larger dispersion coeffi- cient β > 0.6 and for larger variation of the perturbation coefficient ϵ > 0.2, we EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 41 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 Asol = 0.016, Lsol = 30. The vari- able coefficient is q = 1 + ϵ sin(2x) with amplitude ϵ = 0.2. The larger values for the coefficient of nonlinearity and dispersion in- duce 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 am- plitude one-soliton with Asol = 0.055, Lsol = 30 with higher non- linearity 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 42 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 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 Lsol = 5, and the same large amplitude Asol = 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 be- comes 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 co- efficient 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) solu- tion. 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 non- linear 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 EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 43 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 perturba- tion 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 oscilla- tions 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 lin- ear 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 per- turbative coefficient is not compensated by the increase in the soliton amplitude (stronger value for the perturbation is not compensated by higher order of nonlin- earity, 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 Lsol = 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 bot- tom over shallow water waves from studies on non-autonomous Euler equations with 44 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 coefficients depending of the bottom topography [3]. Also the Green-Naghdi equa- tions for nonlinear dispersive gravity waves can be mapped into generalized Boussi- nesq equations in non-autonomous form, with coefficient depending on the bound- ary 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 origi- nal 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]. 5.3. Stability of perturbed soliton solution. In the previous section we dis- cussed how an initial-value problem for the system (2.1)-(2.2) is always locally well posed, for initial conditions provided by smooth Lp(R functions. The numerical study above uses for initial conditions the one-soliton solutions (3.8) of the au- tonomous 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 bi- directional 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 multi- solitons, (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 EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 45 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 Lsol = 20 compare to the perturba- tion 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 dis- persion β 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 signifi- cantly 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 dis- persive Boussinesq-type systems, were obtained when the equations and boundary conditions are formulated in curvilinear coordinates [17]. 6. 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 cor- responding 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 dispersion- less 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 nonlin- ear 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. Acknowledgments. A. Ludu was partial supported by the Program ONR-SFRP- 2021/2023 for the discussions of some of the models presented in this paper. Also he is grateful for enlightening discussions with J. L. Bona and H. Chen. 46 A. LUDU, H. KHANAL, A. S. CARSTEA EJDE-2022/CONF/27 References [1] J. L. Bona, M. Chen, J. -C. Saut; Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I: Derivation and linear theory, J. Nonlin. Sci., 12 (2002) 283-318. [2] J. L. Bona, M. Chen, J. -C. Saut; Global Existence of Smooth Solutions and Stability of Solitary Waves for a Generalized Boussinesq Equation, Comm. Math. Phys., 118 (1988) 15-29. [3] R. Camassa, D. D. Holm, C. D. Levermore; Long-time effects of bottom topography in shallow water, Physica D 98 (1996) 258-286. [4] C.-H. Chang, C.-J. Tang, C. Lin; Vortex generation and flow pattern development after a solitary wave passing over a bottom cavity, Computers and Fluids, 53 (2012) 79-92. [5] C. -L. Chen, X. -Y. Tang, S. -Y. Lou; Solutions of a (2+1)−dimensional dispersive long wave equation, Phys. Rev. E 66, 3 (2002) 036605. [6] S.-J. Chen, X. Lüa, X.-F. Tang; Novel evolutionary behaviors of the mixed solutions to a generalized Burgers equation with variable coefficients, Comm. Nonlin. Sci. Numer Simulat, 95 (2021) 105628. [7] S. -K. Chua; On Weighted Sobolev Spaces, Can. J. Math., 48, 3 (1996) 527-541. [8] D. Clamond, D. Dutykh; Accurate fast computation of steady two-dimensional surface gravity waves in arbitrary depth, J. Fluid Mech., 844 (2018) 491-518. [9] P. A. Clarkson, E. Dowie; Rational solutions of the Boussinesq equation and applications to rogue waves, Trans. Math. Appl., 1 (2017) 1-26. [10] A. Compelli, R. Ivanov, M. Todorov; Hamiltonian models for the propagation of irrotational surface gravity waves over a variable bottom, Phil. Trans. Roy. Soc. A 376, 2111 (2018) 20170091. [11] W. Craig, P. Guyenne, D. P. Nichols, C. Sulem; Hamiltonian long-wave expansions for water waves over a rough bottom, Proc. Roy. Soc. A, 461 (2005) 839-873. [12] C. Dai, J. Zhu, J. Zhang; New exact solutions to the mKdV equation with variable coefficients, Chaos, Solitons and Fractals 27 (2006), 881–886. [13] F. Dias, D. Dutykh; Dynamics of tsunami waves in Extreme man-made and natural hazards in dynamics of structures (Springer, Dordrecht 2007), 201-224. [14] D. Dutykh, D. Clamond; Modified shallow water equations for significantly varying seabeds, Appl. Math. Modelling 40 (2016), 9767–9787. [15] D. Dutykh, O. Goubet; Derivation of dissipative Boussinesq equations using the Dirichlet-to- Neumann operator approach, Math. Comp. Sim., 127 (2016), 80–93. [16] C. A. J. Fletcher; Computational Techniques for Fluid Dynamics 2: Specific Techniques for Different Flow Categories, 2nd Ed, Springer Series in Computational Physics, 1991. [17] A. S. Fokas, A. Nachbin; Water waves over a variable bottom: a non-local formulation and conformal mappings, J. Fluid Mech. 695 (2012), 288-309. [18] M. F. Gobbi, J. T. Kirby, G. E. Wei; A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O(kh)4, J. Fluid Mech., 405 (2000), 181-210. [19] I. M. Gorban; A Numerical Study of Solitary Wave Interactions with a Bottom Step, Con- tinuous and Distributed Systems II: Theory and Applications (2015), 369-387. [20] J. H. Herterich, F. Dias; Extreme long waves over a varying bathymetry numerical, J. Fluid Mech., 878 (2019), 481-501. [21] R. Hirota, J. Satsuma; Nonlinear Evolution Equations Generated from the Bäcklund Trans- formation for the Boussinesq Equation, Prog. Th. Phys. 57, 3 (1977), 797-807. [22] P. A. E. M. Janssen; Nonlinear Four-Wave Interactions and Freak Waves, J. Phys. Oceanog- raphy 33 (2003) 863-884. [23] O. V. Kaptsova, D. O. Kaptsov; Exact solution of Boussinesq equations for propagation of nonlinear waves, Eur. Phys. J. Plus (2020) 135-723. [24] D. J. Kaup; A Higher-Order Water-Wave Equation and the Method for Solving It, Prog. Theor. Phys., 54, 2 (1975) 396-408. [25] B. A. Kupershmidt; Mathematics of dispersive water waves, Commun. Math. Phys., 99, 1 (1985) 51-73. [26] Y. S. Kivshar, B. A. Malomed; Dynamics of solitons in nearly integrable systems, Rev. Mod. Phys., 61, 4 (1989) 763. EJDE-202X/CONF/27 BOUSSINESQ EQUATIONS 47 [27] N. Kolkovska, V. M. Vassilev; Solitary Waves to Boussinesq Equation with Linear Restoring Force. AIP Conf. Proc., 2164 (2019) 110005. [28] Y. Li, W. -X. Ma, J. E. Zhang; Darboux transformations of classical Boussinesq system and its new solutions, Phys. Lett., 275 (2000), 60-66. [29] W.-X. Ma; Integrable couplings of vector AKNS soliton equations, J. Math. Phys., 46 (2005) 033507. [30] L. Molinet, T. Raafat, I. Zaiter; The classical Boussinesq system revisited, arXiv preprint arXiv:2001.11870 (2020) [31] D. E. Mitsotakis; Boussinesq systems in two space dimensions over a variable bottom for the generation and propagation of tsunami waves, Math. Comp. Sim., 80 (2009), 860–873. [32] B. T. Nadiga, L. G. Margolin, P. K. Smolarkiewicz; Different approximations of shallow fluid flow over an obstacle, Phys. Fluids, 8, 8 (1996), 2066-2077. [33] O. Nakoulima et al; Solitary wave dynamics in shallow water over periodic topography, Chaos, 15 (2005) 037107. [34] K. E. Parnell et al.; Ship-induced solitary Riemann waves of depression in Venice Lagoon, Phys. Let. A, 379 (2015), 555-559. [35] S. Pierini, M. Ghil, M. D. Chekroun; Exploring the pullback attractors of a low-order quasi- geostrophic ocean model: The deterministic case, J. Climate 29, 11 (2016), 4185-4202. [36] S. Ponce de León A. R. Osborne; Role of Nonlinear Four-Wave Interactions Source Term on the Spectral Shape, J. Mar. Sci. Eng., 8 (2020) 251. [37] A. A. Saakyan; Convergence of Double Fourier Series after a Change of Variable, Math. Notes, 74, 2 (2003) 255–265; For basic reference to the Lipschitz condition and Fourier series convergence see: R. A. Adams, and J. J. F. Fournier, Sobolev Spaces (Academic Press, 2003). [38] H. H. Sohrab; Basic Real Analysis, vol. 231 (Birkhäuser 2003) p. 142. [39] X. -Y. Tang, S. -Y. Lou, Y. Zhang; Localized excitations in (2 + 1)−dimensional systems, Phys. Rev. E 66, 4 (2002) 046601. [40] B. V. Turesson; Nonlinear Potential Theory and Weighted Sobolev Spaces (Springer 2000) section 1.2.1. [41] R. M. Vargas-Magaña, P. Panayotaros; A non-local Dirichlet-to-Neumann operator: A Whitham–Boussinesq long-wave model for variable topography, Wave Motion 65 (2016) 156–174. [42] V. M. Vassilev, P. A. Djondjorov, M. T. Hadzhilazova, I. M. Mladenov; Traveling Wave Solutions of the One-Dimensional Boussinesq Paradigm Equation, AIP Conf. Proc. 1561 (2013) 327- 332. [43] J. Yu; Revisiting terrain-following Boussinesq equations on a highly variable periodic bed, J. Ocean Eng. Marine Eng., 5 (2019) 403-412. [44] J. Yu; Waveform of gravity and capillary-gravity waves over bathymetry, Phys. Rev. Fluids, 4 (2019) 014806. [45] V. E. Zakharov, E. A. Kuznetsov; Multi-scale expansions in the theory of systems integrable by the inverse scattering transform, Physica 18 D (1986) 455-463. [46] J. E. Zhang, Y. Li; Bidirectional solitons on water, Phys. Rev. E, 67 (2003) 016306. [47] D. Zhao, Y.-J. Zhang, W.-W. Lou, H.-G. Luo; AKNS hierarchy, Darboux transformation and conservation laws of the 1D nonautonomous nonlinear Schrödinger equations, J. Math. Phys., 52 (2011) 043502. Andrei Ludu Embry-Riddle Aeronautical University, Department of Mathematics & Wave Lab, Day- tona Beach, FL, USA Email address: ludua@erau.edu Harihar Khanal Embry-Riddle Aeronautical University, Department of Mathematics , Daytona Beach, FL, USA Email address: harihar.khanal@erau.edu Adrian Stefan Carstea Department of Theoretical Physics, National Institute of Physics and Nuclear Engi- neering, Bucharest-Măgurele 077125, Romania Email address: acarst@theory.nipne.ro 1. Introduction 2. Boussinesq non-autonomous nonlinear system 2.1. Well posed problem for the nonlinear and non-autonomous Boussinesq system 3. Asymptotic approach 3.1. Autonomous nonlinear limit 3.2. Non-autonomous linear limit 4. Multiple-scales method for the Boussinesq system 4.1. Amplitude modulation in the autonomous case: a dispersionless system 4.2. Multi-scale analysis of the non-autonomous nonlinear system 5. Numerical solutions 5.1. Numerical algorithm 5.2. Analysis of numerical results 5.3. Stability of perturbed soliton solution 6. Conclusions Acknowledgments References