Numerical simulation of nonlinear dispersive quantization, preprint 2012

When posed on a periodic domain in one space variable, linear dispersive evolution equations with integral polynomial dispersion relations exhibit strikingly different behaviors depending upon whether the time is rational or irrational relative to the length of the interval, thus producing the Talbot effect of dispersive quantization and fractalization. The goal here is to show that these remarkable phenomena extend to nonlinear dispersive evolution equations. We will present numerical simulations, based on operator splitting methods, of the nonlinear Schrodinger and Korteweg--deVries equations with step function initial data and periodic boundary conditions. For the integrable nonlinear Schrodinger equation, our observations have been rigorously confirmed in a recent paper of Erdoǧan and Tzirakis, [10].


Introduction
In the early 1990's, Michael Berry, [1,2], discovered that the time evolution of rough initial data on periodic domains of the free space Schrödinger equation exhibits radically different behavior depending upon whether the time is a rational or irrational multiple of the length of the space interval. Specifically, given a step function as initial conditions, one finds that, at rational times, the solution is piecewise constant, but discontinuous, whereas at irrational times it is a continuous but nowhere differentiable fractal-like functions. Berry named this the Talbot effect, after an experiment performed by the inventor of the photographic negative, [22]. His discovery was then rigorously analyzed and further investigated in [14,18,19,20,23].
In [17,18], it was shown that the same dispersive quantization and fractalization phenomena appear in a general periodic linear dispersive equation whose dispersion relation is a multiple of a polynomial with integer coefficients (an "integral polynomial"), the prototypical example being the linearized Korteweg-deVries equation. Subsequently, it was found, [6], that the effect persists for more general dispersion relations which are asymptotically polynomial: ω(k) ∼ c k n for large wave numbers k 0, where c ∈ R and 2 ≤ n ∈ N. However, equations having other large wave dispersive asymptotics exhibit a wide variety of very different and surprising behaviors. The paper [6] also displayed some preliminary numerical computations that strongly indicate that the Talbot effect of dispersive quantization and fractalization persists into the nonlinear regime -for both integral and non-integrable evolution equations whose linear part has an integral polynomial dispersion relation.
The goal of the presents the study is to continue this investigationn in the context of two important examples: the nonlinear Schrödinger (nlS) and Korteweg-deVries (KdV) equations, possessing, respectively, elementary second and third order monomial dispersion. The basic numerical tool for this study is the operator splitting method, which also serves to highlight the interplay between the linear and nonlinear parts of the equation. Earlier rigorous results concerning the operator splitting method for the KdV, generalized KdV, and nlS equations can be found in [12,13,16]. See also [9] for additional results comparing the behavior of linear and nonlinear dispersive equations.

Operator Splitting
Let us first summarize the basic ideas behind the operator splitting method for approximating the solutions to nonlinear evolution equations. Details and proofs can be found in various references, for example [13].
Formally, consider the initial value problem for a dynamical differential equation: where K is a differential operator in the spatial variable with no explicit time dependence. We use Φ K (t) to denote the flow induced by the initial value problem, whereby u(t) = Φ K (t)u 0 is the solution to (2.1), which is defined on some time interval t ∈ [0, T ]. The time-independence of the right hand side and uniqueness of solutions to the initial value problem implies the flow condition Φ K (t + s) = Φ K (t) Φ K (s) for any t, s ≥ 0 such that the indicated flow operators are defined. Suppose that we can write K = L + N , where L, N are "simpler" differential operators, meaning that their corresponding flow operators Φ L (t), Φ N (t) can be determined explicitly, or are easily approximated by a numerical scheme. In our applications, L is a linear differential operator, while N is nonlinear. The basic idea of operator splitting is to approximate the flow of their sum K = L + N by repeated application of the individual simpler flows over a sequence of small time steps. For simplicity, we adopt a uniform time step 0 < ∆t 1, and seek to approximate the solution u(t n ) at the successive times t n = n ∆t for n = 0, 1, 2, . . . .
The simplest version is the Godunov splitting scheme, which iteratively defines the approximations In other words, 3) Sspecification of a suitable function space norm allows one to prove that, under appropriate hypotheses on the operators, K, L, N and the initial data, the Godunov approximations converge to the true solution: (2.4) A slightly more sophisticated method, which typically converges eat a faster rate, is the Strang splitting scheme, which replaces the iterative step in (2.2) by (2.5) We refer the reader to [13] for further details concerning operator splitting schemes in the context of the KdV equation.

Dispersive Quantization in Schrödinger equations
We first consider the free space linear Schrödinger equation We are interested in the behavior of rough solutions on a periodic domain, which we take to be 0 < x < 2 π. This is epitomized by considering the step function initial data in analogy with the Riemann problem for first order hyperbolic systems, [21,24]. It is straightforward to write an explicit Fourier series solution to the initial value problem (3.1), (3.2). Michael Berry, [1,2], was the first to discover the Talbot phenomenon of dispersive quantization and fractalization of this solution. Briefly, when the time is rational (relative to π), so t/π ∈ Q, the real part, complex part, and norm of the solution u(t, x) are discontinuous and piecewise constant, whereas when the time is irrational, t/π ∈ Q, the graphs of real part, complex part and norm of the solution are continuous but nowhere differentiable fractal-like functions. Additional details can be found in [14,18,19,20].
Our goal here is to investigate to what extent the dispersive quantization and fractalization phenomena persist into the nonlinear regime. We consider nonlinear Schrödinger equations of the form u t = i (u xx + u |u| p ), (3.3) subject to periodic boundary conditions on [0, 2 π], and use the same step function (3.2) as initial data. In (3.3), the power is assumed to satisfy p ≥ 1, with the particular case p = 2 being an integrable, soliton equation, of importance in fiber optics, [7,24]. In this paper, our numerical experiments will be displayed for the non-integrable case p = 1, although similar phenomena are observed for other powers. We will numerically approximate the solutions to the nonlinear Schrödinger equation (3.3) by use of operator splitting. We decompose its right hand side as the sum of its linear and nonlinear components: The resulting approximations to the solution has been previously analyzed in [13,16]. The flows of the linear dispersive operator L and the nonlinear operator N are easily found, and so we need to understand their interaction in order to work out the behavior of solutions to the full nonlinear initial value problem.
3.1. The linear Schrödinger equation. We begin by plotting some representative graphs of the solution to the periodic initial value problem (3.1), (3.2) for the linear Schrödinger equation. These were obtained by using the Fast Fourier transform (FFT) based on 2048 space nodes. Each row in Figure 3.1 displays the real part, the complex part, and the norm of the solution to the initial value problem at the indicated time. The first three times t = .3, .31, and .314 are irrational (relative to π) and increasingly close to the final rational time t = .1 π, at which time the solution can be rigorously proved to be piecewise constant, [17], although a small but noticeable Gibbs effect appears in the graphs due to the FFT approximation. On the other hand, at each irrational time, the solution is a continuous, but non-differentiable fractallike function, [18]. The sudden appearance of a quantization effect, but modulated by fractal behavior, between times t = .31 and .314 is striking.

The nonlinear equation.
Splitting off the nonlinear part of (3.3) leads to the evolution equation Because no x derivatives are involved, (3.4) is effectively an ordinary differential equation in which x serves as a parameter. Observe that and hence, for each fixed x, the modulus |u(t, x)| is constant as t evolves. This implies that the solution to (3.4) is explicitly given by is the initial condition.

Numerical Results.
In all the plots in Figure 3.2, we take the power p = 1 in (3.3), and use Godunov splitting to approximate the solution to the periodic initial value problem with step function initial data. (Although the Strang splitting method has higher rate of convergence, all the numerical results that it generated turned out to be extremely similar to those obtained via the simpler Godunov scheme, and so we will just display the latter approximations.) The nonlinear flow is computed exactly using (3.6), while the linear flow is computed through use of the FFT based on n = 2048 space nodes. The time step is ∆t = .00005 π and, as before, we display our numerical approximations at the irrational times t = .3, .31, .314 and the rational time t = .1 π. The plots are fairly similar to those from the linear evolution in Figure 3.1, although closer inspection will reveal small differences. The most noticeable is that, at the  rational time t = .1 π the quantized solution is no longer piecewise constant, and its shape exhibits the effect of the nonlinearity.
To analyze the numerical solution, we consider the phenomenon caused by the linear and nonlinear components separately. When time is a rational multiple of π, the linear evolution (3.1) will generate a piecewise constant solution profile. Because the nonlinear equation (3.4) is effectively an ordinary differential equation, with explicit soution formula (3.6), its effect is simply to shift the solution graph vertically by an amount depending upon the initial height, and hence preserves the piecewise constant nature of the solution. On the other hand, when the time is an irrational multiple of π, the linear term will generate a continuous but nowhere  differential profile. The subsequent vertical movement caused by the nonlinear term will be uniform, which means that, on a small time interval, the resulting graph will be very close to nowhere differentiable. The convergence of the operator splitting scheme then implies that the solution of the nonlinear Schrödinger equation will exhibit the same phenomenon. Proof. When u is sufficiently smooth, we have owing to the periodicity of our problem.
When the initial data u(0, x) = f (x) is not smooth, we can smooth it using mollification, and then apply the preceding argument. Specifically, let φ(x) be a standard mollifier, e.g. a Gaussian, and let u (0, x) be the mollified initial data obtained by convolution of u(0, x) and By the preceding calculation, the L 2 norm of the mollified solution u (t, x) is preserved. Moreover, by (3.8), the L 2 norm of the mollified solution converges to the L 2 norm of u(t, x), which completes the proof.
In particular, if the L 2 norm of the initial data of (3.1) is small, so is the L 2 norm of the solution, proving well-posedness.
Next we establish the well-posedness of the nonlinear equation. Proof. Suppose u(0, x) = f (x) and v(0, x) = g(x) are both in L 2 ∩ L ∞ . The corresponding solutions are given by formula (3.6). Thus, where M is a constant depending on the initial data, thus proving well-posedness. Proof. Because the linear term preserves the L 2 norm, we only need to control the growth of the norm due to nonlinear term. Fix a time step ∆t. To compute the approximate solution at the subsequent time t n = n ∆t, we must solve (4.3) n times, and hence (3.9) implies that the growth of the L 2 norm of the numerical solution is bounded by where M < ∞ by our assumption on the initial data. This clearly implies the well-posedness of the Godunov scheme.
As a consequence of Theorem 3 and the results in [3,13], we have the following: Let u and u ∆ denote, respectively, the exact and numerical solutions having the step function as initial data. Similarly, let u and u ∆, be, respectively, the exact and numerical solutions corresponding to the smooth initial data obtained by by mollifying the step function with φ . Then (3.10) When time step ∆t → 0 and → 0, the right hand side of (3.10) tends to 0. This implies that the operator splitting method for the periodic nonlinear Schrödinger equation, with step function as initial data, is convergent.

Dispersive Quantization for the Korteweg-deVries equation
Our second example is the Korteweg-deVries (KdV) equation, a justly famous integrable evolution equation modeling shallow water waves, plasma waves, etc., [7]. We will take the KdV equation in the normalized form u t + u xxx − uu x = 0, (4.1) and consider the periodic initial-boundary value problem on the interval [0, 2 π]. For the linearized KdV equation u t + u xxx = 0, (4.2) (also known as the Airy partial differential equation), the solution to the periodic initial value problem can be easily expressed in Fourier series form. In particular, the dispersive quantization and fractalization phenomenon for the step function initial data were established in [17,18]. 4.1. Numerical solutions. The authors' previous paper, [6], presented the results of numerical experiments on the periodic problem with the step function (3.2) as initial data, using operator splitting to numerically approximate the solution. It was found that the effects of dispersive quantization and fractalization persist into the nonlinear regime. In Figures 4.1 and 4.2 we display some representative graphs of the numerical solution at, respectively, some irrational times, which exhibit a fractal, non-differentiable profile, and at some rational times, whose profiles appear to be quantized and (almost) smoothly varying between the jumps. It is not clear whether the small oscillations appearing between the jumps are due to numerical error or the persistence of a small fractal contribution.

4.2.
Operator splitting. In order to apply operator splitting to the KdV flow, we decompose K[ u ] = − u xxx +uu x into the sum of its linear part L[ u ] = − u xxx , leading to the linearized KdV equation (4.2), and its nonlinear part N [ u ] = uu x , whose corresponding evolution equations is the nonlinear transport equation 3) also known as the Riemann equation or inviscid Burgers' equation, [24]. Unfortunately, we have so far been unable to rigorously establish convergence of the operator splitting method for the periodic KdV equation with rough initial data along the lines used above for the nlS equation.
Let us present what is known about this problem and some discussion of how one might proceed.
First of all, convergence for sufficiently smooth initial data has been established in [12].

4)
where C depends only on s, T , and the initial data u 0 .    The difficulty that is preventing us from extending this result to rough data, such as the step function, is fairly subtle. The key stumbling block is that one can control the L 2 behavior of the linear KdV flow (4.2), and the L 1 behavior of the nonlinear transport flow (4.3). More specifically, the difference between two solutions to the nonlinear component will behave well in L 1 whereas the difference between two solutions of the linear flow behaves well in L 2 . However, there is not an evident method that will allow us to interpolate between these two different Banach spaces.
Nevertheless, the following theorem concerning convergence of the operator splitting scheme to weak solutions has been proved in [13]. Moreover, Zhou, [25], proves that weak solutions to the KdV equation are uniquely determined by their initial data.
Now, if we know that our solution u L (t) has bounded L ∞ norm, then the estimate in Theorem 7 implies that the flow is L 1 contractive. Moreover, Oskolkov, [18], has established a global boundedness result for suitable initial data.
Theorem 8. If the initial data has bounded BV norm, then the resulting solution to the periodic linearized KdV equation is uniformly bounded in L ∞ , the bound depending on the BV and L ∞ norms of the initial data.
Furthermore, Oskolkov proves that when time is irrational (relative to π), the solution of the linearized KdV is nowhere differentiable, and has unbounded BV norm. The difficulty is that the L ∞ bound in Oskolkov's result depends upon the BV norm of the initial data, which can not be easily controlled when the initial data arises from application of the nonlinear transport flow.
On the other hand, consider the nonlinear transport equation (4.3), whose flow is denoted by Φ N (t). Of course, the formation of shock waves implies that Φ N (t) does not preserve smoothness or continuity of the initial data. In [12], the authors establish that the flow is both L 1 contractive and uniformly bounded.
Moreover, the solutions remain uniformly bounded in L ∞ norm: However, when attempting to interpolate these estimates between L 1 and L 2 , at each step in the iterative procedure we will introduce a factor that, unlike we observed in the nonlinear Schrödinger case, we are unable to control as the time step goes to 0 and the number of iterations goes to infinity. For example, although we can bound the nonlinear flow, the known BV norm bound is not easy to apply in the context of Oskolkov's Theorem 8.
Results on the decay of solutions to conservation laws due to Glimm and Lax, [10],imply the decay of the total variation of entropy solutions of nonlinear transport equation: Theorem 10. Suppose that the initial data to the nonlinear transport equation on a periodic domain satisfies u 0 ∈ L ∞ ([0, 2π]). Then the total variation (TV) of the resulting entropy solution u N (t) = Φ N (t)u 0 over the interval will satisfy the decay estimate: T V u N (t) < c/t, where c is a constant that is independent of the initial data.
However, as the time step ∆t → 0, we lose the control of the constant which will appear when we interpolate between L 1 and L 2 . At this point, we do not have any sharper estimate for the BV norm of the nonlinear flow, since, when the time step is irrational, the linear flow will generate the initial data for the nonlinear flow whose BV norm is unbounded. So the similar argument as Schrödinger case does not carry over to the KdV case.
Motivated by Theorem 5, we can try some other ways to show the compactness of the sequence of numerical solutions in L 2 . A standard approach is based on the Kolmogorov-Riesz L p compactness theorem, [11].
(1) F is bounded (2) for every > 0, there is an R > 0 such that |x|>R |f (x)| p dx < p , for every f ∈ F . (4.5) (3) for every > 0 there is a ρ > 0 such that R n |f (x + y) − f (x)| p dx < p , for every f ∈ F and y ∈ R n such that |y| < ρ. Our goal is to prove compactness in L 2 . In our case, condition 2 is redundant. Condition 1 is clear since boundness in L 2 is straightforward: The solution of nonlinear part is stable in L 2 and the L 2 norm of the solution of linear part is preserved. So, no matter many iterations we take, the family of numerical solutions will be uniformly bounded in L 2 . Thus, the remaining task is to verify condition 3. As we can see, when we iterate, the linear flow will behave well in L 2 but we are unable to rigorously control the nonlinear flow as before. On the other hand, our numerical experiments indicate that the numerical solutions are bounded in L ∞ , which will effectively make condition 3 immediate. Motivated by this fact, we can try to prove the uniform boundedness in L ∞ of the solution to the nonlinear part. However, the problem is that we do not have a Maximum Principle, and so cannot uniformly control the L ∞ norm for the linear flow when we perform the iterative step. Thus, we are left tantalizingly short of our goal to prove convergence of our numerical scheme for the KdV evolution of rough initial data. We hope to report further on this problem is a subsequent publication.

4.4.
Final remarks. Although we are as yet unable to prove the convergence of the numerical scheme, we can intuitively understand the behavior of the nonlinear KdV from the behavior of its linearization. The well-known studies of the "zero dispersion limit" initiated by Lax and Levermore, [15], are complemented by the following recent result of Erdoğan and Tzirakis, [8] that compares solutions to the linear and nonlinear equations.
Theorem 12. Consider the Korteweg-de Vries (KdV) equation with periodic boundary conditions. For H s initial data, s > − 1 2 , and for any r < min(3s + 1, s + 1), the difference between the nonlinear and linear evolutions is in H r for all times, with at most polynomially growing H r norm.
This result implies that we can regard the nonlinear KdV flow as a perturbation of the linearized KdV flow. Therefore, the observed persistence of the fractalization and quantization behaviors is to be expected.