Non-stiff methods for Airy flow and the modified Korteweg de Vries equation.

# Non-stiff methods for Airy flow and the modified Korteweg de Vries equation.

Mariano Franco-de-Leon Mariano Franco-de-Leon University of California, Irvine. 540 H Rowland Hall. Irvine, CA 92697-3875
Tel.: 949 307 6787
22email: mfrancod@uci.eduJohn Lowengrub University of California, Irvine. 540 H Rowland Hall. Irvine, CA 92697-3875
Tel.: 949-751-7557
44email: jlowengr@uci.edu
John Lowengrub Mariano Franco-de-Leon University of California, Irvine. 540 H Rowland Hall. Irvine, CA 92697-3875
Tel.: 949 307 6787
22email: mfrancod@uci.eduJohn Lowengrub University of California, Irvine. 540 H Rowland Hall. Irvine, CA 92697-3875
Tel.: 949-751-7557
44email: jlowengr@uci.edu
###### Abstract

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
###### Msc:
65Gxx,65Mxx,65Txx,35Q58,42Axx

\xpatchcmd\@thm

## 1 Introduction

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,

 kt=ksss+32k2ks. (1)

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

 Xt=Vn+Ts,

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 , :

 Xt=(−ks)n+(k22)s, (2)

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:

 Lt=0, (3)
 θt=(2πL)3[θααα+θ3α2]. (4)

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:

 Fˆfm=12π∫2π0f(x)e−imxdx,m=0,±1,±2,... (5)

Then, the Fourier series of is

 If(x)=m=∞∑m=−∞ˆfmeimx. (6)

Denote by , (for periodic functions of zero mean) the spectral derivative and integral operators used in this problem defined in Fourier space by:

 ˆShfm=imˆfm, (7)
 ˆInthfm={ˆfmim,if m≠0,0,if k=0. (8)

Observe that the linear part of equation is diagonalizable by the Fourier transform in the following way

 ∂Fˆθtm∂t+i(2πmL)3Fˆθtm=FˆNLtm, (9)

where is the length of the curve (constant), is the wavenumber, the super index represents time and

 NL(α,t)=((2πL)3θα(α,t)32). (10)

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

 ∂(Ψ(m,t))∂t=rtmFˆNLtm. (11)

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

 ˆftm=1NN/2∑k=−N/2+1f(αk,t)e−imαk;αk=kh, (12)

with inverse Fourier formula given by

 ftk=N/2∑m=−N/2+1ˆftmeimαk. (13)

Non-linear terms are treated in physical space and to avoid convolutions. In other words:

 ˆNLtm=1NN/2∑k=−N/2+1(Shθ)(3)(αk,t)e−imαk, (14)

Implicit time integration methods can now be easily applied.

Based on , the first step is computed using an Euler implementation and the integrating factor method:

 ˆθ1m=ζm(ˆθ0m+ΔtˆNL0m), (15)

where denotes the time step discretization and

 ζm:=exp(−iΔt(m2πL)3). (16)

 ˆθj+1m=ζmˆθjm+Δt2[3ζmˆNLjm−(ζm)2ˆNLj−1m]. (17)

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

 ˆθ1m=ˆθ0m+Δt(ˆL0m+ˆNL0m), (18)

and a Crank-Nicholson-like (CN) method for later steps of the form:

 ˆθj+1m−ˆθj−1m=Δt(ˆLinj+1m+ˆLinj−1m)+2ΔtˆNLjm, (19)

where .

Defining , and

 ζ1m:=1−γ2m1+γ2m+i−2γm1+γ2m=1−iγm1+iγm,ζ2m:=1−iγm1+γ2m, (20)

then, is equivalent to

 ˆθj+1m=ζ1mˆθj−1m+2Δtζ2mˆNLjm, (21)

for each wave number .

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

 ˆθ1m=ˆθ0m12[e−iγm+(1−iγm)]+ˆNL0mΔt2[1+e−iγm]. (22)

## 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:

 Cr[0,2π]:={f:first r % derivatives exist over (0,2π), are of bounded variation  over [0,2π], and whose first r−1 derivatives are 2π-periodic.} (23)

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))

 ˆfm=O(1mr+1), (24)

 |Shf(αi)−fα(αi)|=O(hr−1). (25)

Also, the accuracy of the trapezoidal rule can be estimated (Hammerlin and Hoffmann (1989)) by

 |N/2∑j=−N/2+1|f(αj)|h−∫π−πf(α)dα|=O(1Nr+1). (26)

Approximations are computed with the discrete inner products

 ⟨f,g⟩h:=N/2∑m=−N/2+1hfm¯¯¯¯¯¯gm, ⟨^f,^g⟩:=N/2∑m=−N/2+1^fm¯¯¯¯¯¯¯^gm, (27)

and the associated norms

 ||f||2l2=N/2∑m=−N/2+1|fm|2h, ||^f||2=N/2∑m=−N/2+1|^fm|2. (28)

An immediate consequence of the trapezoidal rule accuracy is

 ||f||2l2−||f||2L2=O(hr+1). (29)

The key ideas for treating stability error besides algebraic manipulation is Plancherel theorem

 12π⟨f,g⟩h=⟨^f,^g⟩⇒||ˆf||=1√2π||f||l2, (30)

that allow us to compute inner products at Fourier or Physical space interchangeably.

The main theoretical results in this section are the following theorems:

###### Theorem 1.

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,

 ||˜θj−θ(⋅,tj)||l2≤C(hr+Δt2). (31)
###### Theorem 2.

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,

 ||˜θj−θ(⋅,tj)||l2≤C(T)(hr+Δt2). (32)
###### Corollary 1.

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,

 ||˜θj−θ(⋅,tj)||l2≤C(T)(hr+Δt2). (33)

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,

 ||Slh(A−l(˙θj))||l2≤C||˙θj||l2,and||A−n(Slh˙θj)||l2≤C||˙θj||l2,for0≤l≤n.

In particular and . For estimates in time we write , for an operator satisfying

 ||A0(Δtn)||l2≤Δtn||f||l2,

where is integrable.

The proofs of theorems (1.) and (2.) are similar. We focus on CN discretization and refer the reader to the appendix for the ADB case.

###### Proof (Theorem 2.)

The error between numerical and exact solution (at a given time ) is denoted by:

 ˙θjm:=˜θjm−θ(αm,tj). (34)

Defining the auxiliary time

 T∗=Sup{t|t≤T,|˙L|

for . We aim to prove that the error of theta at the step also satisfies the estimate and this will imply by induction.

##### Taylor approximations:

for the first step of the induction argument, we calculate upper bounds for Euler step using the Taylor expansion:

 θ1m=θ0m+Δt(L0m+NL0m)+Δt22(θtt)0m+O(Δt3). (36)

Similarly, after the second step. The Crank-Nicholson discretization derived from the Taylor expansion of around and around has the form:

 ˆθj+1m−ˆθj−1m=2(Δtˆ(θt)jm+Δt36ˆ(θttt)jm)+O(Δt4),

where, as usual denotes the temporal derivative of . Using the approximation

 ˆLinj+1m+ˆLinj−1m2=ˆLinjm+Δt2ˆLinttjm+O(Δt5),

we obtain

 ˆθj+1m−ˆθj−1m= (37) 2(Δt(12(ˆLinj+1m+ˆLinj−1m)+ˆNLjm)+Δt36(ˆ(θttt)jm−6ˆLinttjm))+O(Δt4), (38)

where . This is equivalent to

 ˆθj+1m=ζ1mˆθj−1m+2Δtζ2mˆNLjm+Δt33((ˆθttt)j−6ˆLinttjm)+O(Δt4), (39)

where , for each wave number . The numerical solution satisfies

 ˆ˜θ1m=ˆ˜θ0m+Δt(ˆ˜L0m+ˆ˜NL0m), (40)

and

 ˆ˜θj+1m=ζ1mˆ˜θj−1m+2Δtζ2mˆ˜NLjm, (41)

for .

We start simplifying and noticing that

 (θt)t=(θsss+θ3s2)t=θssst+32θ2sθst=[θsss+θ3s2]sss+[θsss+θ3s2]s32θ2s, (42)

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:

 ((θt)t)t=((θsss+12θ3s)sss+32θ2s(θsss+12θ3s)s)t=θsssssst+12(θ3s)ssst+32θ2s(θsss+θ3s2)st+(θsss+12θ3s)s322θsθst= (43)
 θtssssss+32[θ2sθtssss+θts(θ2s)sss]+32θ2s[θtssss+3θsθssθst+32θ2sθtss]+θs[θssss+(θ2s2)s]2, (44)

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,

 ˆ˙θ1m=ˆ˙θ0m+Δt(ˆ˙L0m+ˆ˙NL0m)+A0(Δt2). (45)

Since the error at initial step is zero, Plancherel theorem shows

 ||ˆ˙θ1||2=O(Δt4). (46)

From , and the fact that are integrable, the error evolution after the second step is:

 ˆ˙θj+1m=ζ1mˆ˙θj−1m+2Δtζ2mˆ˙NLjm+A0(Δt3). (47)

To estimate the error consider the inner product:

 ⟨ˆ˙θj+1−ˆ˙θj−1,ˆ˙θj+1+ˆ˙θj−1⟩. (48)

A direct calculation shows how to rewrite this inner product as

 ⟨|ζ1|2(ˆ˙θj−1−ˆ˙θj−3),ˆ˙θj−1+ˆ˙θj−3⟩Jj1+4Δt2⟨|ζ2|2(ˆ˙NLj+ˆ˙NLj−2),ˆ˙NLj−ˆ˙NLj−2)⟩Jj2+⟨ζ1(ˆ˙θj−1−ˆ˙θj−3),2Δtζ2(ˆ˙NLj+ˆ˙NLj−2)⟩Jj3+⟨2Δtζ2(ˆ˙NLj−ˆ˙NLj−2),ζ1(ˆ˙θj−1+ˆ˙θj−3⟩Jj4+⟨ζ1(ˆ˙θj−1−ˆ˙θj−3)+2Δtζ2(ˆ˙NLj−ˆ˙NLj−2),A0(Δt3)⟩Jj5+⟨A0(Δt3),ζ1(ˆ˙θj−1+ˆ˙θj−3)+2Δtζ2(ˆ˙NLj+ˆ˙NLj−2)⟩Jj6. (49)

When taking the sum over time of the left-hand side we obtain a telescopic sum

 n∑j=2⟨ˆ˙θj+1−ˆ˙θj−1,ˆ˙θj+1+ˆ˙θj−1⟩=n∑j=2(||ˆ˙θj+1||2−||ˆ˙θj−1||2+2iIm(⟨ˆ˙θj+1,ˆ˙θj−1⟩))=||ˆ˙θn+1||2+||ˆ˙θn||2−(||ˆ˙θ2||2+||ˆ˙θ1||2)+I1, (50)

where is an imaginary term.

Now we analyze the sum over time of the right-hand side terms .

##### J1 contribution:

a direct calculation shows that

 Jj1=||ˆ˙θj−1||2−||ˆ˙θj−3||2+2iIm⟨ˆ˙θj−1,ˆ˙θj−3⟩.

Therefore, the sum over time is telescopic too

 n∑j=2Jj1=||ˆ˙θn−1||2+||ˆ˙θn−2||2−||ˆ˙θ1||2−||ˆ˙θ0||2+I2, (51)

where is a purely imaginary term.

##### J2 contribution:
similarly, and the sum over time is where is a purely imaginary term. And the following inequality holds (52)
##### J3+j4 contribution:
consider the sum (53)

Therefore, the sum over time is also a telescopic sum

 n∑j=2(J3+J4)j=2Re(⟨ζ1ˆ˙θn−1,2Δtζ2ˆ˙NLn⟩+⟨ζ1ˆ˙θn−2,2Δtζ2ˆ˙NLn−1⟩)+2Re(−⟨ζ1ˆ˙θ1,2Δtζ2ˆ˙NL2⟩−⟨ζ1ˆ˙θ0,2Δtζ2ˆ˙NL