- Email: [email protected]
Approximations of the KdV equation by least squares finite elements
Computer Methods in Applied Mechanics and Engineering North-Holland
93 (1991) l-11
Approximations of the KdV equation by least squares finite elements G.F. Carey and Y. Shen Department of Aerospace
Engineering
and Engineering Mechanics, Austin, TX 78712, USA
The University of Texas at Austin,
Received 24 July 1988 Revised manuscript received July 1989
The problem of approximating solutions to the Korteweg-de Vries (KdV) equation is investigated using a least squares finite element method. The third order KdV equation is recast as a first-order system and a least-squares finite element approach is introduced for the semidiscrete time-differenced form of the resulting equations. Of particular interest are the approximation properties for solitary wave solutions (solitons). We examine the amplitude and phase error for a representative test problem as well as other examples including the passage of one soliton through another.
1. Introduction The solitary wave (soliton) consists of a single ‘hump’ of constant height and speed, and was first observed by Scott Russell in 1844 [l] for waves in channels. Shortly thereafter, Stokes [2] developed approximate expressions for nonlinear periodic waves in deep water and Boussinesq [3] and Rayleigh [4] obtained approximate expressions for the solitary wave. The solitary wave phenomenon is most easily characterized as a special solution of a nonlinear wave equation studied by Korteweg and de Vries [5]. Using perturbation theory for a low long wave in a shallow, narrow channel, they demonstrated solitary wave solutions and also periodic solutions. There is now a very extensive literature on the KdV equation and the mathematical properties of its solutions (see e.g. [6-91). For example, it can be shown that an arbitrary initial disturbance breaks up eventually to a set of solitons. Similarly, despite the nonlinearity of the equation, two or more solitons can pass through one another and emerge unaltered in shape. There have been few previous studies of finite element approximation of solitary waves. Schoombie [lo] employed a Galerkin finite element method and used Hermite cubic splines to accommodate the third derivative in the KdV equation. Argyris and Haase [ll] used cubic Schoenberg splines in their Galerkin finite element formulation and also examined a PetrovGalerkin method using piecewise linear test functions with the Hermite cubic trial basis. Other studies were made by Winther [12], Mitchell and Schoombie [13], and Sanz-Serna and Christie
P41. Least-squares finite element methods (see e.g. [15] have been the subject of recent exploratory (see e.g. [15]) investigations for first order systems ([16, 171) and have proven 00457825/91/$03.50
0
1991 Elsevier Science Publishers
B.V. All rights reserved
G.F. Carey, Y. Shen, Approximations
2
of de KdV equation
effective in the analysis of convection-diffusion, Burgers’ equation and compressible Euler equation problems. The method introduces a numerical dissipation that is similar to that for the Petrov-Galerkin and Taylor-Galerkin schemes for hyperbolic problems of convectiondominated transport processes. Our interest in the present problem stems from two main factors: first, the problem is of third order and has not previously been studied using the least-squares approach; secondly, the soliton solutions are special in that they are solutions to a nonlinear problem that have constant amplitude and velocity, and we wish to assess the effect of amplitude and phase errors introduced by the numerical method.
2. KdV equation The assumption that the wave amplitude is small compared with depth and that channel depth is much smaller than length leads to the KdV equation for long waves, which in normalized form becomes u, +
uu,
+
FU,,,
=
0
(1)
)
where F models the dispersion. The exact solution of (1) for a single soliton of amplitude constant velocity c is
3c moving to the right with
u,(x, t) = 3c sech*[k(x - ct) + x0] ,
(2)
with k = 1 /~(c/E)~“, x,, = k(2 - 1) + d where parameter d provides a reference for the initial position of the soliton. Solitons can pass through one another and emerge unchanged. For 2 or 3 solitons i with speed c,, initially centered about x0( the exact solution is [18]
u(x, t> = B(ln Q
,
(3)
where for two solitons F = 1 +
ez] +
ez2 + A ezl+‘2 ,
2, =
A = [(k, - k,)/(k, + k,)]* >
k; =
q&(x -
c,f> ++I
{(cJE)~‘*,
x0,
=
f d, ,
and for three F = det(l;l,) , 2; =
Fii = 6,. + 2(kjki)“2 e”/(k, + k,) ,
2[ k;(x - c;t) + X”,J )
x0; =
k(3c; - 1) + d; ,
with dj again a reference parameter for the respective reference value for the amplitude of the solitons.
soliton positions,
and B providing
a
G.F.
Carey, Y. Shen, Approximations
3
of de KdV equation
3. First-order system and least-squares FEM The KdV equation 8)t,, 06061,
(1) may be differenced
u (n+l) - $) At
with respect to time to give at t = r3tn+, + (1 -
+ 8(uu, + “Uxxx)(~+l)+ (1 - O)(uu, + EU,,,)(~) = 0.
(4)
Next, we recast this as a first-order system by setting U, = u, u, = w and linearize within each timestep At to obtain n+l) _ J4 Ll( At
+ @&$(*+l)
(n+l) - uf+‘))
@
O(W
Transposing conveniently
(n+l) _
+
+ EW;+‘)) + (1 - @(U@+J@)+ EWF)) = 0,
(1 - B)(u’“’ -
@+I)) + (I_ ux
fj)(w’“’
u_y )= 0) - $‘)
= 0.
(5)
terms at time level n and simplifying, we have a first order system which can be expressed in matrix forms as
&)zr+l)
+ &“)z(“+l’ =f’“’ ,
(6)
where zf = [u u w], A’“’ = A($“‘), B’“’ = B(z’“‘) and f’“’ = f(z’“‘). For an admissible solution 6 satisfying the boundary conditions, the corresponding residual vector is r = A(“)gT+l) + &“)[@+l) _f’“’ , and the least-squares +L 2 sR
residual functional
trial
(7)
on domain R is
rtt dx .
(8)
Setting the first variation FZ = 0 in (8), R [A’“’ @+” +
+ &“) @n+l)]f[A(n)J~+l)
The form of (9) corresponds to a Petrov-Galerkin the operators A’“‘, II’“’ as indicated here. Expanding the integrand expression in (9),
+ &n){(n+l) _f’“‘] dx = 0.
(9)
method where the test functions involve
(n+l)‘A(n)~(“)J;+l) + a{r+l)’ (n)’ @) @+I) + ~5(n+l)(B(n)lA(n)5lo+l) ADS I R [MI + 85 Introducing
(n+l)‘~(n)‘B(n)5!n+l,]dx = R (@‘+l)LA(“)’ + S{(n+l)L&#)j+) I
finite element expansions for the components
dx .
(10)
of 5 on a partition of the domain R
G.F.
4
Carey, Y. Shen, Approximations
of de KdV equation
and setting the variations 85 = (&, O,O), (0, $, 0), (O,O, x), where &, ei, yk are the associated component basis functions for U, u, W, yields the finite element system. It follows from the form of (10) that the discrete system is symmetric. For example, if the same basis is chosen for each solution component, then the element contributions K,, F, from (10) to KV= F can be written conveniently in matrix form as
and F, =
Y
(!&‘A’ + @B’)fdx
,
Ulb)
where, for notational convenience, we have omitted the timestep index n and W is a matrix of element basis functions. The element contributions (11) are accumulated to the global system. Essential boundary conditions are imposed to yield a sparse global system of the form KV=F
(14
to be solved at each timestep. Here K and F depend on the solution at the previous time level. The solution is initiated from the prescribed initial vector V” determined by approximation of U, u and w at t = 0 from specified u(x, 0) = g(x). The full form of the system (12) is given in Appendix A. In the preceding formulation we have treated the nonlinear term in (6) by writing uu, as uu and then linearizing to u (“)u(~+‘). This implies that u, u and w appear explicitly in the KdV equation. Alternatively we can use u (n)~(nf ‘) directly. Numerical comparisons for a single soliton indicate that the difference in Computed solutions with these methods is small (approximately one percent for the large soliton c = 0.3 computed in the numerical studies). Since the least-squares variational problem is based on a first-order system, C” Lagrange basis functions may be used. In particular, we may use piecewise-linear basis functions as compared with the cubic splines utilized in the previous studies cited. Thus, the additional complexity of the first-order system (more unknowns) is offset by the simplicity of the basis. The element integrations are carried out using Gaussian quadrature. Since the additional equations u, = u, u., = w are basically constraints, one might anticipate that similar difficulties to those encountered in penalized variational minimization problems might be encountered. These penalty functionals require reduced numerical integration for the method to be stable and the same situation appears to apply for the least-square method. The form of reduced integration rule used here is examined later in the numerical experiments. In other applications we have shown that the least-squares formulation introduces an artificial dissipation similar to that for the Petrov-Galerkin and Taylor-Galerkin methods. To investigate this further, for the present problem it suffices to integrate the leading term in (10) by parts. The corresponding contribution in the associated Euler-Lagrange equation associated with the weak form in (10) is s 6&‘A’Aj,, where the linearized form of the problem has been assumed. Note that the solution vector J is now differentiated twice. This implies that the least-squares method is equivalent to a Petrov-Galerkin form involving derivatives u,, , u,,,
G.F.
Carey, Y. Shen, Approximations
of de KdV equation
5
W XX;that
is, uXX,u,,, and u,,,,. Clearly the terms involving u,, and u,,,, will effect numerical dissipation and, therefore, the amplitude. Those involving u,,, will reinforce the terms EU,,, and influence the phase. Recall that in the soliton problem this term balances the nonlinear inertial term. Hence we anticipate that there will be some influence on the amplitude and phase properties of the propagating solitary waves. We explore these issues in the following numerical studies.
4. Numerical experiments In the calculations that follow, the approximate problem is posed on a finite domain R with essential boundary data u = 0 at the end points. For the soliton problems of interest here and the short time intervals considered, this provides a good approximation to the actual infinite domain problem for which the analytic solution is known.- _ 4.1. Single soliton As the first case, consider a propagating single soliton of the form given in (2) with E = 4.84 x lo-‘, c = 0.1 and d = 0. For the short time interval 0 < t < 4 considered here, it suffices to compute the approximation on a finite domain 0
x,J X pl7
Ax,, % clPan A% Au,iu,
(x,, up)
0
1
2
3
4
0.8 0.8 0 0.3 0.3000 0.0000 0.000
0.9 0.9 0 0.3 0.2990 0.0010 0.003
1.0 1.0 0 0.3 0.2969 0.0031 0.010
1.1 1.1 0 0.3 0.2955 0.0045 0.015
1.2 1.2 0 0.3 0.2940 0.0060 0.020
G.F.
6
1o
Carey, Y. Shen, Approximations
of de KdV equation
u(x,t)
. t
Fig. 1. Propagating
single solition in (x, t) domain for 0 s t G 4 (100 linear elements,
At = 0.01)
very dissipative for large amplitude solitons. Part of the problem can be attributed to the fact that the curvature w is approximated piecewise linearly and this approximation will deteriorate as the amplitude increases. In the present calculation the maximum error in w = u,,~ is 44%. It should be emphasized, however, that physically the KdV equation is derived under the assumption of long waves of low amplitude, so that the previous high-amplitude waves are not realistic for actual shallow water flow problems. For small amplitude waves that are of practical significance here, the present least-squares method is satisfactory. To examine the effect of mesh size let us consider the first case c = 0.1 with 50 elements. At this coarser mesh level there is a small secondary wave (oscillation) that propagates to the left. Similar observations have been made in other studies by Zabusky [193 and Leibovich and Seebass [20] and this effect is attributed to the negative group velocity associated with a linear wave. Pearson [21] has made an interesting analysis of a similar situation. Further numerical experiments indicate that the linearized problem is stable for l/2 d 8 d 1 and, as 0 increases in this range, dissipation increases (as one might expect in view of the well-known dissipative nature of the backward Euler method, 0 = 1.0). We note that dispersion increases with 8 and At. For instance, increasing the timestep from 0.005 to 0.01 increases the dispersion and dissipation errors: In the previous example with c = 0.3, the numerical dispersion at t = 1 is approximately given by Ax, = 0.02 and dissipation by Au, = 0.087 for calculations with 0 = 0.5, At = 0.01. If we use 0 = 1, the error increases to Ax,, = 0.04, Au, = 0.187. For 0 = 0.5, At = 0.005, the error decreases to Ax,) = 0.01, Au, = 0.061. The influence of At and 8 on dissipation can be seen in the presence of the term At du(‘)u~~+‘) identified in (A.4). The effect of linearization within a timestep is negligible for the step sizes At considered here. (In numerical experiments at At = 0.01 with 5 successive approximation corrections for the nonlinearity, the peak value changed at t = 0.2 from 0.8228 to 0.8286.) The quadrature rule in the previous calculations on linear elements was chosen to under-integrate the integral contribution using a l-point Gauss quadrature. If the order of integration is increased to a 2-point Gauss rule, the method fails. This situation is reminiscent of the type of difficulties that arise when penalty functionals are used to imbed constraints in a
G.F.
Carey, Y. Shen, Approximations
7
of de KdV equation
variational finite element model. (For instance, see [22,23].) An analogous situation occurs here since the additional equations U, = u, u, = w can be viewed as constraints which are included in the corresponding least-squares functionals. Higher-degree elements were also examined. The method does not fail with full integration but reduced integration gave results with better phase and dissipation properties. To examine this situation further we computed the eigenvalues of the element matrices for these terms when different quadrature rules were used and determined the rank of the matrix in each case. When linear elements are used the element matrix is of size 6 x 6 and has full rank if reduced integration is not applied. The desired reduced rank (3) is produced by underintegration. Similarly for quadratic and cubic element matrix contributions, from the ‘constraints’ we should get full rank if reduced integration is not used and the desired reduced rank of 6 and 8, respectively, if underintegration is applied. The reduced rank implies that the matrix contributions associated with the constraints are singular and therefore permits a nontrivial solution. However, we find that full integration yields slightly reduced rank, apparently due to the effect of roundoff in finite precision arithmetic.
4.2. Soliton interaction This second part of the numerical study deals with calculations for two and three interacting solitons. The analytic solutions for these respective cases are given in (3). As indicated previously, solitons of different amplitudes travel at different speeds and pass through each other to emerge unchanged. Numerical solutions for two solitons centered at 0.4 and 0.88, and travelling with speeds c, = 0.3 and c2 = 0.1 are shown in Fig. 2 during 0 G t s 4 for calculations with 100 linear elements on 0 < x < 2 using l-point integration, and 8 = 0.5, At = 0.01. In (3) B = 5.84 x 10-j and d, = d, = -5. Both soliton approximations exhibit numerical dissipation and dispersion as noted previously for the corresponding single wave calculations. Since the approximations are not true solitons, their interaction is not ‘pure’. For example, the small soliton (c = 0.1) now has 4% dispersion and 1% dissipation compared with essentially no
Fig. 2. Interaction
of two solitons in (x, t) domain for 0 s t G 4 (100 linear elements,
At = 0.01).
G.F.
Fig. 3. Interaction
Carey, Y. Shen, Approximations
of de KdV equation
of three solitons in (x, t) domain for 0 G t c 8 (200 linear elements,
At = 0.01).
dispersion or dissipation previously; the larger soliton now has 12% dispersion and 18% dissipation which (fortuitously) are slightly better than 20% and 23%, respectively, for the single soli ton propagation. The interaction of three solitons is shown in Fig. 3 where c, = 0.3, c2 = 0.2, c, = 0.1 and d, = d, = d, = -5, B = 5.84 x lo-’ in (3). The numerical results are for 200 linear elements with l-point integration on 0
Fig. 4. Separation
of initial profile c exp { - b(x - 2)‘) into propagating
waves.
9
G.F. Carey, Y. Shen, Approximations of de KdV equation 5.
Concluding remarks
Soliton solutions to the KdV equation provide an interesting class of test problems since the solitons are themselves of physical significance and they provide the opportunity to study numerical dispersion and numerical dissipation for the nonlinear wave problem. Here we formulate and apply a least-squares finite element method based on a representation of the KdV equation as a first-order system. For solitons of small amplitude the method is seen to exhibit little numerical dissipation or dispersion, but at higher amplitudes these influences are significant. Since the dissipation in this method is proportional to the timestep, reducing the timestep will improve the situation but at the cost of significantly increasing computation time. The results suggest that the least-squares method can be applied (conditionally as noted) to this type of nonlinear wave problem with success. An important point is the need to use reduced integration to produce a viable method.
Appendix A Equations (6)-(11) offer a compact matrix formulation. It is instructive also to develop the contributions explicitly. To do so we transpose terms at time level y1 to rewrite (5) for admissible U, u , w as Yl = uCn+‘) r2 = -&l,“+”
+ At &&n)u(n+U + At &wj;n+‘) - fi , + &j”+‘)
-f,
)
T) = -&J(““) x
+ &Jn+‘)
- f3 )
(A4
where f, = u@) - At(1 - ~~)(u(~)LI@) + EWE)), f2 = (1 - B)(u’“’ - ur’) and f, = (1 - 8)(w’“’ u(n) X )* Then I = 4 Jn (Y: + ri + r:) dx and setting 61= 0 implies Jf2 (Y, 6r, + y2 6r, -t r3 6r,) dx = 0 for arbitrary variations kCnt’), 6u(“+‘), ?3w(““). Introducing (A.l) into 61 = 0, collecting terms and simplifying we get
I
R (r, c%@+‘)- r,8 au?+‘)) dx = 0, with 8u(“+‘) arbitrary,
I
D (q At Ou(‘I 8u@+‘) + r,B &I(“+‘) - r,8 6~;““) with I~u(“~‘) arbitrary,
I
c%(“+‘) = 8~~~“) = 0 ;
&@+I) = 8w(“+‘) = 0 ;
dx = 0,
(A.21
R (I-~At C)& 6w$‘+” + ~$9 SW@+‘))dx = 0, with 8w(“+‘) arbitrary,
ILL@+‘)= ?3wCn+‘)= 0.
Introducing finite element expansions u, = C u~c$~,uh = C uk ekI,, wh = c w[‘yrand setting the variations as the respective basis functions &, I,!J~,y, leads to the finite element system
G.F.
10
In the numerical REMARK.
results we use piecewise-linear
Integrating $& J^
Carey, Y. Shetz, A~~roxi~at~otz.~ of de KdV eyzmtiotz
basis functions for uh, ulr and wk.
by parts in (AL?) we obtain
+(Y?_),8)sU(“+“dx=O,
i
dx = 0, R (Y, At8u (‘I)-I-r,B + (Y~)~B)&_J~‘~+‘) (A.4)
i
n[(-~‘)XA~6t:ir,B]Fw’“+“dx=0.
The third equation in (A.4) contains numerical dissipation.
a term AttSu’“‘~~‘~’ which can be interpreted
as a
References [l] J.S. Russell, Report on waves, Rept. 14th Meeting of British Association for the Advancement (John Murray, London, 1844) 311-390. [Z] G.G. Stokes, On the theory of oscillatory waves, Cambridge Trans. 8 (1847) 441-472. [3] J. Boussinesq, Thtorie des ondes et des remous qui se propagent le long d’un canal rectangulaire
of Science
horizontal,
G.F.
[4] [S] [6] [7] [8] (91 [lo] [ll] [12] [13] [14] [15] [16] [17] [18] [19] (201 [21] [22]
[23]
Carey, Y. Shen, Approximations
of de KdV equation
11
en communiquant au liquide contenu darts ce canal des vitesses sensiblement parailles de la surface au fond, J. Math. Pures Appl. 17 (2) (1872) 55-108. Rayleigh, On waves, Phil. Mag. 5 (1) (1876) 257-279. D.J. Korteweg and G. de Vries, On the change of form of long waves advancing in a rectangular canal and on a new type of long stationary wave, Phil. Mag. 39 (1895) 422-433. C.S. Gardner, J.M. Greene, M.D. Kruskal and R.M. Miura. Method for solving the Kortweg-de Vries equation, Phys. Rev. Lett. 1.5 (1967) 1095-1097. R.M. Miura, Korteweg-de Vries equation and generalizations I. A remarkable explicit nonlinear transformation, J. Math. Phys. 9 (1968) 1202-1204. G.B. Witham, Linear and Nonlinear Waves (Wiley, New York, 1984). A.C. Newell, Solitons in mathematics and physics, Philadelphia (198.5) 23-57. S.W. Schoombie, Finite element methods for the Korteweg-de Vries equation I. Galerkin method with Hermite cubits, Rept. NA/43, Department of Mathematics, University of Dundee, 1980. J. Argyris and M. Haase, An engineer’s guide to soliton phenomena: Application of the finite element method, Comput. Methods Appl. Mech. Engrg. 61 (1) (1987) 71-122. R. Winther. A conservative finite element method for the Korteweg-de Vries equation, Math. Comp. (1981). A.R. Mitchell and S.W. Schoombie, Finite element studies of solitons, in: R.W. Lewis, P. Bettes and E. Hinton, eds., Numerical Methods in Coupled Systems (Wiley, New York, 1984) 465-488. J.M. Sanz-Serna and I. Christie, A simple adaptive technique for nonlinear wave problems, J. Comput. Phys. 67 (2) (1986) 348-360. C.A.J. Fletcher, Computational Galerkin Methods (Springer, New York, 1984). G.F. Carey and B.N. Jiang. Least-squares finite elements for first-order hyperbolic problems, Internat. J. Numer. Methods Engrg. 26 (1988) 81-93. G.F. Carey and B.N. Jiang, Least-squares finite element method for compressible Euler equations, 5th Internat. Conf. Numerical Methods in Laminar and Turbulent Flow, Montreal, Canada, July 1987. R. Hirota, Exact solution of the Korteweg-de Vries equation for multiple collisions of solitons, Phys. Rev. Lett. 27 (1972) 1192-1194. N.J. Zabusky, Solitons and bound states of the time-independent Schrodingei equation, Phys. Rev. 168 (1968) 124-128. S. Leibovich and A.R. Seebass, Nonlinear Waves (Cornell Univ. Press, Ithaca, London, 1974). C.E. Pearson, Note on upwind and downwind effects of numerical approximation to wave equations, Comm. Appl. Numer. Methods (1988). R.L. Sani, P.M. Gresho, R.L. Lee, D. Griffiths and M. Engleman, The cause and cure (?) of the spurious pressure generated by certain FEM solutions of the incompressible Navier Stokes equations, Internat. J. Numer. Methods Fluids 1 (1981) 17-43, 171-204. G.F. Carey and R. Krishnan, Penalty approximation of Stokes flow, Comput. Methods Appl. Mech. Engrg. 35 (2) (1982) 169-206.
93 (1991) l-11
Approximations of the KdV equation by least squares finite elements G.F. Carey and Y. Shen Department of Aerospace
Engineering
and Engineering Mechanics, Austin, TX 78712, USA
The University of Texas at Austin,
Received 24 July 1988 Revised manuscript received July 1989
The problem of approximating solutions to the Korteweg-de Vries (KdV) equation is investigated using a least squares finite element method. The third order KdV equation is recast as a first-order system and a least-squares finite element approach is introduced for the semidiscrete time-differenced form of the resulting equations. Of particular interest are the approximation properties for solitary wave solutions (solitons). We examine the amplitude and phase error for a representative test problem as well as other examples including the passage of one soliton through another.
1. Introduction The solitary wave (soliton) consists of a single ‘hump’ of constant height and speed, and was first observed by Scott Russell in 1844 [l] for waves in channels. Shortly thereafter, Stokes [2] developed approximate expressions for nonlinear periodic waves in deep water and Boussinesq [3] and Rayleigh [4] obtained approximate expressions for the solitary wave. The solitary wave phenomenon is most easily characterized as a special solution of a nonlinear wave equation studied by Korteweg and de Vries [5]. Using perturbation theory for a low long wave in a shallow, narrow channel, they demonstrated solitary wave solutions and also periodic solutions. There is now a very extensive literature on the KdV equation and the mathematical properties of its solutions (see e.g. [6-91). For example, it can be shown that an arbitrary initial disturbance breaks up eventually to a set of solitons. Similarly, despite the nonlinearity of the equation, two or more solitons can pass through one another and emerge unaltered in shape. There have been few previous studies of finite element approximation of solitary waves. Schoombie [lo] employed a Galerkin finite element method and used Hermite cubic splines to accommodate the third derivative in the KdV equation. Argyris and Haase [ll] used cubic Schoenberg splines in their Galerkin finite element formulation and also examined a PetrovGalerkin method using piecewise linear test functions with the Hermite cubic trial basis. Other studies were made by Winther [12], Mitchell and Schoombie [13], and Sanz-Serna and Christie
P41. Least-squares finite element methods (see e.g. [15] have been the subject of recent exploratory (see e.g. [15]) investigations for first order systems ([16, 171) and have proven 00457825/91/$03.50
0
1991 Elsevier Science Publishers
B.V. All rights reserved
G.F. Carey, Y. Shen, Approximations
2
of de KdV equation
effective in the analysis of convection-diffusion, Burgers’ equation and compressible Euler equation problems. The method introduces a numerical dissipation that is similar to that for the Petrov-Galerkin and Taylor-Galerkin schemes for hyperbolic problems of convectiondominated transport processes. Our interest in the present problem stems from two main factors: first, the problem is of third order and has not previously been studied using the least-squares approach; secondly, the soliton solutions are special in that they are solutions to a nonlinear problem that have constant amplitude and velocity, and we wish to assess the effect of amplitude and phase errors introduced by the numerical method.
2. KdV equation The assumption that the wave amplitude is small compared with depth and that channel depth is much smaller than length leads to the KdV equation for long waves, which in normalized form becomes u, +
uu,
+
FU,,,
=
0
(1)
)
where F models the dispersion. The exact solution of (1) for a single soliton of amplitude constant velocity c is
3c moving to the right with
u,(x, t) = 3c sech*[k(x - ct) + x0] ,
(2)
with k = 1 /~(c/E)~“, x,, = k(2 - 1) + d where parameter d provides a reference for the initial position of the soliton. Solitons can pass through one another and emerge unchanged. For 2 or 3 solitons i with speed c,, initially centered about x0( the exact solution is [18]
u(x, t> = B(ln Q
,
(3)
where for two solitons F = 1 +
ez] +
ez2 + A ezl+‘2 ,
2, =
A = [(k, - k,)/(k, + k,)]* >
k; =
q&(x -
c,f> ++I
{(cJE)~‘*,
x0,
=
f d, ,
and for three F = det(l;l,) , 2; =
Fii = 6,. + 2(kjki)“2 e”/(k, + k,) ,
2[ k;(x - c;t) + X”,J )
x0; =
k(3c; - 1) + d; ,
with dj again a reference parameter for the respective reference value for the amplitude of the solitons.
soliton positions,
and B providing
a
G.F.
Carey, Y. Shen, Approximations
3
of de KdV equation
3. First-order system and least-squares FEM The KdV equation 8)t,, 06061,
(1) may be differenced
u (n+l) - $) At
with respect to time to give at t = r3tn+, + (1 -
+ 8(uu, + “Uxxx)(~+l)+ (1 - O)(uu, + EU,,,)(~) = 0.
(4)
Next, we recast this as a first-order system by setting U, = u, u, = w and linearize within each timestep At to obtain n+l) _ J4 Ll( At
+ @&$(*+l)
(n+l) - uf+‘))
@
O(W
Transposing conveniently
(n+l) _
+
+ EW;+‘)) + (1 - @(U@+J@)+ EWF)) = 0,
(1 - B)(u’“’ -
@+I)) + (I_ ux
fj)(w’“’
u_y )= 0) - $‘)
= 0.
(5)
terms at time level n and simplifying, we have a first order system which can be expressed in matrix forms as
&)zr+l)
+ &“)z(“+l’ =f’“’ ,
(6)
where zf = [u u w], A’“’ = A($“‘), B’“’ = B(z’“‘) and f’“’ = f(z’“‘). For an admissible solution 6 satisfying the boundary conditions, the corresponding residual vector is r = A(“)gT+l) + &“)[@+l) _f’“’ , and the least-squares +L 2 sR
residual functional
trial
(7)
on domain R is
rtt dx .
(8)
Setting the first variation FZ = 0 in (8), R [A’“’ @+” +
+ &“) @n+l)]f[A(n)J~+l)
The form of (9) corresponds to a Petrov-Galerkin the operators A’“‘, II’“’ as indicated here. Expanding the integrand expression in (9),
+ &n){(n+l) _f’“‘] dx = 0.
(9)
method where the test functions involve
(n+l)‘A(n)~(“)J;+l) + a{r+l)’ (n)’ @) @+I) + ~5(n+l)(B(n)lA(n)5lo+l) ADS I R [MI + 85 Introducing
(n+l)‘~(n)‘B(n)5!n+l,]dx = R (@‘+l)LA(“)’ + S{(n+l)L&#)j+) I
finite element expansions for the components
dx .
(10)
of 5 on a partition of the domain R
G.F.
4
Carey, Y. Shen, Approximations
of de KdV equation
and setting the variations 85 = (&, O,O), (0, $, 0), (O,O, x), where &, ei, yk are the associated component basis functions for U, u, W, yields the finite element system. It follows from the form of (10) that the discrete system is symmetric. For example, if the same basis is chosen for each solution component, then the element contributions K,, F, from (10) to KV= F can be written conveniently in matrix form as
and F, =
Y
(!&‘A’ + @B’)fdx
,
Ulb)
where, for notational convenience, we have omitted the timestep index n and W is a matrix of element basis functions. The element contributions (11) are accumulated to the global system. Essential boundary conditions are imposed to yield a sparse global system of the form KV=F
(14
to be solved at each timestep. Here K and F depend on the solution at the previous time level. The solution is initiated from the prescribed initial vector V” determined by approximation of U, u and w at t = 0 from specified u(x, 0) = g(x). The full form of the system (12) is given in Appendix A. In the preceding formulation we have treated the nonlinear term in (6) by writing uu, as uu and then linearizing to u (“)u(~+‘). This implies that u, u and w appear explicitly in the KdV equation. Alternatively we can use u (n)~(nf ‘) directly. Numerical comparisons for a single soliton indicate that the difference in Computed solutions with these methods is small (approximately one percent for the large soliton c = 0.3 computed in the numerical studies). Since the least-squares variational problem is based on a first-order system, C” Lagrange basis functions may be used. In particular, we may use piecewise-linear basis functions as compared with the cubic splines utilized in the previous studies cited. Thus, the additional complexity of the first-order system (more unknowns) is offset by the simplicity of the basis. The element integrations are carried out using Gaussian quadrature. Since the additional equations u, = u, u., = w are basically constraints, one might anticipate that similar difficulties to those encountered in penalized variational minimization problems might be encountered. These penalty functionals require reduced numerical integration for the method to be stable and the same situation appears to apply for the least-square method. The form of reduced integration rule used here is examined later in the numerical experiments. In other applications we have shown that the least-squares formulation introduces an artificial dissipation similar to that for the Petrov-Galerkin and Taylor-Galerkin methods. To investigate this further, for the present problem it suffices to integrate the leading term in (10) by parts. The corresponding contribution in the associated Euler-Lagrange equation associated with the weak form in (10) is s 6&‘A’Aj,, where the linearized form of the problem has been assumed. Note that the solution vector J is now differentiated twice. This implies that the least-squares method is equivalent to a Petrov-Galerkin form involving derivatives u,, , u,,,
G.F.
Carey, Y. Shen, Approximations
of de KdV equation
5
W XX;that
is, uXX,u,,, and u,,,,. Clearly the terms involving u,, and u,,,, will effect numerical dissipation and, therefore, the amplitude. Those involving u,,, will reinforce the terms EU,,, and influence the phase. Recall that in the soliton problem this term balances the nonlinear inertial term. Hence we anticipate that there will be some influence on the amplitude and phase properties of the propagating solitary waves. We explore these issues in the following numerical studies.
4. Numerical experiments In the calculations that follow, the approximate problem is posed on a finite domain R with essential boundary data u = 0 at the end points. For the soliton problems of interest here and the short time intervals considered, this provides a good approximation to the actual infinite domain problem for which the analytic solution is known.- _ 4.1. Single soliton As the first case, consider a propagating single soliton of the form given in (2) with E = 4.84 x lo-‘, c = 0.1 and d = 0. For the short time interval 0 < t < 4 considered here, it suffices to compute the approximation on a finite domain 0
x,J X pl7
Ax,, % clPan A% Au,iu,
(x,, up)
0
1
2
3
4
0.8 0.8 0 0.3 0.3000 0.0000 0.000
0.9 0.9 0 0.3 0.2990 0.0010 0.003
1.0 1.0 0 0.3 0.2969 0.0031 0.010
1.1 1.1 0 0.3 0.2955 0.0045 0.015
1.2 1.2 0 0.3 0.2940 0.0060 0.020
G.F.
6
1o
Carey, Y. Shen, Approximations
of de KdV equation
u(x,t)
. t
Fig. 1. Propagating
single solition in (x, t) domain for 0 s t G 4 (100 linear elements,
At = 0.01)
very dissipative for large amplitude solitons. Part of the problem can be attributed to the fact that the curvature w is approximated piecewise linearly and this approximation will deteriorate as the amplitude increases. In the present calculation the maximum error in w = u,,~ is 44%. It should be emphasized, however, that physically the KdV equation is derived under the assumption of long waves of low amplitude, so that the previous high-amplitude waves are not realistic for actual shallow water flow problems. For small amplitude waves that are of practical significance here, the present least-squares method is satisfactory. To examine the effect of mesh size let us consider the first case c = 0.1 with 50 elements. At this coarser mesh level there is a small secondary wave (oscillation) that propagates to the left. Similar observations have been made in other studies by Zabusky [193 and Leibovich and Seebass [20] and this effect is attributed to the negative group velocity associated with a linear wave. Pearson [21] has made an interesting analysis of a similar situation. Further numerical experiments indicate that the linearized problem is stable for l/2 d 8 d 1 and, as 0 increases in this range, dissipation increases (as one might expect in view of the well-known dissipative nature of the backward Euler method, 0 = 1.0). We note that dispersion increases with 8 and At. For instance, increasing the timestep from 0.005 to 0.01 increases the dispersion and dissipation errors: In the previous example with c = 0.3, the numerical dispersion at t = 1 is approximately given by Ax, = 0.02 and dissipation by Au, = 0.087 for calculations with 0 = 0.5, At = 0.01. If we use 0 = 1, the error increases to Ax,, = 0.04, Au, = 0.187. For 0 = 0.5, At = 0.005, the error decreases to Ax,) = 0.01, Au, = 0.061. The influence of At and 8 on dissipation can be seen in the presence of the term At du(‘)u~~+‘) identified in (A.4). The effect of linearization within a timestep is negligible for the step sizes At considered here. (In numerical experiments at At = 0.01 with 5 successive approximation corrections for the nonlinearity, the peak value changed at t = 0.2 from 0.8228 to 0.8286.) The quadrature rule in the previous calculations on linear elements was chosen to under-integrate the integral contribution using a l-point Gauss quadrature. If the order of integration is increased to a 2-point Gauss rule, the method fails. This situation is reminiscent of the type of difficulties that arise when penalty functionals are used to imbed constraints in a
G.F.
Carey, Y. Shen, Approximations
7
of de KdV equation
variational finite element model. (For instance, see [22,23].) An analogous situation occurs here since the additional equations U, = u, u, = w can be viewed as constraints which are included in the corresponding least-squares functionals. Higher-degree elements were also examined. The method does not fail with full integration but reduced integration gave results with better phase and dissipation properties. To examine this situation further we computed the eigenvalues of the element matrices for these terms when different quadrature rules were used and determined the rank of the matrix in each case. When linear elements are used the element matrix is of size 6 x 6 and has full rank if reduced integration is not applied. The desired reduced rank (3) is produced by underintegration. Similarly for quadratic and cubic element matrix contributions, from the ‘constraints’ we should get full rank if reduced integration is not used and the desired reduced rank of 6 and 8, respectively, if underintegration is applied. The reduced rank implies that the matrix contributions associated with the constraints are singular and therefore permits a nontrivial solution. However, we find that full integration yields slightly reduced rank, apparently due to the effect of roundoff in finite precision arithmetic.
4.2. Soliton interaction This second part of the numerical study deals with calculations for two and three interacting solitons. The analytic solutions for these respective cases are given in (3). As indicated previously, solitons of different amplitudes travel at different speeds and pass through each other to emerge unchanged. Numerical solutions for two solitons centered at 0.4 and 0.88, and travelling with speeds c, = 0.3 and c2 = 0.1 are shown in Fig. 2 during 0 G t s 4 for calculations with 100 linear elements on 0 < x < 2 using l-point integration, and 8 = 0.5, At = 0.01. In (3) B = 5.84 x 10-j and d, = d, = -5. Both soliton approximations exhibit numerical dissipation and dispersion as noted previously for the corresponding single wave calculations. Since the approximations are not true solitons, their interaction is not ‘pure’. For example, the small soliton (c = 0.1) now has 4% dispersion and 1% dissipation compared with essentially no
Fig. 2. Interaction
of two solitons in (x, t) domain for 0 s t G 4 (100 linear elements,
At = 0.01).
G.F.
Fig. 3. Interaction
Carey, Y. Shen, Approximations
of de KdV equation
of three solitons in (x, t) domain for 0 G t c 8 (200 linear elements,
At = 0.01).
dispersion or dissipation previously; the larger soliton now has 12% dispersion and 18% dissipation which (fortuitously) are slightly better than 20% and 23%, respectively, for the single soli ton propagation. The interaction of three solitons is shown in Fig. 3 where c, = 0.3, c2 = 0.2, c, = 0.1 and d, = d, = d, = -5, B = 5.84 x lo-’ in (3). The numerical results are for 200 linear elements with l-point integration on 0
Fig. 4. Separation
of initial profile c exp { - b(x - 2)‘) into propagating
waves.
9
G.F. Carey, Y. Shen, Approximations of de KdV equation 5.
Concluding remarks
Soliton solutions to the KdV equation provide an interesting class of test problems since the solitons are themselves of physical significance and they provide the opportunity to study numerical dispersion and numerical dissipation for the nonlinear wave problem. Here we formulate and apply a least-squares finite element method based on a representation of the KdV equation as a first-order system. For solitons of small amplitude the method is seen to exhibit little numerical dissipation or dispersion, but at higher amplitudes these influences are significant. Since the dissipation in this method is proportional to the timestep, reducing the timestep will improve the situation but at the cost of significantly increasing computation time. The results suggest that the least-squares method can be applied (conditionally as noted) to this type of nonlinear wave problem with success. An important point is the need to use reduced integration to produce a viable method.
Appendix A Equations (6)-(11) offer a compact matrix formulation. It is instructive also to develop the contributions explicitly. To do so we transpose terms at time level y1 to rewrite (5) for admissible U, u , w as Yl = uCn+‘) r2 = -&l,“+”
+ At &&n)u(n+U + At &wj;n+‘) - fi , + &j”+‘)
-f,
)
T) = -&J(““) x
+ &Jn+‘)
- f3 )
(A4
where f, = u@) - At(1 - ~~)(u(~)LI@) + EWE)), f2 = (1 - B)(u’“’ - ur’) and f, = (1 - 8)(w’“’ u(n) X )* Then I = 4 Jn (Y: + ri + r:) dx and setting 61= 0 implies Jf2 (Y, 6r, + y2 6r, -t r3 6r,) dx = 0 for arbitrary variations kCnt’), 6u(“+‘), ?3w(““). Introducing (A.l) into 61 = 0, collecting terms and simplifying we get
I
R (r, c%@+‘)- r,8 au?+‘)) dx = 0, with 8u(“+‘) arbitrary,
I
D (q At Ou(‘I 8u@+‘) + r,B &I(“+‘) - r,8 6~;““) with I~u(“~‘) arbitrary,
I
c%(“+‘) = 8~~~“) = 0 ;
&@+I) = 8w(“+‘) = 0 ;
dx = 0,
(A.21
R (I-~At C)& 6w$‘+” + ~$9 SW@+‘))dx = 0, with 8w(“+‘) arbitrary,
ILL@+‘)= ?3wCn+‘)= 0.
Introducing finite element expansions u, = C u~c$~,uh = C uk ekI,, wh = c w[‘yrand setting the variations as the respective basis functions &, I,!J~,y, leads to the finite element system
G.F.
10
In the numerical REMARK.
results we use piecewise-linear
Integrating $& J^
Carey, Y. Shetz, A~~roxi~at~otz.~ of de KdV eyzmtiotz
basis functions for uh, ulr and wk.
by parts in (AL?) we obtain
+(Y?_),8)sU(“+“dx=O,
i
dx = 0, R (Y, At8u (‘I)-I-r,B + (Y~)~B)&_J~‘~+‘) (A.4)
i
n[(-~‘)XA~6t:ir,B]Fw’“+“dx=0.
The third equation in (A.4) contains numerical dissipation.
a term AttSu’“‘~~‘~’ which can be interpreted
as a
References [l] J.S. Russell, Report on waves, Rept. 14th Meeting of British Association for the Advancement (John Murray, London, 1844) 311-390. [Z] G.G. Stokes, On the theory of oscillatory waves, Cambridge Trans. 8 (1847) 441-472. [3] J. Boussinesq, Thtorie des ondes et des remous qui se propagent le long d’un canal rectangulaire
of Science
horizontal,
G.F.
[4] [S] [6] [7] [8] (91 [lo] [ll] [12] [13] [14] [15] [16] [17] [18] [19] (201 [21] [22]
[23]
Carey, Y. Shen, Approximations
of de KdV equation
11
en communiquant au liquide contenu darts ce canal des vitesses sensiblement parailles de la surface au fond, J. Math. Pures Appl. 17 (2) (1872) 55-108. Rayleigh, On waves, Phil. Mag. 5 (1) (1876) 257-279. D.J. Korteweg and G. de Vries, On the change of form of long waves advancing in a rectangular canal and on a new type of long stationary wave, Phil. Mag. 39 (1895) 422-433. C.S. Gardner, J.M. Greene, M.D. Kruskal and R.M. Miura. Method for solving the Kortweg-de Vries equation, Phys. Rev. Lett. 1.5 (1967) 1095-1097. R.M. Miura, Korteweg-de Vries equation and generalizations I. A remarkable explicit nonlinear transformation, J. Math. Phys. 9 (1968) 1202-1204. G.B. Witham, Linear and Nonlinear Waves (Wiley, New York, 1984). A.C. Newell, Solitons in mathematics and physics, Philadelphia (198.5) 23-57. S.W. Schoombie, Finite element methods for the Korteweg-de Vries equation I. Galerkin method with Hermite cubits, Rept. NA/43, Department of Mathematics, University of Dundee, 1980. J. Argyris and M. Haase, An engineer’s guide to soliton phenomena: Application of the finite element method, Comput. Methods Appl. Mech. Engrg. 61 (1) (1987) 71-122. R. Winther. A conservative finite element method for the Korteweg-de Vries equation, Math. Comp. (1981). A.R. Mitchell and S.W. Schoombie, Finite element studies of solitons, in: R.W. Lewis, P. Bettes and E. Hinton, eds., Numerical Methods in Coupled Systems (Wiley, New York, 1984) 465-488. J.M. Sanz-Serna and I. Christie, A simple adaptive technique for nonlinear wave problems, J. Comput. Phys. 67 (2) (1986) 348-360. C.A.J. Fletcher, Computational Galerkin Methods (Springer, New York, 1984). G.F. Carey and B.N. Jiang. Least-squares finite elements for first-order hyperbolic problems, Internat. J. Numer. Methods Engrg. 26 (1988) 81-93. G.F. Carey and B.N. Jiang, Least-squares finite element method for compressible Euler equations, 5th Internat. Conf. Numerical Methods in Laminar and Turbulent Flow, Montreal, Canada, July 1987. R. Hirota, Exact solution of the Korteweg-de Vries equation for multiple collisions of solitons, Phys. Rev. Lett. 27 (1972) 1192-1194. N.J. Zabusky, Solitons and bound states of the time-independent Schrodingei equation, Phys. Rev. 168 (1968) 124-128. S. Leibovich and A.R. Seebass, Nonlinear Waves (Cornell Univ. Press, Ithaca, London, 1974). C.E. Pearson, Note on upwind and downwind effects of numerical approximation to wave equations, Comm. Appl. Numer. Methods (1988). R.L. Sani, P.M. Gresho, R.L. Lee, D. Griffiths and M. Engleman, The cause and cure (?) of the spurious pressure generated by certain FEM solutions of the incompressible Navier Stokes equations, Internat. J. Numer. Methods Fluids 1 (1981) 17-43, 171-204. G.F. Carey and R. Krishnan, Penalty approximation of Stokes flow, Comput. Methods Appl. Mech. Engrg. 35 (2) (1982) 169-206.