- Email: [email protected]
The inverse scattering transform: Tools for the nonlinear fourier analysis and filtering of ocean surface waves
Chats.
Solitom
& Fracfals Vol. 5, No. 12, pp. 2623-2637, 1995 Elsevier Science Ltd Printed in Great Britain 096&0779/95$9.50 + 0.00
0960-0779(94)E0118-9
The Inverse Scattering Transform: Tools for the Nonlinear Fourier Analysis and Filtering of Ocean Surface Waves A. R. OSBORNE Istituto di Fisica Generale dell’Universita, Via Pietro Giuria 1, 10125Torino, Italy (Received 22 November
1993)
Abstract-Surface wave data from the Adriatic Sea are analysed in the light of new data analysis techniques which may be viewed as a nonlinear generalization of the ordinary Fourier transform. Nonlinear Fourier analysis as applied herein arises from the exact spectral solution to large classes of nonlinear wave equations which are integrable by the inverse scattering transform (IST). Numerical methods are discussed which allow for implementation of the approach as a tool for the time series analysis of oceanic wave data. The case for unidirectional propagation in shallow water, where integrable nonlinear wave motion is governed by the Korteweg-deVries (KdV) equation with periodic/quasi-periodic boundary conditions, is considered. Numerical procedures given herein can be used to compute a nonlinear Fourier representation for a given measured time series. The nonlinear oscillation modes (the IST ‘basis functions’) of KdV obey a linear superposition law, just as do the sine waves of a linear Fourier series, However, the KdV basis functions themselves are highly nonlinear, undergo nonlinear interactions with each other and are distinctly non-sinusoidal. Numerical IST is used to analyse Adriatic Sea data and the concept of nonlinear filtering is applied to improve understanding of the dominant nonlinear interactions in the measured wavetrains.
1. INTRODUCTION
This paper summarizes a new numerical approach for the nonlinear Fourier analysis of space and time series of broadband, nonlinear wavetrains. The method, based upon the (periodic/quasi-periodic) inverse scattering transform (IST), is a kind of nonlinear generalization of the ordinary linear Fourier transform or Fourier series. The focus is on nonlinear wave motion for shallow-water waves as governed by the Korteweg-de Vries (KdV) equation. IST may be exploited to determine the numerical inverse scattering transform (NIST) spectrum of a measured or computed wave train which is assumed to be periodic (or quasi-periodic) in space or in time [l]. The approach may also be applied to numerically construct simple or (spectrally) broadband solutions to the KdV equation [2]. The results given herein are based upon previous applications of the periodic scattering transform to the analysis of experimentally measured data and of computer generated signals [3-121. The scope of this paper is to use NIST to analyse measured wave data obtained in the Adriatic Sea on the offshore platform of the Italian National Research Council in 16.5 m of water, about 10 km from Venice, Italy (see also [9, 131). 2. OCEANIC
SURFACE WAVES
The linear Fourier analysis of ocean surface waves is a well-established procedure which has been successfully applied innumerable times since the invention of modern computing machinery about 50 years ago. The technique manifests itself most commonly in the form of power spectral analysis and, as a result, has been responsible for an enhanced 2623
2624
A. R. OSBORNE
understanding of oceanic wind waves [14-161. Spectral analysis, cross-spectral analysis and bi-spectral analysis, together with many other more recently developed techniques, all based upon the linear Fourier transform, are now standard data analysis procedures [17]. Most of these methods are founded upon the assumption that wind waves are ergodic random processes, and thus the early work by Paley, Weiner and Rice [18] has been of prime importance, as was clearly demonstrated by Pierson and coworkers during the 1950s and 1960s (see ref. [15] for an overview); the methods include the characterization of complex surface wave-power spectra and have been extended to include the forecasting and hindcasting of wind waves from real-time and historically available meteorological data [19]. Serious attempts to deal with nonlinear effects in wind waves began with Tick [14] and others [16] who, within a stochastic framework, developed theoretical foundations for computing second and higher order nonlinear effects in random wave trains. Recently, nonlinear random waves for a completely integrable Hamiltonian system (the KdV equation) have been considered for shallow-water wave motion [20]. A review of all of the above work is beyond the scope of the present paper, but recent important results give a hint of the successesof modern approaches [16, 211. 3. THE KdV EQUATION
AND PERIODIC
INVERSE
SCATTERING
THEORY
The KdV equation is a leading-order, singular-perturbative expansion of the free-surface wave amplitude for the Euler equations in shallow water; the fluid is assumed to be irrotational, incompressible, inviscid, unidirectional and to lie on water of constant depth. KdV, in effect, describes the motion of small, finite-amplitude nonlinear wavetrains in shallow water and was the first of many nonlinear wave equations to be integrated by the IST [22-261. In this sense KdV is a completely integrable Hamiltonian system. The (space-like) KdV equation is given in dimensional form by [27, 281: Ilt + Co% + ~Wx + PL
= 0.
(1) Here ~(x, t) is the wave amplitude as a function of space x and time t. The constant coefficients of KdV are given by co = (gh)““, LX= 3c,/2h and fi = c0h2/6. Equation (1) has the linear dispersion relation (set a = 0 in (1)) w = c,-,k- /?k3; g is the acceleration of gravity, co is the linear phase speed and h is the water depth. The subscripts refer to partial derivatives with respect to x and t. For present purposes we consider the Cauchy problem for KdV, i.e. given the spatial behaviour of the wavetrain at t = 0, q(x, 0), (1) determines the motion for all space and time thereafter, n(x, t). Periodic boundary conditions in space, n(x, t) = n(x + L, t), are assumed throughout; L is the spatial period of the wavetrain. Normally one experimentally measures time series (where the wave amplitude is recorded as a function of time at a single spatial location) rather than space series (where the amplitude is recorded as a function the spatial variable at a fixed time). The reasons for preferring time series over space series are often strictly economical, since the cost of instruments to measure the former (say a single wavestaff or pressure recorder) is substantially less than the cost to measure the latter (which requires remote sensing capability). One is thereby motivated to develop procedures for determining the scattering transform of time series, which notationally has the form ~(0, t). To this end one may apply the time-like KdV equation (TKdV) 123, 291: % + Co’% + ~‘w?t + P’%t = 0, (2) where co’ = l/co, (Y’ = -cx/c,,~ and /3’ = -p/co; (2) has the linearized dispersion relation k = w/co + (#co4)&. TKdV is here viewed as solving a boundary value problem: given the
Fourier analysis of ocean surface waves
2625
temporal evolution, ~(0, t), at a fixed spatial location 1: = 0, (2) specifies the motion over all space as a function of time, n(x) t). The temporal variable is assumed to obey periodic boundary conditions (~(x, t) = q(x) t + T)). The scattering transform approach (SKdV) is thereby consistent with the same periodicity assumption used for linear Fourier algorithms (discrete and fast Fourier transforms). Due to recent advances in numerical methods TKdV may now be routinely applied to the time series analysis of experimental data [9-11, 131. All solutions of SKdV (1) may be easily mapped to all solutions of TKdV (2) by a simple transformation given elsewhere [lo, 141. Hence the scattering transform of (2) may be easily expressed in terms of the scattering transform of (1). For present purposes it is only necessary to note that given the IST for (l), the IST for (2) may be easily determined [lo]. Therefore, only the mathematical development of IST for (1) is given here. The spirit of the data analysis approaches discussed here are based upon the fact that the solution to (1) may be described by a generalization of ordinary Fourier series. The inverse scattering transform solution to the periodic KdV equation (1) may be written as a hear superposition of nonlinearly interacting, nonlinear wave.s (hyperelliptic functions), Pj(x; xO9O): N
~Y(x, t) = -El + 2[2Pj(x; x0, t) ‘- &j - Ezj+l]*
(3)
j=l
The constant parameter 3L= a/6/3. Equation (3) is the first of the so-called trace formulae for the KdV equation [30-331 and may be interpreted as a kind of nonlinear Fourier series. The parameters E,j, Ezj+l are constant eigenvalues of the main spectrum of the periodic inverse scattering transform of (1) (see discussion in the next section); x0 is referred to an arbitrary base point in the periodic interval 0 G x c L. The pj are the nonlinear Fourier (oscillation) modes of the periodic KdV equation and they are thereby analogous to the sine waves of linear Fourier analysis. The ,uj evolve in space by the following system of coupled, nolinear, ordinary differential equations: N
(4) j+k
where 2Nfl R&i)
=
n k=l
(Pj
-
E/c).
(5)
Here the aj = f 1 are the signs of the square root of H(Cli). The ,Ujdynamically evolve on two-sheeted Riemann surfaces, each sheet is associated with its index aj; the branch points connecting the surfaces are referred to as band edges and are denoted by the E,j and E2j+l* The spatially and temporally varying pj evolve inside an open band, e.g. in the interval E,j G ~j s Ezj+l, and oscillate between these eigenvalues as a function of x and t, as shown in Fig. 1. Once a particular /Jj reaches a band edge (either E2j or Ezj+l) the sign aj moves to the other Riemann sheet. The fact that the motion lies on two Riemann sheets, together with the strong, stiff nonlinear coupling which occurs among the pj, presents considerable difficulties in numerical integrations of (4) [6]. These difficulties have been largely circumvented by the methods used herein for the time series analysis of nonlinear wavetrains [l]. The temporal evolution of the pj is governed by the following differential equations: 53 = -2[ilq(x, dt
t) - 2&,
(6)
A. R. OSBORNE;
Solution to KdV Equation _ SolutiontoKdVissUm of Hype&Q& Functims
Hyperelliptic Functions -4 h -5 &3 -6 CLq -7 -8 -9
CL3 p2 Pl
Fig. 1. Synthesis of a wave train solution to the KdV equation. In (a) a solution to the KdV equation is shown. In (b) are graphed the associated six nonlinear hyperelliptic functions as a function of X. The exact linear superposition of these modes gives the solution to KdV shown in (a).
where Ln(x, t) is given by the trace formula (3). The ordinary differential equations (ODES) for space (4) and time (6) evolve the ,ui(x, r) (the nonlinear oscillation modes of KdV) in space and time; and the nonlinear Fourier series (3) constructs general solutions to the KdV equation. In the Appendix are described methods for numerically computing the oscillation modes ~j(X, 0) at a particular instant of time, I = 0 (or alternatively pj(O, t) at a particular spatial point, x = 0). The requisite numerical methods are then referred to as nonlinear Fourier analysis procedures for space or time series [l, 21. Generally speaking reference is made to the numerical determination of the main spectrum (E;; 1 < i c 2N + 1) and the auxiliary spectrum (pj(O, 0), ai = fl; 1 c j c N) as the direct scattering transform (see details in Section 4). The computation of the hyperelliptic functions ~j(X; t) as solutions of the nonlinear ODES (4-6) and the construction of solutions of the KdV equation by the nonlinear Fourier series (3) constitutes the inverse scattering transform. 4. CONVENIENT
BASIS FUNCTIONS INVERSE
FOR NUMERICAL COMPUTATION SCATTERING TRANSFORM
OF THE PERIODIC
The spectral eigenvalue problem for KdV (1) is the second-order Schrodinger problem, normally associated with quantum mechanics:
Fourier analysis of ocean surface waves
Yxx + [Aq(x) + k2]W = 0,
2627
i(k2 = E)
(7)
where the potential function, n(x) = &x, 0), is the solution to the KdV equation (1) at an arbitrary time t = 0; k is the spectral wavenumber and E = k2 is the energy. Periodic boundary conditions are assumed so that we take 71(x, 0) = 7(x + L, 0) for L the period. Formally speaking, given 17(x, 0), (7) is to be solved for Y(x, E) and its properties. Extensive details of periodic inverse scattering theory are not treated here; these may be found elsewhere [30-331. For numerical purposes it is convenient to consider a basis of solutions (c, s) of (7) such that: 4x0) 4x0)
03)
The Wronskian W(c, 8) = 1 and therefore (c, s) is a basis set of (7). The matrix o, referred to as the monodromy matrix, carries the solution of (7) from the point x to x + L: c(x + L) + L)
i s(x
::i::
+’ LL;j = (:;:
2:j(:::1
z:::lj.
(9)
a may be viewed as the fundamental matrix of periodic spectral theory for KdV because it contains all spectral information about KdV in the wavenumber domain. The main spectrum of KdV consists of the Ei (1 6 i < 2N + 1) which are the eigenvalues of the Bloch eigenfunction solutions of the Schrodinger equation (7). The auxiliary spectrum is defined as the eigenvalues for which the eigenfunctions s(x) have the boundary conditions s(xo + L) = $(x0) = 0. The specific spectral definition may be summarized: Main spectrum {Ei; 1 s i s 2N + l}: Auxiliary spectrum {cli; 1 c j s N}: {Oj} = {sgn{au(E)
;(q, ~221@)
+ CLAN) = +1
=0
(10)
- a22(E)]EZflj; I c j s N}.
The eigenvalues { Ei; ~j; oj} constitute the direct scattering transform of a discrete wavetrain of N spatial points (or degrees of freedom), such that 1 G i =Z2N + 1; 1 c j G N. The inverse scattering transform, (3-6), then allows for the construction of broadband (in the spectral sense) wavetrain solutions of the KdV equation. See the Appendix for numerical details.
5. NUMERICAL
EXAMPLE
OF NONLINEAR
FOURIER
ANALYSIS
To illustrate the use of the numerical inverse scattering transform for the synthesis of nonlinear wave trains, Fig. 1 shows the numerical construction of a six degree-of-freedom wavetrain [2]. In panel (a) is the wavetrain itself, q(x) 0)) as a function of x . The nonlinear Fourier decomposition of the wavetrain is shown in Fig. l(b) where the hyperelliptic functions pji, j = l-6 are also shown as a function of x ; in the present case the pj are constructed from equation (4) for a rather arbitrary selection of the eigenvalues E,j, E2j+l. The linear superposition of the six oscillation modes gives the solution to KdV as shown in the panel (a). Note that the hyperelliptic oscillation modes are highly non-sinusoidal in appearance due to nonlinear effects. The shaded regions in Fig. l(b) are the ‘open bands’ of the spectrum which define the range (maxima and minima) of the associated rui. The two lowermost hyperelliptic functions are seen to be so close together that there is effectively no ‘gap’ between them; this is the signature of the soliton in the scattering transform, where in the present case there are two solitons in the spectrum (see [2] for further details).
2628
A. R. OSBORNE 6. NONLINEAR
FOURIER
ANALYSIS
OF MEASURED
ADRIATIC
SEA WAVETRAINS
Data from the Adriatic Sea measurement programme is now analysed (see [9, 13, 341). In particular the linear Fourier and scattering transforms are used to analyse nonlinear wave data obtained on the offshore platform of the Italian National Research Council in 16.5 m of water about 10 km from Venice, Italy, in a region where the bottom slope is rather small, e.g. -l/1000 [35]. A typical measured wavetrain, a 2000 point time series with temporal discretization At = 0.25 s, is shown in Fig. 2(a). In order to insure that the inverse scattering transform is an appropriate descriptor of the wave motion, a data set is considered for which most of the wave energy is in the dominant direction of propagation; only 5% of the wave energy is perpendicular to this direction (see [lo, 131 for discussion).
I
0
I
I
100
I
I
I
1
200 300 Time (set)
0.15 0.10 Frequency (Hz)
I
I
400
I
I
500
0.25
Fig. 2. Time series measured in the Adriatic Sea in 16.5 m water depth (a). In (b) the linear Fourier transform of the measured time series is shown.
Fourier analysis of ocean surface waves
2629
This ensures that the waves are, for practical purposes, unidirectional, a requirement of the KdV equation and consequently of the inverse scattering transform. The significant wave height (average of the highest one-third waves) is H, = 2.6 m and the dominant period is Td = 10.1 s. The linear Fourier spectrum of the measured time series in Fig. 2(a) is shown in Fig. 2(b); the results are quite typical of measured ocean wave spectra, e.g. a central peak (around the dominant period, 7’,,) which decays with a power law spectrum at high frequencies. Before continuing the analysis a brief discussion of how appropriate KdV is for describing a particular measured wavetrain is relevant. It goes without saying that if the physics of the wave motion under study is not, in some sense, close to that of the KdV equation, then the results of an IST analysis are of dubious value. Three important tests for assessing the applicability of KdV for a particular data set are [lo, 13, 361: (1) Determine whether the data lie in the KdV region of the Ursell number diagram [36]. (2) Determine if most of the wave energy lies to the left of a KdV ‘cutoff frequency,’ fKdv = 1.36c0/2~rh, in the frequency domain. (3) Determine if there is little directional spreading in the wave field. For the data analysed herein all three criteria are met rather well. The results of the first test are discussed in detail in [13], e.g. the time-like Ursell number, U = 3g H, Td2/321j.h2 - 1. On this basis the Adriatic Sea waves studied herein may be judged to be mildly nonlinear. The second test is verified in Figs 3 (most of the wave energy lies to the left of E kdv = (a-fkdv)2 in panel (a)) and 4 (most of the energy lies to the left of the frequency fKdv in panel (c)). Since only 5% of the wave energy is normal to the dominant wave direction, the last criterion is also satisfied to good accuracy. The above constraints for the selection of experimental data may in some instances be rather conservative; efforts are underway to test the applicability of present approaches to a wider range of data sets. The nonlinear Fourier analysis of the measured wavetrain in Fig. 2(a) is now given in detail. The Floquet discriminant of the measured wavetrain, shown in Fig. 3(a), consists of a graph of the half-trace of the monodromy matrix (A = (a/r1 + ~~~~22)/2, the first of equations (10)) as a function of frequency squared, E = (nf)“. The fluctuations in A(E) for small E are quite large, so that a logarithmic scale has been used to graph the function outside the vertical range -1 < A < 1 (inside the interval (-1 G A c l)A is instead graphed linearly). Where A(E) crosses f 1 defines the main spectrum eigenvalues Ei (1 G i c 2N + 1). The spectrum divides itself into two widely-separated regions (Fig. 3a) of important nonlinear activity corresponding to solitons (on the left) and radiation components (on the right). The soliton part of the Floquet diagram is not easily visible (it is too dense in the E - f2 domain). Therefore, this part of the spectrum has been expanded in scale in Fig. 3(b). The large oscillations to the left represent the soliton modes in the spectrum. The vertical dotted line represents the reference level [5], which is the level upon which the solitons propagate in physical (configuration) space (see discussion below and Fig. 5(b)). The scattering transform spectrum is given in Fig. 4(a); the hyperelliptic function (nonlinear Fourier) amplitudes are graphed as a function of frequency, just as for the linear Fourier transform in Fig. 2(b). The nonlinear Fourier amplitudes (hyperelliptic function amplitudes), Uj = (EZj+l - E,j)/2n, are the widths of the ‘open bands’ in the Floquet spectrum of Fig. 3(b). The radiation spectrum is shown as a solid curve to the right of the reference level, while the solitons are displayed on the left as vertical arrows. About 8% of the wave energy lies in the soliton part of the spectrum. One may contrast the amplitudes of the nonlinear spectrum in Fig. 4(a) with those of the linear spectrum in Fig. 2(b). Note that the radiation components in the scattering transform spectrum are smaller that those for the linear Fourier spectrum. Physically this occurs because part of the energy has been
2630
A. R. OSBORNE:
15
I I
Solitons
I
-15 -0.1
0.0
-15 -0.015
I -0.010
(a)-I
I I
I/
Open Bands I I I 1 0.1 0.2 (xfJ2”;x~ E=
I 0.4
I 0.5
0.6
0.010
0.015
0.020
I -0.005
‘0.000 0.005; E = (I&-)~(Hxf
Fig. 3. Floquet diagram of the measured Adriatic Sea time series in Fig. 2(a) is shown in (a). The Floquet diagram in (a) has been expanded in the soliton (low-frequency) part of the spectrum in (b) to reveal the presence of the reference level and the soliton Floquet oscillations.
transferred from the radiation spectrum to the soliton spectrum due to the presence of nonlinear effects. The modulus of the hyperelliptic functions, an important indicator of nonlinear behaviour [2], is shown in Fig. 4(b). This parameter indicates just how nonlinear the scattering transform spectral components are at a particular frequency (in complete analogy with the modulus of the cnoidal wave solution to KdV for a single spectral component) [2, 51. The modulus indicates strong nonlinear behaviour for values near 1. According to Fig. 4(b) there are therefore two frequency ranges which are of interest in this regard. The first is at low frequency, indicating the presence of solitons in the spectrum. The second is near the
2631
Fourier analysis of ocean surface waves 0.20 0. 18 0. 16 0. 14 0. 12 0.10 0. 08 0. 06 0.04 0.02 0
0.05
0
0.15
0.10
0.20
0.25
Frequency (Hz)
_ _
1 High Nonlinearity 1 I, , in the Soliton
I Hleh Nonlinearitv e-----------d , in the Radiation ’ S&ctrum 1 f KdV /
I I
!A
Y
Filter Band: _ 0.085 - 0.135 Hz
01 0
'
I
0.05
I
I
0.10
I
(
0.15
.'
I I
0.20
0.25
Frequency (Hz) Fig. 4. The scattering transform spectrum of the measured Adriatic Sea time series in Fig. 2(a). Shown are the solitons (vertical arrows) and the radiation spectrum (continuous curve), In (b) is the nonlinear spectral modulus. Values of the modulus near one indicate strong nonlinear behaviour.
peak of the radiative part of the wave train. Nonlinear interactions are quite strong in these two regions. It is of interest to explore these particular cases using nonlinear filtering. To this end the determination of the hyperelliptic functions (nonlinear oscillation modes) of the data as determined by the technique known as base point iteration are now discussed (see Appendix). The first 128 nonlinear modes are computed and their amplitudes are given in Fig. 4(a). The first 14 of these mode amplitudes, computed as a function of time, are shown in Fig. 5(a). The horizontal lines separate each mode from its neighbour on the vertical scale, which has units of squared frequency (or amplitude, these are the units of
2632
A. R. OSBORNE
-0.008 300 200 Time (xx)
2.0
-1.0 Reference Level
-2.0 0
100
300 200 Time (set)
400
500
Fig. 5. The nonlinear oscillation modes in the soliton part of the spectrum (a). In (b) the solitons are constructed by a linear superposition of these nonlinear modes. The solitons are seen to lie beneath the measured wavetrain and to propagate on a ‘reference level’ which lies below the mean water level.
the horizontal coordinate of the Floquet diagram in Fig. 3 or the amplitude of the measured wavetrain in Fig. 2a). It is easily seen that the nonlinear basis functions are distinctly non-sinusoidal, especially for the larger, long-wave modes. The solitons (there are eight of them) are shown in Fig. 5(b). Two applications of nonlinear filtering using the KdV oscillation modes are now discussed. The first is with regard to the soliton part of the spectrum, the second is with regard to the most nonlinear part of the radiation spectrum. In Fig. 5(a) are shown the hyperelliptic modes in the soliton part of the spectrum. Those modes below the reference level are solitonic in origin and can be easily summed to give the soliton part of the measured wavetrain. The reconstructed soliton train is shown in Fig. 5(b), where the
2633
Fourier analysis of ocean surface waves
original wavetrain is also superposed on the figure. The soliton contribution to the wavetrain consists of a long, low-amplitude signal lying beneath the overlying, narrowbanded measured wavetrain, the latter of which is dominated by the radiation modes. The solitons tend to lie beneath the maxima of the local wavegroups; this topic is discussed in detail elsewhere [13]. It is impossible to stress how important the nonlinear filtering process is to the understanding of the soliton dynamics; the author knows of no other method for extracting them from an arbitrary oceanic wavetrain of the type studied here. The most nonlinear components in the radiation spectrum are studied in Fig. 6. Figure 6(a) shows the most nonlinear hyperelliptic modes centered near the peak of the radiation
Filter
I
0.06 0
I
I
100
Range0.085- 0.135Hz
I
200
I 1 I I Sum of Most Nonlinear
2.0
RadiationModes
0
100
.
200
I
L
I
300 Time (set) I
I
300
I
I I So&on Modes
/
rl
I
I
500
400 1
, (b) I .
400
500
Time(set)
Fig. 6. The twenty hyperelliptic functions in the most nonlinear part of the radiation spectrum (see Figure4(b) and the text) (a). The frequency range is 0.085-0.135 Hz. The linear superposition of the radiation modes in (a) is shown in (b) together with the soliton contribution to the wavetrain.
2634
A. R. OSBORNE
spectrum, where the nonlinear spectral modulus is nearly 1, in the frequency range 0.085-0.135 Hz. This range is shown on the graph of the spectral modulus in Fig. 4(b). It is clear that the nonlinear modes in Fig. 6(a) are clearly not sinusoidal and that phase locking plays an important role in their dynamics ([l, 2, 13, 201). There are at least a couple of differences between the nonlinear filtering procedure applied in the present paper and the usual one for linear Fourier analysis: (1) here the spectral modulus is used to select the most nonlinear (and therefore in some sense the most interesting) parts of the spectrum to study and (2j the filtering process is fully nonlinear and often requires an iterative process [2]. The spectral regions which have a large modulus are used to invert (i.e. to reconstruct) and study the wavetrains as a function of time for the most nonlinear components of the spectrum. The linear superposition of certain of these modes gives the most nonlinear part of the radiation contribution to the wavetrain, which is shown in Fig. 6(b). These wavetrains are highly nonlinear, and are generally not represented by the linear Fourier transform, as is clearly evident from Fig. 6(a). For reference the soliton part of the wavetrain is also shown, superposed on the most nonlinear radiation modes in Fig. 6(b). This figure therefore represents the most nonlinear contributions (solitons plus the radiation modes as seen in the time domain) to the measured wavetrain in Fig. 2(a). 7. CONCLUSIONS
There are several significant results which can be highlighted with regard to the present work. Recent numerical progress allows IST to be applied as a nonlinear time series analysis procedure. The data must be screened appropriately to ensure that they correspond to unidirectional, small-but-finite amplitude, long waves in shallow water. Nonlinear filtering may be used to isolate certain ranges of spectral activity. With regard to the latter point, it is worth noting that in the measured wavetrain (Fig. 2a), the soliton components are obscured by the radiation modes. Thus, while soliton components occur in the nonlinear (IST) spectrum (Fig. 4a), they are not directly visible in real space due to the presence of the radiation modes, which energetically dominate the measured wave motion. While the soliton dynamics are physically significant they are not directly visible by an arbitrary observer of the measured wave train. The inverse scattering transform provides a tool for localizing the solitons and for exploring their dynamics by application of the concept of nonlinear filtering. Consider the spectrum in Fig. 4(a) and think of each frequency component as contributing to the nonlinear Fourier series (3). By deleting the terms corresponding to the radiation modes (those to the right of the reference level), and then summing the remaining hyperelliptic terms for the soliton part of the spectrum, one obtains only the contributions that the solitons make to the measured nonlinear time series. In this way one finds a long, low-amplitude wavetrain, which consists of eight nonlinearly interacting solitons. By using the numerical inverse scattering transform as a tool for nonlinear filtering one easily finds the solitons hidden in a sea of background radiation. An associated physical result is that the solitons tend to be phase locked beneath the maxima of the wave packets [9, 131. Complete theoretical understanding of this latter phenomena is still lacking. A further surprising result is that, while the nonlinear spectral modulus is near one for the solitons (as one might a priori suspect), it is also near one for the dominant part of the radiation spectrum as well (Fig. 4b). A formal interpretation of this effect can be given in terms of high-frequency soliton trains [2, 131. The application of nonlinear filtering to the analysis of the Adriatic Sea data in this highly nonlinear region of the radiation spectrum is therefore intimately related to the construction of nonlinear, narrowband wavetrains, an
Fourier analysis of ocean surface waves
2635
example of which is shown Fig. 6(b). Since the nonlinear modes are clearly not sinusoidal (see Fig. 6a), the effect of nonlinear interactions amongst these densely-spaced components is evidently important. It would therefore be inappropriate to interpret the radiation spectrum as being ‘nearly linear’ (historically a convenient designation) as contrasted to solitons as being ‘strongly nonlinear’. Clearly, large nonlinear interactions in the radiation components are an important topic of future research. Additional important consideration rests with understanding the interactions occurring among the radiation modes and the solitons, particularly with regard to the experimentally observed phase locking phenomena. Acknowledgemenrs-This work was supported in part by the Office of Naval Research of the United States of America (Grant NOOO14-92-J-1330)and by the Progetto Salvaguardia di Venezia de1 Consiglio Nazionale delle Ricerche , Italy.
REFERENCES 1. A. R. Osborne, Automatic algorithm for the numerical inverse scattering transform of the Korteweg-deVries equation, submitted for publication (1993). 2. A. R. Osborne, Numerical construction of nonlinear wavetrain solutions of the periodic Korteweg-deVries equation, Phys. Rev. E. 48, 296-309 (1993). 3. A. R. Bishop, M. G. Forest, D. W. McLaughlin and E. A. Overman II, A quasi-periodic route to chaos in a near integrable PDE, Physica D18,293-312 (1986). 4. A. R. Osborne and L. Bergamasco, The small-amplitude limit of the spectral transform for the periodic Korteweg-deVries equation, Nuovo Cimento B85, 229-243 (1985). 5. A. R. Osborne and L. Bergamasco, The solitons of Zabusky and Kruskal revisited: Perspective in terms of the periodic scattering transform, Physica D18, 26-243 (1986). 6. A. R. Osborne and E. Segre. Numerical solutions of the Korteweg-deVries equation using the periodic scattering transform p-representation, Physica D44, 575-604 (1990). 7. G. Terrones, D. W. McLaughlin, E. A. Overman II and A. Pearlstein, Stability and bifurcation of spatially coherent solutions of the damped driven nonlinear Schroedinger equation, SIAM J. Appl. Math. 50, 791-818 (1990). 8. R. Flesch, M. G. Forest and A. Sinha. Physicu D48, 169-208 (1991). 9. A. R. Osborne, E. Segre, G. Boffetta and L. Cavaleri, Soliton basis states in shallow-water ocean surface waves, Phys. Rev. Lett. 64, 1733-595 (1991). 10. A. R. Osborne, Nonlinear Fourier analysis, in Nonlinear Topics in Ocean Physics, edited by A. R. Osborne, pp. 669-700. North-Holland, Amsterdam (1991). 11. A. R. Osborne, Nonlinear Fourier analysis for the infinit-interval Korteweg-de Vries Equation: I. An algorithm for the direct scattering transform, J. Camp. Phys. 94(2), 284-313 (1991). 12. D. W. McLaughlin, C. M. Schober. Chaotic and homoclinic behavior for numerical discretizations of the nonlinear Schroedinger equation, Physicu D57, 447-465 (1992). 13. A. R. Osborne, M. Serio; L. Bergamasco and L. Caveleii, in preparation (1993). 14. Ocean Wave SDectra. Prentice-Hall. Enelewood Cliffs. N.J. (19611. 15. B. Kinsman, &nd Waves, PrenticelHaE, Englewood kliffs, ‘N.J. (1965). 16. A. Brandt, S. E. Ramberg and M. F. Shlesinger, Nonlinear Dynamics of Ocean Waves. World Scientific, Singapore (1992). 17. J. S. Bendat and A. G. Piersol, Random Data: Analysis and Measurement Procedures. Wiley-Interscience, New York (1971). 18. Nelson Wax, editor, Noise and Stochastic Processes. Dover, New York (1954). 19. L. Cavaleri, WAM, a third-generation wave model, in Nonlinear Topics in Ocean Physics, edited by A. R. Osborne, pp. 925-944. North-Holland, Amsterdam (1991). 20. A. R. Osborne, The behavior of solitons in random function solutions of the periodic Korteweg-deVries equation. Phys. Rev. Leti., to appear (1993). 21. S. Elgar and R. T. Guza, J. Fluid Mech. 167, 1 (1986). 22. V. E. Zakharov, S. V. Manakov, S. P. Novikov and M. P. Pitayevsky, Theory of Solitons. The Method of the Inverse Scattering Problem. Nauka, Moscow (1980), in Russian. 23. M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform. SIAM, Philadelphia (1981). 24. R. K. Dodd, J. E. Eilbeck, J. D. Gibbon and H. C. Morris, Solitons and Nonlinear Wave Equations, Academic Press, London (1982). 25. A. C. Newell, Solitons in Mathematics and Physics, SIAM, Philadelphia (1985). 26. A. Degasperis, Nonlinear wave equations solvable by the spectral transform, in Nonlinear Topics in Ocean Physics, edited by A. R. Osborne. Elsevier, Amsterdam (1991). 27. G. B. Whitham, Linear and Nonlinear Waves. Wiley-Interscience, New York (1974). 28. J. W. Miles, Solitary waves, Ann. Rev. Fluid Mech. 12, 11-43 (1980). 29. V. I. Karpman, Non-Linear Waves in Dispersive Media, Pergamon, Oxford (1975).
A. R. OSBORNE
2636
30. B. A. Dubrovin and S. P. Novikov, Periodic and conditionally periodic analogues of the many-soliton solutions of the Korteweg-de Vries equation, Sov. Phys. JEEP 40, 1058-1063 (1974). 31. B. A. Dubrovin, V. B. Matveev, and S. P. Novikov, Nonlinear equation of Korteweg-deVries type, finite-zone linear operators, and Abelian varieties, Russ. Math. Surv. 31, 59-146 (1976). 32. H. Flaschka, and D. W. McLaughlin, Canonically conjugate variables for KdV and Toda lattice under periodic boundary conditions, Prog. Theor. Phys. 55, 438-456 (1976). 33. H. P. McKean and E. Trubowitz, Hill’s operator and hyperelliptic function theory in the presence of infinitely many branch points, Commun. Pure Appl. Math. 29, 143-226 (1976). 34. A. R. Osborne, The numerical inverse scattering transform: nonlinear Fourier analysis and nonlinear filtering of oceanic surface waves, Proc. Aha Huliko’a Winter Workshop, edited by P. Mtiller and D. Henderson. University of Hawaii Press (1993). 35. L. Cavaleri, Experimental characteristics of wind waves, In Topics in Ocean Physics, edited by A. R. Osborne and P. Malanotte Rizzoli. Elsevier, Amsterdam (1982). 36. A. R. Osborne and M. Petti, Laboratory-generated, shallow water surface waves: Analysis using the periodic, inverse scattering transform. Submitted for publication (1993).
APPENDIX The numerical algorithm
Recent work [l, 21 is briefly summarized. The numerical search for the scattering eigenvalues {E,; pi; uj} suggests the need for computing the derivatives of the matrix Cuijwith respect to the energy E. This is because a Newtonian numerical root-finding algorithm is normally used to determine the eigenvalues. To achieve this goal, a matrix method has been developed for obtaining the evolution of the eigenfunction r# as a function of x and E for a particular wavetrain r](x, 0). The analytical determination of derivatives of the matrix elements with respect to E is a feature of the algorithm. Additional considerations with regard to automation of the computer codes are discussed in detail elsewhere [l]. The spectral equations are easily written: VL = -4* (A.11
The subscripts refer to differentiation with respect to x and E: q(x) = An(x) + E. Writing (A.l) in four-vector notation and using a Taylor series expansion for the solution to the scattering equations (A. 1) one obtains:
64.2)
) where
The elements of H are two-by-two matrices. The matrix 0 has ‘zero for its elements and the other matrices are given by: sin (KAx) --
cos @Ax)
T=
K
-resin @An)
(A.4)
cos (KAx) i
and Ax cos (KAx)
_ Ax sin (rcAx) 2K
T E =dr=
aE
-
Ax cos @Ax) + sin(rcAx) 2
-- sin (KAx)
2K2
2K
2K3
-~ Ax ( KAX) 2K
for K = (q)@ = (An(x) + E)@. While K may be either real or imaginary, the matrix TE is always real with determinant 1. This property is exploited in the numerical algorithm below. As in previous numerical problems of this type I assume the wavetrain q(x) has the form of a piecewise constant function with M partitions on the periodic interval (0, L), where the discretization interval is Ax = L/M [lo]. which is associated with a discrete value of the spatial Each partition has wave amplitude n, (0 G n c M-l) variable x, = nAx. It is important to note that the number of degrees of freedom of measured time series of M points is given by the simple relation, N = M/2 [4, 51. The four-by-four scattering matrix M can then be defined:
Fourier analysis of ocean surface waves M =
2637
64.6)
fl H(qnr Ax) n=M-1
The initial conditions of the basis (c, s) at the base point x0 are given by:
(A.7)
From the definition of the matrix (Yone finds:
(A.8) Thus at x = x0 ;(a%,
+
a2222) = (yzl
=
#fll
+
(A.9)
M22)
(A.lO)
Ml2
while the derivatives are given by (A.ll)
Because K = (kn(x, 0) + k2)@ can be either real or imaginary, but Implementation of the numerical algorithm. not complex, the matrix H is always real. This result allows implementation of an algorithm which is entirely real. The following relations have been used in the computer code: if K2 cos (K’Ax) (A.13) TI1 = T12 = fir K2 < 0; cash (K’Ax)
20
(A. 14)
T21
=
-K’ sin (K’Ax) sinh (K’Ax)
K’
if rc23 0 if K2 < 0,
(A.15)
where K’
=
VI/1117 +
k21 =
v/IK’I
(A.16)
and analogously for the matrix Ts. The reconstruction of complex solutions of the KdV equation by (3) ( as well as nonlinear filtering) are accomplished by computing the auxiliary spectra pj(xc = x,) for the M different base points x0 =x0, xi, The approach is formally called base-point iteration and is carried out by computing M different x2,. . .X&f-l. monodromy matrices (A.6) which differ from each other by a horizontal shift Ax in the wave train q,. This procedure arises from the following similarity transformation which arises from (A.6): (A. 17) M(x,+1, E) = H(%, E)M(x,, E)H(h E)-1 The latter expression relates the matrix M(x,+i, E) at a base point x,+i to the previously computed matrix M(x,, E) at the base point x, for a particular value of squared wavenumber E = k:. Values of the auxiliary spectra {p,(x,J} for each x, are computed from the matrices M(x”, E). Knowledge of the auxiliary spectra at every point x, allows reconstruction of the wave train n(x,,) via a discrete version of (3):
Mxn) = -El + :fWj(xrt) - E2j - E2,+11 ,=l
(A.18)
for n = 0, 1, 2,. . . M - 1. This is a finite-term nonlinear generalization of a Fourier series for the discrete wavetrain q(xn). As indicated by the notation, each nonlinear oscillation mode {nj} implicitly depends upon the associated wavenumber ki of the mode. The k, are theoretically given by the simple relation k, = jAk, Ak = 27r/L; surprisingly this is exactly the same as that for the linear Fourier transform, provided that periodicity is assumed. The IST spectrum then consists of the widths of the open bands of the Floquet discriminant, ai = (E2j+l - E2j),/21, graphed as a function of kj (or fj = jAf, Af = l/ir, for T the period of a time series).
Solitom
& Fracfals Vol. 5, No. 12, pp. 2623-2637, 1995 Elsevier Science Ltd Printed in Great Britain 096&0779/95$9.50 + 0.00
0960-0779(94)E0118-9
The Inverse Scattering Transform: Tools for the Nonlinear Fourier Analysis and Filtering of Ocean Surface Waves A. R. OSBORNE Istituto di Fisica Generale dell’Universita, Via Pietro Giuria 1, 10125Torino, Italy (Received 22 November
1993)
Abstract-Surface wave data from the Adriatic Sea are analysed in the light of new data analysis techniques which may be viewed as a nonlinear generalization of the ordinary Fourier transform. Nonlinear Fourier analysis as applied herein arises from the exact spectral solution to large classes of nonlinear wave equations which are integrable by the inverse scattering transform (IST). Numerical methods are discussed which allow for implementation of the approach as a tool for the time series analysis of oceanic wave data. The case for unidirectional propagation in shallow water, where integrable nonlinear wave motion is governed by the Korteweg-deVries (KdV) equation with periodic/quasi-periodic boundary conditions, is considered. Numerical procedures given herein can be used to compute a nonlinear Fourier representation for a given measured time series. The nonlinear oscillation modes (the IST ‘basis functions’) of KdV obey a linear superposition law, just as do the sine waves of a linear Fourier series, However, the KdV basis functions themselves are highly nonlinear, undergo nonlinear interactions with each other and are distinctly non-sinusoidal. Numerical IST is used to analyse Adriatic Sea data and the concept of nonlinear filtering is applied to improve understanding of the dominant nonlinear interactions in the measured wavetrains.
1. INTRODUCTION
This paper summarizes a new numerical approach for the nonlinear Fourier analysis of space and time series of broadband, nonlinear wavetrains. The method, based upon the (periodic/quasi-periodic) inverse scattering transform (IST), is a kind of nonlinear generalization of the ordinary linear Fourier transform or Fourier series. The focus is on nonlinear wave motion for shallow-water waves as governed by the Korteweg-de Vries (KdV) equation. IST may be exploited to determine the numerical inverse scattering transform (NIST) spectrum of a measured or computed wave train which is assumed to be periodic (or quasi-periodic) in space or in time [l]. The approach may also be applied to numerically construct simple or (spectrally) broadband solutions to the KdV equation [2]. The results given herein are based upon previous applications of the periodic scattering transform to the analysis of experimentally measured data and of computer generated signals [3-121. The scope of this paper is to use NIST to analyse measured wave data obtained in the Adriatic Sea on the offshore platform of the Italian National Research Council in 16.5 m of water, about 10 km from Venice, Italy (see also [9, 131). 2. OCEANIC
SURFACE WAVES
The linear Fourier analysis of ocean surface waves is a well-established procedure which has been successfully applied innumerable times since the invention of modern computing machinery about 50 years ago. The technique manifests itself most commonly in the form of power spectral analysis and, as a result, has been responsible for an enhanced 2623
2624
A. R. OSBORNE
understanding of oceanic wind waves [14-161. Spectral analysis, cross-spectral analysis and bi-spectral analysis, together with many other more recently developed techniques, all based upon the linear Fourier transform, are now standard data analysis procedures [17]. Most of these methods are founded upon the assumption that wind waves are ergodic random processes, and thus the early work by Paley, Weiner and Rice [18] has been of prime importance, as was clearly demonstrated by Pierson and coworkers during the 1950s and 1960s (see ref. [15] for an overview); the methods include the characterization of complex surface wave-power spectra and have been extended to include the forecasting and hindcasting of wind waves from real-time and historically available meteorological data [19]. Serious attempts to deal with nonlinear effects in wind waves began with Tick [14] and others [16] who, within a stochastic framework, developed theoretical foundations for computing second and higher order nonlinear effects in random wave trains. Recently, nonlinear random waves for a completely integrable Hamiltonian system (the KdV equation) have been considered for shallow-water wave motion [20]. A review of all of the above work is beyond the scope of the present paper, but recent important results give a hint of the successesof modern approaches [16, 211. 3. THE KdV EQUATION
AND PERIODIC
INVERSE
SCATTERING
THEORY
The KdV equation is a leading-order, singular-perturbative expansion of the free-surface wave amplitude for the Euler equations in shallow water; the fluid is assumed to be irrotational, incompressible, inviscid, unidirectional and to lie on water of constant depth. KdV, in effect, describes the motion of small, finite-amplitude nonlinear wavetrains in shallow water and was the first of many nonlinear wave equations to be integrated by the IST [22-261. In this sense KdV is a completely integrable Hamiltonian system. The (space-like) KdV equation is given in dimensional form by [27, 281: Ilt + Co% + ~Wx + PL
= 0.
(1) Here ~(x, t) is the wave amplitude as a function of space x and time t. The constant coefficients of KdV are given by co = (gh)““, LX= 3c,/2h and fi = c0h2/6. Equation (1) has the linear dispersion relation (set a = 0 in (1)) w = c,-,k- /?k3; g is the acceleration of gravity, co is the linear phase speed and h is the water depth. The subscripts refer to partial derivatives with respect to x and t. For present purposes we consider the Cauchy problem for KdV, i.e. given the spatial behaviour of the wavetrain at t = 0, q(x, 0), (1) determines the motion for all space and time thereafter, n(x, t). Periodic boundary conditions in space, n(x, t) = n(x + L, t), are assumed throughout; L is the spatial period of the wavetrain. Normally one experimentally measures time series (where the wave amplitude is recorded as a function of time at a single spatial location) rather than space series (where the amplitude is recorded as a function the spatial variable at a fixed time). The reasons for preferring time series over space series are often strictly economical, since the cost of instruments to measure the former (say a single wavestaff or pressure recorder) is substantially less than the cost to measure the latter (which requires remote sensing capability). One is thereby motivated to develop procedures for determining the scattering transform of time series, which notationally has the form ~(0, t). To this end one may apply the time-like KdV equation (TKdV) 123, 291: % + Co’% + ~‘w?t + P’%t = 0, (2) where co’ = l/co, (Y’ = -cx/c,,~ and /3’ = -p/co; (2) has the linearized dispersion relation k = w/co + (#co4)&. TKdV is here viewed as solving a boundary value problem: given the
Fourier analysis of ocean surface waves
2625
temporal evolution, ~(0, t), at a fixed spatial location 1: = 0, (2) specifies the motion over all space as a function of time, n(x) t). The temporal variable is assumed to obey periodic boundary conditions (~(x, t) = q(x) t + T)). The scattering transform approach (SKdV) is thereby consistent with the same periodicity assumption used for linear Fourier algorithms (discrete and fast Fourier transforms). Due to recent advances in numerical methods TKdV may now be routinely applied to the time series analysis of experimental data [9-11, 131. All solutions of SKdV (1) may be easily mapped to all solutions of TKdV (2) by a simple transformation given elsewhere [lo, 141. Hence the scattering transform of (2) may be easily expressed in terms of the scattering transform of (1). For present purposes it is only necessary to note that given the IST for (l), the IST for (2) may be easily determined [lo]. Therefore, only the mathematical development of IST for (1) is given here. The spirit of the data analysis approaches discussed here are based upon the fact that the solution to (1) may be described by a generalization of ordinary Fourier series. The inverse scattering transform solution to the periodic KdV equation (1) may be written as a hear superposition of nonlinearly interacting, nonlinear wave.s (hyperelliptic functions), Pj(x; xO9O): N
~Y(x, t) = -El + 2[2Pj(x; x0, t) ‘- &j - Ezj+l]*
(3)
j=l
The constant parameter 3L= a/6/3. Equation (3) is the first of the so-called trace formulae for the KdV equation [30-331 and may be interpreted as a kind of nonlinear Fourier series. The parameters E,j, Ezj+l are constant eigenvalues of the main spectrum of the periodic inverse scattering transform of (1) (see discussion in the next section); x0 is referred to an arbitrary base point in the periodic interval 0 G x c L. The pj are the nonlinear Fourier (oscillation) modes of the periodic KdV equation and they are thereby analogous to the sine waves of linear Fourier analysis. The ,uj evolve in space by the following system of coupled, nolinear, ordinary differential equations: N
(4) j+k
where 2Nfl R&i)
=
n k=l
(Pj
-
E/c).
(5)
Here the aj = f 1 are the signs of the square root of H(Cli). The ,Ujdynamically evolve on two-sheeted Riemann surfaces, each sheet is associated with its index aj; the branch points connecting the surfaces are referred to as band edges and are denoted by the E,j and E2j+l* The spatially and temporally varying pj evolve inside an open band, e.g. in the interval E,j G ~j s Ezj+l, and oscillate between these eigenvalues as a function of x and t, as shown in Fig. 1. Once a particular /Jj reaches a band edge (either E2j or Ezj+l) the sign aj moves to the other Riemann sheet. The fact that the motion lies on two Riemann sheets, together with the strong, stiff nonlinear coupling which occurs among the pj, presents considerable difficulties in numerical integrations of (4) [6]. These difficulties have been largely circumvented by the methods used herein for the time series analysis of nonlinear wavetrains [l]. The temporal evolution of the pj is governed by the following differential equations: 53 = -2[ilq(x, dt
t) - 2&,
(6)
A. R. OSBORNE;
Solution to KdV Equation _ SolutiontoKdVissUm of Hype&Q& Functims
Hyperelliptic Functions -4 h -5 &3 -6 CLq -7 -8 -9
CL3 p2 Pl
Fig. 1. Synthesis of a wave train solution to the KdV equation. In (a) a solution to the KdV equation is shown. In (b) are graphed the associated six nonlinear hyperelliptic functions as a function of X. The exact linear superposition of these modes gives the solution to KdV shown in (a).
where Ln(x, t) is given by the trace formula (3). The ordinary differential equations (ODES) for space (4) and time (6) evolve the ,ui(x, r) (the nonlinear oscillation modes of KdV) in space and time; and the nonlinear Fourier series (3) constructs general solutions to the KdV equation. In the Appendix are described methods for numerically computing the oscillation modes ~j(X, 0) at a particular instant of time, I = 0 (or alternatively pj(O, t) at a particular spatial point, x = 0). The requisite numerical methods are then referred to as nonlinear Fourier analysis procedures for space or time series [l, 21. Generally speaking reference is made to the numerical determination of the main spectrum (E;; 1 < i c 2N + 1) and the auxiliary spectrum (pj(O, 0), ai = fl; 1 c j c N) as the direct scattering transform (see details in Section 4). The computation of the hyperelliptic functions ~j(X; t) as solutions of the nonlinear ODES (4-6) and the construction of solutions of the KdV equation by the nonlinear Fourier series (3) constitutes the inverse scattering transform. 4. CONVENIENT
BASIS FUNCTIONS INVERSE
FOR NUMERICAL COMPUTATION SCATTERING TRANSFORM
OF THE PERIODIC
The spectral eigenvalue problem for KdV (1) is the second-order Schrodinger problem, normally associated with quantum mechanics:
Fourier analysis of ocean surface waves
Yxx + [Aq(x) + k2]W = 0,
2627
i(k2 = E)
(7)
where the potential function, n(x) = &x, 0), is the solution to the KdV equation (1) at an arbitrary time t = 0; k is the spectral wavenumber and E = k2 is the energy. Periodic boundary conditions are assumed so that we take 71(x, 0) = 7(x + L, 0) for L the period. Formally speaking, given 17(x, 0), (7) is to be solved for Y(x, E) and its properties. Extensive details of periodic inverse scattering theory are not treated here; these may be found elsewhere [30-331. For numerical purposes it is convenient to consider a basis of solutions (c, s) of (7) such that: 4x0) 4x0)
03)
The Wronskian W(c, 8) = 1 and therefore (c, s) is a basis set of (7). The matrix o, referred to as the monodromy matrix, carries the solution of (7) from the point x to x + L: c(x + L) + L)
i s(x
::i::
+’ LL;j = (:;:
2:j(:::1
z:::lj.
(9)
a may be viewed as the fundamental matrix of periodic spectral theory for KdV because it contains all spectral information about KdV in the wavenumber domain. The main spectrum of KdV consists of the Ei (1 6 i < 2N + 1) which are the eigenvalues of the Bloch eigenfunction solutions of the Schrodinger equation (7). The auxiliary spectrum is defined as the eigenvalues for which the eigenfunctions s(x) have the boundary conditions s(xo + L) = $(x0) = 0. The specific spectral definition may be summarized: Main spectrum {Ei; 1 s i s 2N + l}: Auxiliary spectrum {cli; 1 c j s N}: {Oj} = {sgn{au(E)
;(q, ~221@)
+ CLAN) = +1
=0
(10)
- a22(E)]EZflj; I c j s N}.
The eigenvalues { Ei; ~j; oj} constitute the direct scattering transform of a discrete wavetrain of N spatial points (or degrees of freedom), such that 1 G i =Z2N + 1; 1 c j G N. The inverse scattering transform, (3-6), then allows for the construction of broadband (in the spectral sense) wavetrain solutions of the KdV equation. See the Appendix for numerical details.
5. NUMERICAL
EXAMPLE
OF NONLINEAR
FOURIER
ANALYSIS
To illustrate the use of the numerical inverse scattering transform for the synthesis of nonlinear wave trains, Fig. 1 shows the numerical construction of a six degree-of-freedom wavetrain [2]. In panel (a) is the wavetrain itself, q(x) 0)) as a function of x . The nonlinear Fourier decomposition of the wavetrain is shown in Fig. l(b) where the hyperelliptic functions pji, j = l-6 are also shown as a function of x ; in the present case the pj are constructed from equation (4) for a rather arbitrary selection of the eigenvalues E,j, E2j+l. The linear superposition of the six oscillation modes gives the solution to KdV as shown in the panel (a). Note that the hyperelliptic oscillation modes are highly non-sinusoidal in appearance due to nonlinear effects. The shaded regions in Fig. l(b) are the ‘open bands’ of the spectrum which define the range (maxima and minima) of the associated rui. The two lowermost hyperelliptic functions are seen to be so close together that there is effectively no ‘gap’ between them; this is the signature of the soliton in the scattering transform, where in the present case there are two solitons in the spectrum (see [2] for further details).
2628
A. R. OSBORNE 6. NONLINEAR
FOURIER
ANALYSIS
OF MEASURED
ADRIATIC
SEA WAVETRAINS
Data from the Adriatic Sea measurement programme is now analysed (see [9, 13, 341). In particular the linear Fourier and scattering transforms are used to analyse nonlinear wave data obtained on the offshore platform of the Italian National Research Council in 16.5 m of water about 10 km from Venice, Italy, in a region where the bottom slope is rather small, e.g. -l/1000 [35]. A typical measured wavetrain, a 2000 point time series with temporal discretization At = 0.25 s, is shown in Fig. 2(a). In order to insure that the inverse scattering transform is an appropriate descriptor of the wave motion, a data set is considered for which most of the wave energy is in the dominant direction of propagation; only 5% of the wave energy is perpendicular to this direction (see [lo, 131 for discussion).
I
0
I
I
100
I
I
I
1
200 300 Time (set)
0.15 0.10 Frequency (Hz)
I
I
400
I
I
500
0.25
Fig. 2. Time series measured in the Adriatic Sea in 16.5 m water depth (a). In (b) the linear Fourier transform of the measured time series is shown.
Fourier analysis of ocean surface waves
2629
This ensures that the waves are, for practical purposes, unidirectional, a requirement of the KdV equation and consequently of the inverse scattering transform. The significant wave height (average of the highest one-third waves) is H, = 2.6 m and the dominant period is Td = 10.1 s. The linear Fourier spectrum of the measured time series in Fig. 2(a) is shown in Fig. 2(b); the results are quite typical of measured ocean wave spectra, e.g. a central peak (around the dominant period, 7’,,) which decays with a power law spectrum at high frequencies. Before continuing the analysis a brief discussion of how appropriate KdV is for describing a particular measured wavetrain is relevant. It goes without saying that if the physics of the wave motion under study is not, in some sense, close to that of the KdV equation, then the results of an IST analysis are of dubious value. Three important tests for assessing the applicability of KdV for a particular data set are [lo, 13, 361: (1) Determine whether the data lie in the KdV region of the Ursell number diagram [36]. (2) Determine if most of the wave energy lies to the left of a KdV ‘cutoff frequency,’ fKdv = 1.36c0/2~rh, in the frequency domain. (3) Determine if there is little directional spreading in the wave field. For the data analysed herein all three criteria are met rather well. The results of the first test are discussed in detail in [13], e.g. the time-like Ursell number, U = 3g H, Td2/321j.h2 - 1. On this basis the Adriatic Sea waves studied herein may be judged to be mildly nonlinear. The second test is verified in Figs 3 (most of the wave energy lies to the left of E kdv = (a-fkdv)2 in panel (a)) and 4 (most of the energy lies to the left of the frequency fKdv in panel (c)). Since only 5% of the wave energy is normal to the dominant wave direction, the last criterion is also satisfied to good accuracy. The above constraints for the selection of experimental data may in some instances be rather conservative; efforts are underway to test the applicability of present approaches to a wider range of data sets. The nonlinear Fourier analysis of the measured wavetrain in Fig. 2(a) is now given in detail. The Floquet discriminant of the measured wavetrain, shown in Fig. 3(a), consists of a graph of the half-trace of the monodromy matrix (A = (a/r1 + ~~~~22)/2, the first of equations (10)) as a function of frequency squared, E = (nf)“. The fluctuations in A(E) for small E are quite large, so that a logarithmic scale has been used to graph the function outside the vertical range -1 < A < 1 (inside the interval (-1 G A c l)A is instead graphed linearly). Where A(E) crosses f 1 defines the main spectrum eigenvalues Ei (1 G i c 2N + 1). The spectrum divides itself into two widely-separated regions (Fig. 3a) of important nonlinear activity corresponding to solitons (on the left) and radiation components (on the right). The soliton part of the Floquet diagram is not easily visible (it is too dense in the E - f2 domain). Therefore, this part of the spectrum has been expanded in scale in Fig. 3(b). The large oscillations to the left represent the soliton modes in the spectrum. The vertical dotted line represents the reference level [5], which is the level upon which the solitons propagate in physical (configuration) space (see discussion below and Fig. 5(b)). The scattering transform spectrum is given in Fig. 4(a); the hyperelliptic function (nonlinear Fourier) amplitudes are graphed as a function of frequency, just as for the linear Fourier transform in Fig. 2(b). The nonlinear Fourier amplitudes (hyperelliptic function amplitudes), Uj = (EZj+l - E,j)/2n, are the widths of the ‘open bands’ in the Floquet spectrum of Fig. 3(b). The radiation spectrum is shown as a solid curve to the right of the reference level, while the solitons are displayed on the left as vertical arrows. About 8% of the wave energy lies in the soliton part of the spectrum. One may contrast the amplitudes of the nonlinear spectrum in Fig. 4(a) with those of the linear spectrum in Fig. 2(b). Note that the radiation components in the scattering transform spectrum are smaller that those for the linear Fourier spectrum. Physically this occurs because part of the energy has been
2630
A. R. OSBORNE:
15
I I
Solitons
I
-15 -0.1
0.0
-15 -0.015
I -0.010
(a)-I
I I
I/
Open Bands I I I 1 0.1 0.2 (xfJ2”;x~ E=
I 0.4
I 0.5
0.6
0.010
0.015
0.020
I -0.005
‘0.000 0.005; E = (I&-)~(Hxf
Fig. 3. Floquet diagram of the measured Adriatic Sea time series in Fig. 2(a) is shown in (a). The Floquet diagram in (a) has been expanded in the soliton (low-frequency) part of the spectrum in (b) to reveal the presence of the reference level and the soliton Floquet oscillations.
transferred from the radiation spectrum to the soliton spectrum due to the presence of nonlinear effects. The modulus of the hyperelliptic functions, an important indicator of nonlinear behaviour [2], is shown in Fig. 4(b). This parameter indicates just how nonlinear the scattering transform spectral components are at a particular frequency (in complete analogy with the modulus of the cnoidal wave solution to KdV for a single spectral component) [2, 51. The modulus indicates strong nonlinear behaviour for values near 1. According to Fig. 4(b) there are therefore two frequency ranges which are of interest in this regard. The first is at low frequency, indicating the presence of solitons in the spectrum. The second is near the
2631
Fourier analysis of ocean surface waves 0.20 0. 18 0. 16 0. 14 0. 12 0.10 0. 08 0. 06 0.04 0.02 0
0.05
0
0.15
0.10
0.20
0.25
Frequency (Hz)
_ _
1 High Nonlinearity 1 I, , in the Soliton
I Hleh Nonlinearitv e-----------d , in the Radiation ’ S&ctrum 1 f KdV /
I I
!A
Y
Filter Band: _ 0.085 - 0.135 Hz
01 0
'
I
0.05
I
I
0.10
I
(
0.15
.'
I I
0.20
0.25
Frequency (Hz) Fig. 4. The scattering transform spectrum of the measured Adriatic Sea time series in Fig. 2(a). Shown are the solitons (vertical arrows) and the radiation spectrum (continuous curve), In (b) is the nonlinear spectral modulus. Values of the modulus near one indicate strong nonlinear behaviour.
peak of the radiative part of the wave train. Nonlinear interactions are quite strong in these two regions. It is of interest to explore these particular cases using nonlinear filtering. To this end the determination of the hyperelliptic functions (nonlinear oscillation modes) of the data as determined by the technique known as base point iteration are now discussed (see Appendix). The first 128 nonlinear modes are computed and their amplitudes are given in Fig. 4(a). The first 14 of these mode amplitudes, computed as a function of time, are shown in Fig. 5(a). The horizontal lines separate each mode from its neighbour on the vertical scale, which has units of squared frequency (or amplitude, these are the units of
2632
A. R. OSBORNE
-0.008 300 200 Time (xx)
2.0
-1.0 Reference Level
-2.0 0
100
300 200 Time (set)
400
500
Fig. 5. The nonlinear oscillation modes in the soliton part of the spectrum (a). In (b) the solitons are constructed by a linear superposition of these nonlinear modes. The solitons are seen to lie beneath the measured wavetrain and to propagate on a ‘reference level’ which lies below the mean water level.
the horizontal coordinate of the Floquet diagram in Fig. 3 or the amplitude of the measured wavetrain in Fig. 2a). It is easily seen that the nonlinear basis functions are distinctly non-sinusoidal, especially for the larger, long-wave modes. The solitons (there are eight of them) are shown in Fig. 5(b). Two applications of nonlinear filtering using the KdV oscillation modes are now discussed. The first is with regard to the soliton part of the spectrum, the second is with regard to the most nonlinear part of the radiation spectrum. In Fig. 5(a) are shown the hyperelliptic modes in the soliton part of the spectrum. Those modes below the reference level are solitonic in origin and can be easily summed to give the soliton part of the measured wavetrain. The reconstructed soliton train is shown in Fig. 5(b), where the
2633
Fourier analysis of ocean surface waves
original wavetrain is also superposed on the figure. The soliton contribution to the wavetrain consists of a long, low-amplitude signal lying beneath the overlying, narrowbanded measured wavetrain, the latter of which is dominated by the radiation modes. The solitons tend to lie beneath the maxima of the local wavegroups; this topic is discussed in detail elsewhere [13]. It is impossible to stress how important the nonlinear filtering process is to the understanding of the soliton dynamics; the author knows of no other method for extracting them from an arbitrary oceanic wavetrain of the type studied here. The most nonlinear components in the radiation spectrum are studied in Fig. 6. Figure 6(a) shows the most nonlinear hyperelliptic modes centered near the peak of the radiation
Filter
I
0.06 0
I
I
100
Range0.085- 0.135Hz
I
200
I 1 I I Sum of Most Nonlinear
2.0
RadiationModes
0
100
.
200
I
L
I
300 Time (set) I
I
300
I
I I So&on Modes
/
rl
I
I
500
400 1
, (b) I .
400
500
Time(set)
Fig. 6. The twenty hyperelliptic functions in the most nonlinear part of the radiation spectrum (see Figure4(b) and the text) (a). The frequency range is 0.085-0.135 Hz. The linear superposition of the radiation modes in (a) is shown in (b) together with the soliton contribution to the wavetrain.
2634
A. R. OSBORNE
spectrum, where the nonlinear spectral modulus is nearly 1, in the frequency range 0.085-0.135 Hz. This range is shown on the graph of the spectral modulus in Fig. 4(b). It is clear that the nonlinear modes in Fig. 6(a) are clearly not sinusoidal and that phase locking plays an important role in their dynamics ([l, 2, 13, 201). There are at least a couple of differences between the nonlinear filtering procedure applied in the present paper and the usual one for linear Fourier analysis: (1) here the spectral modulus is used to select the most nonlinear (and therefore in some sense the most interesting) parts of the spectrum to study and (2j the filtering process is fully nonlinear and often requires an iterative process [2]. The spectral regions which have a large modulus are used to invert (i.e. to reconstruct) and study the wavetrains as a function of time for the most nonlinear components of the spectrum. The linear superposition of certain of these modes gives the most nonlinear part of the radiation contribution to the wavetrain, which is shown in Fig. 6(b). These wavetrains are highly nonlinear, and are generally not represented by the linear Fourier transform, as is clearly evident from Fig. 6(a). For reference the soliton part of the wavetrain is also shown, superposed on the most nonlinear radiation modes in Fig. 6(b). This figure therefore represents the most nonlinear contributions (solitons plus the radiation modes as seen in the time domain) to the measured wavetrain in Fig. 2(a). 7. CONCLUSIONS
There are several significant results which can be highlighted with regard to the present work. Recent numerical progress allows IST to be applied as a nonlinear time series analysis procedure. The data must be screened appropriately to ensure that they correspond to unidirectional, small-but-finite amplitude, long waves in shallow water. Nonlinear filtering may be used to isolate certain ranges of spectral activity. With regard to the latter point, it is worth noting that in the measured wavetrain (Fig. 2a), the soliton components are obscured by the radiation modes. Thus, while soliton components occur in the nonlinear (IST) spectrum (Fig. 4a), they are not directly visible in real space due to the presence of the radiation modes, which energetically dominate the measured wave motion. While the soliton dynamics are physically significant they are not directly visible by an arbitrary observer of the measured wave train. The inverse scattering transform provides a tool for localizing the solitons and for exploring their dynamics by application of the concept of nonlinear filtering. Consider the spectrum in Fig. 4(a) and think of each frequency component as contributing to the nonlinear Fourier series (3). By deleting the terms corresponding to the radiation modes (those to the right of the reference level), and then summing the remaining hyperelliptic terms for the soliton part of the spectrum, one obtains only the contributions that the solitons make to the measured nonlinear time series. In this way one finds a long, low-amplitude wavetrain, which consists of eight nonlinearly interacting solitons. By using the numerical inverse scattering transform as a tool for nonlinear filtering one easily finds the solitons hidden in a sea of background radiation. An associated physical result is that the solitons tend to be phase locked beneath the maxima of the wave packets [9, 131. Complete theoretical understanding of this latter phenomena is still lacking. A further surprising result is that, while the nonlinear spectral modulus is near one for the solitons (as one might a priori suspect), it is also near one for the dominant part of the radiation spectrum as well (Fig. 4b). A formal interpretation of this effect can be given in terms of high-frequency soliton trains [2, 131. The application of nonlinear filtering to the analysis of the Adriatic Sea data in this highly nonlinear region of the radiation spectrum is therefore intimately related to the construction of nonlinear, narrowband wavetrains, an
Fourier analysis of ocean surface waves
2635
example of which is shown Fig. 6(b). Since the nonlinear modes are clearly not sinusoidal (see Fig. 6a), the effect of nonlinear interactions amongst these densely-spaced components is evidently important. It would therefore be inappropriate to interpret the radiation spectrum as being ‘nearly linear’ (historically a convenient designation) as contrasted to solitons as being ‘strongly nonlinear’. Clearly, large nonlinear interactions in the radiation components are an important topic of future research. Additional important consideration rests with understanding the interactions occurring among the radiation modes and the solitons, particularly with regard to the experimentally observed phase locking phenomena. Acknowledgemenrs-This work was supported in part by the Office of Naval Research of the United States of America (Grant NOOO14-92-J-1330)and by the Progetto Salvaguardia di Venezia de1 Consiglio Nazionale delle Ricerche , Italy.
REFERENCES 1. A. R. Osborne, Automatic algorithm for the numerical inverse scattering transform of the Korteweg-deVries equation, submitted for publication (1993). 2. A. R. Osborne, Numerical construction of nonlinear wavetrain solutions of the periodic Korteweg-deVries equation, Phys. Rev. E. 48, 296-309 (1993). 3. A. R. Bishop, M. G. Forest, D. W. McLaughlin and E. A. Overman II, A quasi-periodic route to chaos in a near integrable PDE, Physica D18,293-312 (1986). 4. A. R. Osborne and L. Bergamasco, The small-amplitude limit of the spectral transform for the periodic Korteweg-deVries equation, Nuovo Cimento B85, 229-243 (1985). 5. A. R. Osborne and L. Bergamasco, The solitons of Zabusky and Kruskal revisited: Perspective in terms of the periodic scattering transform, Physica D18, 26-243 (1986). 6. A. R. Osborne and E. Segre. Numerical solutions of the Korteweg-deVries equation using the periodic scattering transform p-representation, Physica D44, 575-604 (1990). 7. G. Terrones, D. W. McLaughlin, E. A. Overman II and A. Pearlstein, Stability and bifurcation of spatially coherent solutions of the damped driven nonlinear Schroedinger equation, SIAM J. Appl. Math. 50, 791-818 (1990). 8. R. Flesch, M. G. Forest and A. Sinha. Physicu D48, 169-208 (1991). 9. A. R. Osborne, E. Segre, G. Boffetta and L. Cavaleri, Soliton basis states in shallow-water ocean surface waves, Phys. Rev. Lett. 64, 1733-595 (1991). 10. A. R. Osborne, Nonlinear Fourier analysis, in Nonlinear Topics in Ocean Physics, edited by A. R. Osborne, pp. 669-700. North-Holland, Amsterdam (1991). 11. A. R. Osborne, Nonlinear Fourier analysis for the infinit-interval Korteweg-de Vries Equation: I. An algorithm for the direct scattering transform, J. Camp. Phys. 94(2), 284-313 (1991). 12. D. W. McLaughlin, C. M. Schober. Chaotic and homoclinic behavior for numerical discretizations of the nonlinear Schroedinger equation, Physicu D57, 447-465 (1992). 13. A. R. Osborne, M. Serio; L. Bergamasco and L. Caveleii, in preparation (1993). 14. Ocean Wave SDectra. Prentice-Hall. Enelewood Cliffs. N.J. (19611. 15. B. Kinsman, &nd Waves, PrenticelHaE, Englewood kliffs, ‘N.J. (1965). 16. A. Brandt, S. E. Ramberg and M. F. Shlesinger, Nonlinear Dynamics of Ocean Waves. World Scientific, Singapore (1992). 17. J. S. Bendat and A. G. Piersol, Random Data: Analysis and Measurement Procedures. Wiley-Interscience, New York (1971). 18. Nelson Wax, editor, Noise and Stochastic Processes. Dover, New York (1954). 19. L. Cavaleri, WAM, a third-generation wave model, in Nonlinear Topics in Ocean Physics, edited by A. R. Osborne, pp. 925-944. North-Holland, Amsterdam (1991). 20. A. R. Osborne, The behavior of solitons in random function solutions of the periodic Korteweg-deVries equation. Phys. Rev. Leti., to appear (1993). 21. S. Elgar and R. T. Guza, J. Fluid Mech. 167, 1 (1986). 22. V. E. Zakharov, S. V. Manakov, S. P. Novikov and M. P. Pitayevsky, Theory of Solitons. The Method of the Inverse Scattering Problem. Nauka, Moscow (1980), in Russian. 23. M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform. SIAM, Philadelphia (1981). 24. R. K. Dodd, J. E. Eilbeck, J. D. Gibbon and H. C. Morris, Solitons and Nonlinear Wave Equations, Academic Press, London (1982). 25. A. C. Newell, Solitons in Mathematics and Physics, SIAM, Philadelphia (1985). 26. A. Degasperis, Nonlinear wave equations solvable by the spectral transform, in Nonlinear Topics in Ocean Physics, edited by A. R. Osborne. Elsevier, Amsterdam (1991). 27. G. B. Whitham, Linear and Nonlinear Waves. Wiley-Interscience, New York (1974). 28. J. W. Miles, Solitary waves, Ann. Rev. Fluid Mech. 12, 11-43 (1980). 29. V. I. Karpman, Non-Linear Waves in Dispersive Media, Pergamon, Oxford (1975).
A. R. OSBORNE
2636
30. B. A. Dubrovin and S. P. Novikov, Periodic and conditionally periodic analogues of the many-soliton solutions of the Korteweg-de Vries equation, Sov. Phys. JEEP 40, 1058-1063 (1974). 31. B. A. Dubrovin, V. B. Matveev, and S. P. Novikov, Nonlinear equation of Korteweg-deVries type, finite-zone linear operators, and Abelian varieties, Russ. Math. Surv. 31, 59-146 (1976). 32. H. Flaschka, and D. W. McLaughlin, Canonically conjugate variables for KdV and Toda lattice under periodic boundary conditions, Prog. Theor. Phys. 55, 438-456 (1976). 33. H. P. McKean and E. Trubowitz, Hill’s operator and hyperelliptic function theory in the presence of infinitely many branch points, Commun. Pure Appl. Math. 29, 143-226 (1976). 34. A. R. Osborne, The numerical inverse scattering transform: nonlinear Fourier analysis and nonlinear filtering of oceanic surface waves, Proc. Aha Huliko’a Winter Workshop, edited by P. Mtiller and D. Henderson. University of Hawaii Press (1993). 35. L. Cavaleri, Experimental characteristics of wind waves, In Topics in Ocean Physics, edited by A. R. Osborne and P. Malanotte Rizzoli. Elsevier, Amsterdam (1982). 36. A. R. Osborne and M. Petti, Laboratory-generated, shallow water surface waves: Analysis using the periodic, inverse scattering transform. Submitted for publication (1993).
APPENDIX The numerical algorithm
Recent work [l, 21 is briefly summarized. The numerical search for the scattering eigenvalues {E,; pi; uj} suggests the need for computing the derivatives of the matrix Cuijwith respect to the energy E. This is because a Newtonian numerical root-finding algorithm is normally used to determine the eigenvalues. To achieve this goal, a matrix method has been developed for obtaining the evolution of the eigenfunction r# as a function of x and E for a particular wavetrain r](x, 0). The analytical determination of derivatives of the matrix elements with respect to E is a feature of the algorithm. Additional considerations with regard to automation of the computer codes are discussed in detail elsewhere [l]. The spectral equations are easily written: VL = -4* (A.11
The subscripts refer to differentiation with respect to x and E: q(x) = An(x) + E. Writing (A.l) in four-vector notation and using a Taylor series expansion for the solution to the scattering equations (A. 1) one obtains:
64.2)
) where
The elements of H are two-by-two matrices. The matrix 0 has ‘zero for its elements and the other matrices are given by: sin (KAx) --
cos @Ax)
T=
K
-resin @An)
(A.4)
cos (KAx) i
and Ax cos (KAx)
_ Ax sin (rcAx) 2K
T E =dr=
aE
-
Ax cos @Ax) + sin(rcAx) 2
-- sin (KAx)
2K2
2K
2K3
-~ Ax ( KAX) 2K
for K = (q)@ = (An(x) + E)@. While K may be either real or imaginary, the matrix TE is always real with determinant 1. This property is exploited in the numerical algorithm below. As in previous numerical problems of this type I assume the wavetrain q(x) has the form of a piecewise constant function with M partitions on the periodic interval (0, L), where the discretization interval is Ax = L/M [lo]. which is associated with a discrete value of the spatial Each partition has wave amplitude n, (0 G n c M-l) variable x, = nAx. It is important to note that the number of degrees of freedom of measured time series of M points is given by the simple relation, N = M/2 [4, 51. The four-by-four scattering matrix M can then be defined:
Fourier analysis of ocean surface waves M =
2637
64.6)
fl H(qnr Ax) n=M-1
The initial conditions of the basis (c, s) at the base point x0 are given by:
(A.7)
From the definition of the matrix (Yone finds:
(A.8) Thus at x = x0 ;(a%,
+
a2222) = (yzl
=
#fll
+
(A.9)
M22)
(A.lO)
Ml2
while the derivatives are given by (A.ll)
Because K = (kn(x, 0) + k2)@ can be either real or imaginary, but Implementation of the numerical algorithm. not complex, the matrix H is always real. This result allows implementation of an algorithm which is entirely real. The following relations have been used in the computer code: if K2 cos (K’Ax) (A.13) TI1 = T12 = fir K2 < 0; cash (K’Ax)
20
(A. 14)
T21
=
-K’ sin (K’Ax) sinh (K’Ax)
K’
if rc23 0 if K2 < 0,
(A.15)
where K’
=
VI/1117 +
k21 =
v/IK’I
(A.16)
and analogously for the matrix Ts. The reconstruction of complex solutions of the KdV equation by (3) ( as well as nonlinear filtering) are accomplished by computing the auxiliary spectra pj(xc = x,) for the M different base points x0 =x0, xi, The approach is formally called base-point iteration and is carried out by computing M different x2,. . .X&f-l. monodromy matrices (A.6) which differ from each other by a horizontal shift Ax in the wave train q,. This procedure arises from the following similarity transformation which arises from (A.6): (A. 17) M(x,+1, E) = H(%, E)M(x,, E)H(h E)-1 The latter expression relates the matrix M(x,+i, E) at a base point x,+i to the previously computed matrix M(x,, E) at the base point x, for a particular value of squared wavenumber E = k:. Values of the auxiliary spectra {p,(x,J} for each x, are computed from the matrices M(x”, E). Knowledge of the auxiliary spectra at every point x, allows reconstruction of the wave train n(x,,) via a discrete version of (3):
Mxn) = -El + :fWj(xrt) - E2j - E2,+11 ,=l
(A.18)
for n = 0, 1, 2,. . . M - 1. This is a finite-term nonlinear generalization of a Fourier series for the discrete wavetrain q(xn). As indicated by the notation, each nonlinear oscillation mode {nj} implicitly depends upon the associated wavenumber ki of the mode. The k, are theoretically given by the simple relation k, = jAk, Ak = 27r/L; surprisingly this is exactly the same as that for the linear Fourier transform, provided that periodicity is assumed. The IST spectrum then consists of the widths of the open bands of the Floquet discriminant, ai = (E2j+l - E2j),/21, graphed as a function of kj (or fj = jAf, Af = l/ir, for T the period of a time series).