- Email: [email protected]
Solitary wave evolution for mKdV equations
ELSEVIER
Wave Motion 21 (1995)
263-275
Solitary wave evolution for mKdV equations N.F. Smyth”T* A.L. Worthy b a Deparfmenf of Mathematics and Statisrics, The King’s Buildings, University of Edinburgh, Mayjeld Road, Edinburgh EH9 332 UK h Departmen! of Mathematics, University of Wollongong, Northfields Avenue, Wollongong, N.S. W 2522, Australia Received 9 May 1994; revised 6 September
1994
Abstract The evolution of an initial condition into soliton(s) is the classic problem for the Korteweg-de Vries (KdV) equation. While this evolution is theoretically given by the inverse scattering solution of the KdV equation, in practice only the final steady state can be easily obtained from inverse scattering. However, an approximate method based on the conservation laws for the KdV equation has been found to give very accurately the evolution of an initial condition into soliton( s) . This approximate method also gives a criterion for the number of solitons formed. In the present work, this method is extended to describe the evolution of an initial condition into solitary wave(s) for mKdV equations, these equations having the same dispersive term as the KdV equation, but a nonlinear term of the form U”Ux, where n 2 1 is a positive integer. It is found that for n < 4, the behaviour of the mKdV equation is similar to the KdV equation in that solitary wave(s) evolve from an arbitrary initial condition. However, for n 2 4, it is found that an initial condition of sufficiently small amplitude decays into dispersive radiation with no solitary wave being formed. For an initial amplitude exceeding a threshold, it is found that the amplitude blows-up. The solutions of the approximate equations are compared with full numerical solutions of the mKdV equation and good agreement is found.
1. Introduction The Korteweg-de
Vries
equation
(1) arises in a large number of physical applications, for example long waves on the surface of a fluid, long internal waves in a stratified fluid and ion plasma waves [ 11. Typically, IQ. ( 1) is solved on the infinite line -CO < x < 00 together with the initial condition u(x,O)
=(/o(x),
-ccl < X < co.
(2)
For initial conditions cl0 which decay sufficiently rapidly as x + 4~~0, the problem ( 1) and (2) can be solved using the inverse scattering method (see Refs. [2,3]). Inverse scattering shows that the initial condition (2) * Corresponding
author. Fax +44-3 l-6506553,
e-mail [email protected].
0165-2125/95/$09.50 @ 1995 Elsevier Science B.V. All rights reserved SSDlOl65-2125(94)00053-O
264
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
breaks up into a finite number of solitons plus dispersive radiation. The number and details of the final steady solitons are given by solving an eigenvalue problem (Schrodinger’s equation plus zero boundary conditions at x = &too) and the time evolution of the solitons to their final steady states and the dispersive radiation are given by the solution of a linear integral equation (the Marchenko equation). While the solution of the eigenvalue problem to find the final steady solitons is reasonably straightforward, the solution of the integral equation is non-trivial. Hence it is not easy to determine the details of the evolution of the initial condition into the final steady solitons. The KdV equation has the mass and energy conservation equations &(n)
+ -$ (3U2fU,,)
=o
and
&(;u2) +$ (Zu’fuu,,g:) =o,
(4)
respectively. To determine the details of the evolution of an initial condition into solitons, an approximate method based on these conservation laws has been developed by Kath and Smyth [4]. For simplicity, the initial condition
Uo(x> = A sech2 g
(5)
was chosen. The soliton was then assumed to evolve as u = a sech2
x - 5(t) P
1
where
6’(t) = V(t),
(7)
V being the soliton’s velocity. The amplitude a and the width /3 of the soliton were both assumed to be functions of time t. Ordinary differential equations for the amplitude a and width p of the soliton were then obtained from the mass and energy conservation equations (3) and (4). These ordinary differential equations include terms representing the effect of the dispersive radiation and it is these terms which drive the solution to the steady state. The key result needed to obtain these ordinary differential equations was that the soliton moves faster than any shed radiation. This is because the phase and group velocities for linear dispersive waves for the KdV equation are negative and the soliton velocity is positive. It was then shown that this result implied that the mass in the radiation could be obtained by integrating the mass conservation equation from x = -cc to the soliton position. This result was further confirmed from the inverse scattering solution. The approximate equations derived from the mass and energy conservation equations for the KdV equation were found to give results in very good agreement with full numerical solutions of the KdV equation and with results from the inverse scattering solution. In the present work, the approximate method of Kath and Smyth [4] will be extended to determine the evolution of the solution of the initial value problem for the so-called mKdV equation
together with the initial condition U(X,O) =&l(x).
(9)
N.E Smyth, A.L. Worthy/Wave
Here n is a positive integer so that u” is always defined. The mKdV equation
and the solitary wave velocity v=
265
Motion 21 (1995) 263-275
has the solitary wave solution
V is
6 ~ a” II + 2
(12)
The mKdV equation (8) arises for tl = 2 in nonlinear optics [ 51 and in the propagation of long internal waves in a fluid when the coefficient of the usual nonlinear term in the KdV equation UU, is zero and the higher order nonlinear term U*U, dominates over higher order dispersive terms [ 61. For n < 4, it is found that the evolution of the initial condition (9) is similar to that for the KdV equation. That is, the initial condition (9) evolves into solitary wave(s) of the form (10) plus dispersive radiation. The solutions of the approximate equations for n < 4 are found to be in good agreement with full numerical solutions of the mKdV equation. However, for n > 4, the approximate equations give that for small enough initial amplitude, the initial condition decays into dispersive radiation only and no solitary waves form. For initial amplitude above a critical value, the approximate equations give that the solitary wave amplitude becomes infinite. The decay of the initial condition into dispersive radiation alone for initial amplitude below a critical value is found to agree with full numerical solutions of the mKdV equation for n 2 4. The full numerical solutions are found to give less conclusive evidence of the blow-up above a critical threshold amplitude, but do indicate that such blow-up can occur (and is not an artifice of the numerical method). These results agree with Benjamin [ 71, Bona [ 81, and Weinstein [9], who found that the solitary wave solution of the mKdV equation is stable for 1 5 n 2 3, and Bona et al. [ lo], who found that the solitary wave solution is unstable for n > 4. The new result of the present work is the existence of a threshold for the initial amplitude below which the initial condition decays into radiation alone and above which it blows-up. The mKdV equation (8) can be solved by the inverse scattering method in the cases n = 1 and tz = 2 only. In the case n = 2, the mKdV equation can be transformed to the KdV equation using the Miura transformation [ 31. So for n > 2, there are no exact methods for finding the solution of (8) and (9). The present work then shows that the method of Kath and Smyth [4] can be successfully extended to equations for which there are no exact solutions.
2. Approximate For simplicity, U,(x)
equations we shall take the initial condition
= A sech”” g,
(9) to be
-cc<<<<.
This initial condition is chosen as it has the same sech*/’ want to study the evolution of the initial condition (13) and energy conservation equations for the mKdV equation energy conservation equations for the mKdV equation (8) i(U)
+ & (3u”+’ + U,,) = 0
(13) profile as the solitary wave solution ( 10). We then for the mKdV equation (8). To do this, the mass are used, as in Kath and Smyth [4]. The mass and are (14)
266
N.E Smyth, A.L. Worthy/Wave
Motion 21 (1995) 263-275
and
(15) respectively. It is expected that the initial condition (13) will evolve into a solitary wave (or solitary waves). seek a solution of the mKdV equation of the form
We then
u=uo+u1,
(16)
where uo = a sech21n :,
(17)
5’(t) = V(t).
e=x-5(t),
(18)
The term ua is the time varying solitary wave which evolves from the initial condition (13) and the term ui accounts for the dispersive radiation which is shed as the solitary wave evolves. It will be assumed that the amplitude of the dispersive radiation is much smaller than the solitary wave amplitude, so that Iui 1< UO.Here t(r) is the position of the solitary wave maximum, V(t) is the solitary wave velocity and the amplitude a and width p of the solitary wave are functions of time t. From the initial condition ( 13), we see that
a(0) = A and p(O) = b.
(19)
To determine the evolution of the solitary wave, equations for a, p and V need to be found. Substituting the expansion (16) into the mKdV equation (8), we obtain the following equation
=-a’ sech21n % + -$ + 2(tGi)a
for UI
- n2Vp2) sech2/” $ tanh $
(4 - n2@p’
[3n2ap2 - 2(n + 2)] sech2-t2’n B tanh 8, P P
(20)
on neglecting terms of order u:. The phase and group velocities for linear dispersive waves for the mKdV equation (8) are -k2 and -3k2 respectively, where k is the wavenumber. Therefore, since the solitary wave velocity is positive, radiation is quickly left behind by the evolving solitary wave. It will hence be assumed that ahead of and in the neighbourhood of the solitary wave, u1 rapidly decays to zero. For this to be the case, the solitary wave position must be correct. This can be seen from Eq. (20) for ut when ( 11) holds, so that un is an exact solitary wave. In this case, a homogeneous solution of (20) is
u1 = uox = --
2a
x - vt
sech21n -
nP
P
x - vt
tanh -
P
.
(21)
Adding a small amount of this solution to ua then amounts to just a phase shift in ua, which results in a shift in the position of the solitary wave. If the solitary wave position has been accurately determined however, then ui must rapidly decay to zero ahead of and in the vicinity of the solitary wave. Assuming this to be the case, then integrating Eq. (20) (or the mass conservation equation ( 14)) from x = ,$ to x = co, we obtain
$ap)=3a’+’
- Va - iapp2,
(22)
N.F: Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
where B( p, q) is the Beta function.
M=
s UI
26-l
Let us set
dx,
(23)
so that M is then the mass shed in the radiation. the solitary wave as
Integrating
(17)
from x = -cc
to x = cc gives the mass in
m
u. dx = B
(24)
Hence from (23) and (24)) mass conservation B
becomes
( >-&3)+~=0. ;,;
(25)
The consideration of mass conservation is now completed and we have two equations for the four unknowns a, /?, M and V. To obtain the two equations needed to complete the system, energy conservation is now considered. From the energy conservation equation ( 15) and the expansion ( 16)) we see that the energy density is
(26) Now it has been assumed above that ~1 is zero in the neighbourhood of the solitary wave. Hence the cross term ~0~1 is zero and the contribution to the energy density from the radiation is second order in ui. If this second order term is neglected, then integrating the energy conservation equation ( 15) from n = --00 to x = cc gives ;B
( >;(a’P) ;,;
=0
(27)
on using ( 17). To complete the system of equations describing the evolution of the solitary wave, an expression velocity V must be found. This expression is found from the moment of energy equation
for the
for the mKdV equation (8). Moment of energy, rather than moment of mass, was chosen to calculate the velocity since, as explained above, any contribution to the moment of energy equation from ui is second order. Then integrating the moment of energy equation from x = --oo to x = cc and neglecting these second order contributions gives
!LV= dt
Wn
+ 1)
(n+2)(n+4ja”-
n(n1t4)p-2
on using the energy conservation equation (27). Note that for large times, the relation ( 11) holds and this velocity approaches the solitary wave velocity (12), as required. The system of equations (22), (25), (27) and (29) gives a description of the evolving solitary wave which is in good agreement with numerical solutions of the mKdV equation. However, as in Kath and Smyth [ 41, it is found that the accuracy of the equations can be improved by including an estimate of the energy associated with the shed radiation. In perturbed soliton problems, mass and energy leak away from the sohton due to the perturbed soliton approaching a non-zero constant as x -+ --cc [ 111, resulting in the formation of a so-called
268
N.E Smyth. A.L. Worthy/Wave Motion 21 (1995) 263-275
shelf behind the soliton. It is this leakage of mass and energy that gives rise to the dispersive radiation behind the solitary wave. As in Kath and Smyth [4] for the KdV equation, it can easily be shown that the same effect occurs for the mKdV equation. Let us then suppose that the mass loss from the solitary wave is due to the solitary wave approaching a small non-zero value u, in the shelf, which we will take to occur near x = -L. The shelf is assumed to follow at a fixed distance behind the solitary wave, so that it moves with the solitary wave velocity V. Then integrating the mass conservation equation (14) from x = -L to x = 0;) (i.e. over the solitary wave) gives
u dx = -Vu,
(30)
-L
on neglecting
terms quadratic
and higher in ui. Comparing
(23),
(25) and (30),
we see that (31)
In a similar manner,
integrating
the energy conservation
equation
(15) gives
(32)
on neglecting
cubic and higher terms in ui. Using (31), this becomes
(33) -L
This is an expression for the energy shed from the evolving solitary wave. Combining energy conservation result (27) yields the energy conservation equation B
(;,a>i
(a’p)
=
-$ (T>‘.
it with the first order
(34)
It can be seen that the velocity expression (29) can be zero at non-zero values of a, which can introduce a non-physical singularity into the energy conservation equation (34). The above derivation of the energy loss, based upon the formation of a shelf behind the solitary wave, is strictly speaking valid in the limit of a slowly varying solitary wave. It is therefore reasonable to use the solitary wave velocity (12) for the velocity V in (34), this velocity expression not having the artificial singularity of (29). Furthermore, using a different expression for V in (34) will produce errors for only a finite time. So, since the rate of energy loss is small, the total integrated difference will also be small. The system of equations governing the evolution of the solitary wave is then (22), (25), (29) and (34). These equations are the same as those of Kath and Smyth [4] in the KdV case n = 1. While these equations were found to give results in good agreement with numerical solutions of the KdV equation, Kath and Smyth [ 41 found that better agreement could be obtained by generalising the form of the velocity expression (29) for n = 1 in the following manner V = 2a + ua - 2vpe2,
(35)
where Y is a constant. In the limit of a soliton, so that ap2 = 2, this velocity approaches the correct soliton velocity, independent of the value of v. Kath and Smyth [4] found that the value v = 3/2 gave the best
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
269
agreement between solutions of the approximate conservation equations and full numerical solutions of the KdV equation, independent of the initial values of a and p. The parameter v in the velocity expression (35) is a decay constant controlling the rate at which the solution approaches the steady state, but it does not have a large effect upon the value of this steady state. The same method for determining an improved expression for V will be used in the present work. We therefore set v=
6 -----un n+2
+ zX2”- Y
2(fi + 2) P_2 3n2 ’
(36)
noting that in the limit of a solitary wave, in which the amplitude-width relation (11) holds, this velocity expression is the same as the solitary wave velocity (12). Again, the parameter v will be chosen to give the best agreement with numerical solutions of the mKdV equation. We note that the moment of energy result (29) gives 18n u=
(37)
(n+2)(n+4)‘
While the major role of Y is to determine the rate at which the solution approaches the steady state, it also has an effect on the stability of the approximate conservation equations. To study this, let us ignore the small loss of energy in (34), so that a2P = A2b, on noting (36) is a04-n
(38)
the initial
=
values of a and /?. Then the steady state of the approximate
3 (A’b)’ n2 2 (n+2>
and
PO = ~u;‘d2 &in
equations
(22),
(27)
and
’ (39)
Linearising
about this steady state by setting a = an + 2, we obtain
it’=--&
(s-v)
(n-4)u;ff5i?,
(40)
where E=A2b
and
The approximate
V>-&
p=B
equations
(41) are then stable for (42)
for n < 4 and v<-
3n n+2
(43)
for n > 4. Comparing the value (37) of v given by the moment of energy equation with the stability constraints (42) and (43), we see that the moment of energy velocity expression (29) gives unstable equations for n = 3 and stable equations for all other values of n. Furthermore, the equations for n = 2 and n = 4 are neutrally
270
N.F: Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
t Fig. 1. Solitary wave amplitude a as a function of time t for n = 3, initial amplitude and h = 1.O. Conservation laws: --; numerical solution of mKdV equation: - - -.
(a) A = 0.7 and initial width h = 1.O and (b) A = 0.9
stable as the growth rate in (40) is zero for these values of n. This result for n = 4 is not surprising as it will be found in the next section that for n = 4 the solitary wave solution of the mKdV equation is unstable. However, it is very surprising that for n = 3 moment of energy gives an unstable set of equations. Numerical solutions of the mKdV equation for IZ = 3 show no qualitative differences from numerical solutions for other values of n (n < 4). In the next section, it will be found that by adjusting the value of v in (36), very good agreement can be obtained with numerical solutions of the mKdV equation. The same also holds true for the n = 2 case. For n = 2, it is surprising that moment of energy yields only a neutrally stable set of equations since the mKdV equation for n = 2 can be transformed to the KdV equation by the Miura transformation and moment of energy for the KdV equation yields a stable set of equations.
3. Results and comparison
with numerical
solutions
In this section, solutions of the approximate equations (22), (25) and (34) together with the velocity expression (36) will be compared with full numerical solutions of the mKdV equation (8) with the initial condition ( 13). The approximate equations are nonlinear and cannot be solved exactly in general. For the results presented in this section, these equations are solved numerically using a fourth order Runge-Kutta method. For the numerical integration of the mKdV equation, an adaption of the pseudo-spectral method of Fornberg and Whitham [ 121 was used. This involves using FIT to evaluate the space derivatives and a fourth order Runge-Kutta method to propagate the scheme forward in time. The cases n = 2, n = 3, n = 4 and n = 5 will be dealt with as these give the full range of possible behaviours. The case n = 1 (the KdV equation) was dealt with in Kath and Smyth [4 J. For n = 2, the mKdV equation is exactly integrable, but not for n = 3, n = 4 and n = 5. It will be found that for n < 4, the behaviour of the solution is similar to the n = 1 KdV case. For n 2 4, it will be found that the initial condition decays into dispersive radiation with no solitary wave being formed if the initial amplitude is small enough. For initial amplitude above a critical value, the amplitude of the solitary wave blows-up. Let us consider the case n = 3 first since there is no inverse scattering solution of the mKdV equation with n = 3. The best fit of the solutions of the approximate conservation equations with full numerical solutions was found for the parameter value Y = 2.2 in the velocity expression (36). This value of Y was found to give the best fit independent of the initial values of a and p (A and 6). Fig. la shows a comparison between the
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
271
solitary wave amplitude as given by the solution of the approximate conservation equations and as given by the full numerical solution for an initial amplitude of A = 0.7 and an initial width b = 1. It can be seen that the agreement is very good for the final steady state amplitude. However, the agreement between the transient behaviours is not quite so good, but is adequate. Fig. lb shows the same comparison for A = 0.9 and b = 1. The agreement between the solutions is similar to that in Fig. la, showing that the value of v is independent of the initial condition chosen. Fig. 2 shows a comparison between the final steady amplitudes of the solitary wave produced as given by the approximate conservation equations and by the full numerical solution of the mKdV equation for a range of initial amplitudes and initial width b = 1.O. It can be seen that the comparison is excellent, as can be noted from Figs. la and lb for special cases. When the area of the initial condition (13) is large enough, more than one solitary wave will eventually form. The initial amplitude and width at which this change occurs can be determined from the approximate conservation equations (22)) (25)) (34) and (36). The point at which the change from producing one final solitary wave to producing two final solitary waves occurs is given by the values of A and b for which M’ = 0 initially. This is because if M’ = 0 initially, then there is exactly enough mass in the initial condition to produce one solitary wave and any more mass than this will result in another solitary wave being produced. Setting M’ = 0 in (25)) we see from (22) that a’ = 0 and p’ = 0 and that a and p are related by the solitary wave relation ( 11) . Then using the expression (36) for the solitary wave velocity, we find from (22) that
For a initially greater than this value, 2 solitary waves will form from the initial condition. waves will then form when A” >
g&2
For n = 3, 2 solitary
(45)
This result was found to be in good agreement with numerical solutions, although it was difficult to verify the exact point at which 2 solitary waves first formed due to the formation time for the second solitary wave becoming infinite as the change over point was approached. For b = 1, 2 solitary waves are formed when A > 0.71814. In Fig. 2, the amplitude plotted when A > 0.71814 is the amplitude of the hugest solitary wave. We then see that the approximate conservation equations give good results even when 2 solitary waves are formed. The mKdV equation for n = 2 was also considered, the mKdV equation being related to the KdV equation by the Miura transformation in this case [ 31. The value of v in the velocity expression (36) which gave the best agreement with numerical solutions was found to be v = 2.0. Fig. 3 shows a comparison between the final steady solitary wave amplitude as given by the approximate conservation equations and that given by the numerical solution of the mKdV equation as a function of the initial amplitude A for b = 1. Again it can be seen that the agreement is excellent. Furthermore, the relation (44) shows that 2 solitons will form for A2 > ;b-2,
(46)
which is A > 0.81650 for b = 1. Hence Fig. 3 shows that the approximate conservation equations again give good agreement for the largest soliton even when 2 solitons are formed. The point at which two solitons first form, given by (46), was again found to be in good agreement with numerical solutions. The mKdV equation for n = 2 and n = 3 has similar behaviour to the KdV equation. However, the behaviour for n 2 4 is completely different. This can be seen by solving the approximate conservation equations (22), (25), (34) and (36) when the energy leakage is ignored in Eq. (34). The energy leakage is usually small, so the results obtained by neglecting it are a good approximation to the solution of the full equations. On ignoring the energy leakage, Eq. (34) gives
272
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275 1.6
l
1.3
t
1.4 -
1.2 -
t
.
1.1 -
1.2 _ P
l-
.
1.
a .
.
0.9 _
.
a 0.8 -
.
0.8 -
.
0.7 . .
0.6 -
.
0.6 -
.
0.5 _ '+I.6
0.65 t
0.7
0.75
0.6
0.65
0.9
0.95
1
1.05
t
0.1 0.6
0.65
0.7
03s
A
0.8
0.81
0.9
0.95
I
,.oi
A
Fig. 2. Final steady solitary wave amplitude a as a function of initial amplitude A for n = 3 and initial width b = 1. Conservation laws: 0 ; numerical solution of mKdV equation: +. Fig. 3. Final steady solitary wave amplitude a as a function of initial amplitude A for n = 2 and initial width b = 1. Conservation laws: 0 ; numerical solution of mKdV equation: +. a2P =
E,
(47)
where E is a constant.
Eq. (22) then becomes (48)
Integrating
this equation,
$ (a-4
_ A-“)
we obtain
+ ;E-2
(a-2
_ A-2)
+ 2E-4
log (2a2 - 3E2)A2
= 1
(3
_ v) t,
(49)
(2A2-3.!?)a2 ITE 2
9
('- i?E-2a)A+ ;(a-5 -A-5)+ GE-2 $fEe2A)a
(a-4
-A-4)+
;($,'E-4
(a-3
_ A-3)
(l-
+i(g)3E-6
(a-‘-
A-=)
and
a6=A6
l+ [
3 -A6(2
- v)(4
CLUE
+ (i)4Ee8
(a-’
-I 1
- E-=)t
for the cases n = 2, n = 3 and n = 4, respectively. and
~4 = B
- A-‘)
= s(g
- .>,
(51)
Here
(52)
It can be seen from (49) and (50) that a approaches a constant value as t + cc for n = 2 and n = 3. However, for the case n = 4, (51) gives that a either approaches 0 as t + 00 or a approaches oc in finite time, depending on the sign of (2 - Y) (4 - Ed=). For A small (i.e. E small), numerical results show that a decreases as t increases. We see from (48) that we then require Y > 2 for this to be the case. The best fit with numerical solutions was found for Y = 4.2. Thus from the solution (51) for a it is seen that for A below a
213
N.F: Smyth. A.L. Worthy/Wave Motion 21 (1995) 263-275
t
t
Fig. 4. Solitary wave amplitude ll as a function of time t for n = 4. (a) Initial amplitude A = 0.4 and initial width h = 2.5. Approximate conservation equations: -; numerical solution: - - -. (b) Numerical solution for initial amplitude A = 0.8 and initial width h = I .O.
threshold (E < l/2), a -+ 0 as t + co and for A above this threshold, a + 0;) in finite time. This result is confirmed by full numerical solutions of the mKdV equation for n = 4. Fig. 4a shows the numerical solution of the mKdV equation for n = 4, initial amplitude A = 0.4 and initial width b = 2.5. It can be seen that the solitary wave amplitude is approaching 0 rather than approaching a steady non-zero value. The solitary wave decays to zero and all of its initial mass and energy is shed as dispersive radiation. This is in accord with the predictions of (5 I), for which 4 - Ee2 is negative for these values of A and b. Also shown in Fig. 4a is the solution of the approximate conservation equations. The agreement between the approximate and numerical solutions is reasonable. The agreement is less close than for the previous values of n as in this case the solitary wave is decaying to zero. The assumption of small radiation made in the derivation of the approximate conservation equations is then less valid. In choosing the value Y = 4.2, a compromise was made between agreement with the numerical solution for small and large times. Other values of 1/ will give better agreement at either small times or large times, but not both. The numerical evidence for the blow-up of the solitary wave amplitude for E > l/2 is less clear. This behaviour of the solution for n = 4 is consistent with the instability of the solitary wave solution of the mKdV equation for n = 4 found by Bona et al. [ lo]. However, their stability analysis does not give the threshold behaviour for the evolution of the initial condition as predicted by the approximate conservation equations. It was found numerically that the amplitude did become very large in finite time for E sufficiently large. Furthermore, the numerically determined initial amplitude at which the solitary wave evolution behaviour changed from decay to blow-up agreed very well with the approximate relation E = l/2. For example, taking b = 1, it was found numerically that the value of the initial amplitude A at which the amplitude evolution changed from decay to growth was A = 0.707, which is in excellent agreement with the approximate value from E = I /2, A = fi. Indeed the degree of agreement is surprising given the approximations made in deriving the governing equations (22)) (25), (34) and (36). To obtain the best agreement between the full numerical and approximate solutions, the value of v in the velocity expression (36) was determined by fitting with the numerical solution. However, the value of the threshold (E = l/2) at which the solitary wave evolution changes from decay to blow-up is independent of the value of v, as can be seen from the solution (5 1) . The prediction of this change in behaviour then does not depend on any fit with the full numerical solution. Fig. 4b shows the numerical solitary wave amplitude as a function of time for initial amplitude A = 0.8 and initial width b = 1.O (for which E > l/2). It can be seen that the amplitude is rapidly becoming large. Indeed, the amplitude was found to become beyond the computer capacity for t not much larger than 6. This
214
N.E Smyrh, A.L. Worthy/Wave Motion 21 (1995) 263-275
numerical blow-up was not due to numerical error as it was found to occur even if the space and time steps were decreased. It is of course problematical to say from numerical results that the solitary wave amplitude is becoming infinite since as the amplitude increases, the width decreases, which results in the loss of spatial resolution over the solitary wave. However, there is numerical evidence for the predicted behaviour of the n = 4 case, in particular the blow-up for sufficiently large values of E. As the value of n increases, the evolution of the initial condition (13) becomes slower. Let us now consider the mKdV equation for n 2 5. For low initial values of A, numerical solutions show that the amplitude a must decrease with time t. From (48) we see that this requires U>-
3n n+2’
We then see from the stability condition (43) that the solitary wave is unstable according to the approximate conservation equations, which agrees with the results of Bona et al. [lo]. This is also borne out by numerical solutions of the mKdV equation. These numerical solutions show that the behaviour of the mKdV equation for n > 4 is similar to its behaviour for n = 4 in that for values of the initial amplitude below a critical value, the amplitude of the solitary wave decays to zero as t + co, while for values of the initial amplitude above this critical value, the amplitude approaches infinity as t -+ co. This change in the behaviour of the solution depending upon the value of the initial amplitude can be seen from (48). Since Y satisfies (53), we see that for
a decreases A”
>
as r increases,
while for
2(n + 2) b_2 3n2
’
a increases as t increases. The instability then manifests itself as the solitary wave decaying to zero as c -+ cc when (54) holds and the solitary wave amplitude becoming infinite as t + co when (55) holds. As for the case n = 4 discussed above, the change in the behaviour of the solution at A” = 2(n
+ 2, b-2
3n2 is independent of the value of v and so is independent of any fitting of the value of v from the full numerical solution. To verify these results, we shall consider the case n = 5 in detail. For n = 5, the best agreement with numerical solutions was obtained for Y = 2.46. Fig. 5a shows the numerical solution of the mKdV equation for A = 0.7 and b = 1.0, these values satisfying (54). It can be seen that the solitary wave amplitude is approaching zero, although on a slow time scale. Also shown is the solution of the approximate conservation equations and it can be seen that the agreement between the two solutions is reasonable. As for the n = 4 case, due to the solitary wave decaying to zero, the dispersive radiation is eventually not small and hence good agreement is not to be expected. Fig. 5b shows the numerical solution of the mKdV equation for A = 0.9 and b = 1.O, these values satisfying (55). It can be seen that the solitary wave amplitude is increasing and seems to be blowing-up. There is thus numerical confirmation of the two different types of behaviour depending on the values of A and b. As for the case n = 4 discussed above, the value of the threshold (56) at which the solution behaviour changes from decay to growth can be compared with that given by the full numerical solution of the mKdV equation. When b = 1, the relation (56) gives that the solution behaviour changes when the initial amplitude is A = 0.7148, which is in excellent agreement with the numerical value of A = 0.714. As for the case n = 4, this degree of agreement is surprising given the approximations made in deriving the governing equations (22), (25), (34) and (36), but it does however confirm the basic validity of these assumptions.
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
275
a
t Fig. 5. Solitary wave amplitude a as a function of time t for n = 5. (a) Initial amplitude A = 0.7 and initial width h = 1.0. Approximate conservation equations: --; numerical solution: - - -. (b) Numerical solution for initial amplitude A = 0.9 and initial width h = I .O.
Based on the predictions of the approximate conservation equations, the behaviour of the mKdV equation for II > 5 is expected to be similar to the II = 5 case, which is in accord with the results of Bona et al. [ IO]. While Bona et al. [lo] show that the solitary wave solution of the mKdV equation is unstable for n 2 4, their stability analysis does not predict the ultimate behaviour of the solitary wave due to small perturbations. The existence of a critical threshold for n 2 4 below which the initial condition decays into radiation and above which the initial condition blows-up is the new result predicted by the approximate conservation equations.
Acknowledgement The authors would like to thank Professor W.L. Kath for help with the numerical presented in this paper.
routines
used for the work
References 1 I ] G.L. Lamb, Elemenrs of Soliton Theory, Wiley-Interscience,
New York ( 1980).
[ 21 A.C. Newell, Solitons in Mathematics and Physics, SIAM, Philadelphia ( 1985). [ 3 1 G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York ( 1974). [ 4 1 W.L. Kath and N.F. Smyth, “Soliton evolution and radiation loss for the Korteweg-de Vries equation”, Phys. Rev. E 51 ( 1)) 66 l-670 (1995).
[ 5 1 G.P Agrawal, Nonlinear Fibre Optics, Academic Press, Boston ( 1989). [ 6 ] J.A. Gear and R. Grimshaw, “A second-order theory for solitary waves in shallow fluids”, Phys. Fluids 26(l), 14-29 (1983). [ 7 1 T.B. Benjamin, “The stability of solitary waves”, Proc. R. Sot. London, Sect. A 328, 153-183 ( 1972). [ 8 1 J.L. Bona, “On the stability of solitary waves”, Proc. R. Sot. London. Sect. A 344, 363-374 ( 1975). 19 ] M. Weinstein, “Lyapunov stability of ground states of nonlinear dispersive evolution equations”, Commun. Pure Appl. Mafh. 39, 51-68 (1986). 1 101 J.L. Bona, PE. Souganidis and W.A. Strauss, “Stability and instability of solitary waves of Korteweg-de Vries type”, Pm,. R. Sot. London, Sect. A 411, 395-412 (1987). [ I I] C.J. Knickerbocker and A.C. Newell, “Shelves and the Korteweg-de Vries equation”, J. Fluid Mech. 98, 803-8 I8 ( 1980). [ 12 1 B. Fomberg and G.B. Whitham, “A numerical and theoretical study of certain nonlinear wave phenomena”, Phil. Trans. R. Sac. London A 289, 373-403 ( 1978).
Wave Motion 21 (1995)
263-275
Solitary wave evolution for mKdV equations N.F. Smyth”T* A.L. Worthy b a Deparfmenf of Mathematics and Statisrics, The King’s Buildings, University of Edinburgh, Mayjeld Road, Edinburgh EH9 332 UK h Departmen! of Mathematics, University of Wollongong, Northfields Avenue, Wollongong, N.S. W 2522, Australia Received 9 May 1994; revised 6 September
1994
Abstract The evolution of an initial condition into soliton(s) is the classic problem for the Korteweg-de Vries (KdV) equation. While this evolution is theoretically given by the inverse scattering solution of the KdV equation, in practice only the final steady state can be easily obtained from inverse scattering. However, an approximate method based on the conservation laws for the KdV equation has been found to give very accurately the evolution of an initial condition into soliton( s) . This approximate method also gives a criterion for the number of solitons formed. In the present work, this method is extended to describe the evolution of an initial condition into solitary wave(s) for mKdV equations, these equations having the same dispersive term as the KdV equation, but a nonlinear term of the form U”Ux, where n 2 1 is a positive integer. It is found that for n < 4, the behaviour of the mKdV equation is similar to the KdV equation in that solitary wave(s) evolve from an arbitrary initial condition. However, for n 2 4, it is found that an initial condition of sufficiently small amplitude decays into dispersive radiation with no solitary wave being formed. For an initial amplitude exceeding a threshold, it is found that the amplitude blows-up. The solutions of the approximate equations are compared with full numerical solutions of the mKdV equation and good agreement is found.
1. Introduction The Korteweg-de
Vries
equation
(1) arises in a large number of physical applications, for example long waves on the surface of a fluid, long internal waves in a stratified fluid and ion plasma waves [ 11. Typically, IQ. ( 1) is solved on the infinite line -CO < x < 00 together with the initial condition u(x,O)
=(/o(x),
-ccl < X < co.
(2)
For initial conditions cl0 which decay sufficiently rapidly as x + 4~~0, the problem ( 1) and (2) can be solved using the inverse scattering method (see Refs. [2,3]). Inverse scattering shows that the initial condition (2) * Corresponding
author. Fax +44-3 l-6506553,
e-mail [email protected].
0165-2125/95/$09.50 @ 1995 Elsevier Science B.V. All rights reserved SSDlOl65-2125(94)00053-O
264
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
breaks up into a finite number of solitons plus dispersive radiation. The number and details of the final steady solitons are given by solving an eigenvalue problem (Schrodinger’s equation plus zero boundary conditions at x = &too) and the time evolution of the solitons to their final steady states and the dispersive radiation are given by the solution of a linear integral equation (the Marchenko equation). While the solution of the eigenvalue problem to find the final steady solitons is reasonably straightforward, the solution of the integral equation is non-trivial. Hence it is not easy to determine the details of the evolution of the initial condition into the final steady solitons. The KdV equation has the mass and energy conservation equations &(n)
+ -$ (3U2fU,,)
=o
and
&(;u2) +$ (Zu’fuu,,g:) =o,
(4)
respectively. To determine the details of the evolution of an initial condition into solitons, an approximate method based on these conservation laws has been developed by Kath and Smyth [4]. For simplicity, the initial condition
Uo(x> = A sech2 g
(5)
was chosen. The soliton was then assumed to evolve as u = a sech2
x - 5(t) P
1
where
6’(t) = V(t),
(7)
V being the soliton’s velocity. The amplitude a and the width /3 of the soliton were both assumed to be functions of time t. Ordinary differential equations for the amplitude a and width p of the soliton were then obtained from the mass and energy conservation equations (3) and (4). These ordinary differential equations include terms representing the effect of the dispersive radiation and it is these terms which drive the solution to the steady state. The key result needed to obtain these ordinary differential equations was that the soliton moves faster than any shed radiation. This is because the phase and group velocities for linear dispersive waves for the KdV equation are negative and the soliton velocity is positive. It was then shown that this result implied that the mass in the radiation could be obtained by integrating the mass conservation equation from x = -cc to the soliton position. This result was further confirmed from the inverse scattering solution. The approximate equations derived from the mass and energy conservation equations for the KdV equation were found to give results in very good agreement with full numerical solutions of the KdV equation and with results from the inverse scattering solution. In the present work, the approximate method of Kath and Smyth [4] will be extended to determine the evolution of the solution of the initial value problem for the so-called mKdV equation
together with the initial condition U(X,O) =&l(x).
(9)
N.E Smyth, A.L. Worthy/Wave
Here n is a positive integer so that u” is always defined. The mKdV equation
and the solitary wave velocity v=
265
Motion 21 (1995) 263-275
has the solitary wave solution
V is
6 ~ a” II + 2
(12)
The mKdV equation (8) arises for tl = 2 in nonlinear optics [ 51 and in the propagation of long internal waves in a fluid when the coefficient of the usual nonlinear term in the KdV equation UU, is zero and the higher order nonlinear term U*U, dominates over higher order dispersive terms [ 61. For n < 4, it is found that the evolution of the initial condition (9) is similar to that for the KdV equation. That is, the initial condition (9) evolves into solitary wave(s) of the form (10) plus dispersive radiation. The solutions of the approximate equations for n < 4 are found to be in good agreement with full numerical solutions of the mKdV equation. However, for n > 4, the approximate equations give that for small enough initial amplitude, the initial condition decays into dispersive radiation only and no solitary waves form. For initial amplitude above a critical value, the approximate equations give that the solitary wave amplitude becomes infinite. The decay of the initial condition into dispersive radiation alone for initial amplitude below a critical value is found to agree with full numerical solutions of the mKdV equation for n 2 4. The full numerical solutions are found to give less conclusive evidence of the blow-up above a critical threshold amplitude, but do indicate that such blow-up can occur (and is not an artifice of the numerical method). These results agree with Benjamin [ 71, Bona [ 81, and Weinstein [9], who found that the solitary wave solution of the mKdV equation is stable for 1 5 n 2 3, and Bona et al. [ lo], who found that the solitary wave solution is unstable for n > 4. The new result of the present work is the existence of a threshold for the initial amplitude below which the initial condition decays into radiation alone and above which it blows-up. The mKdV equation (8) can be solved by the inverse scattering method in the cases n = 1 and tz = 2 only. In the case n = 2, the mKdV equation can be transformed to the KdV equation using the Miura transformation [ 31. So for n > 2, there are no exact methods for finding the solution of (8) and (9). The present work then shows that the method of Kath and Smyth [4] can be successfully extended to equations for which there are no exact solutions.
2. Approximate For simplicity, U,(x)
equations we shall take the initial condition
= A sech”” g,
(9) to be
-cc<<<<.
This initial condition is chosen as it has the same sech*/’ want to study the evolution of the initial condition (13) and energy conservation equations for the mKdV equation energy conservation equations for the mKdV equation (8) i(U)
+ & (3u”+’ + U,,) = 0
(13) profile as the solitary wave solution ( 10). We then for the mKdV equation (8). To do this, the mass are used, as in Kath and Smyth [4]. The mass and are (14)
266
N.E Smyth, A.L. Worthy/Wave
Motion 21 (1995) 263-275
and
(15) respectively. It is expected that the initial condition (13) will evolve into a solitary wave (or solitary waves). seek a solution of the mKdV equation of the form
We then
u=uo+u1,
(16)
where uo = a sech21n :,
(17)
5’(t) = V(t).
e=x-5(t),
(18)
The term ua is the time varying solitary wave which evolves from the initial condition (13) and the term ui accounts for the dispersive radiation which is shed as the solitary wave evolves. It will be assumed that the amplitude of the dispersive radiation is much smaller than the solitary wave amplitude, so that Iui 1< UO.Here t(r) is the position of the solitary wave maximum, V(t) is the solitary wave velocity and the amplitude a and width p of the solitary wave are functions of time t. From the initial condition ( 13), we see that
a(0) = A and p(O) = b.
(19)
To determine the evolution of the solitary wave, equations for a, p and V need to be found. Substituting the expansion (16) into the mKdV equation (8), we obtain the following equation
=-a’ sech21n % + -$ + 2(tGi)a
for UI
- n2Vp2) sech2/” $ tanh $
(4 - n2@p’
[3n2ap2 - 2(n + 2)] sech2-t2’n B tanh 8, P P
(20)
on neglecting terms of order u:. The phase and group velocities for linear dispersive waves for the mKdV equation (8) are -k2 and -3k2 respectively, where k is the wavenumber. Therefore, since the solitary wave velocity is positive, radiation is quickly left behind by the evolving solitary wave. It will hence be assumed that ahead of and in the neighbourhood of the solitary wave, u1 rapidly decays to zero. For this to be the case, the solitary wave position must be correct. This can be seen from Eq. (20) for ut when ( 11) holds, so that un is an exact solitary wave. In this case, a homogeneous solution of (20) is
u1 = uox = --
2a
x - vt
sech21n -
nP
P
x - vt
tanh -
P
.
(21)
Adding a small amount of this solution to ua then amounts to just a phase shift in ua, which results in a shift in the position of the solitary wave. If the solitary wave position has been accurately determined however, then ui must rapidly decay to zero ahead of and in the vicinity of the solitary wave. Assuming this to be the case, then integrating Eq. (20) (or the mass conservation equation ( 14)) from x = ,$ to x = co, we obtain
$ap)=3a’+’
- Va - iapp2,
(22)
N.F: Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
where B( p, q) is the Beta function.
M=
s UI
26-l
Let us set
dx,
(23)
so that M is then the mass shed in the radiation. the solitary wave as
Integrating
(17)
from x = -cc
to x = cc gives the mass in
m
u. dx = B
(24)
Hence from (23) and (24)) mass conservation B
becomes
( >-&3)+~=0. ;,;
(25)
The consideration of mass conservation is now completed and we have two equations for the four unknowns a, /?, M and V. To obtain the two equations needed to complete the system, energy conservation is now considered. From the energy conservation equation ( 15) and the expansion ( 16)) we see that the energy density is
(26) Now it has been assumed above that ~1 is zero in the neighbourhood of the solitary wave. Hence the cross term ~0~1 is zero and the contribution to the energy density from the radiation is second order in ui. If this second order term is neglected, then integrating the energy conservation equation ( 15) from n = --00 to x = cc gives ;B
( >;(a’P) ;,;
=0
(27)
on using ( 17). To complete the system of equations describing the evolution of the solitary wave, an expression velocity V must be found. This expression is found from the moment of energy equation
for the
for the mKdV equation (8). Moment of energy, rather than moment of mass, was chosen to calculate the velocity since, as explained above, any contribution to the moment of energy equation from ui is second order. Then integrating the moment of energy equation from x = --oo to x = cc and neglecting these second order contributions gives
!LV= dt
Wn
+ 1)
(n+2)(n+4ja”-
n(n1t4)p-2
on using the energy conservation equation (27). Note that for large times, the relation ( 11) holds and this velocity approaches the solitary wave velocity (12), as required. The system of equations (22), (25), (27) and (29) gives a description of the evolving solitary wave which is in good agreement with numerical solutions of the mKdV equation. However, as in Kath and Smyth [ 41, it is found that the accuracy of the equations can be improved by including an estimate of the energy associated with the shed radiation. In perturbed soliton problems, mass and energy leak away from the sohton due to the perturbed soliton approaching a non-zero constant as x -+ --cc [ 111, resulting in the formation of a so-called
268
N.E Smyth. A.L. Worthy/Wave Motion 21 (1995) 263-275
shelf behind the soliton. It is this leakage of mass and energy that gives rise to the dispersive radiation behind the solitary wave. As in Kath and Smyth [4] for the KdV equation, it can easily be shown that the same effect occurs for the mKdV equation. Let us then suppose that the mass loss from the solitary wave is due to the solitary wave approaching a small non-zero value u, in the shelf, which we will take to occur near x = -L. The shelf is assumed to follow at a fixed distance behind the solitary wave, so that it moves with the solitary wave velocity V. Then integrating the mass conservation equation (14) from x = -L to x = 0;) (i.e. over the solitary wave) gives
u dx = -Vu,
(30)
-L
on neglecting
terms quadratic
and higher in ui. Comparing
(23),
(25) and (30),
we see that (31)
In a similar manner,
integrating
the energy conservation
equation
(15) gives
(32)
on neglecting
cubic and higher terms in ui. Using (31), this becomes
(33) -L
This is an expression for the energy shed from the evolving solitary wave. Combining energy conservation result (27) yields the energy conservation equation B
(;,a>i
(a’p)
=
-$ (T>‘.
it with the first order
(34)
It can be seen that the velocity expression (29) can be zero at non-zero values of a, which can introduce a non-physical singularity into the energy conservation equation (34). The above derivation of the energy loss, based upon the formation of a shelf behind the solitary wave, is strictly speaking valid in the limit of a slowly varying solitary wave. It is therefore reasonable to use the solitary wave velocity (12) for the velocity V in (34), this velocity expression not having the artificial singularity of (29). Furthermore, using a different expression for V in (34) will produce errors for only a finite time. So, since the rate of energy loss is small, the total integrated difference will also be small. The system of equations governing the evolution of the solitary wave is then (22), (25), (29) and (34). These equations are the same as those of Kath and Smyth [4] in the KdV case n = 1. While these equations were found to give results in good agreement with numerical solutions of the KdV equation, Kath and Smyth [ 41 found that better agreement could be obtained by generalising the form of the velocity expression (29) for n = 1 in the following manner V = 2a + ua - 2vpe2,
(35)
where Y is a constant. In the limit of a soliton, so that ap2 = 2, this velocity approaches the correct soliton velocity, independent of the value of v. Kath and Smyth [4] found that the value v = 3/2 gave the best
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
269
agreement between solutions of the approximate conservation equations and full numerical solutions of the KdV equation, independent of the initial values of a and p. The parameter v in the velocity expression (35) is a decay constant controlling the rate at which the solution approaches the steady state, but it does not have a large effect upon the value of this steady state. The same method for determining an improved expression for V will be used in the present work. We therefore set v=
6 -----un n+2
+ zX2”- Y
2(fi + 2) P_2 3n2 ’
(36)
noting that in the limit of a solitary wave, in which the amplitude-width relation (11) holds, this velocity expression is the same as the solitary wave velocity (12). Again, the parameter v will be chosen to give the best agreement with numerical solutions of the mKdV equation. We note that the moment of energy result (29) gives 18n u=
(37)
(n+2)(n+4)‘
While the major role of Y is to determine the rate at which the solution approaches the steady state, it also has an effect on the stability of the approximate conservation equations. To study this, let us ignore the small loss of energy in (34), so that a2P = A2b, on noting (36) is a04-n
(38)
the initial
=
values of a and /?. Then the steady state of the approximate
3 (A’b)’ n2 2 (n+2>
and
PO = ~u;‘d2 &in
equations
(22),
(27)
and
’ (39)
Linearising
about this steady state by setting a = an + 2, we obtain
it’=--&
(s-v)
(n-4)u;ff5i?,
(40)
where E=A2b
and
The approximate
V>-&
p=B
equations
(41) are then stable for (42)
for n < 4 and v<-
3n n+2
(43)
for n > 4. Comparing the value (37) of v given by the moment of energy equation with the stability constraints (42) and (43), we see that the moment of energy velocity expression (29) gives unstable equations for n = 3 and stable equations for all other values of n. Furthermore, the equations for n = 2 and n = 4 are neutrally
270
N.F: Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
t Fig. 1. Solitary wave amplitude a as a function of time t for n = 3, initial amplitude and h = 1.O. Conservation laws: --; numerical solution of mKdV equation: - - -.
(a) A = 0.7 and initial width h = 1.O and (b) A = 0.9
stable as the growth rate in (40) is zero for these values of n. This result for n = 4 is not surprising as it will be found in the next section that for n = 4 the solitary wave solution of the mKdV equation is unstable. However, it is very surprising that for n = 3 moment of energy gives an unstable set of equations. Numerical solutions of the mKdV equation for IZ = 3 show no qualitative differences from numerical solutions for other values of n (n < 4). In the next section, it will be found that by adjusting the value of v in (36), very good agreement can be obtained with numerical solutions of the mKdV equation. The same also holds true for the n = 2 case. For n = 2, it is surprising that moment of energy yields only a neutrally stable set of equations since the mKdV equation for n = 2 can be transformed to the KdV equation by the Miura transformation and moment of energy for the KdV equation yields a stable set of equations.
3. Results and comparison
with numerical
solutions
In this section, solutions of the approximate equations (22), (25) and (34) together with the velocity expression (36) will be compared with full numerical solutions of the mKdV equation (8) with the initial condition ( 13). The approximate equations are nonlinear and cannot be solved exactly in general. For the results presented in this section, these equations are solved numerically using a fourth order Runge-Kutta method. For the numerical integration of the mKdV equation, an adaption of the pseudo-spectral method of Fornberg and Whitham [ 121 was used. This involves using FIT to evaluate the space derivatives and a fourth order Runge-Kutta method to propagate the scheme forward in time. The cases n = 2, n = 3, n = 4 and n = 5 will be dealt with as these give the full range of possible behaviours. The case n = 1 (the KdV equation) was dealt with in Kath and Smyth [4 J. For n = 2, the mKdV equation is exactly integrable, but not for n = 3, n = 4 and n = 5. It will be found that for n < 4, the behaviour of the solution is similar to the n = 1 KdV case. For n 2 4, it will be found that the initial condition decays into dispersive radiation with no solitary wave being formed if the initial amplitude is small enough. For initial amplitude above a critical value, the amplitude of the solitary wave blows-up. Let us consider the case n = 3 first since there is no inverse scattering solution of the mKdV equation with n = 3. The best fit of the solutions of the approximate conservation equations with full numerical solutions was found for the parameter value Y = 2.2 in the velocity expression (36). This value of Y was found to give the best fit independent of the initial values of a and p (A and 6). Fig. la shows a comparison between the
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
271
solitary wave amplitude as given by the solution of the approximate conservation equations and as given by the full numerical solution for an initial amplitude of A = 0.7 and an initial width b = 1. It can be seen that the agreement is very good for the final steady state amplitude. However, the agreement between the transient behaviours is not quite so good, but is adequate. Fig. lb shows the same comparison for A = 0.9 and b = 1. The agreement between the solutions is similar to that in Fig. la, showing that the value of v is independent of the initial condition chosen. Fig. 2 shows a comparison between the final steady amplitudes of the solitary wave produced as given by the approximate conservation equations and by the full numerical solution of the mKdV equation for a range of initial amplitudes and initial width b = 1.O. It can be seen that the comparison is excellent, as can be noted from Figs. la and lb for special cases. When the area of the initial condition (13) is large enough, more than one solitary wave will eventually form. The initial amplitude and width at which this change occurs can be determined from the approximate conservation equations (22)) (25)) (34) and (36). The point at which the change from producing one final solitary wave to producing two final solitary waves occurs is given by the values of A and b for which M’ = 0 initially. This is because if M’ = 0 initially, then there is exactly enough mass in the initial condition to produce one solitary wave and any more mass than this will result in another solitary wave being produced. Setting M’ = 0 in (25)) we see from (22) that a’ = 0 and p’ = 0 and that a and p are related by the solitary wave relation ( 11) . Then using the expression (36) for the solitary wave velocity, we find from (22) that
For a initially greater than this value, 2 solitary waves will form from the initial condition. waves will then form when A” >
g&2
For n = 3, 2 solitary
(45)
This result was found to be in good agreement with numerical solutions, although it was difficult to verify the exact point at which 2 solitary waves first formed due to the formation time for the second solitary wave becoming infinite as the change over point was approached. For b = 1, 2 solitary waves are formed when A > 0.71814. In Fig. 2, the amplitude plotted when A > 0.71814 is the amplitude of the hugest solitary wave. We then see that the approximate conservation equations give good results even when 2 solitary waves are formed. The mKdV equation for n = 2 was also considered, the mKdV equation being related to the KdV equation by the Miura transformation in this case [ 31. The value of v in the velocity expression (36) which gave the best agreement with numerical solutions was found to be v = 2.0. Fig. 3 shows a comparison between the final steady solitary wave amplitude as given by the approximate conservation equations and that given by the numerical solution of the mKdV equation as a function of the initial amplitude A for b = 1. Again it can be seen that the agreement is excellent. Furthermore, the relation (44) shows that 2 solitons will form for A2 > ;b-2,
(46)
which is A > 0.81650 for b = 1. Hence Fig. 3 shows that the approximate conservation equations again give good agreement for the largest soliton even when 2 solitons are formed. The point at which two solitons first form, given by (46), was again found to be in good agreement with numerical solutions. The mKdV equation for n = 2 and n = 3 has similar behaviour to the KdV equation. However, the behaviour for n 2 4 is completely different. This can be seen by solving the approximate conservation equations (22), (25), (34) and (36) when the energy leakage is ignored in Eq. (34). The energy leakage is usually small, so the results obtained by neglecting it are a good approximation to the solution of the full equations. On ignoring the energy leakage, Eq. (34) gives
272
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275 1.6
l
1.3
t
1.4 -
1.2 -
t
.
1.1 -
1.2 _ P
l-
.
1.
a .
.
0.9 _
.
a 0.8 -
.
0.8 -
.
0.7 . .
0.6 -
.
0.6 -
.
0.5 _ '+I.6
0.65 t
0.7
0.75
0.6
0.65
0.9
0.95
1
1.05
t
0.1 0.6
0.65
0.7
03s
A
0.8
0.81
0.9
0.95
I
,.oi
A
Fig. 2. Final steady solitary wave amplitude a as a function of initial amplitude A for n = 3 and initial width b = 1. Conservation laws: 0 ; numerical solution of mKdV equation: +. Fig. 3. Final steady solitary wave amplitude a as a function of initial amplitude A for n = 2 and initial width b = 1. Conservation laws: 0 ; numerical solution of mKdV equation: +. a2P =
E,
(47)
where E is a constant.
Eq. (22) then becomes (48)
Integrating
this equation,
$ (a-4
_ A-“)
we obtain
+ ;E-2
(a-2
_ A-2)
+ 2E-4
log (2a2 - 3E2)A2
= 1
(3
_ v) t,
(49)
(2A2-3.!?)a2 ITE 2
9
('- i?E-2a)A+ ;(a-5 -A-5)+ GE-2 $fEe2A)a
(a-4
-A-4)+
;($,'E-4
(a-3
_ A-3)
(l-
+i(g)3E-6
(a-‘-
A-=)
and
a6=A6
l+ [
3 -A6(2
- v)(4
CLUE
+ (i)4Ee8
(a-’
-I 1
- E-=)t
for the cases n = 2, n = 3 and n = 4, respectively. and
~4 = B
- A-‘)
= s(g
- .>,
(51)
Here
(52)
It can be seen from (49) and (50) that a approaches a constant value as t + cc for n = 2 and n = 3. However, for the case n = 4, (51) gives that a either approaches 0 as t + 00 or a approaches oc in finite time, depending on the sign of (2 - Y) (4 - Ed=). For A small (i.e. E small), numerical results show that a decreases as t increases. We see from (48) that we then require Y > 2 for this to be the case. The best fit with numerical solutions was found for Y = 4.2. Thus from the solution (51) for a it is seen that for A below a
213
N.F: Smyth. A.L. Worthy/Wave Motion 21 (1995) 263-275
t
t
Fig. 4. Solitary wave amplitude ll as a function of time t for n = 4. (a) Initial amplitude A = 0.4 and initial width h = 2.5. Approximate conservation equations: -; numerical solution: - - -. (b) Numerical solution for initial amplitude A = 0.8 and initial width h = I .O.
threshold (E < l/2), a -+ 0 as t + co and for A above this threshold, a + 0;) in finite time. This result is confirmed by full numerical solutions of the mKdV equation for n = 4. Fig. 4a shows the numerical solution of the mKdV equation for n = 4, initial amplitude A = 0.4 and initial width b = 2.5. It can be seen that the solitary wave amplitude is approaching 0 rather than approaching a steady non-zero value. The solitary wave decays to zero and all of its initial mass and energy is shed as dispersive radiation. This is in accord with the predictions of (5 I), for which 4 - Ee2 is negative for these values of A and b. Also shown in Fig. 4a is the solution of the approximate conservation equations. The agreement between the approximate and numerical solutions is reasonable. The agreement is less close than for the previous values of n as in this case the solitary wave is decaying to zero. The assumption of small radiation made in the derivation of the approximate conservation equations is then less valid. In choosing the value Y = 4.2, a compromise was made between agreement with the numerical solution for small and large times. Other values of 1/ will give better agreement at either small times or large times, but not both. The numerical evidence for the blow-up of the solitary wave amplitude for E > l/2 is less clear. This behaviour of the solution for n = 4 is consistent with the instability of the solitary wave solution of the mKdV equation for n = 4 found by Bona et al. [ lo]. However, their stability analysis does not give the threshold behaviour for the evolution of the initial condition as predicted by the approximate conservation equations. It was found numerically that the amplitude did become very large in finite time for E sufficiently large. Furthermore, the numerically determined initial amplitude at which the solitary wave evolution behaviour changed from decay to blow-up agreed very well with the approximate relation E = l/2. For example, taking b = 1, it was found numerically that the value of the initial amplitude A at which the amplitude evolution changed from decay to growth was A = 0.707, which is in excellent agreement with the approximate value from E = I /2, A = fi. Indeed the degree of agreement is surprising given the approximations made in deriving the governing equations (22)) (25), (34) and (36). To obtain the best agreement between the full numerical and approximate solutions, the value of v in the velocity expression (36) was determined by fitting with the numerical solution. However, the value of the threshold (E = l/2) at which the solitary wave evolution changes from decay to blow-up is independent of the value of v, as can be seen from the solution (5 1) . The prediction of this change in behaviour then does not depend on any fit with the full numerical solution. Fig. 4b shows the numerical solitary wave amplitude as a function of time for initial amplitude A = 0.8 and initial width b = 1.O (for which E > l/2). It can be seen that the amplitude is rapidly becoming large. Indeed, the amplitude was found to become beyond the computer capacity for t not much larger than 6. This
214
N.E Smyrh, A.L. Worthy/Wave Motion 21 (1995) 263-275
numerical blow-up was not due to numerical error as it was found to occur even if the space and time steps were decreased. It is of course problematical to say from numerical results that the solitary wave amplitude is becoming infinite since as the amplitude increases, the width decreases, which results in the loss of spatial resolution over the solitary wave. However, there is numerical evidence for the predicted behaviour of the n = 4 case, in particular the blow-up for sufficiently large values of E. As the value of n increases, the evolution of the initial condition (13) becomes slower. Let us now consider the mKdV equation for n 2 5. For low initial values of A, numerical solutions show that the amplitude a must decrease with time t. From (48) we see that this requires U>-
3n n+2’
We then see from the stability condition (43) that the solitary wave is unstable according to the approximate conservation equations, which agrees with the results of Bona et al. [lo]. This is also borne out by numerical solutions of the mKdV equation. These numerical solutions show that the behaviour of the mKdV equation for n > 4 is similar to its behaviour for n = 4 in that for values of the initial amplitude below a critical value, the amplitude of the solitary wave decays to zero as t + co, while for values of the initial amplitude above this critical value, the amplitude approaches infinity as t -+ co. This change in the behaviour of the solution depending upon the value of the initial amplitude can be seen from (48). Since Y satisfies (53), we see that for
a decreases A”
>
as r increases,
while for
2(n + 2) b_2 3n2
’
a increases as t increases. The instability then manifests itself as the solitary wave decaying to zero as c -+ cc when (54) holds and the solitary wave amplitude becoming infinite as t + co when (55) holds. As for the case n = 4 discussed above, the change in the behaviour of the solution at A” = 2(n
+ 2, b-2
3n2 is independent of the value of v and so is independent of any fitting of the value of v from the full numerical solution. To verify these results, we shall consider the case n = 5 in detail. For n = 5, the best agreement with numerical solutions was obtained for Y = 2.46. Fig. 5a shows the numerical solution of the mKdV equation for A = 0.7 and b = 1.0, these values satisfying (54). It can be seen that the solitary wave amplitude is approaching zero, although on a slow time scale. Also shown is the solution of the approximate conservation equations and it can be seen that the agreement between the two solutions is reasonable. As for the n = 4 case, due to the solitary wave decaying to zero, the dispersive radiation is eventually not small and hence good agreement is not to be expected. Fig. 5b shows the numerical solution of the mKdV equation for A = 0.9 and b = 1.O, these values satisfying (55). It can be seen that the solitary wave amplitude is increasing and seems to be blowing-up. There is thus numerical confirmation of the two different types of behaviour depending on the values of A and b. As for the case n = 4 discussed above, the value of the threshold (56) at which the solution behaviour changes from decay to growth can be compared with that given by the full numerical solution of the mKdV equation. When b = 1, the relation (56) gives that the solution behaviour changes when the initial amplitude is A = 0.7148, which is in excellent agreement with the numerical value of A = 0.714. As for the case n = 4, this degree of agreement is surprising given the approximations made in deriving the governing equations (22), (25), (34) and (36), but it does however confirm the basic validity of these assumptions.
N.E Smyth, A.L. Worthy/Wave Motion 21 (1995) 263-275
275
a
t Fig. 5. Solitary wave amplitude a as a function of time t for n = 5. (a) Initial amplitude A = 0.7 and initial width h = 1.0. Approximate conservation equations: --; numerical solution: - - -. (b) Numerical solution for initial amplitude A = 0.9 and initial width h = I .O.
Based on the predictions of the approximate conservation equations, the behaviour of the mKdV equation for II > 5 is expected to be similar to the II = 5 case, which is in accord with the results of Bona et al. [ IO]. While Bona et al. [lo] show that the solitary wave solution of the mKdV equation is unstable for n 2 4, their stability analysis does not predict the ultimate behaviour of the solitary wave due to small perturbations. The existence of a critical threshold for n 2 4 below which the initial condition decays into radiation and above which the initial condition blows-up is the new result predicted by the approximate conservation equations.
Acknowledgement The authors would like to thank Professor W.L. Kath for help with the numerical presented in this paper.
routines
used for the work
References 1 I ] G.L. Lamb, Elemenrs of Soliton Theory, Wiley-Interscience,
New York ( 1980).
[ 21 A.C. Newell, Solitons in Mathematics and Physics, SIAM, Philadelphia ( 1985). [ 3 1 G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York ( 1974). [ 4 1 W.L. Kath and N.F. Smyth, “Soliton evolution and radiation loss for the Korteweg-de Vries equation”, Phys. Rev. E 51 ( 1)) 66 l-670 (1995).
[ 5 1 G.P Agrawal, Nonlinear Fibre Optics, Academic Press, Boston ( 1989). [ 6 ] J.A. Gear and R. Grimshaw, “A second-order theory for solitary waves in shallow fluids”, Phys. Fluids 26(l), 14-29 (1983). [ 7 1 T.B. Benjamin, “The stability of solitary waves”, Proc. R. Sot. London, Sect. A 328, 153-183 ( 1972). [ 8 1 J.L. Bona, “On the stability of solitary waves”, Proc. R. Sot. London. Sect. A 344, 363-374 ( 1975). 19 ] M. Weinstein, “Lyapunov stability of ground states of nonlinear dispersive evolution equations”, Commun. Pure Appl. Mafh. 39, 51-68 (1986). 1 101 J.L. Bona, PE. Souganidis and W.A. Strauss, “Stability and instability of solitary waves of Korteweg-de Vries type”, Pm,. R. Sot. London, Sect. A 411, 395-412 (1987). [ I I] C.J. Knickerbocker and A.C. Newell, “Shelves and the Korteweg-de Vries equation”, J. Fluid Mech. 98, 803-8 I8 ( 1980). [ 12 1 B. Fomberg and G.B. Whitham, “A numerical and theoretical study of certain nonlinear wave phenomena”, Phil. Trans. R. Sac. London A 289, 373-403 ( 1978).