Eventual periodicity of the forced oscillations for a Korteweg–de Vries type equation on a bounded domain using a sinc collocation method

Accepted Manuscript Eventual periodicity of the forced oscillations for a Korteweg–de Vries type equation on a bounded domain using a sinc collocation method Kamel Al-Khaled, Nicholas Haynes, William Schiesser, Muhammad Usman

PII: DOI: Reference:

S0377-0427(17)30412-0 http://dx.doi.org/10.1016/j.cam.2017.08.023 CAM 11275

To appear in:

Journal of Computational and Applied Mathematics

Received date : 27 August 2015 Revised date : 25 May 2017 Please cite this article as: K. Al-Khaled, N. Haynes, W. Schiesser, M. Usman, Eventual periodicity of the forced oscillations for a Korteweg–de Vries type equation on a bounded domain using a sinc collocation method, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.08.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Eventual periodicity of the forced oscillations for a Korteweg-de Vries type equation on a bounded domain using a sinc collocation method Kamel Al-Khaled1 , Nicholas Haynes2 , William Schiesser3 , and Muhammad Usman∗4 1

Department of Mathematics and Statistics, Jordan University of Science and Technology, Irbid, 22110, Jordan 2 Department of Physics, Duke University, Durham, NC 27708 3 Lehigh University, Bethlehem, PA 18015 4 University of Dayton, 300 College Park, Dayton OH 45469-2316 September 10, 2017

Abstract We demonstrate numerically the eventual time periodicity of solutions u(., t) to the Korteweg-de Vries type equation with periodic forcing at one end using the sinc-collocation method. This method approximates the space dimension of the solution with a cardinal expansion of sinc functions, thus allowing the avoidance of a costly finite difference grid for a third order boundary value problem. The first order time derivative is approximated with a θ-weighted finite difference method. The sinc-collocation method was found to be more robust and more efficient than other numerical schemes when applied to this problem.

Keywords: KdV equation, eventual periodicity, Sinc collocation method.

1

Introduction

The Korteweg-de Vries (KdV) equation was first derived by Boussinesq [10] in 1870, and then by Korteweg and de Vries in 1895 [17], as a model for long-crested small-amplitude long waves propagating on the surface of water. The same partial differential equation serves as a model for unidirectional ∗

Corresponding author: [email protected]. Authors names are in alphabetical order.

1

propagation of waves in a variety of physical systems. In the last century, the study of these types of waves, sometimes called solitons, has developed into a rich area of research in physics, applied mathematics, and related disciplines [4, 5, 6, 13, 16, 19]. Besides fluid dynamics, the KdV equation has found applications describing nonlinear heat transfer in crystal lattices [12], ion-acoustic waves in plasmas [33], blood flow in arteries [15], and even cosmology [18]. In [4], the authors modeled a wave-maker situation using the initial-boundary-value problem (IBVP) of KdV type equation ηt + ηx + ηηx + ηxxx = 0, x, t ≥ 0, η(x, 0) = 0, η(0, t) = h(t),

(1)

where x = 0 corresponds to the position where the wave-maker is mounted. There were several practical assumptions including perfect fluid, channel has been infinitely away from the wave-maker hence ignoring the wave reflections at the end of the channel. In practice, this is similar to a situation in which a beach is very gently sloping so that little or no energy comes back, or, in a channel, the experiment takes place over a time scale such that the wave motion does not reach the end of the channel. If there is significant reflected wave motion moving back toward the wave maker, this model is not appropriate, as this is model for unidirectional wave propagation. To take account of two-way propagation at the KdV level of approximation, a Boussinesq system of equations would be needed as in [9]. Bona, Sun and Zhang [8] studied a damped KdV type equation on a quarter plane domain  ut + ux + uux + uxxx − αuxx + γu = 0, 0 < x < ∞, (2) u(0, t) = h(t), t > 0, where γu and αuxx are damping terms. They proved that for γ > 0 and under the assumption of smallness of h(t) there exists a unique time periodic solution that is globally exponentially stable. Further, the eventual periodicity of u(., t) was also shown, that is ||u(., t) − u(., t + T )|| → 0 as t → ∞, where T is the period of boundary forcing h(t). It was shown in [29, 30] that the solution of the IBVP   ut + ux + uux + uxxx = 0, 0 < x < 1, t ≥ 0 u(x, 0) = φ(t),  u(0, t) = h(t), u(1, t) = 0, ux (1, t) = 0,

(3)

with h(t) periodic such that h(t) = h(t+T ) exhibits asymptotic time periodic behavior at any fixed spatial point, assuming that the amplitude of the boundary forcing term h(t) is small. So the wavemaker is putting energy in 2

a finite channel from the left boundary (x = 0) while the right end (x = 1) of the channel is free. One may question on the need of dissipation in the model. In fact, if h(t) = 0 then by multiplying the above equation with u and integration by parts yields Z d 1 2 u (x, t) dx = −u2x (0, t), dt 0 for any t ≥ 0. This shows that there is a weak dissipation from the left end of the boundary, although there is no explicit dissipation in the system. Recently, in [31] the authors have shown eventual periodicity of (3) using the Adomian decomposition method. In this paper, we demonstrate this behavior numerically using the sinc-collocation method. We consider two similar equations, the KdV equation with periodic forcing at the boundary [29, 30],   ut + αux + βuux + γuxxx = 0, 0 < x < 1, t > 0, u(x, 0) = 0, 0 < x < 1, (4)  u(0, t) = h(t), u(1, t) = 0, ux (1, t) = 0, t > 0, and the KdV equation with periodic external forcing,   ut + αux + βuux + γuxxx = h(t), 0 < x < 1, t > 0, u(x, 0) = 0, 0 < x < 1,  u(0, t) = 0, u(1, t) = 0, ux (1, t) = 0, t > 0.

(5)

Among the computational methods for the IBVPs of nonlinear partial differential equations (NPDEs) the popular ones are finite difference (FD), finite element (FE), spectral methods, and collocation methods. Spectral methods are complicated to implement and discretization of NPDEs leads to the solution of large system of linear or nonlinear equations involving a full matrix [11]. Spectral methods cannot represent the physical process in spectral space, and are also limited to IBVPs with periodic boundary conditions. Recently, there has been considerable work using radial basis functions [14] and sinc-collocation methods to solve nonlinear iIBVPs. Papers [28, 20] consider versions of the regularized long wave equation for which analytical solutions are not known. A good comparison of the performance of the sinccollocation method to the exact solutions of nonlinear PDEs can be found in [3] and [21]. The sinc-collocation method is used to examine the behavior of a classic version of the KdV equation with no forcing in [22]. Finally, periodic behavior of the KdV equation with small amplitude forcing at the boundary was observed numerically using the dual Petrov-Galerkin method in [24]. Other significant work on similar IBVPs of the KdV type equation is discussed in [32] and in [25] on bounded domain with zero boundary conditions. Most of the numerical techniques in the literature are either for the 3

initial value problem of KdV equation or on the unbounded domain. In this work we discuss an initial and boundary problem of KdV type problem. We present below a use of the sinc-collocation method to observe timeperiodic behavior of the KdV with forcing of arbitrary amplitude, and a detailed stability analysis of the problem. In section 2, we provide preliminary information about the sinc function and how it is used to approximate the solution of a partial differential equation. In section 3, we present the stability analysis, deriving necessary and sufficient bounds on the weight θ between the discretized time steps and confirming these bounds numerically. Finally, section 4 contains numerical results, including an implementation of the sinc method to a version of the KdV equation with a known exact solution and the calculation of conserved quantities to demonstrate the accuracy of the sinc method and as a way to compare it to other methods. In this section we also apply the technique to our problem (4), in order to demonstrate the eventual periodicity. We conclude with a brief discussion in section 5.

2

The sinc function and sinc-collocation method

The sinc function is commonly defined as  1 if z = 0, sinc(z) = sin(πz) if z = 6 0, πz

(6)

where z ∈ C. For a series of nodes equally spaced δ apart, the sinc function can be written as:   x − jδ S(j, δ)(x) = sinc . (7) δ Following the work done by Stenger [26, 27], the Whittaker cardinal function C(f, δ) of a function f is defined as C(f, δ)(x) =

∞ X

f (jδ)S(j, δ)(x),

(8)

j=−∞

granted the sum converges. Note that as C(f, δ) was invented as way of describing entire functions, it is exact if f ∈ B2 (Dd ). For practical purposes we choose the truncated sum for some j = −N, . . . , N , where N is the number of sinc grid points. Definition 1. Let d > 0 and let Dd denote the strip of width 2d about the real axis: Dd = {z ∈ C I : |=(z)| < d}. 4

In addition, for  ∈ (0, 1), let Dd () denote the rectangle in the complex plane: Dd () = {z ∈ C : |<(z)| < 1/, |=(z)| < d/(1 − )}. Let B2 (Dd ) denote the family of all functions g that are analytic in Dd , such that Z d |g(x + iy)|dy → 0, x → ∓∞, −d

and such that N2 (g, Dd ) = lim

→0

Z

∂Dd ()

1/2 |g(z)|2 |dz| < ∞.

Now, in order to approximate the solution to (4), the derivatives of S(j, δ)(x) must be calculated. Referring to [27], any derivative can be calculated according to (2r)

Sj

(xi ) =

and (2r+1)

Sj

d2r S(j, δ)(x) dx2r x=xi 
π 2r (−1)r if j = i,   δ 2r+1 ,   =     (−1)i−j Pr−1 (−1)l+1 2r! π 2l (i − j)2l , if j 6= i, l=0 (2l+1)! δ 2r (i−j)2r

(xi ) =

=

    

d2r+1 S(j, δ)(x) 2r+1 dx x=xi

0,

(−1)(i−j) δ 2r+1 (i−j)2r+1

(9)

if j = i, (10)

Pr

l (2r+1)! 2l 2l l=0 (−1) (2l+1)! π (i − j) , if j 6= i,

where r = 0, 1, 2, . . . and xi are the collocation points such that xj = a + (j − 1)δ (assuming x1 = a). For the purpose of applying this theory to the solution of the KdV equation, the derivatives of interest are  1 if j = i, (0) Sj (xi ) = (11) 0 if j 6= i,   0, if j = i,  (1) (12) Sj (xi ) =   (−1)i−j if j 6= i, δ(i−j) 5

(3) Sj (xi )

=

    

0, (−1)i−j δ 3 (i−j)3

if j = i,

  6 − π 2 (i − j)2 . if j 6= i.

(13)

To implement the sinc-collocation method, the time derivative is discretized with a θ-weighted scheme. This turns (4) into h i uk+1 − uk k+1 k+1 k+1 + θ (α + βu )ux + γuxxx δt h i + (1 − θ) (α + βuk )ukx + γukxxx = 0, (14)

where uk = u(x, t = tk ) is the value of the solution at the k th time step. A total of N collocation points are used in the spatial dimension such that |b − a| x0 = a, xN = b, and xi = a + (i − 1)δx , where δx = . The solution N can then be interpolated as k

u(xi , tk ) ≡ u (xi ) ≈

N X

ukj Sj (xi ).

(15)

j=1

The term uk+1 uk+1 on the left-hand side of (14) must be linearized before x continuing. This is accomplished by multiplying the Taylor expansions of uk+1 and uk+1 and disregarding terms O(δt2 ) and above. Namely, x uk+1 = u(tk + δt ), = u(tk ) + δt ut (tk ) + O(δt2 ),

(16)

and uk+1 = ux (tk + δt ), x = ux (tk ) + δt uxt (tk ) + O(δt2 ),

(17)

so h ih i 2 2 uk+1 uk+1 = u(t ) + δ u (t ) + O(δ ) u (t ) + δ u (t ) + O(δ ), t t k x k t xt k k x t t ! ! uk+1 − uk uk+1 − ukx x = uk ukx + δt ukx + δt uk + O(δt2 ), δt δt ≈ uk+1 ukx + uk uk+1 − uk ukx . x

(18)

Substituting (15) and (18) into (14) yields the following system of equa-

6

tions; N X

(0) uk+1 Sj (xi ) j

j=0

+ θδt β

=

"

N X

+ δt θα

N X

(1) uk+1 Sj (xi ) j

+ δt θγ

j=0

N X

(0) ukl Sl (xi )

N X

(0) ukj Sj (xi )

j=0

(1) uk+1 Sj (xi ) j

+

N X

(0) uk+1 Sj (xi ) j

j=0

− δt (1 − θ)α

N X

(1) ukj Sj (xi )

j=0

+ δt (2θ − 1)β for i = 1, . . . , N − 1, and

(3)

uk+1 Sj (xi ) j

j=0

j=0

l=0

N X

PN

N X

#

(1) ukl Sl (xi )

l=0

− δt (1 − θ)γ (0)

ukj Sj (xi )

j=0

k (0) j=0 uj Sj (x0 )

N X

N X

N X

(3)

ukj Sj (xi )

j=0

(1)

ukl Sl (xi ), (19)

l=0

= h(tk ),

PN

k (0) j=0 uj Sj (xN )

(20)

= 0,

PN

k (1) j=0 uj Sj (xN )

= 0.

In the interest of succinctness, we now switch to a matrix representation of (19) and (20), where U n = [uk1 , uk2 , ..., nkN ]T , S (0) = (Sj (xi )), i, j = 1, ..., N, S (1) = (Sj0 (xi )), i, j = 1, ..., N, (21) S (3) = (Sj000 (xi )), i, j = 1, ..., N, N1k = Uk ∗ S (1) , N2k = Ukx ∗ S (0) , where ∗ stands for component by component multiplication. Equations (19) and (20) now become the following matrix equation " #   S (0) + δt θαS (1) + δt θβ N1k + N2k + δt θγS (3) uk+1 "

= S

(0)

− (1 − θ)δt αS

(1)

+ (2θ − 7

1)δt βN1k

− (1 − θ)δt γS

(3)

#

uk ,

or Auk+1 = Buk ,

(22)

where, h i A = S (0) + αδt θS (1) + βδt θ[Un ∗ S (1) + Unx ∗ S (0) ] + δt θγS (3) ,

h i B = S (0) − α(1 − θ)δt S (1) + (2θ − 1)δt βUn ∗ S (1) − (1 − θ)δt γS (3) .

Using the above equations, we arrive at

Un+1 = I (0) A−1 Bun .

(23)

We can obtain the coefficients of the approximate solution by solving the system in equation (23) using any iterative technique. For the convergence of the method, we state the following two theorems, and the proofs are very similar to the proofs of the results in [1, 2]. Theorem 1. Let the function u(x, t) be as in equation (1) or (2), and let the matrix U be defined as in (23). Then for sufficiently large N , there exists a constant C independent of N such that √ sup k u(xi , tk ) − U k≤ CN 2 exp(− πN ). (xi ,tk )

Theorem 2. Given a constant R > 0, there is a constant T0 > 0 such that if k U1 − U0 k≤ R/2, then the iterative scheme (23) converges to the unique solution.

3

Stability analysis

We present in this section an analysis of the stability of the sinc method for solving the initial and boundary value problem of KdV equation (14). Using the method outlined in Section 2, the evolution of error can be written as " #   S (0) + δt θαS (1) + δt θβ N1k + N2k + δt θγS (3) ek+1 "

#

= S (0) − (1 − θ)δt αS (1) + (2θ − 1)δt βN1k − (1 − θ)δt γS (3) ek , (24) where S (0) , S (1) , S (3) , and N2k are defined by (21), and the error is defined as: ek = |ukexact − ukapprox |, (25)

where ukexact and ukapprox are the exact and sinc-approximated solutions at time tk , respectively. 8

Equation (24) can be written as ek+1 = C ek ,

(26)

where C≡A

−1

"

B= S "

× S

(0)

(0)

+ δt θαS

(1)

− (1 − θ)δt αS

+ δt θβ

(1)



N1k

+ (2θ −

+

N2k



1)δt βN1k

+ δt θγS

(3)

#−1

− (1 − θ)δt γS

(3)

#

. (27)

The scheme is considered numerically stable if ρ(C) ≤ 1, where ρ(·) denotes the spectral radius. Stability is assured if 1 − (1 − θ)δ (αλ + γλ ) + δ (2θ − 1)βλ t 1 3 t N1 (28) ≤ 1, 1 + θδt [αλ1 + γλ3 + β(λN1 + λN2 )]

for all possible values of λ0 , λ1 , λ3 , λN1 , and λN2 , where λ0 , λ1 , λ3 , λN1 , and λN2 are eigenvalues of S (0) , S (1) , S (3) , N1 , and N2 , respectively. Referring to Section 2, one can see that S (0) is just the identity matrix, so λ0 = 1. Furthermore, {S (1) }j,i =

(−1)i−j (−1)j−i =− = −{S (1) }i,j , δx (j − i) δx (i − j)

(29)

and {S (3) }j,i =

(−1)j−i [6 − π 2 (j − i)2 ] δx3 (j − i)3 (−1)i−j =− 3 [6 − π 2 (i − j)2 ] = −{S (3) }i,j , (30) δx (i − j)3

and {S (1) }i,i = {S (3) }j,j = 0. S (1) and S (3) are therefore skew-symmetric and thus have pure imaginary eigenvalues, i.e. λ1 = i|λ1 | and λ3 = i|λ3 |. In general, λN1 and λN2 are complex. To reduce notational complexity, we deI fine λ ≡ δt (α|λ1 | + γ|λ3 |), as well as λR N1 ≡ δt β Re(λN1 ), λN1 ≡ δt β Im(λN1 ), R I λN2 ≡ δt β Re(λN2 ), and λN2 ≡ δt β Im(λN2 ). Now, (28) can be rewritten 1 + (2θ − 1)λR + i(θ − 1)λ + (2θ − 1)λI  N1 N1 (31) ≤ 1, R I I 1 + θ(λR N1 + λN2 ) + iθ(λ + λN1 + λN2 ) which implies that   I + i (θ − 1)λ + (2θ − 1)λ 1 + (2θ − 1)λR N1 N1 R I I ≤ 1 + θ(λR + λ ) + iθ(λ + λ + λ ) N1 N2 N1 N2 ; (32) 9

i.e., i2 h i2 h + (θ − 1)λ + (2θ − 1)λIN1 1 + (2θ − 1)λR N1 h i2 h i2 R I I . (33) ≤ 1 + θ(λR + λ ) + θ(λ + λ + λ ) N1 N2 N1 N2

After algebraic manipulation, this condition becomes I 2(θ − 1)λR N1 + (2θ − 1)(θ − 1)λλN1

2 I R R I I − 2θλR N2 − 2θ (λλN2 + λN1 λN2 + λN1 λN2 )

≤ (3θ − 1)(1 − θ)|λN1 |2 + θ2 |λN2 |2 + (2θ − 1)λ2 . (34)

This condition must hold for all eigenvalues of the respective matrices for the method to be stable. Notice that if 1/2 ≤ θ ≤ 1, the right-hand side of (34) is nonnegative, whereas the lefthand side may run negative, depending on the choice of eigenvalues. We conclude that the condition 1/2 ≤ θ ≤ 1 is necessary, but not sufficient, for stability of the sinc collocation method for (4). To find a sufficient condition, we examine the method when θ = 1. In this case, (31) becomes I 1 + λR N1 + iλN1 (35) ≤ 1; R R I I 1 + λN1 + λN2 + i(λ + λN1 + λN2 ) i.e.,

R I R I I 1 + λR N1 + iλN1 ≤ 1 + λN1 + λN2 + i(λ + λN1 + λN2 ) R I I ≤ 1 + λR N1 + iλN1 + λN2 + i(λ + λN2 ) , (36)

where the right-most side of (36) is obtained by the triangle inequality. I Eliminating common terms, stability requires that |λR N2 + i(λ + λN2 )| ≥ 0, which is trivially true for any choice of eigenvalues. Thus we conclude that the condition θ = 1 is sufficient for stability of the method. To demonstrate the stability conditions, we compute the spectral radius of the iteration matrix C, defined by (27), for values of θ between 0 and 1. Recall that stability of the method is guaranteed if ρ(C) ≤ 1. Figure 1 shows that the condition 1/2 ≤ θ ≤ 1 is very nearly sufficient. Indeed, numerical calculations show that ρ(C) − 1 < 10−3 for θ > 0.5001, and that ρ(C) = 1 to within machine precision for θ = 1. For θ = 0 the method is conditionally stable with the condition on time step: δt ≤

β 2 Re2 (λ

2βRe(λN1 ) 2 N1 ) + [(α|λ1 | + γ|λ3 |) + βIm(λN1 )] 10

1.5

ρ(C)

1

0.5

0

0

0.2

0.4

0.6

0.8

1

θ

Figure 1: Demonstration of convergence dependence on θ

4

Numerical results

Before presenting the numerical solution to (14), we test the sinc method described in the previous section on a version of the KdV equation with a known analytic solution. Namely, the problem ut + uux + uxxx = 0, u(x = 0) = u(x = L) = ux (x = L) = 0, u(t = 0) = Asech2 (κx − x0 ),

(37)

where A = 12κ2 , has the solution u(t, x) = Asech2 (κx − ωt − x0 ), where ω = 4κ3 .

(38)

In the spirit of [21], we choose κ = 0.3 and x0 = 0, and calculate the solution on the space and time intervals [−30, 30] and [0, 10], respectively, with N = 100 and nt = 1000. A surface plot of this numerical solution is shown in figure (2). The processing took approximately 3.49 seconds on an Intel Core i7 CPU running at 2.8 GHz. The accuracy of the approximation shown in figure (2) is evaluated two different ways. First, because the exact solution to (37) is known (namely, (38)), the error at each grid point can be calculated directly. A surface plot of the error, e(x, t) = uexact (x, t) − uapprox (x, t), is shown in figure(3). The error at any grid point was less than 8 × 10−7 , and ||e||∞ = 2.60 × 10−4 . In addition, the condition number of the inverted matrix, cond(A), was less than 1.33 for every time step, indicating that A is exceptionally wellconditioned. Following [21], another way of evaluating the accuracy of the solution is by computing the so-called conserved quantities, which result from general-

11

1.5 1 10 0.5 8 0 −0.5 30

6 4 20

10

2 0

−10 x

−20

−30

0

t

Figure 2: Sinc-collocation approximation of (37)

Figure 3: Error in solving problem (37) using the sinc method ized conservation laws. The first five invariant quantities are given by Z ∞ I1 = udx, −∞ Z ∞ 1 2 I2 = u dx, −∞ 2 Z ∞ 1  I3 = u3 + u2x dx, (39) 2 −∞ Z ∞  I4 = 5u4 + 10uu2x + u2xx dx, and Z−∞  ∞  21u5 + 105u2 u2x + 21uu2x + u2xxx dx . I5 = −∞ 12

Ideally, there would be no variation in these quantities from the initial value. Table (1) provides the values of the maximum normalized absolute variation of the j th observed conserved quantity, IjObs , from the theoretical value, IjT h , defined as I T h − I Obs j j V = max (40) , t IjT h

for the five conserved quantities. Also included in Table (1) are the variations in conserved quantities for the solution to problem (37) using the method of lines (c.f. [23]). The spatial dimension was discretized using 100 grid points, which required 3261 grid points in the time dimension. The CPU time was 967 seconds, the approximate solution deviated from the exact solution by less than 8 × 10−4 at each grid point, and ||e||∞ = 5.18 × 10−2 . Table 1: Max normalized absolute variation of the conserved quantities. Vsinc VM oL I1 3.27 × 10−7 4.13 × 10−6 I2 1.85 × 10−11 1.96 × 10−8 I3 1.89 × 10−5 1.88 × 10−5 I4 6.94 × 10−5 6.93 × 10−5 I5 1.54 × 10−3 1.54 × 10−3 In order to judge how the number of collocation points, N , affected the accuracy and efficiency of the method, the running time and maximum error were recorded while varying the number of collocation points, as seen in figure (4). Note the interesting behavior as N reaches 80; the running time continues to increase exponentially without making the solution any more accurate.

13

3

0

10

10 Running time Maximum error

−1

10 2

10

−2

10 1

−3

10

−4

Error

Time (s)

10

10

0

10

−5

10 −1

10

−6

10 −2

10

−7

1

10

2

10 N

10 3 10

Figure 4: Accuracy and efficiency of the solution to (37) A finite difference algorithm was also attempted, but due to the sensitivity of the stability on parameters, the solution would not converge for (37). In contrast, the stability of the sinc method is independent of all parameters except θ. We conclude that the sinc method is both significantly more efficient and accurate than the method of lines, and much more robust than the finite difference method. Note that because calculating the conserved quantities in (39) requires numerically reconstructing derivatives of the solution, creating the values in table (1) is computationally expensive. Therefore, even though the sinc method is relatively efficient for calculating u(x, t) using a large number of grid points, evaluating the accuracy of the solution using conserved quantities quickly becomes infeasible. As an example, though constructing the solution shown in figure (2) took only about 3.5 seconds, generating the values of Vsinc in table (1) required approximately 10 minutes. With the accuracy of the method established, we now turn to solving the original problems, (4) and (5). Figure (5) shows the computed solution of (4), with parameters N = 100, nt = 1000, L = 1, α = 0.9, β = 10−3 , γ = 10−5 , T = 1, θ = 1, and h(t) = sin(20πt). Of interest here is a numerical demonstration that the results of [29] – namely that the solution of (4) exhibits eventual time-periodic behavior – apply for forcing of arbitrary amplitude, and that the period of the solution is the same as the period of the forcing. By overlaying the solutions in figures (5) in (6), it becomes readily apparent that the solutions eventually settle into periodic behavior, and except for a phase shift, this periodic behavior exactly matches the period of forcing at the boundary.

14

x = L/4 1

0.5

0.5

0

0

u

u

x=0 1

−0.5

−0.5

−1

−1 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

t x = L/2

0.8

1

0.8

1

x=L

1

1

0.5

0.5

0

0

u

u

0.6 t

−0.5

−0.5

−1

−1 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

t

0.6 t

Figure 5: Solution of (4) at several points in the tank

1 0.8 0.6 0.4

u

0.2 0 −0.2 −0.4 Forcing u(x=L/4) u(x=L/2)

−0.6 −0.8 −1 0

0.2

0.4

0.6

0.8

1

t

Figure 6: Overlaid solutions to (4) Figure (7) shows the solution to (5) at the midpoint in the tank overlaid with the external forcing applied to the tank, with parameters N = 100, nt = 1000, L = 10, α = 1, β = 10−3 , γ = 10−5 , T = 1, θ = 1, and h(t) = cos(20πt). If the tank is taken to be semi-infinite, then the solution is identical to that shown in figure (7) for every value of x. If the tank is of 15

finite length, the solution shown is not stable – after some time, reflections from the ends of the tank disturb the periodicity. Notice that the solution shown is phase-shifted −π/2 from the forcing. Thus, each point in the tank behaves as if it was a mass on a damped spring being driven by a sinusoidal force at its natural frequency. 1 Forcing u(x=L/2)

0.8 0.6 0.4

u

0.2 0 −0.2 −0.4 −0.6 −0.8 −1

0

0.2

0.4

0.6

0.8

1

t

Figure 7: Solution to (5) overlaid with external forcing Figure (8) represents a simulations of periods calculated from computed time series data u(0.5, t) with cross correlation with the boundary data with known periods, both seem to be in agreement. Finally figure (9) is 3D plot of solution of (4) with boundary forcing 0.5 sin(10πt).

16

1.2 Period of boundary data Period of numerical data u(0.5,t) obtained from cross correlation

1

Period

0.8

0.6

0.4

0.2

0 1

2

3

4

5

6

7

8

9

n

Figure 8: Periods computed using cross correlation between the boundary data and numerical data u(0.5, t) for problem (4), where n denotes the number of trials.

17

Figure 9: 3D plot of solution to (4).

5

Conclusion

The sinc-collocation method was described in detail and implemented to compute a numerical solution to the nonlinear third order Korteweg-de Vries equation with periodic forcing at the boundary. A detailed stability analysis was provided, which produced a necessary and a sufficient bound on an input parameter, θ, and this bound was confirmed numerically and shown to be very nearly sufficient. The method was first implemented on a test problem with a known analytic solution to evaluate the accuracy of the solution. The test problem also allowed the sinc method to be compared to two other schemes, the method of lines and the finite difference method, and it was found that the sinc method enjoys a great deal of robustness, stability, and computational efficiency that the other methods do not. Finally, the sinc method was used to numerically solve the problem in question, and to show that the solution becomes asymptotically time periodic (with the same period as the boundary forcing), and that this effect could be observed with any small amplitude of boundary forcing.

18

References [1] Kamel Al-Khaled, Numerical study of Fisher’s reaction-diffusion equation by the Sinc collocation method. J. Comput. Appl. Math., 137, (2001), Pages 245-255. [2] Kamel Al-Khaled, Sinc numerical solution for solitons and solitary Waves, J. Comput. Appl. Math. Vol. 130, No. 1-2, pp. 283-292 (2001). [3] N. Bellomo and L. Rodolfi. Solution of nonlinear initial-boundary value problems by sinc collocation-interpolation methods. Computers Math. Applic., 29(4):15–28, 1995. [4] J.L. Bona, G. Pritchard, and L.R. Scott. An evaluation of a model equation for waver waves. Philos. Trans. Roy. Soc. London A, 302:457– 510, 1981. [5] J.L. Bona and R. Smith. The initial value problem for the Korteweg-de Vries equation. Philos. Trans. Roy. Soc. London A, 278:555–601, 1978. [6] J.L. Bona and R. Winther. The Korteweg de-Vries equation in a quarter plane, continuous dependence results. Differential Integral Equations, 2:228–250, 1989. [7] J. L. Bona, S. M. Sun and B.-Y. Zhang, A non-homogeneous boundaryvalue problem for the Korteweg-de Vries equation in a quarter plane, Trans. American Math. Soc. 354 (2002), 427-490. [8] J. L. Bona, S. M. Sun and B.-Y. Zhang, Forced Oscillations of a Damped Korteweg-de Vries Equation in a quarter plane, Comm. in PDEs, 5 (2003), 369-400. [9] J. L. Bona, S. M. Sun and B.-Y. Zhang, A non-homogeneous boundaryvalue problem for the Korteweg-de Vries Equation on a finite domain, Comm. in PDEs, 28 (2003), 1391-1436. [10] J. V. Boussinesq, Th´eorie de l’intumescence liquide appel´ee onde solitaire ou de translation se propageant dans un canal rectangulaire, C. R. Acad. Sci. Paris, 72(1871), pp. 755–759. [11] A. Bueno-Orovio, V. M. Prez-Garca, and F. H. Fenton Spectral Methods for Partial Differential Equations in Irregular Domains: The Spectral Smoothed Boundary Method, SIAM J. Sci. Comput., 28(3), 886?900, 2006. [12] E. Fermi, J. Pasta, and S. Ulam. Studies in nonlinear problems, I. Los Alamos Report LA 1940, 1955.

19

[13] B. Fornberg and G.B. Whitham. A numerical and theoretical study of certain nonlinear wave phenomena. Philos. Trans. Roy. Soc. London A, 289:373–404, 1978. [14] S. Haq, N. Bibi, S.I.A. Tirmizi, and M. Usman. Meshless method of lines for the numerical solution of generalized Kuramoto-Shivashinsky equation. Applied Mathematics and Computation, (217)6: 2404–2413, 2010. [15] Y. Hashizume. Nonlinear pressure waves in a fluid-filled elastic tube. J. Phys. Soc. Jpn., 54:3305–3312, 1985. [16] T. Kawahara. Oscillatory solitary waves in dispersive media. J. Phys. Soc. Jpn., 33:260–264, 1972. [17] D.J. Korteweg and F. de Vries. On the change of form of long waves advancing in a rectangular canal. Philos. Mag., 39:422–433, 1895. [18] James E. Lidsay. Cosmology and the Korteweg-de Vries equation. arXiv:1205.5641v2 [astro-ph.CO], 2012. [19] A.R. Mitchell and S.W. Schoombie. Finite element studies of solitons, pages 465–488. Wiley, 1984. [20] Reza Mokhtari and Maryam Mohammadi. Numerical solution of GRLW equation using sinc-collocation method. Comp. Phys. Sim., 181:1266– 1274, 2010. [21] R. Revelli and L. Ridolfi. Sinc collocation-interpolation method for the simulation of nonlinear waves. Computers Math. Applic., 46:1443–1453, 2003. [22] Mehri Sajjadian. Numerical solutions of Korteweg de Vries and Korteweg de Vries-Burger’s equations using computer programming. arXiv:1209.1782v1 [math.NA], 2012. [23] W.E. Schiesser. The Numerical Method of Lines. Academic Press, 1991. [24] Jie Shen, Jiahong Wu, and Juan-Ming Yuan. Eventual periodicity for the KdV equation on a half-line. Physica D, 227:105–119, 2007. [25] J. O. Skogestad and H. Kalisch A boundary value problem for the KdV equation:comparison of finite difference and Chebyshev methods. Math. Comput. Simulation, 80(1):151–163, 2009. [26] Frank Stenger. Numerical methods based on Whittaker cardinal, or sinc functions. SIAM Review, 23(2):165–224, 1981.

20

[27] Frank Stenger. Numerical Methods Based on Sinc and Analytic Functions. Springer, 1993. [28] Siraj ul Islam, Sirajul Haq, and Arshed Ali. A meshfree method for the numerical solution of the RLW equation. J. Comput. Appl. Math., 223:997–1012, 2009. [29] Muhammad Usman and Bing-Yu Zhang. Forced oscillations of a class of nonlinear dispersive wave equations and their stability. J. Syst. Sci. and Complexity, 20:284–292, 2007. [30] Muhammad Usman and Bing-Yu Zhang. Forced Oscillations of the Korteweg-de Vries Equation on a Bounded Domain and their Stability. Discrete and Continuous Dynamical Systems-Series A., Vol. 26, No.4, p 1509–1523, 2010. [31] Rahmat Ali Khan and Muhammad Usman. Eventual periodicity of forced oscillations of the Korteweg-de Vries type equation. Applied Mathematical Modelling, Volume 36, Issue 2, Pages 736-742, February 2012. [32] N. Yi, Y. Huang and H. Liu. A direct discontinuous Galerkin method for the generalized Korteweg-de Vries equation: energy conservation and boundary effect. J. Comput. Phys., 242:351–366, 2013. [33] N.J. Zubusky and M.D. Kruskal. Interactions of “solitons” in a collisionless plasma and the recurrence of initial states. Phys. Rev. Lett., 15:240–243, 1965.

21