Numerical Study of the Generalized Korteweg–de Vries Equations with Oscillating Nonlinearities and Boundary Conditions

The focus here is upon the generalized Korteweg–de Vries equation, ut+ux+1pupx+uxxx=0,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} u_t + u_x + \frac{1}{p} \left( u^p \right) _x +u_{xxx} \, = \, 0, \end{aligned}$$\end{document}where p=2,3,…\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p = 2, 3, \ldots $$\end{document}. When p≥5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p \ge 5$$\end{document}, it is thought that the equation is not globally well posed in time for L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}-based Sobolev class data. Various numerical simulations carried out by multiple research groups indicate that solutions can blowup in finite time for large, smooth initial data. This is known to be the case in the critical case p=5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p = 5$$\end{document}, but remains a conjecture for supercritical values of p. Studied here are methods for controlling this potential blow up. Several candidates are put forward; the addition of dissipation or of higher order dispersion are two obvious candidates. However, these apparently can only work for a limited range of nonlinearities. However, the introduction of high frequency temporal oscillations appear to be more effective. Both temporal oscillation of the nonlinearity and of the boundary condition in an initial-boundary-value configuration are considered. The bulk of the discussion will turn around this prospect in fact.


Background and Motivation
The Korteweg-de Vries equation, u t + u x + uu x + u x x x = 0, written here in scaled, dimensionless variables, was originally derived by Boussinesq in 1877 [27]. It was put forward by him as a model for unidirectional long wavelength surface water waves of small amplitude. Rederived by Korteweg and his student de Vries in 1895 [54], it went unremarked by the European schools of hydrodynamics for decades.
This changed in the 1950s when this equation also appeared as the continuum limit of a mass and string model studied by Fermi, Pasta and Ulam with help from Tsingou [35,41]. It also arose a little later in a plasma physics model derived by Gardner and Morikawa [43]. The advent of the inverse scattering method for solving the Korteweg-de Vries equation in the middle 1960s by Gardner, Greene, Kruskal and Miura, and with later help from Su, (see the review article of Miura [58]) brought this equation and its relatives into a central position in the mathematical firmament. Subsequent laboratory and field studies as well as theoretical work have shown the efficacy of Korteweg-de Vries type models in describing surface water waves (see, for example [18,45,46,76]). Since the time of its ascendancy in the scientific world, the Korteweg-de Vries equation has appeared as a model for a considerable variety of other real-world phenomena.
One of these is in internal wave theory. In the idealized format of a two-fluid system with a lighter layer resting upon a denser lower layer, the so-called extended Korteweg-de Vries equation (also known as the Gardner equation) u t + u x + auu x + u 2 u x + u x x x = 0 arises at a low level of approximation (see [63]). Here, a is a constant dependent upon various details of the physical situation. In particular, at a certain critical ratio of depths and densities, the dominant nonlinear term is no longer quadratic, but is instead cubic. One can view this as an invitation to the study of what is known as the generalized Korteweg-de Vries equation, namely where p = 2, 3, . . .. (It must be confessed that for larger values of p, this equation arose first in mathematical investigations of the interaction between nonlinearity and dispersion and not as a model of a physical situation (see for example [3])).
In the critical case when p = 5, the pure initial-value problem for (1) on the whole line is known to have solutions that blow up in finite time (see the sustained and detailed work of Merle and his several collaborators, starting with [57]). It is an open question whether or not the Cauchy problem for the generalized Korteweg-de Vries equation (gKdV equation henceforth) with a supercritical nonlinearity p ≥ 6 has global solutions for large, smooth initial data. Numerical simulations (e.g. [10,12]) indicate probably not, but rigorous results attesting to this have not yet been forthcoming.
There are several possible ways of controlling the singularity formation that seems to occur when p > 5. One of these might be to append dissipation, so leading to an equation of the form where δ > 0. In fact, as seen in [13], this does not appear to stop the blow up occurring for large initial data. More precisely, for given initial data that blows up in the absence of dissipation, there is a critical positive value δ c of δ such that for δ > δ c , the solution is global and uniformly bounded. However, for δ < δ c , the solution continues to blow up in finite time. One might suppose that the failure of the −u x x term to control the singularity formation is because it is a lower-order term in the equation. Appending a term of the form νu x x x x instead might do the trick. However, this also appears not to work as the remarks in [15] indicate. Another way the blow up might be controlled is to add higher-order dispersion, e.g. something like where δ = 0. This, does indeed manage the blow up as long as p < 8. However, numerical experiments show that for p > 9, blow up reasserts itself. A third possibility is based on recent theoretical studies that have shown that a timedependent oscillation added to the nonlinearity can avert blow-up. Indeed, the solution converges in the limit of large-frequency oscillation to a global solution of a certain limiting problem. It is our purpose here to study this latter phenomenon numerically in the situation where the problem is posed as a periodic initial-value problem. We also study a related problem, namely a boundary-value problem for the supercritical gKdV equation with oscillating boundary conditions. It is found here also that large frequency oscillation can kill potential blow-up phenomena.
For the initial-value problem, (1) is modified by the addition of temporal modulation of its nonlinearity, viz.
where ω 1 and g 1 is a mean-zero periodic function used to manage the nonlinearity. In this situation, posed as an initial-value problem on the entire real line R, the solution is known to be global for large values of ω (see [29,30,60,61] for the case p = 2). In particular, Panthee and Scialom showed in [61] that solutions to the initial-value problem (IVP in what follows) associated with (2) converge to a limit problem as ω → ∞ (see (5) in Section 2). (Note that for the pure initial-value problem, the convective term u x can be dispensed with by moving to traveling coordinates. This is not the case for the boundary-value problem to be discussed presently.) For the pure initial-value problem, a Fourier spectral method is put forward, tested for accuracy and convergence and used in our study. Approximate solutions to the boundary-value problem are obtained via a spectral element method.
The present, quantitative appraisal of the initial-value problem for (2) is placed in a periodic context rather than on the whole real line R. The practice of approximating KdV-type equations on R by associated periodic problems is commonplace and goes back to the work of Zabusky and others in the 1960s. However, it is worth noting that rigorous theory asserting the validity of this procedure on finite, but long time intervals may be found in [32] (and see also [8]).
A related issue arises for the gKdV equation posed as a non-homogeneous problem, viz.
where the solution is forced from the left-hand boundary. Here, g 2 is again a periodic function. While there is no theory showing this problem to be globally well posed in time, the approximate solutions computed here indicate that there is again no blow up for large values of ω whereas the addition of dissipation or higher order dispersion has only limited success in this aspect. The overall conclusion derived from this study is that high frequency oscillation can mitigate the effects of supercritical nonlinearity, thereby resulting in a problem that is globally well posed.
Theory for the initial-value problem for real-valued solutions of the gKdV equations has a long history, starting with the case p = 2 in the 1970s (see [19,48]). In the ensuing years, the theory for these equations has become quite subtle (see for example the monograph [70] of Tao or the nice review lectures of Erdogan and Tzirakis [39], but the reader is cautioned that there are many, many references). However, in the present essay, we will only need the relatively elementary results derived already by Kato in [49] (and see also [3]). These state that the IVP for the gKdV equation (1) is globally well posed in the L 2 -based Sobolev space H k if k ≥ 1 and p = 2, 3, 4. If p = 5, the problem is globally well posed provided the L 2 -norm of the initial data is not too large. The same is true for p > 5, but now under the assumption that the H 1 -norm of the initial data is not too large. These results apply equally to the problem posed on the whole real line and to the periodic initial-value problem.
For initial-boundary-value problems (IBVP) such as (3), much less is available in the literature. Problems such as these for the KdV equation ( p = 2) itself have been studied, but for higher power nonlinearities, not much is known. The articles [22,23] provide a guide to the available literature.
Considerable effort has also been made to develop numerical methods for studying solutions of the KdV and gKdV equations. The extensive literature in this area includes finite difference, finite element, finite volume and spectral methods for the KdV equation itself (see e.g. [2,4,5,37,72] and the references therein). More recently, discontinuous Galerkin methods have been adapted to the gKdV equations and related models [9,59,74,75]. These have the great advantage that local spatial refinement is very easily implemented, which helps especially with solutions that appear to be blowing up. For the gKdV equations, a fully discrete scheme for the (periodic) initialvalue problem using the finite element method for the spatial approximation and a diagonally implicit Runge-Kutta method for the time stepping was introduced in [10][11][12]. These works provided numerical evidence of the existence of solutions that blow up in finite time. More recent numerical contributions to the study of gKdV equations may be found in the works of Klein and his collaborators (see [38,51,52]).
In the present paper, a Fourier spectral method is implemented to solve the periodic initial-value problem (2), whereas a spectral element and a collocation method are used to investigate the initial-boundary-value problem (3). For dispersive equations posed as periodic initial-value problems, the spectral method has the advantage of achieving high accuracy with a relatively coarse spatial discretization. Such methods have been studied in [44,56,67], for example. In the presence of non-periodic and nonhomogeneous boundary conditions, Chebyshev collocation methods for KdV-type equations were proposed in [62,68], and Legendre-Galerkin methods were developed in [55,65]. Both of these methods have been implemented for the spatial discretization of the gKdV equation (3) in the presence of a time-oscillating boundary condition. However, details of this work are only reported here for the spectral element method. The convergence of the two methods has been tested and both are found to achieve spectral accuracy in space and the appropriate algebraic convergence rate in time. As they produced very similar results, only the details of the spectral element method are presented.
The plan of the paper is to first investigate the time-oscillating nonlinearity in Sect. 2. This will include details of our numerical scheme as well as convergence and accuracy tests. Once the computer code is deemed satisfactory, it is used in an exploratory mode to understand in more detail the route to global well posedness despite the presence of a supercritical nonlinearity. A similar pattern is followed in Sect. 3, which deals with the oscillating non-homogeneous boundary-value problem. A short conclusion follows Sect. 3.
A final comment concerns investigations of other equations in the setting of timeoscillating nonlinearities. We point to the works of Abdullaev et al., Konotop and Pacciani, and Zharnitsky et al. [1,53,77], respectively, in the context of Bose-Einstein condensates as well as investigations of the influence of time-oscillation on nonlinear Schrödinger equations by Goubet and collaborators in [34,36], Gabitov and Lushnikov in [42], and that of Cazenave and Scialom [30].

Oscillating Nonlinearities
The focus in this section is the periodic initial-boundary-value problem, which is (2), repeated here for convenience. The function u 0 is a smooth, periodic function with period 2M > 0, g = g(τ ) is a periodic function of τ while the parameter ω is positive and will eventually be taken large. Attention is given to the supercritical cases p = 6 and 7. (The reason for choosing both an even and an odd value of p will appear presently. As mentioned earlier, the initial-value problem for (4) has a satisfactory local existence theory.) In [61], the authors showed that as ω −→ ∞, the solution u of (4) converges to the solution U of with the same initial value, uniformly on bounded time intervals. Here, if is the period of g, then is its temporal average. The concern here is not with issues of how rough u 0 or g can be taken and still have a satisfactory theory of this sort. Indeed, to justify rigorously the convergence of our numerical schemes in reasonably strong spaces, a fair amount of smoothness is required. Notice that regardless of the choice of g, the L 2 -norm of solutions is conserved, which is to say It is straightforward to deduce that as long as the L ∞ -norm of a solution remains bounded on bounded time intervals, then the solution is global, as the L ∞ -bound allows the local existence theory to be successfully iterated, achieving a solution which is bounded on bounded time intervals in whatever space the initial data allows. Hence, singularity formation can be tested by monitoring the L ∞ -norm of a solution. It is at this point that the difference between even and odd values of p presents itself. For even values, the sign in front of the nonlinear term does not matter and solutions blow up with either sign. This is no longer true for odd values of p. Indeed, if p is odd and there is a minus sign in front of the nonlinear term, the solution are global no matter how large the initial data and despite the supercriticality of the nonlinearity. This is the so-called defocusing case. The global well posedness subsists on the conservation law, Together with the invariance of the L 2 -norm, this implies the H 1 -norm to be uniformly bounded and hence that the solution is uniformly bounded. As remarked, this implies the solution is global in time.
The next subsection describes and tests a Fourier-spectral method for approximating solutions of (4). A sequence of numerical experiments is then presented which casts more light on the results in [61].

Fourier-Spectral Scheme
A change of the spatial variable puts (4) on the interval [0, 2π ] and the problem becomes Taking the discrete Fourier transform of this equation yields where u is the Fourier transform in space of u. A standard problem now arises, namely that the linear dispersion relation features high frequencies (large wave numbers k corresponding to very short wavelengths), thereby leading to a very stiff system of ordinary differential equations when a time-stepping is implemented. Indeed, the CFL stability condition would mandate very small time steps indeed. To ameliorate this issue, introduce an integrating factor as in [31,73]. Setting there obtains while the initial data remain unchanged. The stiff term is thereby removed. In numerical computations, Equation (10) is discretized in the form where F is the discrete Fourier transform operator (see e.g. [66,73]) and β = ik 3 π/M 3 . The standard fourth-order Runge-Kutta method is employed for the temporal discretization.
Since exact solutions of (4) are not available, non-homogeneous equations of the form Here, N is the number of Fourier modes being kept in the simulation and the errors are reported at t = 1. Notice the spectral accuracy with a source term f are considered. Choosing the corresponding forcing function f is easily computed. Of course, the solution (13) is not periodic in space. However, since u converges to 0 exponentially rapidly as |x| → ∞, the initial-value problem on the whole line can be considered as a periodic initial-value problem for x ∈ (−M, M). As long as the solution does not have significant amplitude at the boundaries (see again [32]), these two problems yield essentially the same answer on the spatial interval [−M, M].
To test the numerical convergence, define the relative L 2 -error where u ex is the exact solution (13) and u num is the numerical approximation of (12).
Since the L 2 -norm is preserved up to roundoff error in the code, it seems wise to also calculate the relative error in the H 1 -seminorm. For the tests whose resulting errors are recorded below, the following oscillation function and parameters are used: and the simulations were carried out up to a final time t = 1. Tables 1 and 2 display the performance of the numerical method when p = 6 and 7. Notice the spectral convergence in space and the fourth-order convergence in time.

The Supercritical Case p = 6
Reported first is a sequence of numerical simulations of the gKdV equation with time-oscillating nonlinearity as in (4) for the supercritical case p = 6. The auxiliary specifications used in these simulations are As mentioned in (7), the L 2 -norm of solutions is invariant in time. Setting g(t) = sin(200π t), the change in the L 2 -norm of the numerical solution is around 6 × 10 −11 from t = 0 to t = 5, another check on the computer code. Figures 1, 2 and 3 each display spatial traces of the numerical approximations at increasing times of various simulations of (4) with the specifications in (17). The different figures correspond to different choices of time-oscillating functions g and choices of the oscillation frequency ω. As anticipated, with a low frequency parameter (e.g. ω = π ), the numerical approximations appear to blow-up for both supercritical values p = 6, 7. The oscillation is slow enough that it does not have time to effect the solution before it has already gone into blowup mode. Moreover, the structure of the blowing-up peak is very much like what is observed without the oscillation (see our  [12] corresponds to p − 1 in our notation.) However, numerical approximations associated with larger values of ω show a different behavior. For instance, when the associated oscillation is g(ωt) = sin(ωt) and ω is large, there is no indication of singularity formation, though of course one does see the effect of the oscillation-see Figs. 2 and 3. In fact, the L ∞ -norm of the approximate solution appears to decrease more or less monotonically as a function of t in both of these simulations. This result is explained by the theory in [61]. From (6) with g(τ ) = sin(ωτ ), it is clear that m(g) = 0. Thus, for large frequencies, it is expected that the solution will resemble a purely dispersive solution of the linearized On the other hand, if the oscillation takes the form g(ωt) = cos 2 (ωt), the numerical results tell another story. From (6) with this latter specification of g, one calculates that m(g) = 1/2. Thus, blowup is to be expected, at least for large values of ω. When ω = π and the auxiliary data shown in (17) is used, Fig. 5 displays the same sort of blowup seen with the oscillation sin(π t). However, for quite a large middle range, roughly ω in the interval [19.5π, 4292π ], the solutions no longer appear to blow up. In these cases, the time-oscillation parameter ω is apparently not large enough that the solution is modeled well by the solution of But, it does appear that it is large enough to control the blowup. An example is shown in Fig. 6 where ω = 100π . As theory assures, for truly large values of ω, say ω = 10 5 ,  (18) in Fig. 7. At present, we do not have an explanation for this somewhat counterintuitive behavior in the middle range of frequency of oscillation.

The Supercritical Case p = 7
In this section, numerical experiments are reported for p = 7 in (4). Even and odd powers of the nonlinearity are different. Even powers do not depend on the sign in front of the nonlinearity, whereas odd ones do. More precisely, if p is odd and the gKdV equation features a minus sign in front of the nonlinear term, the so-called defocusing case, then solutions of the initial-value problem are global and remain uniformly bounded, are used in what follows in this subsection. Simulations were made with g(t) = sin 2 (ωt) where ω = π , 200π , and 10 5 . Just as for the case p = 6, two different regimes were detected. When the value of ω was small, the solutions blew up quickly, though not so quickly as did the case p   Figure 8 presents the outcome of three simulations, shown in pairs on single graphs. The inputs are identical except for the frequencies, which took the values ω = 18.75π, 19π and 19.25π . As the oscillation is sin(ωt) again, we know already that for small values of ω, blowup is expected whereas for large values the solution should be global. Here, a transition is exhibited. Blowup appears for ω = 18.75π which transitions to a solution at ω = 19.25π where blowup seems about to manifest itself, but can't quite beat the dispersive effect of the oscillation. If p is even and the time oscillation is non-positive so that g(t) ≤ 0, the defocusing gKdV-type equations are presented. The defocusing nature of the equation guarantees the global well-posedness corresponding to quite reasonable classes of initial data u 0 ; for more theoretical details, see e.g. [69,71]. In Fig. 9, the outcome of a numerical simulation using g(t) = − sin 2 (π t) is displayed. The numerical results show slowly decreasing solutions corresponding to the order one value ω = π , on account of the defocusing nature of the forcing.

Boundary Oscillation
A new aspect of control of singularity formation is dealt with in this section. The context is the initial-boudary-value problem (3) with p = 6 and 7. The original mathematical problem that arose when trying to check the accuracy of the KdV approximation is the half-line problem, Boundary-value problems of this sort have arisen when modeling long-crested waves generated by a wavemaker. In such a situation, the function g(t) is determined from a measurement taken at a single spatial point as a function of time at one end of the medium of propagation (see [18,46,76]). Just as for the pure initial-value problem (2), the solution of the two-point boundary-value problem (3) is known to approximate well the solution of (20), when the latter is restricted to the spatial interval [0, L], on a time interval of order L (see [32]). As remarked earlier, the convective term u x cannot be eliminated by a traveling change of variables as the left-hand boundary is then deformed. In this problem, the time-oscillating parameter ω appears in the left-hand boundary condition, u(0, t) = g 1 (ωt) where g 1 is a periodic function.
The theory for either of the problems (3) and (20) is somewhat more complicated than for the initial-value problem. For p = 2, the Korteweg-de Vries equation itself, there is a global theory in [24,25] of smooth solutions for (20) provided the initial-and boundary-data satisfy certain obvious compatibility conditions at (x, t) = (0, 0) (the number of compatibility conditions depends upon the smoothness class in which one is seeking a solution). Such smooth solutions can be approximated by the numerical techniques to be introduced presently. Rigorous theory to this effect for the numerical approximations can be set in order, but this is not pursued here. The theory does depend upon the smoothness of the solutions of course. No smoothness means no approximation in continuous function spaces, for example. Later works slightly generalized these results (see for example [17] and [40] and especially the references to the Russian literature in the paper of Faminski). Using a clever reduction to a forced problem, Colliander and Kenig [33] dealt with general positive integer values of p, but obtained only very weak solutions. So far as we are aware, there is not even a local theory of strong solutions in the literature for non-homogeneous boundary-value problems for the generalized KdV equation (1), let alone for the oscillatory problems under consideration here. In what follows, we shall assume that strong solutions obtain for smooth, compatible initial data, at least locally in time.
Two spectral methods were implemented for the spatial discretization of the gKdV equation in the presence of lateral boundaries; the Chebyshev-collocation method and the Legendre-Galerkin method. Using exact solutions, the convergence of these two methods was tested and both achieve spectral accuracy in space and fourth-order accuracy in time. Only the Legendre-Galerkin results are reported here, as those for the Chebyshev method produced very similar outcomes. A large number of experiments were then initiated, some of which will be described after preliminary discussion of the accuracy of the scheme.
A simple change of the spatial variable puts the problem in the form where p is a positive integer, g(τ ) is a periodic function and ω > 0. This problem has different characteristics from those of the initial-value problem. For a start, energy is not preserved here, as the wavemaker (represented by the left-hand boundary condition) may be introducing or extracting energy from the system.

Description of Spectral Approximations
The spectral element method to approximate solutions of (21) is introduced now. Before describing the method, it is convenient to homogenize the boundary conditions. Define In the new variable v, the initial-boundary-value problem (21) becomes Notice that we are implicitly presuming that u 0 (−1) = g(0), otherwise there is a mismatch at (0, 0). This is in fact the lowest level compatibility condition that was mentioned earlier. Since homogeneous boundary conditions are imposed in (23), the standard framework of spectral methods is available for this problem. The Legendre-Galerkin spectral element method The numerical scheme proposed here is based on the ideas in [44,55,65]. The scheme actually approximates the solution of the homogenized equation (23) and the solution to the oscillating boundary gKdV equation (21) is retrieved by reversing the change of variables (22).
For the temporal discretization, a Crank-Nicolson/leapfrog scheme is employed, where v n (·) is the approximation of v(·, t n ) with t n = n t. The scheme is initiated by using an implicit Euler method with Picard iteration for the first time step. This comprises an unconditionally stable scheme. Let P N denote the linear space spanned by the Legendre polynomials of degree at most N . The subspaces are natural for the spatial situation of interest here. Let T be the final time to which the simulation is aimed. The problem then devolves to finding v k N ∈ V N , k = 0, 1, . . . K = T / t such that for any w ∈ V * N , Here, v n N is the approximation of v(·, n t), I C is the interpolation operator at the Chebyshev-Gauss-Lobatto points, and the operator H (v n N ) := (v n N + rg n ) p has the nonlinearity.
A convenient set of basis functions in this situation is which satisfy the boundary conditions The function L n (x) is the standard Legendre polynomial satisfying with the usual normalization L n (1) = 1. For N ≥ 3, we have Taking w = ψ m , m = 1, . . . N − 3 and writing v N in terms of this basis, viz v N = N −3 k=0ṽ k φ k , the scheme just outlined can be written in matrix form in the usual way. The convergence of this numerical scheme is tested by creating an exact solution of (23). We use v(x, t) = sin(π x)(1 − x) cos(π t), with g(t) = cos(π t).
The function v satisfies the three homogeneous boundary conditions, and while it does not satisfy the equation, a relevant forcing term is easily calculated. Table 3 displays the relative L 2 -and H 1 -seminorm-errors of the Legendre-Galerkin spectral method just outlined for solving (23) versus the time discretization t with p = 2 and p = 6. The relative L 2 -and H 1 -seminorm-errors are defined just as in (14) and (15). Seen clearly is the second-order convergence rate one associates to Crank-Nicholson time stepping. Table 4 shows the relative L 2 -and H 1 -seminormerrors of the Legendre-Galerkin spectral method for solving (23) versus the number N of modes being kept in the simulation, again for p = 2 and p = 6. Note the spectral accuracy that obtains.

The Supercritical Case p = 6
In this section, numerical results for the gKdV equations in a finite domain (21) with oscillating boundary conditions are presented. To begin, take p = 6 and N = 100, t = 10 −4 .  Numerical approximations of (21) with p = 6 and the specifications (32) and several choices of boundary conditions and frequencies ω at the left-hand side of the interval were run. Figure 10 shows a typical outcome when ω = π and g(t) = 3 sin(ωt). With this low temporal frequency, the numerical solution appears to blow up. The singularity formation starts to show itself near the left-hand boundary which is not surprising. For example, the amplitude of the solution shown in Fig. 10 goes from about 4.5 to almost 300 in the last time interval of length 0.005. A very similar result obtained when ω = π and g(t) = 3 sin 2 (ωt). The display is redundant and so is omitted. Notice the blow-up structure for this boundary forced problem is different from that of the pure initial-value problem for p = 6 (see the case p = 5 in [12]). However, when the frequency ω of the boundary data is increased, the numerical approximations no longer show indications of blowing up. Instead their L 2 -norm appears to stabilize around something like a limit cycle. This is true for both the mean-zero boundary data 3 sin(20π t) (see Fig. 11) and the positive-mean boundary oscillation 3 sin 2 (60π t) (see Fig. 13). This stabilization of the L 2 -norm continued far beyond t = 1. Indeed, as seen in Fig. 12, the solution itself gives way to a structure that is periodic in t. This is perhaps not surprising in hindsight, as it is known experimentally [18] that real water waves in a channel forced periodically by a wavemaker become periodic in time at each point down the channel. It is also known rigorously that the KdV and the BBM equation on the half line, forced periodically from the left, yields a solution that is asymptotically periodic in time at each spatial point (see [18,21,26]). Figure 13 displays the stable numerical solution of the gKdV equation with higher time frequency parameter ω = 60π . Since the boundary condition is non-negative, it requires higher time-oscillatory frequencies than the previous case in Fig. 11 to obtain bounded numerical solutions. In Fig. 13A, the stabilization of the L 2 -energy of the numerical solution is shown from t = 0 to t = 1. This stabilization of the L 2 -norm appears to continue indefinitely. Figure 13B displays the numerical approximation at t = 5. The numerical experiments are repeated with an odd power p = 7 of the nonlinearity. The Legendre-Galerkin approximation continues to be used in the investigation, though one finds the same results using the Chebyshev-collocation method. The following specifications were used: N = 100, t = 10 −4 , p = 7, g(t) = 2.5 sin(ωt).
In Fig. 14, the numerical solution appears to be forming a singularity when ω = π . The blow-up pattern is the same as in the case p = 6. However, when ω is increased to 20π , the numerical solution of the gKdV equation shows no signs of unbounded behavior. Figure 15A displays the L 2 -energy of the numerical solutions from t = 0 to t = 1; for t > 1, the same pattern is observed. Figure 15B shows the numerical

Conclusion
Numerical studies of the supercritical gKdV equations adorned with high frequency oscillations are presented. In Sect. 2, the IVP associated to the supercritical gKdV equation (4) with a temporal oscillation attached to the nonlinearity was considered. As the frequency of the time-oscillation increases, the numerical solution resembles more and more the solution of the limit problems (5) in which the oscillation is replaced by its mean value. This is consistent with the theoretical results [61]. When p is odd, the defocusing case can appear with the correct sign of the forcing. When this is simulated, there are no signs of singularities. This provides evidence that the blow-up solutions observed in other circumstances are not caused by instabilities connected with the oscillations. In Sect. 3, the supercritical gKdV equations (21) have been studied in a finite domain with zero boundary conditions at the right and an oscillating boundary condition at the left. Due to the non-homogeneous boundary conditions, mathematical studies of the gKdV equations (21) are technically more difficult and relatively little progress has been made in this direction. For numerical experiments, an accurate spectral collocation method and spectral element method were implemented and tested. Applying the numerical scheme, a blow-up structure, which started to appear near the oscillating boundary, was found. As the temporal frequency of the boundary oscillation was increased, non-explosive numerical solutions of (21) are observed. In the case of boundary oscillation, the mean-zero aspect does not appear to play a role. Several interesting avenues of investigation present themselves as a consequence of what has been observed. The question of what other equations can be managed by way of oscillations in their nonlinearity is an obvious line of development. As far as boundary-value problem are concerned, this is largely virgin territory, both analytically and numerically. This also seems worth further investigation.