- Email: [email protected]
Local discontinuous Galerkin methods for the Kuramoto–Sivashinsky equations and the Ito-type coupled KdV equations
Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447 www.elsevier.com/locate/cma
Local discontinuous Galerkin methods for the Kuramoto–Sivashinsky equations and the Ito-type coupled KdV equations q Yan Xu a
a,*
, Chi-Wang Shu
b
Department of Mathematics, University of Science and Technology of China, Hefei, Anhui 230026, PR China b Division of Applied Mathematics, Brown University, Providence, RI 02912, USA Accepted 17 June 2005
Abstract In this paper we develop a local discontinuous Galerkin method to solve the Kuramoto–Sivashinsky equations and the Ito-type coupled KdV equations. The L2 stability of the schemes is obtained for both of these nonlinear equations. We use both the traditional nonlinearly stable explicit high order Runge–Kutta methods and the explicit exponential time differencing method for the time discretization; the latter can achieve high order accuracy and maintain good stability while avoiding the very restrictive explicit stability limit of the former when the PDE contains higher order spatial derivatives. Numerical examples are shown to demonstrate the accuracy and capability of these methods. Ó 2005 Elsevier B.V. All rights reserved. MSC classification: 65M60; 35Q53 Keywords: Local discontinuous Galerkin method; Kuramoto–Sivashinsky equation; Ito-type coupled KdV equations; Stability; Exponential time differencing method
q
Research partially supported by the Chinese Academy of Sciences grant 2004-1-8 while the author was in residence at the University of Science and Technology of China. Additional support is provided by ARO grant DAAD19-00-1-0405, NSF grant DMS0207451 and AFOSR grant F49620-02-1-0113. * Corresponding author. E-mail addresses: [email protected] (Y. Xu), [email protected] (C.-W. Shu). 0045-7825/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.06.021
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3431
1. Introduction In this paper we develop a local discontinuous Galerkin method to solve two classes of nonlinear wave equations formulated by the Kuramoto–Sivashinsky type equations ut þ f ðuÞx ðaðuÞux Þx þ ðr0 ðuÞgðrðuÞx Þx Þx þ ðsðux Þuxx Þxx ¼ 0;
ð1:1Þ
and the Ito-type coupled KdV equations ut þ auux þ bvvx þ cuxxx ¼ 0; vt þ bðuvÞx ¼ 0;
ð1:2Þ
where f(u), a(u) P 0, r(u), s(p) P 0 and g(p) are arbitrary (smooth) nonlinear functions and a, b and c are arbitrary constants. The schemes we present extend the previous work of Yan and Shu [24,25] on local discontinuous Galerkin methods solving partial differential equations with higher order spatial derivatives. The Kuramoto–Sivashinsky equation ut þ uux þ auxx þ buxxxx ¼ 0;
ð1:3Þ
which is a special case of (1.1), is a canonical evolution equation which has attracted considerable attention over the last decades. When the coefficients a and b are both positive, its linear terms describe a balance between long-wave instability and short-wave stability, with the nonlinear term providing a mechanism for energy transfer between wave modes. It arises in a broad spectrum of contexts and admits various fascinating solutions like traveling waves of permanent form and chaos [14,19]. It is one of the simplest partial differential equations which is capable of exhibiting chaotic behavior. The chaotic behavior typically occurs when (1.3) is integrated over finite x-domain with periodic boundary conditions. The well-posedness of Eq. (1.3) was proved in [21]. For numerical study, Hooper and Grimshaw [13] identified three classes of traveling wave solutions of the Kuramoto–Sivashinsky equation and labeled them as regular shocks, oscillatory shocks and solitary waves. In [26], Yang further studied the traveling-wave solutions of the Kuramoto– Sivashinsky equation. Coupled nonlinear equations in which a KdV structure is embedded occur naturally in shallow water wave problems [22]. For a = 6, b = 2, c = 1, d = 2, Eq. (1.2) represents the ItoÕs equation, which describes the interaction process of two internal long waves and has infinitely many conserved quantities [15]. It is noted that this choice of parameters is not unique [1]. In [17], it is shown that this coupled system can be a member of a bi-Hamiltonian integrable hierarchy. The discontinuous Galerkin (DG) method we discuss in this paper is a class of finite element methods using completely discontinuous piecewise-polynomial space for the numerical solution and the test functions in the spatial variables, usually coupled with explicit and nonlinearly stable high order Runge–Kutta time discretization [20]. It was first developed for hyperbolic conservation laws containing first derivatives by Cockburn et al. in a series of papers [8,7,5,9]. For a detailed description of the method as well as its implementation and applications, we refer the readers to the lecture notes [4], the survey paper [6], other papers in that Springer volume, and the review paper [11]. These discontinuous Galerkin methods were generalized to solve a convection diffusion equation (containing second derivatives) by Cockburn and Shu [10]. Their work was motivated by the successful numerical experiments of Bassi and Rebay [2] for the compressible Navier–Stokes equations. Later, Yan and Shu developed a local discontinuous Galerkin method for a general KdV type equation containing third derivatives in [24], and they generalized the local discontinuous Galerkin method to PDEs with fourth and fifth spatial derivatives in [25]. Levy et al. [18] developed local discontinuous Galerkin methods for solving nonlinear dispersive equations that have compactly supported traveling wave solutions, the so-called
3432
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
‘‘compactons’’. Recently, Xu and Shu [23] further developed the local discontinuous Galerkin method to solve three classes of nonlinear wave equations formulated by the general KdV–Burgers type equations, the general fifth-order KdV type equations and the fully nonlinear K(n, n, n) equations. These discontinuous Galerkin methods have several attractive properties. It can be easily designed for any order of accuracy. In fact, the order of accuracy can be locally determined in each cell, thus allowing for efficient p adaptivity. It can be used on arbitrary triangulations, even those with hanging nodes, thus allowing for efficient h adaptivity. The methods have excellent parallel efficiency. It is extremely local in data communications. The evolution of the solution in each cell needs to communicate only with the immediate neighbors, regardless of the order of accuracy. Finally, it has excellent provable nonlinear stability. One can typically prove a strong L2 stability and a cell entropy inequality for the square entropy, for quite general nonlinear cases, for any orders of accuracy on arbitrary triangulations in any space dimension, without the need for nonlinear limiters. The paper is organized as follows. In Section 2, we present and analyze the local discontinuous Galerkin methods for Eqs. (1.1) and (1.2). In Section 2.1, we present the methods for the Kuramoto–Sivashinsky type equations. We prove a theoretical result of the L2 stability for the nonlinear case. In Section 2.2, we present the local discontinuous Galerkin methods for Ito-type coupled KdV equations and give a theoretical result of the L2 stability. In Section 3 we describe the exponential time differencing method for time discretization, which can achieve high order accuracy and maintain good stability. Section 4 contains numerical results for the nonlinear problems to demonstrate the accuracy and capability of the methods. Concluding remarks are given in Section 5.
2. The local discontinuous Galerkin method 2.1. The Kuramoto–Sivashinsky type equations In this section, we present and analyze a local discontinuous Galerkin method for the following Kuramoto–Sivashinsky type equations ut þ f ðuÞx ðaðuÞux Þx þ ðr0 ðuÞgðrðuÞx Þx Þx þ ðsðux Þuxx Þxx ¼ 0;
ð2:1Þ
with an initial condition uðx; 0Þ ¼ u0 ðxÞ
ð2:2Þ
and periodic boundary conditions uðx; tÞ ¼ uðx þ L; tÞ;
ð2:3Þ
where L is the period in the x direction. Here f(u), a(u) P 0, r(u), s(p) P 0 and g(p) are arbitrary (smooth) nonlinear functions. Notice that the assumption of periodic boundary conditions is for simplicity only and is not essential: the method can be easily designed for non-periodic boundary conditions and such non-periodic boundary conditions are used in the numerical experiments in next section. We denote the mesh by I j ¼ ½xj12 ; xjþ12 , for j = 1, . . . , N. The center of the cell is xj ¼ ðxj12 þ xjþ12 Þ=2 and Dxj ¼ xjþ12 xj12 . We denote by uþ and u the values of u at xjþ12 , from the right cell, Ij + 1, and from the jþ1 jþ1 2
2
left cell, Ij, respectively. We define the piecewise-polynomial space VDx as the space of polynomials of the degree k in each cell Ij, i.e. V Dx ¼ fv :
v 2 P k ðI j Þ for x 2 I j ; j ¼ 1; . . . ; N g:
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3433
To define the local discontinuous Galerkin method, we rewrite Eq. (2.1) as a first-order system: ut þ f ðuÞx ðbðuÞvÞx þ ðr0 ðuÞpÞx þ wx ¼ 0; ð2:4Þ v BðuÞx ¼ 0; p gðqÞx ¼ 0; q rðuÞx ¼ 0; w ðlðzÞhÞx ¼ 0; h QðzÞx ¼ 0; z ux ¼ 0; pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi where bðuÞ ¼ aðuÞ, B(u) = òub(s)ds, lðzÞ ¼ sðzÞ, Q(z) = òzl(s)ds. Now we can use the local discontinuous Galerkin method to Eq. (2.4), resulting in the following scheme: find u, v, p, q, w, h, z 2 VDx such that, for all test functions q, /, u, w, g, n, f 2 VDx, Z Z ut q dx ðf ðuÞ bðuÞv þ r0 ðuÞp þ wÞqx dx Ij
Ij
^ ^v þ r^0 ^p þ w ^ Þjþ1 q ^ Þj1 qþ þ ðf^ ^ b^v þ r^0 ^ pþw 1 ðf b^ 1 ¼ 0; 2 jþ2 2 j2 Z Z ^ jþ1 / 1 þ B ^ j1 /þ 1 ¼ 0; v/ dx þ BðuÞ/x dx B 2 jþ2 2 j2 I I Zj Zj pu dx þ gðqÞux dx g^jþ12 u þ g^j12 uþ ¼ 0; jþ12 j12 Ij Ij Z Z qw dx þ rðuÞwx dx ^rjþ12 w rj12 wþ ¼ 0; jþ12 þ ^ j12 Ij Ij Z Z þ ^^ wg dx þ lðzÞhgx dx ð^l^ hÞjþ1 g 1 þ ðlhÞj1 g 1 ¼ 0; 2 jþ2 2 j2 Ij Ij Z Z ^ jþ1 n 1 þ Q ^ j1 nþ 1 ¼ 0; hn dx þ QðzÞnx dx Q 2 jþ2 2 j2 Ij Ij Z Z zf dx þ ufx dx ^ ujþ12 f uj12 fþ ¼ 0: jþ1 þ ^ j1 Ij
Ij
2
ð2:5Þ
2
The ‘‘hat’’ terms in (2.5) in the cell boundary terms from integration by parts are the so-called ‘‘numerical fluxes’’, which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure stability. It turns out that we can take the simple choices such that f^ ¼ f^ ðu ; uþ Þ; ^r ¼ rðu Þ;
^ ¼ Bðu Þ; ^v ¼ vþ ; g^ ¼ g^ðq ; qþ Þ; B ^ ¼ Qðzþ Þ; ^ ^p ¼ pþ ; Q ^ ¼ wþ ; h ¼ h ; w
Bðuþ Þ Bðu Þ ^ ; b¼ uþ u
rðuþ Þ rðu Þ r^0 ¼ ; uþ u
^u ¼ u ;
ð2:6Þ
þ ^l ¼ Qðz Þ Qðz Þ ; þ z z
where we have omitted the half-integer indices j þ 12 as all quantities in (2.6) are computed at the same points (i.e. the interfaces between the cells). Here f^ ðu ; uþ Þ and ^ gðq ; qþ Þ are monotone fluxes, i.e. Lips^ chitz continuous in both arguments, consistent (i.e. f ðu; uÞ ¼ f ðuÞ), non-decreasing in the first argument and non-increasing in the second argument. Examples of monotone fluxes which are suitable for discontinuous Galerkin methods can be found in, e.g., [8]. We could for example use the simple Lax–Friedrichs flux 1 f^ ðu ; uþ Þ ¼ ðf ðu Þ þ f ðuþ Þ aðuþ u ÞÞ; 2
a ¼ max jf 0 ðuÞj;
where the maximum is taken over a relevant range of u. This flux is used in the numerical experiments in next section. The algorithm is now well defined. ^ and ^v from We remark that the choice for the fluxes (2.6) is not unique. In fact the crucial part is taking B ^ ^ ^ and ^u opposite sides, taking ^r and ^ p from opposite sides, taking Q and h from opposite sides and taking w from opposite sides.
3434
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
With such a choice of fluxes we can get the theoretical results of the L2 stability. ^ jþ1 such that the solution to Proposition 2.1 (cell entropy inequality). There exist numerical entropy fluxes U 2 the scheme (2.5), (2.6) satisfies Z Z 1 d 2 ^ ^ u ðx; tÞdx þ Ujþ12 Uj12 þ ðh2 þ v2 Þdx 6 0: 2 dt I j Ij Proof. Since (2.5) holds for any test functions in VDx, we can choose q ¼ u;
/ ¼ v;
u ¼ q;
w ¼ p;
g ¼ z;
n ¼ h;
f ¼ w:
With these choices of test functions and summing up the seven equations in (2.5) we obtain Z Z ^ jþ1 U ^ j1 þ Hj1 þ ðh2 þ v2 Þdx ¼ 0; ut u dx þ U 2 2 2 Ij
Ij
where the numerical entropy flux is given by ^ ¼ F ðu Þ þ Gðq Þ þ ðf^ ^ U b^v þ r^0 ^ p þ wþ Þu g^q ^l^hz ðrðu Þ ^rÞp ^ ^ þ ðQðz Þ QÞh þ ðBðu Þ BÞv and the extra term H is given by ^ H ¼ ½F ðuÞ f^ ½u ½GðqÞ þ g^½q þ ½rðuÞp ^r½p r^0 ^p½u ½BðuÞv þ ^b^v½u þ B½v ^ ½QðzÞh þ ^l^ h½z þ Q½h: Here F ðuÞ ¼
Z
u
f ðuÞdu;
GðqÞ ¼
Z
q
gðqÞdq
and [v] = v+ v denotes the jump of v. With the definition (2.6) of the numerical fluxes and with simple algebraic manipulations, we easily obtain ^ ¼ 0; ½BðuÞv þ ^ b^v½u þ B½v
p½u ¼ 0; ½rðuÞp ^r½p r^0 ^
^ ¼0 ½QðzÞh þ ^l^h½z þ Q½h
and hence H ¼ ½F ðuÞ f^ ½u ½GðqÞ þ g^½q ¼
Z
uþ
ðf ðsÞ f^ ðu ; uþ ÞÞds
u
Z
qþ
ðgðsÞ g^ðq ; qþ ÞÞds P 0; q
where the last inequality follows from the monotonicity of the fluxes f^ and ^ g. Hence Z Z ^ jþ1 U ^ j1 þ ðh2 þ v2 Þdx 6 0; ut u dx þ U 2 2 Ij
Ij
which is the cell entropy inequality.
h
Summing up the cell entropy inequalities, we obtain Corollary 2.2 (L2 stability). The solution to the scheme (2.5), (2.6) satisfies the L2 stability Z L d u2 ðx; tÞdx 6 0: dt 0
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3435
2.2. Ito-type coupled equations In this section, we present and analyze a local discontinuous Galerkin method for the following Ito-type coupled nonlinear problem ut þ auux þ bvvx þ cuxxx ¼ 0;
ð2:7Þ
vt þ bðuvÞx ¼ 0; with an initial condition uðx; 0Þ ¼ u0 ðxÞ;
vðx; 0Þ ¼ v0 ðxÞ
ð2:8Þ
and periodic boundary conditions uðx; tÞ ¼ uðx þ L; tÞ;
vðx; tÞ ¼ vðx þ L; tÞ;
ð2:9Þ
where L is the period in the x direction. Notice that the assumption of periodic boundary conditions is for simplicity only and is not essential: the method can be easily designed for non-periodic boundary conditions. To define the local discontinuous Galerkin method, we rewrite Eq. (2.7) as a first-order system ut þ f ðu; vÞx þ px ¼ 0; p cqx ¼ 0;
q ux ¼ 0;
ð2:10Þ
vt þ gðu; vÞx ¼ 0; where f ðu; vÞ ¼ a2 u2 þ b2 v2 and g(u, v) = buv. Now we can define the local discontinuous Galerkin method to Eq. (2.10), resulting in the following scheme: find u, p, q, v 2 VDx such that, for all test functions q, /, u, w 2 VDx, Z Z ut q dx ðf ðu; vÞ þ pÞqx dx þ ðf^ þ ^ pÞjþ1 q ðf^ þ ^pÞj1 qþ ¼ 0; jþ1 j1 Z Z Z
Ij
p/ dx þ c Ij
qu dx þ Ij
Ij
Z Z
Ij
2
2
2
q/x dx c^ qjþ12 / qj12 /þ ¼ 0; jþ1 þ c^ j1
uux dx Ij
vt w dx Ij
2
Ij
Z
2
2
ð2:11Þ ^ ujþ12 u jþ12
þ
^ uj12 uþ j12
¼ 0;
^j1 wþ gðu; vÞwx dx þ g^jþ12 w 1 ¼ 0: jþ1 g 2 j 2
2
The ‘‘hat’’ terms in (2.11) are the numerical fluxes. It turns out that we can again take the simple choices such that 1 1 f^ ¼ ðf þ þ f kðuþ u ÞÞ; g^ ¼ ðgþ þ g kðvþ v ÞÞ; 2 2 q if c < 0; þ ^ p¼p ; ^ u¼u ; ^ q¼ qþ if c > 0; where k¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 max ða þ bÞu ða bÞ u2 þ 4b2 v2 ; ða þ bÞu þ ða bÞ u2 þ 4b2 v2 2
ð2:12Þ
3436
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
is the maximum eigenvalue of the Jacobian of the flux (f, g)T. The algorithm is now well defined. Notice that f^ and g^ are the standard Lax–Friedrichs fluxes in this system case. We remark that the choice for the fluxes (2.12) is not unique. In fact the crucial part is taking ^p and ^u from opposite sides. With such a choice of fluxes we can get the theoretical results of the L2 stability. Proposition 2.3 (L2 stability). The solution to the scheme (2.11), (2.12) satisfies the L2 stability Z L d ðu2 þ v2 Þdx 6 0: dt 0 Proof. We prove the case for c < 0 and take c = 1, the case for c > 0 can be proved in a similar way. Since (2.11) holds for any test functions in VDx, we can choose q ¼ u;
/ ¼ q;
u ¼ p;
w ¼ v:
With these choices of test functions and summing up the four equations in (2.11) we obtain Z ^ jþ1 H ^ j1 þ Hj1 ¼ 0; ðut u þ vt vÞdx þ H 2 2 2
ð2:13Þ
Ij
where the numerical entropy flux is given by ^ ¼ pþ u þ 1 ðq2 Þ a ðu3 Þ b ðuv2 Þ þ f^ u þ g^v H 2 6 2 and the extra term H is given by 1 2 1 2 2 H ¼ ½q þ ð½u ða½u þ 6kÞ þ 3½v ðb½u þ 2kÞÞ: 2 12 Here [v] = v+ v denotes the jump of v. Note that 1 k P ðja þ bj þ ja bjÞjuj ¼ maxfjaj; jbjgjuj: 2 Hence we obtain ½u2 ða½u þ 6kÞ þ 3½v2 ðb½u þ 2kÞ P 0: This implies H P 0. Summing up (2.13) over j and taking into account the periodic boundary condition, we obtain the desired L2 stability. h
3. Exponential time differencing methods The Runge–Kutta discontinuous Galerkin method is originally designed with explicit and nonlinearly stable high order Runge–Kutta time discretizations [20]. This is adequate for hyperbolic problems with only first derivatives, but the explicit stability limit for the time step decreases dramatically for problems containing higher order spatial derivatives. It is thus important to consider other more efficient time discretization methods. In this section we describe the exponential time differencing fourth-order Runge–Kutta (ETD4RK) method which was developed by Cox and Matthews in [12], among general formulas for ETD-schemes to arbitrary order (also called ELP-schemes [3]). In a recent paper Kassam and Trefethen [16] compared various fourth-order methods for solving equations of the form (3.2), and conclude that the best by a clear margin was a modification by using the complex integral to evaluate the coefficients
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3437
in the update formula. In essence, for nonlinear time dependent equations, the ETD schemes provide a systematic coupling of the explicit treatment of nonlinearities and the implicit and possibly exact integration of the stiff linear parts of the equations, while achieving high order accuracy and maintaining good stability. We consider partial differential equations of the form ut ¼ Lu þ Nðu; tÞ;
ð3:1Þ
where L is a linear term and N is nonlinear. Many interesting equations are of this form, where typically L represents the stiff part of the equation. When discretizing in space we obtain a system of ODEs of the form u_ ¼ Lu þ N ðu; tÞ:
ð3:2Þ
The exponential time differencing (ETD) methods can be described in the context of solving (3.2). Integrating the equation over a single time step from t = tn to tn+1 = tn + Dt, we obtain Z Dt uðtnþ1 Þ ¼ eLDt uðtn Þ þ eLDt eLs N ðuðtn þ sÞ; tn þ sÞds: ð3:3Þ 0
Eq. (3.3) is exact, and the various ETD schemes come from the approximations to the integral [12]. Cox and Matthews derived a set of ETD methods based on Runge–Kutta time-stepping, which they called ETDRK schemes. For completeness, we give below the formulae for the fourth-order scheme (ETD4RK): an ¼ eLDt=2 un þ L1 ðeLDt=2 IÞN ðun ; tn Þ; bn ¼ eLDt=2 un þ L1 ðeLDt=2 IÞN ðan ; tn þ Dt=2Þ; cn ¼ eLDt=2 an þ L1 ðeLDt=2 IÞð2N ðbn ; tn þ Dt=2Þ N ðun ; tn ÞÞ;
ð3:4Þ
2
unþ1 ¼ eLDt un þ Dt2 L3 f½4 LDt þ eLDt ð4 3LDt þ ðLDtÞ ÞN ðun ; tn Þ þ 2½2 þ LDt þ eLDt ð2 þ LDtÞðN ðan ; tn þ Dt=2Þ þ N ðbn ; tn þ Dt=2ÞÞ þ ½4 3LDt ðLDtÞ2 þ eLDt ð4 LDtÞN ðcn ; tn þ DtÞg: To overcome the vulnerability of the error cancelations in the high order ETD and ETDRK schemes, and to generalize the ETD schemes to non-diagonal problems, modified ETD schemes are proposed in [16] by using the complex contour integrals Z 1 f ðLÞ ¼ f ðzÞðzI LÞ1 dz ð3:5Þ 2p C on a suitable contour C to evaluate the coefficients in the update formula for ETD4RK. The explicit exponential time differencing method for time discretization can achieve high order accuracy and maintain good stability while avoiding the severe explicit stability time step limit of the traditional Runge–Kutta discontinuous Galerkin methods which use explicit and nonlinearly stable high order Runge–Kutta time discretization due to the presence of the high order derivative terms.
4. Numerical results In this section we provide numerical examples to illustrate the accuracy and capability of the method. Time discretization is by the third order explicit Runge–Kutta method in [20] or by the exponential time differencing fourth-order Runge–Kutta method in [16]. The time step of the exponential time differencing fourth-order Runge–Kutta methods is Dt = O(Dx), and the time step of the third-order Runge–Kutta
3438
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Table 4.1 Accuracy test for the Kuramoto–Sivashinsky equation (4.1) with the exact solution (4.2) N
L2 error
Order
L1 error
Order
0
p
40 80 160 320
2.38E01 1.58E01 9.41E02 5.24E02
– 0.60 0.74 0.84
1.37 8.81E01 5.21E01 2.91E01
– 0.64 0.76 0.84
p1
40 80 160 320
6.08E02 1.49E02 3.78E03 9.57E04
– 2.02 1.98 1.98
6.64E01 1.82E01 4.64E02 1.19E02
– 1.87 1.97 1.96
p2
40 80 160 320
9.15E03 1.06E03 1.32E04 1.65E05
3.45 3.11 3.00 3.00
1.49E01 1.73E02 2.43E03 3.05E04
2.58 3.11 2.83 2.99
p3
20 40 80 160
3.45E02 9.42E04 7.87E05 5.01E06
– 5.20 3.58 3.97
2.62E01 1.12E02 1.30E03 8.78E05
– 4.55 3.11 3.89
c = 0.1, k ¼ 12
qffiffiffiffi 11 19, x0 = 10. Exact boundary condition. Uniform meshes with N cells at time t = 1.
Fig. 4.1. The shock profile wave propagation of the Kuramoto–Sivashinsky equation (4.1) with initial condition (4.4), c = 5, k ¼ 12 x0 = 12. Exact boundary condition in [30, 30], P2 elements with 160 cells.
qffiffiffiffi 11 , 19
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3439
Table 4.2 Accuracy test for the Kuramoto–Sivashinsky equation (4.5) with the exact solution (4.6) N
L2 error
Order
L1 error
Order
0
p
40 80 160 320
1.33E02 6.69E03 3.35E03 1.68E03
– 0.98 1.00 1.00
7.56E02 3.87E02 1.95E02 9.80E03
– 0.92 0.98 1.00
p1
40 80 160 320
1.29E03 3.26E04 8.17E05 2.05E05
– 1.99 2.00 2.00
1.04E02 2.59E03 6.54E04 1.64E04
– 2.01 1.98 2.00
p2
40 80 160 320
6.48E05 8.04E06 1.00E06 1.25E07
– 3.01 3.00 3.00
7.93E04 1.06E04 1.34E05 1.67E06
– 2.91 2.99 3.00
p3
40 80 160 320
5.92E06 3.83E07 2.42E08 1.64E09
– 3.95 3.99 3.88
6.07E05 4.00E06 2.56E07 1.61E08
– 3.92 3.97 3.99
ffi, x0 = 10. Exact boundary condition. Uniform meshes with N cells at time t = 1. c = 0.2, k ¼ 2p1ffiffiffi 19
ffi, x0 = 25. Exact Fig. 4.2. The shock profile wave propagation of the KS equation (4.5) with initial condition (4.8), c = 5, k ¼ 2p1ffiffiffi 19 boundary condition in [50, 50], P2 elements with 200 cells.
3440
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Table 4.3 Accuracy test for the Kuramoto–Sivashinsky equation (4.9) with the exact solution (4.10) N
L2 error
Order
L1 error
Order
1
p
80 160 320 640
1.65E01 2.33E02 4.68E03 1.28E03
– 2.82 2.32 1.87
1.01 2.01E01 6.49E02 1.83E02
– 2.32 1.63 1.82
p2
80 160 320 640
5.98E03 7.43E04 9.30E05 1.16E05
– 3.01 3.00 3.00
1.07E01 1.44E02 1.84E03 2.31E04
– 2.89 2.97 3.00
p3
40 80 160 320
1.10E02 6.84E04 4.36E05 2.75E06
– 4.00 3.97 3.99
9.31E02 1.18E02 7.36E04 5.01E05
– 2.98 4.01 3.88
c = 6, r = 4, k ¼ 12, x0 = 10. Periodic boundary condition. Uniform meshes with N cells at time t = 1.
Fig. 4.3. The solitary wave propagation of the Kuramoto–Sivashinsky equation (4.9) with initial condition (4.11), c = 6, r = 4, k ¼ 12, x0 = 10. Periodic boundary condition in [30, 30], P2 elements with 250 cells.
method is Dt = O(Dx4) for the KS equation and Dt = O(Dx3) for the ItoÕs equation. We have verified with the aid of successive mesh refinements, that in all cases, the results shown are numerically convergent.
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3441
Example 4.1. We show an accuracy test for the Kuramoto–Sivashinsky equation ut þ uux þ uxx þ uxxxx ¼ 0 with the exact solution rffiffiffiffiffi 15 11 uðx; tÞ ¼ c þ 9 tanhðkðx ct x0 ÞÞ þ 11tanh3 ðkðx ct x0 ÞÞ 19 19
ð4:1Þ
ð4:2Þ
and the boundary condition uð30; tÞ ¼ g1 ðtÞ;
uð30; tÞ ¼ g2 ðtÞ;
ð4:3Þ
where gi(t) corresponds to the data from the exact solution. Notice that the local discontinuous Galerkin method allows for an easy implementation of such boundary conditions. The L2 and L1 errors and the numerical orders of accuracy are contained in Table 4.1. We can see that the method with Pk elements gives a uniform (k + 1)th order of accuracy in both norms.
Fig. 4.4. Chaotic solution of the Kuramoto–Sivashinsky equation (4.12) with initial condition (4.13). Periodic boundary condition in [0, 32p], P2 elements with two different meshes using N = 300 and N = 600 uniform cells.
3442
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
We also present the regular shock profile wave propagation for the Kuramoto–Sivashinsky equation (4.1) in Fig. 4.1 with the initial condition rffiffiffiffiffi 15 11 uðx; 0Þ ¼ c þ ð4:4Þ 9 tanhðkðx x0 ÞÞ þ 11tanh3 ðkðx x0 ÞÞ : 19 19 We can see clearly that the moving shock profile is resolved very well. Example 4.2. We show an accuracy test for the Kuramoto–Sivashinsky equation ut þ uux uxx þ uxxxx ¼ 0
ð4:5Þ
with the exact solution uðx; tÞ ¼ c þ
15 pffiffiffiffiffi 3 tanhðkðx ct x0 ÞÞ þ tanh3 ðkðx ct x0 ÞÞ 19 19
ð4:6Þ
and the boundary condition uð50; tÞ ¼ h1 ðtÞ;
uð50; tÞ ¼ h2 ðtÞ;
ð4:7Þ
Fig. 4.5. The chaotic solution of the Kuramoto–Sivashinsky equation (4.12) with the Gaussian initial condition (4.14). Periodic boundary condition in [16, 16], P2 elements with two different meshes using N = 160 and N = 320 uniform cells.
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3443
where hi(t) corresponds to the data from the exact solution. The result is contained in Table 4.2. The scheme clearly demonstrates an order of accuracy of k + 1 for Pk elements in both L2 and L1 norms for this problem with non-periodic boundary conditions. Next we show the shock profile wave propagation for the Kuramoto–Sivashinsky equation (4.5) in Fig. 4.2 with the initial condition 15 uðx; 0Þ ¼ c þ pffiffiffiffiffi 3 tanhðkðx x0 ÞÞ þ tanh3 ðkðx x0 ÞÞ : ð4:8Þ 19 19 We can again observe that the moving shock profile is resolved very well. Example 4.3. We show an accuracy test for the Kuramoto–Sivashinsky equation ut þ uux þ uxx þ ruxxx þ uxxxx ¼ 0
ð4:9Þ
with the exact solution uðx; tÞ ¼ c þ 9 15 tanhðkðx ct x0 ÞÞ þ tanh2 ðkðx ct x0 ÞÞ tanh3 ðkðx ct x0 ÞÞ :
ð4:10Þ
Table 4.3 give the errors of numerical solution at t = 1 using the periodic boundary condition. We can see from the table that the orders of accuracy are comparable to that for the non-periodic case. Next we show the solitary wave propagation for the Kuramoto–Sivashinsky equation (4.9) in Fig. 4.3 with the initial condition
Fig. 4.6. Numerical results of u for the ItoÕs equation (4.15) with the initial condition (4.16). Periodic boundary condition in [0, 2p], P2 elements with 80 cells.
3444
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Fig. 4.7. Numerical results of v for the ItoÕs equation (4.15) with the initial condition (4.16). Periodic boundary condition in [0, 2p], P2 elements with 80 cells.
uðx; 0Þ ¼ c þ 9 15 tanhðkðx x0 ÞÞ þ tanh2 ðkðx x0 ÞÞ tanh3 ðkðx x0 ÞÞ :
ð4:11Þ
This moving solitary wave is clearly well resolved. Example 4.4. The Kuramoto–Sivashinsky equation ut þ uux þ uxx þ uxxxx ¼ 0
ð4:12Þ
has received a lot of interest because it is one of the simplest partial differential equations which is capable of exhibiting chaotic behavior. The chaotic behavior typically occurs when (4.12) is integrated over finite xdomain with periodic boundary conditions. First, we take the initial condition x x uðx; 0Þ ¼ cos 1 þ sin ð4:13Þ 16 16 and present the results in Fig. 4.4 for two different meshes with N = 300 and N = 600 uniform cells respectively. We can see clearly that the results are numerically convergent for the quite chaotic behavior. Next we take the Gaussian initial condition 2
uðx; 0Þ ¼ ex
ð4:14Þ
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3445
Fig. 4.8. Numerical results of u for the ItoÕs equation (4.15) with the initial condition (4.17). Periodic boundary condition in [15, 15], P2 elements with 160 cells.
and present the results in Fig. 4.5 for two different meshes with N = 160 and N = 320 uniform cells respectively. We can again see clearly that the results are numerically convergent for this case. Example 4.5. In this example, we show the numerical results for the ItoÕs equation ut ð3u2 þ v2 Þx uxxx ¼ 0; vt 2ðuvÞx ¼ 0;
ð4:15Þ
with an initial condition uðx; 0Þ ¼ cosðxÞ;
vðx; 0Þ ¼ cosðxÞ:
ð4:16Þ
We choose an Gaussian initial condition uðx; 0Þ ¼ expðx2 Þ;
vðx; 0Þ ¼ expðx2 Þ:
ð4:17Þ
Notice that there is a dispersive term in the first equation for u in (4.15), while there is no such dispersive term in the second equation for v. From Figs. 4.6–4.8 and 4.9, we observe the result for u behaving like dispersive wave solutions and the result for v behaving like shock wave solutions.
3446
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Fig. 4.9. Numerical results of v for the ItoÕs equation (4.15) with the initial condition (4.17). Periodic boundary condition in [15, 15], P2 elements with 160 cells.
5. Concluding remarks We have developed local discontinuous Galerkin methods to solve the Kuramoto–Sivashinsky type equations and the Ito-type coupled KdV equations and have proven the L2 stability of these methods. The exponential time differencing fourth-order Runge–Kutta methods have been used for treating the nonlinearity and stiffness of the systems in order to improve computational efficiency. Numerical examples are shown to illustrate the accuracy and capability of the methods. Although not addressed in this paper, these methods are flexible for general geometry, unstructured meshes and h-p adaptivity, and have excellent parallel efficiency. They should provide a useful class of numerical tools for solving nonlinear wave equations. An important issue not addressed in this paper is L2 a priori error estimates. From the cell entropy inequality and approximation results, one can derive L2 a priori error estimates following the lines in [10,24], for linear PDEs. The proof is not completely straightforward and sometimes one would lose half an order or even one order in accuracy, because of a lack of control for certain jump terms at cell boundaries. It is more challenging to perform L2 a priori error estimates for nonlinear PDEs. Such error estimates are left for future work. References [1] M. Antonowicz, A.P. Fordy, Coupled KdV equations with multi-Hamiltonian structures, Physica D 28 (1987) 345–357. [2] F. Bassi, S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations, J. Comput. Phys. 131 (1997) 267–279.
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3447
[3] G. Beylkin, J.M. Keiser, L. Vozovoi, A new class of time discretization schemes for the solution of nonlinear PDEs, J. Comput. Phys. 147 (1998) 362–387. [4] B. Cockburn, Discontinuous Galerkin methods for methods for convection-dominated problems, in: T.J. Barth, H. Deconinck (Eds.), High-Order Methods for Computational Physics, Lecture Notes in Computational Science and Engineering, vol. 9, Springer, 1999, pp. 69–224. [5] B. Cockburn, S. Hou, C.-W. Shu, The Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case, Math. Comput. 54 (1990) 545–581. [6] B. Cockburn, G. Karniadakis, C.-W. Shu, The development of discontinuous Galerkin methods, in: B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Lecture Notes in Computational Science and Engineering, vol. 11, Springer, 2000, pp. 3–50, Part I: Overview. [7] B. Cockburn, S.-Y. Lin, C.-W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one dimensional systems, J. Comput. Phys. 84 (1989) 90–113. [8] B. Cockburn, C.-W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Math. Comput. 52 (1989) 411–435. [9] B. Cockburn, C.-W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, J. Comput. Phys. 141 (1998) 199–224. [10] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463. [11] B. Cockburn, C.-W. Shu, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (2001) 173–261. [12] S.M. Cox, P.C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys. 176 (2002) 430–455. [13] A.P. Hooper, R. Grimshaw, Travelling wave solutions of the Kuramoto–Sivashinsky equation, Wave Motion 10 (1988) 405–420. [14] J.M. Hyman, B. Nicolaenko, The Kuramoto–Sivashinsky equation: a bridge between PDEs and dynamical systems, Physica D 18 (1986) 113–126. [15] M. Ito, Symmetries and conservation laws of a coupled nonlinear wave equation, Phys. Lett. A 91 (1982) 335–338. [16] A.K. Kassam, L.N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput. 26 (2005) 1214–1233. [17] B.A. Kupershmidt, A coupled Korteweg–de Vries equation with dispersion, J. Phys. A: Math. Gen. 18 (1985) L571–L573. [18] D. Levy, C.-W. Shu, J. Yan, Local discontinuous Galerkin methods for nonlinear dispersive equations, J. Comput. Phys. 196 (2004) 751–772. [19] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation, Physica D 19 (1986) 89–111. [20] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [21] E. Tadmor, The well-posedness of the Kuramoto–Sivashinsky equation, SIAM J. Math. Anal. 17 (1986) 884–893. [22] G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York, 1974. [23] Y. Xu, C.-W. Shu, Local discontinuous Galerkin methods for three classes of nonlinear wave equations, J. Comput. Math. 22 (2004) 250–274. [24] J. Yan, C.-W. Shu, A local discontinuous Galerkin method for KdV type equations, SIAM J. Numer. Anal. 40 (2002) 769–791. [25] J. Yan, C.-W. Shu, Local discontinuous Galerkin methods for partial differential equations with higher order derivatives, J. Sci. Comput. 17 (2002) 27–47. [26] T.S. Yang, On traveling-wave solutions of the Kuramoto–Sivashinsky equation, Physica D 110 (1997) 25–42.
Local discontinuous Galerkin methods for the Kuramoto–Sivashinsky equations and the Ito-type coupled KdV equations q Yan Xu a
a,*
, Chi-Wang Shu
b
Department of Mathematics, University of Science and Technology of China, Hefei, Anhui 230026, PR China b Division of Applied Mathematics, Brown University, Providence, RI 02912, USA Accepted 17 June 2005
Abstract In this paper we develop a local discontinuous Galerkin method to solve the Kuramoto–Sivashinsky equations and the Ito-type coupled KdV equations. The L2 stability of the schemes is obtained for both of these nonlinear equations. We use both the traditional nonlinearly stable explicit high order Runge–Kutta methods and the explicit exponential time differencing method for the time discretization; the latter can achieve high order accuracy and maintain good stability while avoiding the very restrictive explicit stability limit of the former when the PDE contains higher order spatial derivatives. Numerical examples are shown to demonstrate the accuracy and capability of these methods. Ó 2005 Elsevier B.V. All rights reserved. MSC classification: 65M60; 35Q53 Keywords: Local discontinuous Galerkin method; Kuramoto–Sivashinsky equation; Ito-type coupled KdV equations; Stability; Exponential time differencing method
q
Research partially supported by the Chinese Academy of Sciences grant 2004-1-8 while the author was in residence at the University of Science and Technology of China. Additional support is provided by ARO grant DAAD19-00-1-0405, NSF grant DMS0207451 and AFOSR grant F49620-02-1-0113. * Corresponding author. E-mail addresses: [email protected] (Y. Xu), [email protected] (C.-W. Shu). 0045-7825/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.06.021
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3431
1. Introduction In this paper we develop a local discontinuous Galerkin method to solve two classes of nonlinear wave equations formulated by the Kuramoto–Sivashinsky type equations ut þ f ðuÞx ðaðuÞux Þx þ ðr0 ðuÞgðrðuÞx Þx Þx þ ðsðux Þuxx Þxx ¼ 0;
ð1:1Þ
and the Ito-type coupled KdV equations ut þ auux þ bvvx þ cuxxx ¼ 0; vt þ bðuvÞx ¼ 0;
ð1:2Þ
where f(u), a(u) P 0, r(u), s(p) P 0 and g(p) are arbitrary (smooth) nonlinear functions and a, b and c are arbitrary constants. The schemes we present extend the previous work of Yan and Shu [24,25] on local discontinuous Galerkin methods solving partial differential equations with higher order spatial derivatives. The Kuramoto–Sivashinsky equation ut þ uux þ auxx þ buxxxx ¼ 0;
ð1:3Þ
which is a special case of (1.1), is a canonical evolution equation which has attracted considerable attention over the last decades. When the coefficients a and b are both positive, its linear terms describe a balance between long-wave instability and short-wave stability, with the nonlinear term providing a mechanism for energy transfer between wave modes. It arises in a broad spectrum of contexts and admits various fascinating solutions like traveling waves of permanent form and chaos [14,19]. It is one of the simplest partial differential equations which is capable of exhibiting chaotic behavior. The chaotic behavior typically occurs when (1.3) is integrated over finite x-domain with periodic boundary conditions. The well-posedness of Eq. (1.3) was proved in [21]. For numerical study, Hooper and Grimshaw [13] identified three classes of traveling wave solutions of the Kuramoto–Sivashinsky equation and labeled them as regular shocks, oscillatory shocks and solitary waves. In [26], Yang further studied the traveling-wave solutions of the Kuramoto– Sivashinsky equation. Coupled nonlinear equations in which a KdV structure is embedded occur naturally in shallow water wave problems [22]. For a = 6, b = 2, c = 1, d = 2, Eq. (1.2) represents the ItoÕs equation, which describes the interaction process of two internal long waves and has infinitely many conserved quantities [15]. It is noted that this choice of parameters is not unique [1]. In [17], it is shown that this coupled system can be a member of a bi-Hamiltonian integrable hierarchy. The discontinuous Galerkin (DG) method we discuss in this paper is a class of finite element methods using completely discontinuous piecewise-polynomial space for the numerical solution and the test functions in the spatial variables, usually coupled with explicit and nonlinearly stable high order Runge–Kutta time discretization [20]. It was first developed for hyperbolic conservation laws containing first derivatives by Cockburn et al. in a series of papers [8,7,5,9]. For a detailed description of the method as well as its implementation and applications, we refer the readers to the lecture notes [4], the survey paper [6], other papers in that Springer volume, and the review paper [11]. These discontinuous Galerkin methods were generalized to solve a convection diffusion equation (containing second derivatives) by Cockburn and Shu [10]. Their work was motivated by the successful numerical experiments of Bassi and Rebay [2] for the compressible Navier–Stokes equations. Later, Yan and Shu developed a local discontinuous Galerkin method for a general KdV type equation containing third derivatives in [24], and they generalized the local discontinuous Galerkin method to PDEs with fourth and fifth spatial derivatives in [25]. Levy et al. [18] developed local discontinuous Galerkin methods for solving nonlinear dispersive equations that have compactly supported traveling wave solutions, the so-called
3432
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
‘‘compactons’’. Recently, Xu and Shu [23] further developed the local discontinuous Galerkin method to solve three classes of nonlinear wave equations formulated by the general KdV–Burgers type equations, the general fifth-order KdV type equations and the fully nonlinear K(n, n, n) equations. These discontinuous Galerkin methods have several attractive properties. It can be easily designed for any order of accuracy. In fact, the order of accuracy can be locally determined in each cell, thus allowing for efficient p adaptivity. It can be used on arbitrary triangulations, even those with hanging nodes, thus allowing for efficient h adaptivity. The methods have excellent parallel efficiency. It is extremely local in data communications. The evolution of the solution in each cell needs to communicate only with the immediate neighbors, regardless of the order of accuracy. Finally, it has excellent provable nonlinear stability. One can typically prove a strong L2 stability and a cell entropy inequality for the square entropy, for quite general nonlinear cases, for any orders of accuracy on arbitrary triangulations in any space dimension, without the need for nonlinear limiters. The paper is organized as follows. In Section 2, we present and analyze the local discontinuous Galerkin methods for Eqs. (1.1) and (1.2). In Section 2.1, we present the methods for the Kuramoto–Sivashinsky type equations. We prove a theoretical result of the L2 stability for the nonlinear case. In Section 2.2, we present the local discontinuous Galerkin methods for Ito-type coupled KdV equations and give a theoretical result of the L2 stability. In Section 3 we describe the exponential time differencing method for time discretization, which can achieve high order accuracy and maintain good stability. Section 4 contains numerical results for the nonlinear problems to demonstrate the accuracy and capability of the methods. Concluding remarks are given in Section 5.
2. The local discontinuous Galerkin method 2.1. The Kuramoto–Sivashinsky type equations In this section, we present and analyze a local discontinuous Galerkin method for the following Kuramoto–Sivashinsky type equations ut þ f ðuÞx ðaðuÞux Þx þ ðr0 ðuÞgðrðuÞx Þx Þx þ ðsðux Þuxx Þxx ¼ 0;
ð2:1Þ
with an initial condition uðx; 0Þ ¼ u0 ðxÞ
ð2:2Þ
and periodic boundary conditions uðx; tÞ ¼ uðx þ L; tÞ;
ð2:3Þ
where L is the period in the x direction. Here f(u), a(u) P 0, r(u), s(p) P 0 and g(p) are arbitrary (smooth) nonlinear functions. Notice that the assumption of periodic boundary conditions is for simplicity only and is not essential: the method can be easily designed for non-periodic boundary conditions and such non-periodic boundary conditions are used in the numerical experiments in next section. We denote the mesh by I j ¼ ½xj12 ; xjþ12 , for j = 1, . . . , N. The center of the cell is xj ¼ ðxj12 þ xjþ12 Þ=2 and Dxj ¼ xjþ12 xj12 . We denote by uþ and u the values of u at xjþ12 , from the right cell, Ij + 1, and from the jþ1 jþ1 2
2
left cell, Ij, respectively. We define the piecewise-polynomial space VDx as the space of polynomials of the degree k in each cell Ij, i.e. V Dx ¼ fv :
v 2 P k ðI j Þ for x 2 I j ; j ¼ 1; . . . ; N g:
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3433
To define the local discontinuous Galerkin method, we rewrite Eq. (2.1) as a first-order system: ut þ f ðuÞx ðbðuÞvÞx þ ðr0 ðuÞpÞx þ wx ¼ 0; ð2:4Þ v BðuÞx ¼ 0; p gðqÞx ¼ 0; q rðuÞx ¼ 0; w ðlðzÞhÞx ¼ 0; h QðzÞx ¼ 0; z ux ¼ 0; pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi where bðuÞ ¼ aðuÞ, B(u) = òub(s)ds, lðzÞ ¼ sðzÞ, Q(z) = òzl(s)ds. Now we can use the local discontinuous Galerkin method to Eq. (2.4), resulting in the following scheme: find u, v, p, q, w, h, z 2 VDx such that, for all test functions q, /, u, w, g, n, f 2 VDx, Z Z ut q dx ðf ðuÞ bðuÞv þ r0 ðuÞp þ wÞqx dx Ij
Ij
^ ^v þ r^0 ^p þ w ^ Þjþ1 q ^ Þj1 qþ þ ðf^ ^ b^v þ r^0 ^ pþw 1 ðf b^ 1 ¼ 0; 2 jþ2 2 j2 Z Z ^ jþ1 / 1 þ B ^ j1 /þ 1 ¼ 0; v/ dx þ BðuÞ/x dx B 2 jþ2 2 j2 I I Zj Zj pu dx þ gðqÞux dx g^jþ12 u þ g^j12 uþ ¼ 0; jþ12 j12 Ij Ij Z Z qw dx þ rðuÞwx dx ^rjþ12 w rj12 wþ ¼ 0; jþ12 þ ^ j12 Ij Ij Z Z þ ^^ wg dx þ lðzÞhgx dx ð^l^ hÞjþ1 g 1 þ ðlhÞj1 g 1 ¼ 0; 2 jþ2 2 j2 Ij Ij Z Z ^ jþ1 n 1 þ Q ^ j1 nþ 1 ¼ 0; hn dx þ QðzÞnx dx Q 2 jþ2 2 j2 Ij Ij Z Z zf dx þ ufx dx ^ ujþ12 f uj12 fþ ¼ 0: jþ1 þ ^ j1 Ij
Ij
2
ð2:5Þ
2
The ‘‘hat’’ terms in (2.5) in the cell boundary terms from integration by parts are the so-called ‘‘numerical fluxes’’, which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure stability. It turns out that we can take the simple choices such that f^ ¼ f^ ðu ; uþ Þ; ^r ¼ rðu Þ;
^ ¼ Bðu Þ; ^v ¼ vþ ; g^ ¼ g^ðq ; qþ Þ; B ^ ¼ Qðzþ Þ; ^ ^p ¼ pþ ; Q ^ ¼ wþ ; h ¼ h ; w
Bðuþ Þ Bðu Þ ^ ; b¼ uþ u
rðuþ Þ rðu Þ r^0 ¼ ; uþ u
^u ¼ u ;
ð2:6Þ
þ ^l ¼ Qðz Þ Qðz Þ ; þ z z
where we have omitted the half-integer indices j þ 12 as all quantities in (2.6) are computed at the same points (i.e. the interfaces between the cells). Here f^ ðu ; uþ Þ and ^ gðq ; qþ Þ are monotone fluxes, i.e. Lips^ chitz continuous in both arguments, consistent (i.e. f ðu; uÞ ¼ f ðuÞ), non-decreasing in the first argument and non-increasing in the second argument. Examples of monotone fluxes which are suitable for discontinuous Galerkin methods can be found in, e.g., [8]. We could for example use the simple Lax–Friedrichs flux 1 f^ ðu ; uþ Þ ¼ ðf ðu Þ þ f ðuþ Þ aðuþ u ÞÞ; 2
a ¼ max jf 0 ðuÞj;
where the maximum is taken over a relevant range of u. This flux is used in the numerical experiments in next section. The algorithm is now well defined. ^ and ^v from We remark that the choice for the fluxes (2.6) is not unique. In fact the crucial part is taking B ^ ^ ^ and ^u opposite sides, taking ^r and ^ p from opposite sides, taking Q and h from opposite sides and taking w from opposite sides.
3434
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
With such a choice of fluxes we can get the theoretical results of the L2 stability. ^ jþ1 such that the solution to Proposition 2.1 (cell entropy inequality). There exist numerical entropy fluxes U 2 the scheme (2.5), (2.6) satisfies Z Z 1 d 2 ^ ^ u ðx; tÞdx þ Ujþ12 Uj12 þ ðh2 þ v2 Þdx 6 0: 2 dt I j Ij Proof. Since (2.5) holds for any test functions in VDx, we can choose q ¼ u;
/ ¼ v;
u ¼ q;
w ¼ p;
g ¼ z;
n ¼ h;
f ¼ w:
With these choices of test functions and summing up the seven equations in (2.5) we obtain Z Z ^ jþ1 U ^ j1 þ Hj1 þ ðh2 þ v2 Þdx ¼ 0; ut u dx þ U 2 2 2 Ij
Ij
where the numerical entropy flux is given by ^ ¼ F ðu Þ þ Gðq Þ þ ðf^ ^ U b^v þ r^0 ^ p þ wþ Þu g^q ^l^hz ðrðu Þ ^rÞp ^ ^ þ ðQðz Þ QÞh þ ðBðu Þ BÞv and the extra term H is given by ^ H ¼ ½F ðuÞ f^ ½u ½GðqÞ þ g^½q þ ½rðuÞp ^r½p r^0 ^p½u ½BðuÞv þ ^b^v½u þ B½v ^ ½QðzÞh þ ^l^ h½z þ Q½h: Here F ðuÞ ¼
Z
u
f ðuÞdu;
GðqÞ ¼
Z
q
gðqÞdq
and [v] = v+ v denotes the jump of v. With the definition (2.6) of the numerical fluxes and with simple algebraic manipulations, we easily obtain ^ ¼ 0; ½BðuÞv þ ^ b^v½u þ B½v
p½u ¼ 0; ½rðuÞp ^r½p r^0 ^
^ ¼0 ½QðzÞh þ ^l^h½z þ Q½h
and hence H ¼ ½F ðuÞ f^ ½u ½GðqÞ þ g^½q ¼
Z
uþ
ðf ðsÞ f^ ðu ; uþ ÞÞds
u
Z
qþ
ðgðsÞ g^ðq ; qþ ÞÞds P 0; q
where the last inequality follows from the monotonicity of the fluxes f^ and ^ g. Hence Z Z ^ jþ1 U ^ j1 þ ðh2 þ v2 Þdx 6 0; ut u dx þ U 2 2 Ij
Ij
which is the cell entropy inequality.
h
Summing up the cell entropy inequalities, we obtain Corollary 2.2 (L2 stability). The solution to the scheme (2.5), (2.6) satisfies the L2 stability Z L d u2 ðx; tÞdx 6 0: dt 0
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3435
2.2. Ito-type coupled equations In this section, we present and analyze a local discontinuous Galerkin method for the following Ito-type coupled nonlinear problem ut þ auux þ bvvx þ cuxxx ¼ 0;
ð2:7Þ
vt þ bðuvÞx ¼ 0; with an initial condition uðx; 0Þ ¼ u0 ðxÞ;
vðx; 0Þ ¼ v0 ðxÞ
ð2:8Þ
and periodic boundary conditions uðx; tÞ ¼ uðx þ L; tÞ;
vðx; tÞ ¼ vðx þ L; tÞ;
ð2:9Þ
where L is the period in the x direction. Notice that the assumption of periodic boundary conditions is for simplicity only and is not essential: the method can be easily designed for non-periodic boundary conditions. To define the local discontinuous Galerkin method, we rewrite Eq. (2.7) as a first-order system ut þ f ðu; vÞx þ px ¼ 0; p cqx ¼ 0;
q ux ¼ 0;
ð2:10Þ
vt þ gðu; vÞx ¼ 0; where f ðu; vÞ ¼ a2 u2 þ b2 v2 and g(u, v) = buv. Now we can define the local discontinuous Galerkin method to Eq. (2.10), resulting in the following scheme: find u, p, q, v 2 VDx such that, for all test functions q, /, u, w 2 VDx, Z Z ut q dx ðf ðu; vÞ þ pÞqx dx þ ðf^ þ ^ pÞjþ1 q ðf^ þ ^pÞj1 qþ ¼ 0; jþ1 j1 Z Z Z
Ij
p/ dx þ c Ij
qu dx þ Ij
Ij
Z Z
Ij
2
2
2
q/x dx c^ qjþ12 / qj12 /þ ¼ 0; jþ1 þ c^ j1
uux dx Ij
vt w dx Ij
2
Ij
Z
2
2
ð2:11Þ ^ ujþ12 u jþ12
þ
^ uj12 uþ j12
¼ 0;
^j1 wþ gðu; vÞwx dx þ g^jþ12 w 1 ¼ 0: jþ1 g 2 j 2
2
The ‘‘hat’’ terms in (2.11) are the numerical fluxes. It turns out that we can again take the simple choices such that 1 1 f^ ¼ ðf þ þ f kðuþ u ÞÞ; g^ ¼ ðgþ þ g kðvþ v ÞÞ; 2 2 q if c < 0; þ ^ p¼p ; ^ u¼u ; ^ q¼ qþ if c > 0; where k¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 max ða þ bÞu ða bÞ u2 þ 4b2 v2 ; ða þ bÞu þ ða bÞ u2 þ 4b2 v2 2
ð2:12Þ
3436
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
is the maximum eigenvalue of the Jacobian of the flux (f, g)T. The algorithm is now well defined. Notice that f^ and g^ are the standard Lax–Friedrichs fluxes in this system case. We remark that the choice for the fluxes (2.12) is not unique. In fact the crucial part is taking ^p and ^u from opposite sides. With such a choice of fluxes we can get the theoretical results of the L2 stability. Proposition 2.3 (L2 stability). The solution to the scheme (2.11), (2.12) satisfies the L2 stability Z L d ðu2 þ v2 Þdx 6 0: dt 0 Proof. We prove the case for c < 0 and take c = 1, the case for c > 0 can be proved in a similar way. Since (2.11) holds for any test functions in VDx, we can choose q ¼ u;
/ ¼ q;
u ¼ p;
w ¼ v:
With these choices of test functions and summing up the four equations in (2.11) we obtain Z ^ jþ1 H ^ j1 þ Hj1 ¼ 0; ðut u þ vt vÞdx þ H 2 2 2
ð2:13Þ
Ij
where the numerical entropy flux is given by ^ ¼ pþ u þ 1 ðq2 Þ a ðu3 Þ b ðuv2 Þ þ f^ u þ g^v H 2 6 2 and the extra term H is given by 1 2 1 2 2 H ¼ ½q þ ð½u ða½u þ 6kÞ þ 3½v ðb½u þ 2kÞÞ: 2 12 Here [v] = v+ v denotes the jump of v. Note that 1 k P ðja þ bj þ ja bjÞjuj ¼ maxfjaj; jbjgjuj: 2 Hence we obtain ½u2 ða½u þ 6kÞ þ 3½v2 ðb½u þ 2kÞ P 0: This implies H P 0. Summing up (2.13) over j and taking into account the periodic boundary condition, we obtain the desired L2 stability. h
3. Exponential time differencing methods The Runge–Kutta discontinuous Galerkin method is originally designed with explicit and nonlinearly stable high order Runge–Kutta time discretizations [20]. This is adequate for hyperbolic problems with only first derivatives, but the explicit stability limit for the time step decreases dramatically for problems containing higher order spatial derivatives. It is thus important to consider other more efficient time discretization methods. In this section we describe the exponential time differencing fourth-order Runge–Kutta (ETD4RK) method which was developed by Cox and Matthews in [12], among general formulas for ETD-schemes to arbitrary order (also called ELP-schemes [3]). In a recent paper Kassam and Trefethen [16] compared various fourth-order methods for solving equations of the form (3.2), and conclude that the best by a clear margin was a modification by using the complex integral to evaluate the coefficients
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3437
in the update formula. In essence, for nonlinear time dependent equations, the ETD schemes provide a systematic coupling of the explicit treatment of nonlinearities and the implicit and possibly exact integration of the stiff linear parts of the equations, while achieving high order accuracy and maintaining good stability. We consider partial differential equations of the form ut ¼ Lu þ Nðu; tÞ;
ð3:1Þ
where L is a linear term and N is nonlinear. Many interesting equations are of this form, where typically L represents the stiff part of the equation. When discretizing in space we obtain a system of ODEs of the form u_ ¼ Lu þ N ðu; tÞ:
ð3:2Þ
The exponential time differencing (ETD) methods can be described in the context of solving (3.2). Integrating the equation over a single time step from t = tn to tn+1 = tn + Dt, we obtain Z Dt uðtnþ1 Þ ¼ eLDt uðtn Þ þ eLDt eLs N ðuðtn þ sÞ; tn þ sÞds: ð3:3Þ 0
Eq. (3.3) is exact, and the various ETD schemes come from the approximations to the integral [12]. Cox and Matthews derived a set of ETD methods based on Runge–Kutta time-stepping, which they called ETDRK schemes. For completeness, we give below the formulae for the fourth-order scheme (ETD4RK): an ¼ eLDt=2 un þ L1 ðeLDt=2 IÞN ðun ; tn Þ; bn ¼ eLDt=2 un þ L1 ðeLDt=2 IÞN ðan ; tn þ Dt=2Þ; cn ¼ eLDt=2 an þ L1 ðeLDt=2 IÞð2N ðbn ; tn þ Dt=2Þ N ðun ; tn ÞÞ;
ð3:4Þ
2
unþ1 ¼ eLDt un þ Dt2 L3 f½4 LDt þ eLDt ð4 3LDt þ ðLDtÞ ÞN ðun ; tn Þ þ 2½2 þ LDt þ eLDt ð2 þ LDtÞðN ðan ; tn þ Dt=2Þ þ N ðbn ; tn þ Dt=2ÞÞ þ ½4 3LDt ðLDtÞ2 þ eLDt ð4 LDtÞN ðcn ; tn þ DtÞg: To overcome the vulnerability of the error cancelations in the high order ETD and ETDRK schemes, and to generalize the ETD schemes to non-diagonal problems, modified ETD schemes are proposed in [16] by using the complex contour integrals Z 1 f ðLÞ ¼ f ðzÞðzI LÞ1 dz ð3:5Þ 2p C on a suitable contour C to evaluate the coefficients in the update formula for ETD4RK. The explicit exponential time differencing method for time discretization can achieve high order accuracy and maintain good stability while avoiding the severe explicit stability time step limit of the traditional Runge–Kutta discontinuous Galerkin methods which use explicit and nonlinearly stable high order Runge–Kutta time discretization due to the presence of the high order derivative terms.
4. Numerical results In this section we provide numerical examples to illustrate the accuracy and capability of the method. Time discretization is by the third order explicit Runge–Kutta method in [20] or by the exponential time differencing fourth-order Runge–Kutta method in [16]. The time step of the exponential time differencing fourth-order Runge–Kutta methods is Dt = O(Dx), and the time step of the third-order Runge–Kutta
3438
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Table 4.1 Accuracy test for the Kuramoto–Sivashinsky equation (4.1) with the exact solution (4.2) N
L2 error
Order
L1 error
Order
0
p
40 80 160 320
2.38E01 1.58E01 9.41E02 5.24E02
– 0.60 0.74 0.84
1.37 8.81E01 5.21E01 2.91E01
– 0.64 0.76 0.84
p1
40 80 160 320
6.08E02 1.49E02 3.78E03 9.57E04
– 2.02 1.98 1.98
6.64E01 1.82E01 4.64E02 1.19E02
– 1.87 1.97 1.96
p2
40 80 160 320
9.15E03 1.06E03 1.32E04 1.65E05
3.45 3.11 3.00 3.00
1.49E01 1.73E02 2.43E03 3.05E04
2.58 3.11 2.83 2.99
p3
20 40 80 160
3.45E02 9.42E04 7.87E05 5.01E06
– 5.20 3.58 3.97
2.62E01 1.12E02 1.30E03 8.78E05
– 4.55 3.11 3.89
c = 0.1, k ¼ 12
qffiffiffiffi 11 19, x0 = 10. Exact boundary condition. Uniform meshes with N cells at time t = 1.
Fig. 4.1. The shock profile wave propagation of the Kuramoto–Sivashinsky equation (4.1) with initial condition (4.4), c = 5, k ¼ 12 x0 = 12. Exact boundary condition in [30, 30], P2 elements with 160 cells.
qffiffiffiffi 11 , 19
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3439
Table 4.2 Accuracy test for the Kuramoto–Sivashinsky equation (4.5) with the exact solution (4.6) N
L2 error
Order
L1 error
Order
0
p
40 80 160 320
1.33E02 6.69E03 3.35E03 1.68E03
– 0.98 1.00 1.00
7.56E02 3.87E02 1.95E02 9.80E03
– 0.92 0.98 1.00
p1
40 80 160 320
1.29E03 3.26E04 8.17E05 2.05E05
– 1.99 2.00 2.00
1.04E02 2.59E03 6.54E04 1.64E04
– 2.01 1.98 2.00
p2
40 80 160 320
6.48E05 8.04E06 1.00E06 1.25E07
– 3.01 3.00 3.00
7.93E04 1.06E04 1.34E05 1.67E06
– 2.91 2.99 3.00
p3
40 80 160 320
5.92E06 3.83E07 2.42E08 1.64E09
– 3.95 3.99 3.88
6.07E05 4.00E06 2.56E07 1.61E08
– 3.92 3.97 3.99
ffi, x0 = 10. Exact boundary condition. Uniform meshes with N cells at time t = 1. c = 0.2, k ¼ 2p1ffiffiffi 19
ffi, x0 = 25. Exact Fig. 4.2. The shock profile wave propagation of the KS equation (4.5) with initial condition (4.8), c = 5, k ¼ 2p1ffiffiffi 19 boundary condition in [50, 50], P2 elements with 200 cells.
3440
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Table 4.3 Accuracy test for the Kuramoto–Sivashinsky equation (4.9) with the exact solution (4.10) N
L2 error
Order
L1 error
Order
1
p
80 160 320 640
1.65E01 2.33E02 4.68E03 1.28E03
– 2.82 2.32 1.87
1.01 2.01E01 6.49E02 1.83E02
– 2.32 1.63 1.82
p2
80 160 320 640
5.98E03 7.43E04 9.30E05 1.16E05
– 3.01 3.00 3.00
1.07E01 1.44E02 1.84E03 2.31E04
– 2.89 2.97 3.00
p3
40 80 160 320
1.10E02 6.84E04 4.36E05 2.75E06
– 4.00 3.97 3.99
9.31E02 1.18E02 7.36E04 5.01E05
– 2.98 4.01 3.88
c = 6, r = 4, k ¼ 12, x0 = 10. Periodic boundary condition. Uniform meshes with N cells at time t = 1.
Fig. 4.3. The solitary wave propagation of the Kuramoto–Sivashinsky equation (4.9) with initial condition (4.11), c = 6, r = 4, k ¼ 12, x0 = 10. Periodic boundary condition in [30, 30], P2 elements with 250 cells.
method is Dt = O(Dx4) for the KS equation and Dt = O(Dx3) for the ItoÕs equation. We have verified with the aid of successive mesh refinements, that in all cases, the results shown are numerically convergent.
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3441
Example 4.1. We show an accuracy test for the Kuramoto–Sivashinsky equation ut þ uux þ uxx þ uxxxx ¼ 0 with the exact solution rffiffiffiffiffi 15 11 uðx; tÞ ¼ c þ 9 tanhðkðx ct x0 ÞÞ þ 11tanh3 ðkðx ct x0 ÞÞ 19 19
ð4:1Þ
ð4:2Þ
and the boundary condition uð30; tÞ ¼ g1 ðtÞ;
uð30; tÞ ¼ g2 ðtÞ;
ð4:3Þ
where gi(t) corresponds to the data from the exact solution. Notice that the local discontinuous Galerkin method allows for an easy implementation of such boundary conditions. The L2 and L1 errors and the numerical orders of accuracy are contained in Table 4.1. We can see that the method with Pk elements gives a uniform (k + 1)th order of accuracy in both norms.
Fig. 4.4. Chaotic solution of the Kuramoto–Sivashinsky equation (4.12) with initial condition (4.13). Periodic boundary condition in [0, 32p], P2 elements with two different meshes using N = 300 and N = 600 uniform cells.
3442
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
We also present the regular shock profile wave propagation for the Kuramoto–Sivashinsky equation (4.1) in Fig. 4.1 with the initial condition rffiffiffiffiffi 15 11 uðx; 0Þ ¼ c þ ð4:4Þ 9 tanhðkðx x0 ÞÞ þ 11tanh3 ðkðx x0 ÞÞ : 19 19 We can see clearly that the moving shock profile is resolved very well. Example 4.2. We show an accuracy test for the Kuramoto–Sivashinsky equation ut þ uux uxx þ uxxxx ¼ 0
ð4:5Þ
with the exact solution uðx; tÞ ¼ c þ
15 pffiffiffiffiffi 3 tanhðkðx ct x0 ÞÞ þ tanh3 ðkðx ct x0 ÞÞ 19 19
ð4:6Þ
and the boundary condition uð50; tÞ ¼ h1 ðtÞ;
uð50; tÞ ¼ h2 ðtÞ;
ð4:7Þ
Fig. 4.5. The chaotic solution of the Kuramoto–Sivashinsky equation (4.12) with the Gaussian initial condition (4.14). Periodic boundary condition in [16, 16], P2 elements with two different meshes using N = 160 and N = 320 uniform cells.
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3443
where hi(t) corresponds to the data from the exact solution. The result is contained in Table 4.2. The scheme clearly demonstrates an order of accuracy of k + 1 for Pk elements in both L2 and L1 norms for this problem with non-periodic boundary conditions. Next we show the shock profile wave propagation for the Kuramoto–Sivashinsky equation (4.5) in Fig. 4.2 with the initial condition 15 uðx; 0Þ ¼ c þ pffiffiffiffiffi 3 tanhðkðx x0 ÞÞ þ tanh3 ðkðx x0 ÞÞ : ð4:8Þ 19 19 We can again observe that the moving shock profile is resolved very well. Example 4.3. We show an accuracy test for the Kuramoto–Sivashinsky equation ut þ uux þ uxx þ ruxxx þ uxxxx ¼ 0
ð4:9Þ
with the exact solution uðx; tÞ ¼ c þ 9 15 tanhðkðx ct x0 ÞÞ þ tanh2 ðkðx ct x0 ÞÞ tanh3 ðkðx ct x0 ÞÞ :
ð4:10Þ
Table 4.3 give the errors of numerical solution at t = 1 using the periodic boundary condition. We can see from the table that the orders of accuracy are comparable to that for the non-periodic case. Next we show the solitary wave propagation for the Kuramoto–Sivashinsky equation (4.9) in Fig. 4.3 with the initial condition
Fig. 4.6. Numerical results of u for the ItoÕs equation (4.15) with the initial condition (4.16). Periodic boundary condition in [0, 2p], P2 elements with 80 cells.
3444
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Fig. 4.7. Numerical results of v for the ItoÕs equation (4.15) with the initial condition (4.16). Periodic boundary condition in [0, 2p], P2 elements with 80 cells.
uðx; 0Þ ¼ c þ 9 15 tanhðkðx x0 ÞÞ þ tanh2 ðkðx x0 ÞÞ tanh3 ðkðx x0 ÞÞ :
ð4:11Þ
This moving solitary wave is clearly well resolved. Example 4.4. The Kuramoto–Sivashinsky equation ut þ uux þ uxx þ uxxxx ¼ 0
ð4:12Þ
has received a lot of interest because it is one of the simplest partial differential equations which is capable of exhibiting chaotic behavior. The chaotic behavior typically occurs when (4.12) is integrated over finite xdomain with periodic boundary conditions. First, we take the initial condition x x uðx; 0Þ ¼ cos 1 þ sin ð4:13Þ 16 16 and present the results in Fig. 4.4 for two different meshes with N = 300 and N = 600 uniform cells respectively. We can see clearly that the results are numerically convergent for the quite chaotic behavior. Next we take the Gaussian initial condition 2
uðx; 0Þ ¼ ex
ð4:14Þ
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3445
Fig. 4.8. Numerical results of u for the ItoÕs equation (4.15) with the initial condition (4.17). Periodic boundary condition in [15, 15], P2 elements with 160 cells.
and present the results in Fig. 4.5 for two different meshes with N = 160 and N = 320 uniform cells respectively. We can again see clearly that the results are numerically convergent for this case. Example 4.5. In this example, we show the numerical results for the ItoÕs equation ut ð3u2 þ v2 Þx uxxx ¼ 0; vt 2ðuvÞx ¼ 0;
ð4:15Þ
with an initial condition uðx; 0Þ ¼ cosðxÞ;
vðx; 0Þ ¼ cosðxÞ:
ð4:16Þ
We choose an Gaussian initial condition uðx; 0Þ ¼ expðx2 Þ;
vðx; 0Þ ¼ expðx2 Þ:
ð4:17Þ
Notice that there is a dispersive term in the first equation for u in (4.15), while there is no such dispersive term in the second equation for v. From Figs. 4.6–4.8 and 4.9, we observe the result for u behaving like dispersive wave solutions and the result for v behaving like shock wave solutions.
3446
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
Fig. 4.9. Numerical results of v for the ItoÕs equation (4.15) with the initial condition (4.17). Periodic boundary condition in [15, 15], P2 elements with 160 cells.
5. Concluding remarks We have developed local discontinuous Galerkin methods to solve the Kuramoto–Sivashinsky type equations and the Ito-type coupled KdV equations and have proven the L2 stability of these methods. The exponential time differencing fourth-order Runge–Kutta methods have been used for treating the nonlinearity and stiffness of the systems in order to improve computational efficiency. Numerical examples are shown to illustrate the accuracy and capability of the methods. Although not addressed in this paper, these methods are flexible for general geometry, unstructured meshes and h-p adaptivity, and have excellent parallel efficiency. They should provide a useful class of numerical tools for solving nonlinear wave equations. An important issue not addressed in this paper is L2 a priori error estimates. From the cell entropy inequality and approximation results, one can derive L2 a priori error estimates following the lines in [10,24], for linear PDEs. The proof is not completely straightforward and sometimes one would lose half an order or even one order in accuracy, because of a lack of control for certain jump terms at cell boundaries. It is more challenging to perform L2 a priori error estimates for nonlinear PDEs. Such error estimates are left for future work. References [1] M. Antonowicz, A.P. Fordy, Coupled KdV equations with multi-Hamiltonian structures, Physica D 28 (1987) 345–357. [2] F. Bassi, S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations, J. Comput. Phys. 131 (1997) 267–279.
Y. Xu, C.-W. Shu / Comput. Methods Appl. Mech. Engrg. 195 (2006) 3430–3447
3447
[3] G. Beylkin, J.M. Keiser, L. Vozovoi, A new class of time discretization schemes for the solution of nonlinear PDEs, J. Comput. Phys. 147 (1998) 362–387. [4] B. Cockburn, Discontinuous Galerkin methods for methods for convection-dominated problems, in: T.J. Barth, H. Deconinck (Eds.), High-Order Methods for Computational Physics, Lecture Notes in Computational Science and Engineering, vol. 9, Springer, 1999, pp. 69–224. [5] B. Cockburn, S. Hou, C.-W. Shu, The Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case, Math. Comput. 54 (1990) 545–581. [6] B. Cockburn, G. Karniadakis, C.-W. Shu, The development of discontinuous Galerkin methods, in: B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Lecture Notes in Computational Science and Engineering, vol. 11, Springer, 2000, pp. 3–50, Part I: Overview. [7] B. Cockburn, S.-Y. Lin, C.-W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one dimensional systems, J. Comput. Phys. 84 (1989) 90–113. [8] B. Cockburn, C.-W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Math. Comput. 52 (1989) 411–435. [9] B. Cockburn, C.-W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, J. Comput. Phys. 141 (1998) 199–224. [10] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463. [11] B. Cockburn, C.-W. Shu, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (2001) 173–261. [12] S.M. Cox, P.C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys. 176 (2002) 430–455. [13] A.P. Hooper, R. Grimshaw, Travelling wave solutions of the Kuramoto–Sivashinsky equation, Wave Motion 10 (1988) 405–420. [14] J.M. Hyman, B. Nicolaenko, The Kuramoto–Sivashinsky equation: a bridge between PDEs and dynamical systems, Physica D 18 (1986) 113–126. [15] M. Ito, Symmetries and conservation laws of a coupled nonlinear wave equation, Phys. Lett. A 91 (1982) 335–338. [16] A.K. Kassam, L.N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput. 26 (2005) 1214–1233. [17] B.A. Kupershmidt, A coupled Korteweg–de Vries equation with dispersion, J. Phys. A: Math. Gen. 18 (1985) L571–L573. [18] D. Levy, C.-W. Shu, J. Yan, Local discontinuous Galerkin methods for nonlinear dispersive equations, J. Comput. Phys. 196 (2004) 751–772. [19] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation, Physica D 19 (1986) 89–111. [20] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [21] E. Tadmor, The well-posedness of the Kuramoto–Sivashinsky equation, SIAM J. Math. Anal. 17 (1986) 884–893. [22] G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York, 1974. [23] Y. Xu, C.-W. Shu, Local discontinuous Galerkin methods for three classes of nonlinear wave equations, J. Comput. Math. 22 (2004) 250–274. [24] J. Yan, C.-W. Shu, A local discontinuous Galerkin method for KdV type equations, SIAM J. Numer. Anal. 40 (2002) 769–791. [25] J. Yan, C.-W. Shu, Local discontinuous Galerkin methods for partial differential equations with higher order derivatives, J. Sci. Comput. 17 (2002) 27–47. [26] T.S. Yang, On traveling-wave solutions of the Kuramoto–Sivashinsky equation, Physica D 110 (1997) 25–42.