- Email: [email protected]
Sinc numerical solution for solitons and solitary waves
Journal of Computational and Applied Mathematics 130 (2001) 283–292 www.elsevier.nl/locate/cam
Sinc numerical solution for solitons and solitary waves Kamel Al-Khaled Department of Mathematical Sciences, Jordan University of Science and Technology, P.O. Box 3030, Irbid 22110, Jordan Received 17 May 1999; received in revised form 27 September 1999
Abstract A numerical scheme using Sinc–Galerkin method is developed to approximate the solution for the Korteweg– de Vries model equation. Sinc approximation to both derivatives and inde0nite integral reduce the integral equation to an explicit system of algebraic equations, then using various properties of Sinc functions, it is shown that the Sinc solution produce an error of order O(exp(−c=h)) for some positive constants c; h. The method is applied to a few test examples to illustrate c 2001 Elsevier Science B.V. All rights reserved. the accuracy and the implementation of the method. MSC: 35Q53; 35B Keywords: Korteweg– de Vries equation; Pulse solution; Sinc function
1. Introduction The theory of nonlinear dispersive wave motion has recently undergone much study, especially by Whitham [14]. It can be shown [14] that the theory of water waves for the case of shallow water and waves of small amplitude can be approximately described by the Korteweg–de Vries equation ut + (c + u)ux + uxxx = 0;
(x; t) ∈ R × (0; T0 );
(1.1)
where c and are given constants, and u gives the height of a wave above some equilibrium level. Since the amplitude of these waves is assumed to be small, it can serve as a perturbation parameter. These problems have been studied by many authors [2,7–9,11,15]. However, they used a formal perturbation technique. Sometimes called multiscale expansion, or, using evens functions techniques, as in [11]. Also 0nite di@erence method for the Korteweg–de Vries equation have been analyzed in [5]. One aspect that has been investigated is the linearized form of Eq. (1.1): ut + cux + uxxx = 0 c 2001 Elsevier Science B.V. All rights reserved. 0377-0427/01/$ - see front matter PII: S 0 3 7 7 - 0 4 2 7 ( 9 9 ) 0 0 3 7 6 - 3
(1.2)
284
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
which has traveling wave solutions u(x; t) = a cos(kx − !t), where a is constant and ! = !(k) = ck − k 3 . The existence of traveling wave solutions to (1.2) already has been studied in [4,12]. If we drop the third derivative term in (1.1), we have ut + (c + u)ux = 0 which is a quasi-linear 0rst-order wave equation whose wave speed depends on the amplitude and has the implicit solutions u(x; t) = a cos[kx − k(c + u)t]. If c = 0; = 1 in Eq. (1.1) we get another form of Korteweg–de Vries equation ut + uux + uxxx = 0:
(1.3)
This nonlinear equation admits traveling wave solutions of di@erent types. One particular type of traveling wave that arises from the Korteweg–de Vries equation is the soliton, or solitary wave [10, p. 38]. The same equation (1.3) has also come up in the theory of plasma and several other branches of physics. In general, it is diFcult to specify all the terms in a perturbation series for problems involving partial di@erential equations. Thus only the 0rst few terms in the series are determined, which is not enough sometimes. Also more that one type of expansion may be necessary to completely describe the perturbation or asymptotic solution of a given problem. For example, an expansion may break down in some region or may be insuFcient to satisfy the data for the problem. These diFculties signify that the given expansion is not uniformly valid over the entire region of interest. There exist relatively few types of procedures for obtaining approximate solutions of ordinary or partial di@erential equations. While these methods are referred to by a variety of names, such as Rayleigh’s method, Galerkin’s method, or the collocation method, they are all, in e@ect, variants of the same method. One 0rst selects a suitable basis {Sj }nj=1 , and then attempts an approximate solution of the forms n
aj Sj
(1.4)
j=1
in which the coeFcients aj are unknown. Form (1.4) is substituted into the di@erential equation to be solved, and the coeFcients aj are then determined in one of the variety of ways, depending on which of the above procedures is used. We rely heavily on collocation method in the sequel. In this paper the method of solving (1.1) subject to the initial condition u(x; 0) = u0 (x);
x∈R
(1.5)
is based on using the Sinc method, which builds an approximate solution valid on the entire spatial domain and on a small interval in the time domain. The main idea is to replace di@erential and integral equations by their Sinc approximations. The ease of implementation coupled with the exponential convergence rate have demonstrated the viability of the method. One avenue that deserves attention is, approximation by Sinc functions handles singularities in the problem. The paper has been organized into four sections. Section 2 contains, de0nitions and some results of the Sinc functions theory. In Section 3 Sinc solution is developed for the Korteweg–de Vries equation, along with a detailed description of the associated errors. The resulting discrete system for (1.1) is then complied, the matrix structure is examined in detail and appropriate bounds are given also in Section 3. In Section 4 numerical examples are presented which illustrate the exponential convergence of the method.
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
285
2. Sinc function preliminaries The goal of this section is to recall notation and de0nitions of Sinc function, state some known results, and derive useful formulas that are important for this paper. First, we denote the set of all integers, the set of all real numbers, the set of all complex numbers by Z; R and C, respectively. Let f be a function de0ned on R and h ¿ 0 a stepsize. Then the Whittaker cardinal function is de0ned by the series ∞
C(f; h; x) =
f(kh)S(k; h)(x);
k=−∞
whenever this series converges and where S(k; h)(x) =
sin[(x − kh)=h] (x − kh)=h
is known as the kth Sinc function. Also for positive integer N , de0ne CN (f; h; x) =
N
f(kh)S(k; h)(x):
(2.1)
k= −N
Denition 1. Let d ¿ 0, and let Dd denote the region {z = x + iy: |y| ¡ d} in the complex plane C, and the conformal map of a simply connected domain D in the complex domain onto Dd such that (a) = −∞ and (b) = ∞, where a and b are boundary points of D, i.e., a; b ∈ @D. Let denote the inverse map of , and let the arc , with endpoints a and b (a; b ∈ ), given by = (−∞; ∞). For h ¿ 0, let the points xk on be given by xk = (kh); k ∈ Z; (z) = exp((z)). For 16p ¡ ∞, let H (D) denote the family of all functions f that are analytic in D, such that |f(z)| |d z| ¡ ∞. @D Corresponding to the number !, let L! (D) denote the family of all functions f analytic for which there exists a constant C1 such that |f(z)|6C1
| (z)|! ; [1 + | (z)|]2!
z ∈ D:
Now if x is on the arc , then by introducing the conformal map , and a “nulli0er” function g the following theorem gives a formula for approximating f(m) on . Let g be an analytic function on D, and for k ∈ Z set
(z) − jh Sj (z) = g(z)Sinc = g(z)S(j; h) ◦ (z); h
z ∈ D:
Theorem 2.1 (Stenger [13, p. 208]). Let f=g ∈ H (D); and let sup
−=h6t6=h
n d −n d x g(x) exp(it(x)) 6C2 h ;
x∈
286
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
for n = 0; 1; 2; : : : ; m; with C2 a constant depending only on m; ; and g. If f=g ∈ L! (D); ! a positive constant; then taking h = d=!N it follows that N √ (n) f(x ) j Sj(n) (x) 6C3 N (n+1)=2 exp(− d!N ) sup f (x) − g(xj ) x∈ j=−N
for n = 0; 1; 2; : : : ; m; with C3 a constant depending only on m; ; g; d; !; and f. The approximation of the mth derivative of f in Theorem 2.1 is simply an mth derivative of CN (f=g; h; x)=g in (2.1). The weight function g is chosen relative to the order of the derivative that is to be approximated. For instance, to approximate the mth derivative the choice g(x) = 1=( (x))m often suFce, [13]. So the approximation of f by Sinc expansion is given by N f(xj )
f (x) ≈
j=−N
g(xj )
Sj (x):
(2.2)
The derivatives of Sinc functions evaluated at the nodes will also be needed and these quantities are delineated by (q) ≡ hq %jk
dq [Sj ◦ (x)]|x=xk : dq
In particular, the following convenient notation will be useful in formulating the discrete system:
%(0) jk
= [S(j; h) ◦ (x)]|x=xk =
%(1) jk = h
1; j = k; 0; j = k; 0;
j = k; d [S(j; h) ◦ (x)]|x=xk = (−1) j−k d ; j = k j−k
and %(3) jk
0;
j = k; d j−k =h [S(j; h) ◦ (x)]|x=xk = (−1) d3 [6 − 2 (j − k)2 ]; j = k: (j − k)3 3
3
So that the approximation in (2.2) at the Sinc nodes xk takes the form
f (xk ) ≈
N j=−N
%(1) f(xj ) jk : + %(0) jk g (xj ) h g(xj )
(2.3)
System (2.3) is more conveniently recorded by de0ning the vector f = (f−N ; : : : ; f0 ; : : : ; fN )T . (q) Then de0ne the m × m; (m = 2N + 1) Toeplitz matrices Im(q) = [%jk ]; q = 0; 1; 2; whose jkth entry is (q) ; q = 0; 1; 2. Also de0ne the diagonal matrix D(g) = diag[g(x−N ); : : : ; g(x0 ); : : : ; g(xN )]T . given by %jk
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
287
The matrix Im(0) is the identity matrix. The matrix Im(1) takes the form (−1)m−1 0 −1 : : : m−1 .. 1 : : : . (1) (2.4) Im = . . . .. .. .. −1 (−1)m ::: ::: 0 m−1 and (−1)m−1 (6 − (m − 1)2 2 ) 2 0 −(6 − ) ::: (m − 1)3 .. 2 6− . (3) .. .. (2.5) Im = : . . −(6 − 2 ) m 2 2 (−1) (6 − (m − 1) ) 2 (6 − ) 0 3 (m − 1) System (2.3) takes the form −1 (1) (0) I D(1=g) + Im D(g =g) f ≡ A1 f : (2.6) f ≈ h m For the present paper the interval in Theorem 2.1 is (−∞; ∞). Therefore, to approximate the 0rst derivative we take (x)=x, and g(x)=1= (x), the square matrix (2.6) becomes A1 =(−1=h)Im(1) . Then the approximation of the 0rst derivative evaluated at the vector nodes xi can be written as (2.7) f (xi ) ∼ = A1 f (xi ): The same way we can approximate the third derivative by f (xi ) ∼ = A3 f (xi ); 3
(2.8)
)Im(3) .
where the matrix A3 is de0ned by A3 = (−1=h Let us describe Sinc de0nite integration over an interval. At the outset, we de0ne the numbers k (−1) (−1) %k = 12 + 0 sin(t)=t dt; k ∈ Z, and we de0ne Toeplitz matrix %k−j of order m by (−1) ] I (−1) = [%k−j
(2.9)
(−1) denoting the kjth element of I (−1) . For h= d=(!N ), de0ne the matrix B=hI (−1) D(1=, ). with %k−j x Let F ∈ H (D); set F˜ = (F(z−N ); : : : ; F(zN ))T ; de0ne (FF)(x) = a F(t) dt, set G˜ = (G−N ; : : : ; G0 ; : : : ; ˜ then de0ne FN F by GN )T = BF,
N (x)GN e kh GN + (FN F)(x) = Gk − 1 + (x) k=−N 1 + e kh
S(k; h) ◦ ,(x)
then for F=, ∈ L! (D) we have [13, p. 219]
√ sup |(FF)(x) − (FN F)(x)|6C4 exp(− d!N ); x∈
where C4 is a constant that is independent of N .
(2.10)
288
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
3. Implementation of the method We determine the Sinc approximation for Korteweg–de Vries equation under the assumption that the initial condition in (1.5) belongs to the class of functions L! (D). Integrating Eq. (1.1) with respect to t, we get u(x; t) = −
t
0
[uxxx (x; /) + (c + u(x; /))ux (x; /)] d/ + u0 (x):
(3.1)
To obtain a direct discretization of Eq. (3.1), and since the domain is R×(0; T0 ), the relevant maps are de0ned as follows: In the space direction, choose the map (x) = x which maps the in0nite strip Dd = {0 = 1 + i2: |2| ¡ d} onto Dd . In the time direction, choose the map ,(t) = log(t=T0 − t) which carries the eye-shaped region D3 = {t = x + iy: |arg(t=T0 − t)| ¡ d6=2} onto the in0nite strip Dd . The compositions S(m; hx ) ◦ (x); m = −Nx ; : : : ; Nx and S(k; ht ) ◦ ,(t); k = −Nt ; : : : ; Nt de0ne the basis elements for (−∞; ∞) and (0; T0 ), respectively, the mesh sizes hx and ht represent the mesh sizes in the in0nite strip Dd for the uniform grid {ihx }; −∞ ¡ i ¡ ∞, and {jht }; −∞ ¡ j ¡ ∞. The Sinc grid points xi ∈ (−∞; ∞) in Dd and tj ∈ (0; T0 ) in D3 are the inverse images of the equispaced grid points; that is, xi = −1 (ihx ) = ihx
and
tj = , −1 (jht ) = (T0 exp(jht ))=(1 + exp(jht )):
In Eq. (3.1) let us carry out the Sinc approximation of ux and uxxx . To proceed use Eqs. (2.7), (2.8) and replace ux by −Im(1)x u(xi ; t)=hx and uxxx by −Im(3)x u(xi ; t)=h3x where mx = 2Nx + 1 and the skew-symmetric matrices Im(1)x and Im(3)x as de0ned in (2.4), (2.5), respectively. Next in Eq. (3.1) evaluate u(x; t) at the x-nodes, and replace ux (x; t); uxxx (x; t) by their approximations, we get the Volterra integral equation u(t) = −
0
t
[A3 u(/) + (c + u(/))A1 u(/)] d/ + u0 ;
(3.2)
where the square matrices A1 ; A3 are given by A1 ∼ = −Im(1)x =hx and A3 ∼ = −Im(3)x =h3x , with u(t) = T T [u−Nx (t); : : : ; uNx (t)] , where ui (t) = u(xi ; t) and u0 = [u0 (z−Nx ); : : : ; u0 (zNx )] . We next collocate with respect to the t-variable via the use of the de0nite integration formula, (see Eq. (2.10)) with the conformal map ,(t) = log(t=T0 − t). Thus, de0ne the matrix B by B = −1 ht Im(−1) D(1=, ) with the nodes t = , (jh ) for j = −N ; : : : ; N , where h = d=(!Nt ), and Im(−1) j t t t t t t 0 0 0 as de0ned in (2.9), with mt = 2Nt + 1. De0ne the matrix U by U = [u (xi ; 0)]. Then the solution of Eq. (3.2) in matrix form is given by the rectangular mx × mt matrix U = [ui; j ]: U = −[A3 U + (C + U ) ◦ A1 U ]BT + U 0 ;
(3.3)
where the notation “◦” denotes the Hadamard matrix multiplication. Note that in our discretization we are taking the time nodes as rows and the space nodes as columns, so the matrix (A3 U + (C + U ) ◦ A1 U ) forms the vector nodes for the integral in (3.2). In (3.3) the vector U 0 has the same dimensions as the vector U and every column of U 0 consists of the same vector u0 . Also the mx × mt matrix C has each entry as the number c. We use the notation U = [u(xi ; tj ] to denote the mx × mt matrix of node values of the function u where mx = 2Nx + 1; mt = 2Nt + 1. If u(x; t) is the exact solution of the integral equation (3.1). Then the approximation of ux (x; t) in matrix form has an exponential error (see Theorem 2.1)
[ux (x; t)] − A1 U 6C3 Nx exp(− d!Nx ):
(3.4)
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
289
Here we assumed that C3 is bounded with respect to t, and the matrix A1 as de0ned above. Also we require that the error in approximating the third derivative uxxx (x; t) is exponentially small, i.e., with A3 as de0ned above we have
[uxxx (x; t)] − A3 U 6C3 Nx2 exp(− d!Nx ):
t
(3.5)
We also require that the error in approximating the integral K(x; t) = − 0 G(/) d/, where G(t)= [uxxx (x; t)+(c +u(x; t))ux (x; t)] is exponentially small. This means, the approximation of the integral in matrix form takes the form
[K(xi ; tj )] − BU 6C4 exp(− d!Nt );
(3.6)
where again C4 is bounded with respect to x, and the matrix B as de0ned above. We shall obtain a small error in our approximation to Eq. (1.1), as noted in the following two theorems. The proof resembles the proof of Theorems 4:1 and 4:2 in [1]. Theorem 3.1. Let u0 ∈ L! (D); let the function u(x; t) be as in Eq. (3:1); and let the matrix U be de8ned as in (3:3). Then for Nx ; Nt ¿ 16=d! there exists a constant C5 independent of Nx ; Nt such that √ sup [u(xi ; tj )] − U 6C5 N 2 exp(− d!N ); (xi ; tj )
where N = min{Nx ; Nt }. Theorem 3.2. Given a constant R ¿ 0; there is a constant T0 ¿ 0 such that if U 1 − U 0 ¡ R=2; then the solution U = G(U ) + U 0 has a unique solution. Moreover; the iteration scheme U n+1 = G(U n ) + U 0 converges to this unique solution. 4. Numerical examples There are many numerical simulation results which seem to support the stability of pulses (for detailed simulations see [8,9]). However, one should notice that all these numerical simulations are performed in an interval of 0nite length. One might expect that the stability=instability problem of pulses can be numerically realized by taking suFcient calculation length. But still another diFculty remains. If one wants to approximate a pulse solution numerically in a wide interval, then the instability will grow [12, p. 490]. We therefore consider it is diFcult to realize the stability=instability problem for pulses numerically. However, in this section, at least, suggest what happens in the idealized in0nite interval situation. The model problem (1.1) was tested on four examples. For the four examples sequences of runs with Nx = Nt = 8; 16; 24; 32 are reported. The problem included here illustrate various features of the method, demonstrating the ease of implementation and assembly of the discrete system, the choice of parameters and the exponential convergence rate. Example 4.1. We solve the model equation ut + 6uux + uxxx = 0;
x∈R
(4.1)
290
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292 Table 1 The error in the approximate solution for Example 4.1 N x = Nt
T0 = 3 supuij − u(xi ; tj )
8 16 24 32
1:2698 × 10−3 5:8980 × 10−4 5:2229 × 10−4 4:0157 × 10−4
Table 2 Results for Example 4.2 using Sinc approach N x = Nt
hx = ht
supuij − u(xi ; tj )
8 16 24 32
0.785 0.555 0.453 0.392
1:5040 × 10−2 2:2252 × 10−3 6:8275 × 10−6 1:0636 × 10−6
with the initial data of a single soliton, say u(x; 0) = 2 sech2 x: Here the exact solution is given by u(x; t) = 2 sech2 (x − 4t). We solve this model equation using our approach when d = =2; ! = 1. Table 1 shows that the method converges for T0 = 3. The second column reports the supremum norm of the error between the exact solution u(xi ; tj ) and the approximate solution uij . Again the exponential convergence rate decreases the error that are initially present for small Nx ; Nt : Example 4.2. Let us consider the same problem in the last example, and solve the Eq. (4.1) with the initial condition u(x; 0) = 6 sech2 x. In this case the exact solution of Eq. (4.1) is given by 12[3 + 4 cosh(2x − 8t) + cosh(4x − 64t)] u(x; t) = : [3 cosh(x − 28t) + cosh(3x − 36t)]2 We compared the performance of the method presented here by the 0nite element scheme. Finite element methods for the Korteweg–de Vries equation have been analyzed by Winther [15], where he proved convergence for a class of these equations. Tables 2 and 3 refer to a single-soliton solution and correspond, respectively, to the Sinc method and to the 0nite element method. The tables display information at T0 = 3. We recall that, even though the tables correspond to T0 = 3, the integration was followed up to T0 = 20 and there was no problems of divergence of our scheme. Comparing the errors for the two schemes, we see that for the 0nite element method we do have an error of order O(h2 ), while for the Sinc method it is O(exp(−c=h)) where h is the step size and c ¿ 0. And as predicted, the method using Sinc basis functions is more accurate for soliton propagation. Example 4.3. The known exact one-soliton solutions of Eq. (4.1) exhibit exponential decay for large x. The inverse scattering transform [3] tells us that the solution of (4.1) with generic initial
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
291
Table 3 Results for Example 4.2 using 0nite elements method hx = ht
supuij − u(xi ; tj )
0.785 0.555 0.453 0.392
2:1102 × 10−2 6:4541 × 10−3 8:5009 × 10−4 5:2001 × 10−5
Table 4 Results for Example 4.3 using Sinc approach N x = Nt
hx = ht
supuij − u(xi ; tj )
8 16 24 32
0.785 0.555 0.453 0.392
8:5423 × 10−3 5:7101 × 10−3 1:3652 × 10−3 8:8602 × 10−5
condition u(x; 0) = u0 (x) (with u0 (x) → 0 suFciently rapidly as |x| → ∞) consists of a train of solitons moving to the right, along with a dispersive wave travelling to the left. As in the 0 above two examples, we can argue that when √ u (x) has compact support, the solution emerge from the rightmost soliton decaying as exp(− cx), where c is the relevant speed, but for suFciently large x the presence of the solitons will not yet be felt, and the behavior of the solution will be determined by the √ Green’s function of the linearized equation ut + uxxx = 0, i.e., u will decay roughly as exp(−2x3=2 =3 3t). Thus there is a transition in the nature of the decay. If there is more than one soliton in the soliton train, say two, with speeds c1 ; c2 (c2 ¿ c1 ), the a problem arises. The √ solution emerges from the faster-moving soliton decaying as exp(− c2 x), but because the tail of the slower-moving soliton falls slower, it is possible that it will return to dominate, i.e., the decay √ will slow to exp(− c1 x). We note that the exact two-soliton solutions (c2 − c1 )(c2 cosh2 !1 + c1 sinh2 !2 ) ; √ 2( c2 cos h!1 cos g!2 − c1 sin h!1 sin h!2 )2 √ √ where !1 = c1 (x − c1 t)=2 and !2 = c2 (x − c2 t)=2 exhibit exactly this phenomenon: when x is large √ u ∼ exp(− c1 x), i.e., the tail of the slow soliton dominates the decay. Table 4 refers to two-soliton solution with c1 = 6; c2 = 10. The table displays information at T0 = 3, with ! = 1; d = =2. u(x; t) =
√
Example 4.4. Vanden-Broeck [16] considered the time-independent surface water of an incompressible and inviscid Puid Pow over a bump in a two-dimensional channel. The model equation is the forced Korteweg–de Vries equation (fkdv): ut + 7ux + 2!uux + uxxx = fx (x);
x ∈ R; t ¿ 0:
The forcing represented by the function f(x) in the fkdv equation is due to the bump on the bottom of the channel. Solitary wave means that the free-surface elevation u(x) has the property u(∓∞) = ux (∓∞) = uxx (∓∞) = 0:
292
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292 Table 5 Results for Example 4.4 using Sinc approach N x = Nt
hx = ht
supuij − u(xi ; tj )
8 16 24 32
0.785 0.555 0.453 0.392
4:1314 × 10−2 1:6575 × 10−2 3:5789 × 10−3 2:4625 × 10−4
To illustrate our numerical results, let us take ! = − 34 ; = − 16 ; 7 = 3, and let the forcing function f be de0ned by the equation
f(x) =
sin(x); |x| ¡ 1; 0; |x| ¿ 1:
We use our approach and compare our result with the solution in [6]. Table 5 displays some results, where the solution u(xi ; tj ) is assumed to be the solution in [6], and uij is our Sinc solution. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
K. Al-Khaled, A Sinc–Galerkin approach on the p-system, Saitama Math. J. 16 (1998) 1–13. R. Courant, D. Hilbert, Methods of Mathematical Physics, Wiley-Interscience, New York, 1962. P.G. Drazin, R.S Johnson, Solitons: An Introduction, Cambridge University Press, Cambridge, 1989. N.M. Ercolani, D.W. Mclaughlin, H. Roitner, Attractors and transients for a perturbed Kdv equation, A nonlinear spectral analysis, J. Nonlinear Sci. 3 (1993) 477–539. Feng Bao-Feng, M. Taketomo, A 0nite di@erence method for the Korteweg– de Vries and the Kadomtsev–Petviashvili equations, J. Comput. Appl. Math. 1 (90) (1998) 97–118. L. Gong, S.S. Shen, Multiple supercritical solitary wave solution of the stationary forced Korteweg– de Vries equation and their stability, SIAM J. Appl. Math. 54 (5) (1994) 1266–1290. R. Grimshaw, H. Mitsudera, Slowly varying solitary wave solutions of the perturbed Korteweg– de Vries equation revisited, Stud. Appl. Math. 90 (1993) 75–86. T. Kawahara, S. Toh, Nonlinear dispersive periodic waves in the presence of instability and damping, Phys. Fluids 28 (1985) 1636–1638. Y. Kodama, M. Ablowitz, Perturbations of solitons and solitary waves, Stud. Appl. Math. 64 (1994) 225–245. D.J. Logan, An Introduction to Nonlinear Partial Di@erential Equations, Wiley, New York, 1994. T. Ogawa, H. Suzuki, On the spectra of pulses in a nearly integrable system, SIAM J. Appl. Math. 57 (2) (1997) 485–500. T. Ogawa, Traveling wave solutions to a perturbed Korteweg– de Vries equation, Hiroshima Math. J. 24 (1994) 401–422. F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer, New York, 1993. G.B. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, New York, 1974. R. Winther, A conservative 0nite element method for the Korteweg– de Vries equation, Math. Comput. 34 (1980) 23–43. J.M. Vanden-Broeck, Free-surface Pow over a semi-circular obstruction in a channel, Phys. Fluids 30 (1987) 2315– 2317.
Sinc numerical solution for solitons and solitary waves Kamel Al-Khaled Department of Mathematical Sciences, Jordan University of Science and Technology, P.O. Box 3030, Irbid 22110, Jordan Received 17 May 1999; received in revised form 27 September 1999
Abstract A numerical scheme using Sinc–Galerkin method is developed to approximate the solution for the Korteweg– de Vries model equation. Sinc approximation to both derivatives and inde0nite integral reduce the integral equation to an explicit system of algebraic equations, then using various properties of Sinc functions, it is shown that the Sinc solution produce an error of order O(exp(−c=h)) for some positive constants c; h. The method is applied to a few test examples to illustrate c 2001 Elsevier Science B.V. All rights reserved. the accuracy and the implementation of the method. MSC: 35Q53; 35B Keywords: Korteweg– de Vries equation; Pulse solution; Sinc function
1. Introduction The theory of nonlinear dispersive wave motion has recently undergone much study, especially by Whitham [14]. It can be shown [14] that the theory of water waves for the case of shallow water and waves of small amplitude can be approximately described by the Korteweg–de Vries equation ut + (c + u)ux + uxxx = 0;
(x; t) ∈ R × (0; T0 );
(1.1)
where c and are given constants, and u gives the height of a wave above some equilibrium level. Since the amplitude of these waves is assumed to be small, it can serve as a perturbation parameter. These problems have been studied by many authors [2,7–9,11,15]. However, they used a formal perturbation technique. Sometimes called multiscale expansion, or, using evens functions techniques, as in [11]. Also 0nite di@erence method for the Korteweg–de Vries equation have been analyzed in [5]. One aspect that has been investigated is the linearized form of Eq. (1.1): ut + cux + uxxx = 0 c 2001 Elsevier Science B.V. All rights reserved. 0377-0427/01/$ - see front matter PII: S 0 3 7 7 - 0 4 2 7 ( 9 9 ) 0 0 3 7 6 - 3
(1.2)
284
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
which has traveling wave solutions u(x; t) = a cos(kx − !t), where a is constant and ! = !(k) = ck − k 3 . The existence of traveling wave solutions to (1.2) already has been studied in [4,12]. If we drop the third derivative term in (1.1), we have ut + (c + u)ux = 0 which is a quasi-linear 0rst-order wave equation whose wave speed depends on the amplitude and has the implicit solutions u(x; t) = a cos[kx − k(c + u)t]. If c = 0; = 1 in Eq. (1.1) we get another form of Korteweg–de Vries equation ut + uux + uxxx = 0:
(1.3)
This nonlinear equation admits traveling wave solutions of di@erent types. One particular type of traveling wave that arises from the Korteweg–de Vries equation is the soliton, or solitary wave [10, p. 38]. The same equation (1.3) has also come up in the theory of plasma and several other branches of physics. In general, it is diFcult to specify all the terms in a perturbation series for problems involving partial di@erential equations. Thus only the 0rst few terms in the series are determined, which is not enough sometimes. Also more that one type of expansion may be necessary to completely describe the perturbation or asymptotic solution of a given problem. For example, an expansion may break down in some region or may be insuFcient to satisfy the data for the problem. These diFculties signify that the given expansion is not uniformly valid over the entire region of interest. There exist relatively few types of procedures for obtaining approximate solutions of ordinary or partial di@erential equations. While these methods are referred to by a variety of names, such as Rayleigh’s method, Galerkin’s method, or the collocation method, they are all, in e@ect, variants of the same method. One 0rst selects a suitable basis {Sj }nj=1 , and then attempts an approximate solution of the forms n
aj Sj
(1.4)
j=1
in which the coeFcients aj are unknown. Form (1.4) is substituted into the di@erential equation to be solved, and the coeFcients aj are then determined in one of the variety of ways, depending on which of the above procedures is used. We rely heavily on collocation method in the sequel. In this paper the method of solving (1.1) subject to the initial condition u(x; 0) = u0 (x);
x∈R
(1.5)
is based on using the Sinc method, which builds an approximate solution valid on the entire spatial domain and on a small interval in the time domain. The main idea is to replace di@erential and integral equations by their Sinc approximations. The ease of implementation coupled with the exponential convergence rate have demonstrated the viability of the method. One avenue that deserves attention is, approximation by Sinc functions handles singularities in the problem. The paper has been organized into four sections. Section 2 contains, de0nitions and some results of the Sinc functions theory. In Section 3 Sinc solution is developed for the Korteweg–de Vries equation, along with a detailed description of the associated errors. The resulting discrete system for (1.1) is then complied, the matrix structure is examined in detail and appropriate bounds are given also in Section 3. In Section 4 numerical examples are presented which illustrate the exponential convergence of the method.
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
285
2. Sinc function preliminaries The goal of this section is to recall notation and de0nitions of Sinc function, state some known results, and derive useful formulas that are important for this paper. First, we denote the set of all integers, the set of all real numbers, the set of all complex numbers by Z; R and C, respectively. Let f be a function de0ned on R and h ¿ 0 a stepsize. Then the Whittaker cardinal function is de0ned by the series ∞
C(f; h; x) =
f(kh)S(k; h)(x);
k=−∞
whenever this series converges and where S(k; h)(x) =
sin[(x − kh)=h] (x − kh)=h
is known as the kth Sinc function. Also for positive integer N , de0ne CN (f; h; x) =
N
f(kh)S(k; h)(x):
(2.1)
k= −N
Denition 1. Let d ¿ 0, and let Dd denote the region {z = x + iy: |y| ¡ d} in the complex plane C, and the conformal map of a simply connected domain D in the complex domain onto Dd such that (a) = −∞ and (b) = ∞, where a and b are boundary points of D, i.e., a; b ∈ @D. Let denote the inverse map of , and let the arc , with endpoints a and b (a; b ∈ ), given by = (−∞; ∞). For h ¿ 0, let the points xk on be given by xk = (kh); k ∈ Z; (z) = exp((z)). For 16p ¡ ∞, let H (D) denote the family of all functions f that are analytic in D, such that |f(z)| |d z| ¡ ∞. @D Corresponding to the number !, let L! (D) denote the family of all functions f analytic for which there exists a constant C1 such that |f(z)|6C1
| (z)|! ; [1 + | (z)|]2!
z ∈ D:
Now if x is on the arc , then by introducing the conformal map , and a “nulli0er” function g the following theorem gives a formula for approximating f(m) on . Let g be an analytic function on D, and for k ∈ Z set
(z) − jh Sj (z) = g(z)Sinc = g(z)S(j; h) ◦ (z); h
z ∈ D:
Theorem 2.1 (Stenger [13, p. 208]). Let f=g ∈ H (D); and let sup
−=h6t6=h
n d −n d x g(x) exp(it(x)) 6C2 h ;
x∈
286
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
for n = 0; 1; 2; : : : ; m; with C2 a constant depending only on m; ; and g. If f=g ∈ L! (D); ! a positive constant; then taking h = d=!N it follows that N √ (n) f(x ) j Sj(n) (x) 6C3 N (n+1)=2 exp(− d!N ) sup f (x) − g(xj ) x∈ j=−N
for n = 0; 1; 2; : : : ; m; with C3 a constant depending only on m; ; g; d; !; and f. The approximation of the mth derivative of f in Theorem 2.1 is simply an mth derivative of CN (f=g; h; x)=g in (2.1). The weight function g is chosen relative to the order of the derivative that is to be approximated. For instance, to approximate the mth derivative the choice g(x) = 1=( (x))m often suFce, [13]. So the approximation of f by Sinc expansion is given by N f(xj )
f (x) ≈
j=−N
g(xj )
Sj (x):
(2.2)
The derivatives of Sinc functions evaluated at the nodes will also be needed and these quantities are delineated by (q) ≡ hq %jk
dq [Sj ◦ (x)]|x=xk : dq
In particular, the following convenient notation will be useful in formulating the discrete system:
%(0) jk
= [S(j; h) ◦ (x)]|x=xk =
%(1) jk = h
1; j = k; 0; j = k; 0;
j = k; d [S(j; h) ◦ (x)]|x=xk = (−1) j−k d ; j = k j−k
and %(3) jk
0;
j = k; d j−k =h [S(j; h) ◦ (x)]|x=xk = (−1) d3 [6 − 2 (j − k)2 ]; j = k: (j − k)3 3
3
So that the approximation in (2.2) at the Sinc nodes xk takes the form
f (xk ) ≈
N j=−N
%(1) f(xj ) jk : + %(0) jk g (xj ) h g(xj )
(2.3)
System (2.3) is more conveniently recorded by de0ning the vector f = (f−N ; : : : ; f0 ; : : : ; fN )T . (q) Then de0ne the m × m; (m = 2N + 1) Toeplitz matrices Im(q) = [%jk ]; q = 0; 1; 2; whose jkth entry is (q) ; q = 0; 1; 2. Also de0ne the diagonal matrix D(g) = diag[g(x−N ); : : : ; g(x0 ); : : : ; g(xN )]T . given by %jk
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
287
The matrix Im(0) is the identity matrix. The matrix Im(1) takes the form (−1)m−1 0 −1 : : : m−1 .. 1 : : : . (1) (2.4) Im = . . . .. .. .. −1 (−1)m ::: ::: 0 m−1 and (−1)m−1 (6 − (m − 1)2 2 ) 2 0 −(6 − ) ::: (m − 1)3 .. 2 6− . (3) .. .. (2.5) Im = : . . −(6 − 2 ) m 2 2 (−1) (6 − (m − 1) ) 2 (6 − ) 0 3 (m − 1) System (2.3) takes the form −1 (1) (0) I D(1=g) + Im D(g =g) f ≡ A1 f : (2.6) f ≈ h m For the present paper the interval in Theorem 2.1 is (−∞; ∞). Therefore, to approximate the 0rst derivative we take (x)=x, and g(x)=1= (x), the square matrix (2.6) becomes A1 =(−1=h)Im(1) . Then the approximation of the 0rst derivative evaluated at the vector nodes xi can be written as (2.7) f (xi ) ∼ = A1 f (xi ): The same way we can approximate the third derivative by f (xi ) ∼ = A3 f (xi ); 3
(2.8)
)Im(3) .
where the matrix A3 is de0ned by A3 = (−1=h Let us describe Sinc de0nite integration over an interval. At the outset, we de0ne the numbers k (−1) (−1) %k = 12 + 0 sin(t)=t dt; k ∈ Z, and we de0ne Toeplitz matrix %k−j of order m by (−1) ] I (−1) = [%k−j
(2.9)
(−1) denoting the kjth element of I (−1) . For h= d=(!N ), de0ne the matrix B=hI (−1) D(1=, ). with %k−j x Let F ∈ H (D); set F˜ = (F(z−N ); : : : ; F(zN ))T ; de0ne (FF)(x) = a F(t) dt, set G˜ = (G−N ; : : : ; G0 ; : : : ; ˜ then de0ne FN F by GN )T = BF,
N (x)GN e kh GN + (FN F)(x) = Gk − 1 + (x) k=−N 1 + e kh
S(k; h) ◦ ,(x)
then for F=, ∈ L! (D) we have [13, p. 219]
√ sup |(FF)(x) − (FN F)(x)|6C4 exp(− d!N ); x∈
where C4 is a constant that is independent of N .
(2.10)
288
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
3. Implementation of the method We determine the Sinc approximation for Korteweg–de Vries equation under the assumption that the initial condition in (1.5) belongs to the class of functions L! (D). Integrating Eq. (1.1) with respect to t, we get u(x; t) = −
t
0
[uxxx (x; /) + (c + u(x; /))ux (x; /)] d/ + u0 (x):
(3.1)
To obtain a direct discretization of Eq. (3.1), and since the domain is R×(0; T0 ), the relevant maps are de0ned as follows: In the space direction, choose the map (x) = x which maps the in0nite strip Dd = {0 = 1 + i2: |2| ¡ d} onto Dd . In the time direction, choose the map ,(t) = log(t=T0 − t) which carries the eye-shaped region D3 = {t = x + iy: |arg(t=T0 − t)| ¡ d6=2} onto the in0nite strip Dd . The compositions S(m; hx ) ◦ (x); m = −Nx ; : : : ; Nx and S(k; ht ) ◦ ,(t); k = −Nt ; : : : ; Nt de0ne the basis elements for (−∞; ∞) and (0; T0 ), respectively, the mesh sizes hx and ht represent the mesh sizes in the in0nite strip Dd for the uniform grid {ihx }; −∞ ¡ i ¡ ∞, and {jht }; −∞ ¡ j ¡ ∞. The Sinc grid points xi ∈ (−∞; ∞) in Dd and tj ∈ (0; T0 ) in D3 are the inverse images of the equispaced grid points; that is, xi = −1 (ihx ) = ihx
and
tj = , −1 (jht ) = (T0 exp(jht ))=(1 + exp(jht )):
In Eq. (3.1) let us carry out the Sinc approximation of ux and uxxx . To proceed use Eqs. (2.7), (2.8) and replace ux by −Im(1)x u(xi ; t)=hx and uxxx by −Im(3)x u(xi ; t)=h3x where mx = 2Nx + 1 and the skew-symmetric matrices Im(1)x and Im(3)x as de0ned in (2.4), (2.5), respectively. Next in Eq. (3.1) evaluate u(x; t) at the x-nodes, and replace ux (x; t); uxxx (x; t) by their approximations, we get the Volterra integral equation u(t) = −
0
t
[A3 u(/) + (c + u(/))A1 u(/)] d/ + u0 ;
(3.2)
where the square matrices A1 ; A3 are given by A1 ∼ = −Im(1)x =hx and A3 ∼ = −Im(3)x =h3x , with u(t) = T T [u−Nx (t); : : : ; uNx (t)] , where ui (t) = u(xi ; t) and u0 = [u0 (z−Nx ); : : : ; u0 (zNx )] . We next collocate with respect to the t-variable via the use of the de0nite integration formula, (see Eq. (2.10)) with the conformal map ,(t) = log(t=T0 − t). Thus, de0ne the matrix B by B = −1 ht Im(−1) D(1=, ) with the nodes t = , (jh ) for j = −N ; : : : ; N , where h = d=(!Nt ), and Im(−1) j t t t t t t 0 0 0 as de0ned in (2.9), with mt = 2Nt + 1. De0ne the matrix U by U = [u (xi ; 0)]. Then the solution of Eq. (3.2) in matrix form is given by the rectangular mx × mt matrix U = [ui; j ]: U = −[A3 U + (C + U ) ◦ A1 U ]BT + U 0 ;
(3.3)
where the notation “◦” denotes the Hadamard matrix multiplication. Note that in our discretization we are taking the time nodes as rows and the space nodes as columns, so the matrix (A3 U + (C + U ) ◦ A1 U ) forms the vector nodes for the integral in (3.2). In (3.3) the vector U 0 has the same dimensions as the vector U and every column of U 0 consists of the same vector u0 . Also the mx × mt matrix C has each entry as the number c. We use the notation U = [u(xi ; tj ] to denote the mx × mt matrix of node values of the function u where mx = 2Nx + 1; mt = 2Nt + 1. If u(x; t) is the exact solution of the integral equation (3.1). Then the approximation of ux (x; t) in matrix form has an exponential error (see Theorem 2.1)
[ux (x; t)] − A1 U 6C3 Nx exp(− d!Nx ):
(3.4)
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
289
Here we assumed that C3 is bounded with respect to t, and the matrix A1 as de0ned above. Also we require that the error in approximating the third derivative uxxx (x; t) is exponentially small, i.e., with A3 as de0ned above we have
[uxxx (x; t)] − A3 U 6C3 Nx2 exp(− d!Nx ):
t
(3.5)
We also require that the error in approximating the integral K(x; t) = − 0 G(/) d/, where G(t)= [uxxx (x; t)+(c +u(x; t))ux (x; t)] is exponentially small. This means, the approximation of the integral in matrix form takes the form
[K(xi ; tj )] − BU 6C4 exp(− d!Nt );
(3.6)
where again C4 is bounded with respect to x, and the matrix B as de0ned above. We shall obtain a small error in our approximation to Eq. (1.1), as noted in the following two theorems. The proof resembles the proof of Theorems 4:1 and 4:2 in [1]. Theorem 3.1. Let u0 ∈ L! (D); let the function u(x; t) be as in Eq. (3:1); and let the matrix U be de8ned as in (3:3). Then for Nx ; Nt ¿ 16=d! there exists a constant C5 independent of Nx ; Nt such that √ sup [u(xi ; tj )] − U 6C5 N 2 exp(− d!N ); (xi ; tj )
where N = min{Nx ; Nt }. Theorem 3.2. Given a constant R ¿ 0; there is a constant T0 ¿ 0 such that if U 1 − U 0 ¡ R=2; then the solution U = G(U ) + U 0 has a unique solution. Moreover; the iteration scheme U n+1 = G(U n ) + U 0 converges to this unique solution. 4. Numerical examples There are many numerical simulation results which seem to support the stability of pulses (for detailed simulations see [8,9]). However, one should notice that all these numerical simulations are performed in an interval of 0nite length. One might expect that the stability=instability problem of pulses can be numerically realized by taking suFcient calculation length. But still another diFculty remains. If one wants to approximate a pulse solution numerically in a wide interval, then the instability will grow [12, p. 490]. We therefore consider it is diFcult to realize the stability=instability problem for pulses numerically. However, in this section, at least, suggest what happens in the idealized in0nite interval situation. The model problem (1.1) was tested on four examples. For the four examples sequences of runs with Nx = Nt = 8; 16; 24; 32 are reported. The problem included here illustrate various features of the method, demonstrating the ease of implementation and assembly of the discrete system, the choice of parameters and the exponential convergence rate. Example 4.1. We solve the model equation ut + 6uux + uxxx = 0;
x∈R
(4.1)
290
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292 Table 1 The error in the approximate solution for Example 4.1 N x = Nt
T0 = 3 supuij − u(xi ; tj )
8 16 24 32
1:2698 × 10−3 5:8980 × 10−4 5:2229 × 10−4 4:0157 × 10−4
Table 2 Results for Example 4.2 using Sinc approach N x = Nt
hx = ht
supuij − u(xi ; tj )
8 16 24 32
0.785 0.555 0.453 0.392
1:5040 × 10−2 2:2252 × 10−3 6:8275 × 10−6 1:0636 × 10−6
with the initial data of a single soliton, say u(x; 0) = 2 sech2 x: Here the exact solution is given by u(x; t) = 2 sech2 (x − 4t). We solve this model equation using our approach when d = =2; ! = 1. Table 1 shows that the method converges for T0 = 3. The second column reports the supremum norm of the error between the exact solution u(xi ; tj ) and the approximate solution uij . Again the exponential convergence rate decreases the error that are initially present for small Nx ; Nt : Example 4.2. Let us consider the same problem in the last example, and solve the Eq. (4.1) with the initial condition u(x; 0) = 6 sech2 x. In this case the exact solution of Eq. (4.1) is given by 12[3 + 4 cosh(2x − 8t) + cosh(4x − 64t)] u(x; t) = : [3 cosh(x − 28t) + cosh(3x − 36t)]2 We compared the performance of the method presented here by the 0nite element scheme. Finite element methods for the Korteweg–de Vries equation have been analyzed by Winther [15], where he proved convergence for a class of these equations. Tables 2 and 3 refer to a single-soliton solution and correspond, respectively, to the Sinc method and to the 0nite element method. The tables display information at T0 = 3. We recall that, even though the tables correspond to T0 = 3, the integration was followed up to T0 = 20 and there was no problems of divergence of our scheme. Comparing the errors for the two schemes, we see that for the 0nite element method we do have an error of order O(h2 ), while for the Sinc method it is O(exp(−c=h)) where h is the step size and c ¿ 0. And as predicted, the method using Sinc basis functions is more accurate for soliton propagation. Example 4.3. The known exact one-soliton solutions of Eq. (4.1) exhibit exponential decay for large x. The inverse scattering transform [3] tells us that the solution of (4.1) with generic initial
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292
291
Table 3 Results for Example 4.2 using 0nite elements method hx = ht
supuij − u(xi ; tj )
0.785 0.555 0.453 0.392
2:1102 × 10−2 6:4541 × 10−3 8:5009 × 10−4 5:2001 × 10−5
Table 4 Results for Example 4.3 using Sinc approach N x = Nt
hx = ht
supuij − u(xi ; tj )
8 16 24 32
0.785 0.555 0.453 0.392
8:5423 × 10−3 5:7101 × 10−3 1:3652 × 10−3 8:8602 × 10−5
condition u(x; 0) = u0 (x) (with u0 (x) → 0 suFciently rapidly as |x| → ∞) consists of a train of solitons moving to the right, along with a dispersive wave travelling to the left. As in the 0 above two examples, we can argue that when √ u (x) has compact support, the solution emerge from the rightmost soliton decaying as exp(− cx), where c is the relevant speed, but for suFciently large x the presence of the solitons will not yet be felt, and the behavior of the solution will be determined by the √ Green’s function of the linearized equation ut + uxxx = 0, i.e., u will decay roughly as exp(−2x3=2 =3 3t). Thus there is a transition in the nature of the decay. If there is more than one soliton in the soliton train, say two, with speeds c1 ; c2 (c2 ¿ c1 ), the a problem arises. The √ solution emerges from the faster-moving soliton decaying as exp(− c2 x), but because the tail of the slower-moving soliton falls slower, it is possible that it will return to dominate, i.e., the decay √ will slow to exp(− c1 x). We note that the exact two-soliton solutions (c2 − c1 )(c2 cosh2 !1 + c1 sinh2 !2 ) ; √ 2( c2 cos h!1 cos g!2 − c1 sin h!1 sin h!2 )2 √ √ where !1 = c1 (x − c1 t)=2 and !2 = c2 (x − c2 t)=2 exhibit exactly this phenomenon: when x is large √ u ∼ exp(− c1 x), i.e., the tail of the slow soliton dominates the decay. Table 4 refers to two-soliton solution with c1 = 6; c2 = 10. The table displays information at T0 = 3, with ! = 1; d = =2. u(x; t) =
√
Example 4.4. Vanden-Broeck [16] considered the time-independent surface water of an incompressible and inviscid Puid Pow over a bump in a two-dimensional channel. The model equation is the forced Korteweg–de Vries equation (fkdv): ut + 7ux + 2!uux + uxxx = fx (x);
x ∈ R; t ¿ 0:
The forcing represented by the function f(x) in the fkdv equation is due to the bump on the bottom of the channel. Solitary wave means that the free-surface elevation u(x) has the property u(∓∞) = ux (∓∞) = uxx (∓∞) = 0:
292
K. Al-Khaled / Journal of Computational and Applied Mathematics 130 (2001) 283–292 Table 5 Results for Example 4.4 using Sinc approach N x = Nt
hx = ht
supuij − u(xi ; tj )
8 16 24 32
0.785 0.555 0.453 0.392
4:1314 × 10−2 1:6575 × 10−2 3:5789 × 10−3 2:4625 × 10−4
To illustrate our numerical results, let us take ! = − 34 ; = − 16 ; 7 = 3, and let the forcing function f be de0ned by the equation
f(x) =
sin(x); |x| ¡ 1; 0; |x| ¿ 1:
We use our approach and compare our result with the solution in [6]. Table 5 displays some results, where the solution u(xi ; tj ) is assumed to be the solution in [6], and uij is our Sinc solution. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
K. Al-Khaled, A Sinc–Galerkin approach on the p-system, Saitama Math. J. 16 (1998) 1–13. R. Courant, D. Hilbert, Methods of Mathematical Physics, Wiley-Interscience, New York, 1962. P.G. Drazin, R.S Johnson, Solitons: An Introduction, Cambridge University Press, Cambridge, 1989. N.M. Ercolani, D.W. Mclaughlin, H. Roitner, Attractors and transients for a perturbed Kdv equation, A nonlinear spectral analysis, J. Nonlinear Sci. 3 (1993) 477–539. Feng Bao-Feng, M. Taketomo, A 0nite di@erence method for the Korteweg– de Vries and the Kadomtsev–Petviashvili equations, J. Comput. Appl. Math. 1 (90) (1998) 97–118. L. Gong, S.S. Shen, Multiple supercritical solitary wave solution of the stationary forced Korteweg– de Vries equation and their stability, SIAM J. Appl. Math. 54 (5) (1994) 1266–1290. R. Grimshaw, H. Mitsudera, Slowly varying solitary wave solutions of the perturbed Korteweg– de Vries equation revisited, Stud. Appl. Math. 90 (1993) 75–86. T. Kawahara, S. Toh, Nonlinear dispersive periodic waves in the presence of instability and damping, Phys. Fluids 28 (1985) 1636–1638. Y. Kodama, M. Ablowitz, Perturbations of solitons and solitary waves, Stud. Appl. Math. 64 (1994) 225–245. D.J. Logan, An Introduction to Nonlinear Partial Di@erential Equations, Wiley, New York, 1994. T. Ogawa, H. Suzuki, On the spectra of pulses in a nearly integrable system, SIAM J. Appl. Math. 57 (2) (1997) 485–500. T. Ogawa, Traveling wave solutions to a perturbed Korteweg– de Vries equation, Hiroshima Math. J. 24 (1994) 401–422. F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer, New York, 1993. G.B. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, New York, 1974. R. Winther, A conservative 0nite element method for the Korteweg– de Vries equation, Math. Comput. 34 (1980) 23–43. J.M. Vanden-Broeck, Free-surface Pow over a semi-circular obstruction in a channel, Phys. Fluids 30 (1987) 2315– 2317.