Non-stiff methods for Airy flow and the modified Korteweg de Vries equation.
In this paper, we implement non-stiff interface tracking methods for the evolution of 2-D curves that follow Airy flow, a curvature-dependent dispersive geometric evolution law. The curvature of the curve satisfies the modified Korteweg de Vries equation, a dispersive non-linear soliton equation. We present a fully discrete space-time analysis of the equations (proof of convergence) and numerical evidence that confirms the accuracy, convergence, efficiency, and stability of the methods.
Keywords:dispersivesolitons numerical solution mKdV Airy flow
Geometric curve flow models are an important class of methods for interface motion; where we understand an interface as a geometrical one-dimensional surface with no thickness. Under these laws, curves evolve according to local functionals of their geometrical properties. A classical example is the mean curvature flow (Helfrich (1973), Hou et al. (2001), Tsai and Yue (1996)). The governing equations are parabolic partial differential equations. In the materials science context, mean curvature flows are related to the motion of grain boundaries that separate crystallites (grains) with different crystallographic symmetries. Another type of geometric evolution, where the governing equations of evolution are dispersive rather than parabolic, has been garnering increased attention. Dispersive equations arise in a variety of applications (collision-free hydromagnetic waves, ion-acoustic waves in cold plasma, electrostatic fields of graphene, human arm movement, computer vision (Smyth and Worthy (1995), Ho and Roy (2015), Miura (1976), Flash and Handzel (2007), Giblin and Sapiro (1998))), and their mathematical theories have revealed strong relations with differential geometry, geometrical analysis, soliton theory, and integrable systems (Miura (1976), Palais (1997), Terng (2014), Chou and Qu (2002), Colliander et al. (2004)).
In this article, we present the development, implementation, and analysis of schemes to obtain numerical (periodic in space) solutions for the modified Korteweg de Vries (mKdV) equation,
The mKdV equation is the first non-linear generalization of the KdV soliton model for water waves (1895). The contributions of Zabusky, Gardner, Green, Kruskal, Miura, Lax (1968) displayed their striking properties including: the preservation of form through non-linear interactions, decomposition of waves into smaller solitons, different families of solutions, infinite number of conservation laws, its relations with the Schrödinger operator and the eigenvalue problem (Miura (1976), Guan and Kuksin (2014)), the Miura transform to obtain solutions (and well-posedness) of KdV from solutions (and well-posedness) of mKdV (Miura (1976), Colliander et al. (2003), Gardner et al. (1974), Kevrekidis et al. (2004)), and the inverse scattering transform (IST). For other type of mKdV solutions (e.g. kinks, breathers), or periodic domains (also well-posed Colliander et al. (2004)), this approach is not plausible since decay at is a crucial hypothesis for IST. Other analytical techniques to find solutions of (periodic and non-periodic) mKdV-like equations include the use of Jacobi, Weierstrass functions, Hamiltonian structures, Bäcklund-Darboux transforms, the tangent hyperbolic method (Terng et al. (1997), Wang and Xiang (2013), Malfliet (2004), Deconinck and Nivala (2011), Lax (2005), Zheng et al. (2013)). Nevertheless, there is a lot of work to develop regarding the orbital stability of these periodic (and non-periodic) waves under perturbations of the underlying solution, soliton resolution conjecture, collisions, multisolitons, compactons, generalizations (gKdV), nonlinear Schrodinger-Airy system, new solutions, and its relations with other equations (Colliander et al. (2003), Rosenau and Hyman (1993), Tao (2008), Bonanno (2015), Mousavian et al. (2011), Song (2012), Pava (2009), Kodama and Hasegawa (1987), Guo (2009)).
Over the past 30 years, the numerical study of initial value problems of free surface flows has been increasingly important in representing systems of partial differential equations, not just for physical modeling, but also as an empirical tool to analyze theoretical aspects of the underling system. The primary classes of algorithms (Lagrangian and Eulerian), as well as mixed approaches, had been focus on the solution of parabolic (dissipative case) partial differential equations (Smyth and Worthy (1995), Palais (2014), Burchard et al. (2001), Liu (2014), Helal and Mehanna (2007), Benson (1992), Chang et al. (1996), Tsai and Yue (1996), Li et al. (2009), Rossman and Boulos (1996), Leung et al. (2011)). There are far fewer methods (Adomain decomposition, finite differences, radial basis functions, pseudo-spectral methods) developed for simulating dispersive geometric evolution equations (Smyth and Worthy (1995), Palais (2014), Helal and Mehanna (2007), Yagmurlu et al. (2016), Dağ and Dereli (2008)).
In this paper, we exploit the theory behind dispersive equations and geometric curve flows by evolving solutions of a closed curve under Airy flow. Then we recover mKdV solutions from the curvature of the curve (See ), instead of solving mKdV directly, and gaining one degree of smoothness in the numerical implementation. The evolution of any -D closed smooth planar curve with spatial ( periodic) parameter , time variable , can be described as
where s, n denote the tangent and outward-normal unit vectors respectively, and are the corresponding normal and tangential velocities. In Airy flow the normal and tangential velocities are , :
where is the curvature along the curve, denotes the arc-length parameter, and subscripts represent partial differentiation.
The high number of spatial derivatives, nonlinearity, and dispersive effects represent particular challenges when solving these equations numerically. Explicit time stepping methods undergo severe time constraints. In addition, certain spatial discretizations may lead to numerical instabilities. As observed previously (Ceniceros and Hou (1998), Robertson and Sherwin (1999), Beale et al. (1996), Beale et al. (1994)), even spectral accuracy does not guarantee stability. Further, time-step constraints may be amplified during the evolution due to clustering of points at the interface. The tangential velocity for Airy flow enforces equal arc-length parametrization at all times provided it is satisfied at the initial step. In this way, is everywhere equal to its mean and evolves according to the length of the curve, a uniform discretization in is then uniform in (i.e. ). Numerically, this choice of frame avoids the time-step restrictions for stability due to clustering of grid points at the interface. Another feature for curvature-dependent problems is the relation , between the curvature and , the angle that makes the tangent vector and the -axis ()). Using the arc-length parameter and , , (- formulation), as dynamical variables (Baker and Shelley (1990)) instead of coordinates, equation becomes:
We can then obtain by integrating the expression , Hou et al. (1994), and recover solutions of the mKdV equation from .
The linear term of the equation displays the reason of stiffness whose stability constraint for an explicit method has the form where and is the grid spacing in . A stable and accurate discretization must guarantee a perfect balance between nonlinear and dispersive effects. We use the small-scale decomposition (SSD) of the equations, developed by Lowengrub, and Shelley (HLS) Hou et al. (1994) to examine the source of stiffness at small scales at which curvature acts as a linear operator.
Linear analysis and numerical conservation of first integrals of motion (conservation of mass, momentum, and energy for the problem of the real line (Miura (1976), Miura et al. (1968), Dingemans (1968))) for mKdV equation are used to test the accuracy of the numerical methods. Semi-discrete (continuous time) analysis (e.g. Beale et al. (1996), Ceniceros and Hou (1998)) suggested that numerical filters need to be used to overcome instabilities generated by truncation and aliasing errors arising when computing spatial derivatives (Krasny (1986)). In contrast, our fully discrete space-time analysis of convergence demonstrates that the use of the filter is not related to convergence, but may enhance stability.
The paper is organized as follows: In section 2 we describe the numerical schemes used to treat the nonlinear dispersive equation (4) needed to evolve Airy flow. Our most important theoretical (proof of) convergence results are given in section 3. As a first accuracy test, linear versions of the solution for Airy flow and mKdV are derived in section 4 and compared against the numerical solutions. Additionally, numerical results including accuracy, convergence, stability, dynamics and the use of filters, is covered on section 5. Concluding remarks are given in section 6, and technical computations in the appendix A.
2 Numerical Methods
Next, we introduce the notation to describe the schemes and the convergence analysis (Beale et al. (1996), Ceniceros and Hou (1998), Canuto et al. (1988)): arbitrary smooth functions are expressed by , and constants (independents of discretization) are written generically as . For a complex valued function defined over , the (continuous) Fourier coefficients of are:
Then, the Fourier series of is
Denote by , (for periodic functions of zero mean) the spectral derivative and integral operators used in this problem defined in Fourier space by:
Observe that the linear part of equation is diagonalizable by the Fourier transform in the following way
where is the length of the curve (constant), is the wavenumber, the super index represents time and
Consider a linear propagator method to absorb the leading order (linear term) prior to discretization. Several researchers in different contexts have used linear propagator schemes, e.g. simulations for Navier-Stokes equations, Hele-Shaw flows, reaction-diffusion systems, multicomponent fluids, multiphase materials (Hou et al. (2001), Hou et al. (1994), Crapper (1970), Leo et al. (2000), Nie et al. (2006)) to name a few. Consider the integrating factor and the function , thus equation is equivalent to
These formulation motivates the use of Discrete Fourier Transform (DFT). In parallel with the continuous case ,, given a periodic function , whose values are known on a uniform grid of mesh size ( is a power of two), the -th discrete Fourier coefficients of are defined as
with inverse Fourier formula given by
Non-linear terms are treated in physical space and to avoid convolutions. In other words:
Implicit time integration methods can now be easily applied.
2.1 Linear propagator method and Adams-Bashforth (ADB)
Based on , the first step is computed using an Euler implementation and the integrating factor method:
where denotes the time step discretization and
Subsequent steps are calculated with the second order Adams-Bashforth (ADB) method:
Notice how at the th time-step is propagated forward to the next step at the exact exponential rate associated with the linear term. If , this yields to the exact solution of the linear problem. In the case of Airy Flow, the length of the curve is constant, thus is constant over time.
2.2 Crank-Nicholson (CN):
It is also possible to discretize using an Euler discretization for the first step
and a Crank-Nicholson-like (CN) method for later steps of the form:
Defining , and
then, is equivalent to
for each wave number .
2.3 Crank-Nicholson and Adams-Bashforth (CNADB)
The scheme CNADB is a modification of CN, where the first step after initialization is the average of the schemes used for CN and ADB discretizations, that is
3 Analytical convergence
To prove the convergence of the presented schemes, we denote to the exact continuous solution evaluated at the grid points by , and we use for the discrete approximations. Purely imaginary terms are denoted by . For simplicity, we omit the time notation where the specific time is not relevant for the computation.
We work with the following space of functions:
The existence of the first derivatives is understood in the almost everywhere Riemann-Stieltjes sense (Canuto et al. (1988)).
The main tool to handle truncation error is the spectral accuracy of the method. In other words, the Fourier coefficients of any satisfies the decay condition (Canuto et al. (1988))
which implies (Tadmor (1987))
Also, the accuracy of the trapezoidal rule can be estimated (Hammerlin and Hoffmann (1989)) by
Approximations are computed with the discrete inner products
and the associated norms
An immediate consequence of the trapezoidal rule accuracy is
The key ideas for treating stability error besides algebraic manipulation is Plancherel theorem
that allow us to compute inner products at Fourier or Physical space interchangeably.
The main theoretical results in this section are the following theorems:
Assume that for there exists a regular solution of the system of evolution equations and (for Airy flow and the mKdV equation) with belonging to for and whose second derivative is continuous with respect to time. If denotes the numerical solution obtained with the scheme ,, then for and we have,
Assume that for there exists a regular solution of the system of evolution equations and (for Airy flow and the mKdV equation) with belonging to for and whose third derivative is continuous with respect to time. If denotes the numerical solution obtained with the scheme , then for and we have,
Assume that for there exists a regular solution of the system of evolution equations and (for Airy flow and the mKdV equation) with belonging to for and whose third derivative is continuous with respect to time. If denotes the numerical solution obtained with the scheme ,, then for and we have,
At this point, we introduce the notation that will be used in the error analysis. We define the discrete -th order smoothing operator, written generically , as an operator satisfying,
In particular and . For estimates in time we write , for an operator satisfying
where is integrable.
Proof (Theorem 2.)
The error between numerical and exact solution (at a given time ) is denoted by:
Defining the auxiliary time
for . We aim to prove that the error of theta at the step also satisfies the estimate and this will imply by induction.
for the first step of the induction argument, we calculate upper bounds for Euler step using the Taylor expansion:
Similarly, after the second step. The Crank-Nicholson discretization derived from the Taylor expansion of around and around has the form:
where, as usual denotes the temporal derivative of . Using the approximation
where . This is equivalent to
where , for each wave number . The numerical solution satisfies
We start simplifying and noticing that
involves spatial derivatives of order for theta, and we have used the fact that is at least two times continuously differentiable with respect time to commute derivatives. Similarly, we compute:
note that involves derivatives of order in space. In the expression the term contains the most (9 to be precise) derivatives for , which by hypothesis , we know these are integrable. Using computation , observe that also involves spatial derivatives of order 9 for .
If is substracted from and using the fact that is integrable we obtain the equation for the first step,
Since the error at initial step is zero, Plancherel theorem shows
From , and the fact that are integrable, the error evolution after the second step is:
To estimate the error consider the inner product:
A direct calculation shows how to rewrite this inner product as
When taking the sum over time of the left-hand side we obtain a telescopic sum
where is an imaginary term.
Now we analyze the sum over time of the right-hand side terms .
a direct calculation shows that
Therefore, the sum over time is telescopic too
where is a purely imaginary term.
contribution:similarly, and the sum over time is where is a purely imaginary term. And the following inequality holds (52)
contribution:consider the sum (53)
Therefore, the sum over time is also a telescopic sum