Linear and nonlinear waves - Scholarpedia (original) (raw)

The study of waves can be traced back to antiquity where philosophers, such as Pythagoras (c. 560-480 BC), studied the relation of pitch and length of string in musical instruments. However, it was not until the work of Giovani Benedetti (1530-90), Isaac Beeckman (1588-1637) and Galileo (1564-1642) that the relationship between pitch and frequency was discovered. This started the science of acoustics, a term coined by Joseph Sauveur (1653-1716) who showed that strings can vibrate simultaneously at a fundamental frequency and at _integral multiples_that he called harmonics. Isaac Newton (1642-1727) was the first to calculate the speed of sound in his Principia. However, he assumed isothermal conditions so his value was too low compared with measured values. This discrepancy was resolved by Laplace (1749-1827) when he included adiabatic heating and cooling effects. The first analytical solution for a vibrating string was given by Brook Taylor (1685-1731). After this, advances were made by Daniel Bernoulli (1700-82), Leonard Euler (1707-83) and Jean d'Alembert (1717-83) who found the first solution to the linear wave equation, see section (The linear wave equation). Whilst others had shown that a wave can be represented as a sum of simple harmonic oscillations, it was Joseph Fourier (1768-1830) who conjectured that arbitrary functions can be represented by the superposition of an infinite sum of sines and cosines - now known as the Fourier series. However, whilst his conjecture was controversial and not widely accepted at the time, Dirichlet subsequently provided a proof, in 1828, that all functions satisfying Dirichlet's conditions(i.e. non-pathological piecewise continuous) could be represented by a convergent Fourier series. Finally, the subject of _classical acoustics_was laid down and presented as a coherent whole by John William Strutt (Lord Rayleigh, 1832-1901) in his treatise Theory of Sound. The science of modern acoustics has now moved into such diverse areas as sonar, auditoria, electronic amplifiers, etc.

The study of hydrostatics and hydrodynamics was being pursued in parallel with the study of acoustics. Everyone is familiar with Archimedes (c. 287-212 BC) eureka moment; however he also discovered many principles of hydrostatics and can be considered to be the father of this subject. The theory of fluids in motion began in the 17th century with the help of practical experiments of flow from reservoirs and aqueducts, most notably by Galileo's student Benedetto Castelli. Newton also made contributions in the Principia with regard to resistance to motion, also that the minimum cross-section of a stream issuing from a hole in a reservoir is reached just outside the wall (the vena contracta). Rapid developments using advanced calculus methods by Siméon-Denis Poisson (1781-1840), Claude Louis Marie Henri Navier (1785-1836), Augustin Louis Cauchy (1789-1857), Sir George Gabriel Stokes (1819-1903), Sir George Biddell Airy (1801-92), and others established a rigorous basis for hydrodynamics, including vortices and water waves, see section (Physical wave types). This subject now goes under the name of fluid dynamics and has many branches such as multi-phase flow, turbulent flow, inviscid flow, aerodynamics, meteorology, etc.

The study of electromagnetism was again started in antiquity, but very few advances were made until a proper scientific basis was finally initiated by William Gilbert (1544-1603) in his_De Magnete_. However, it was only late in the 18th century that real progress was achieved when Franz Ulrich Theodor Aepinus (1724-1802), Henry Cavendish (1731-1810), Charles-Augustin de Coulomb (1736-1806) and Alessandro Volta (1745-1827) introduced the concepts of charge, capacity and potential. Additional discoveries by Hans Christian Ørsted (1777-1851), André-Marie Ampère (1775-1836) and Michael Faraday (1791-1867) found the connection between electricity and magnetism and a full unified theory in rigorous mathematical terms was finally set out by James Clerk Maxwell (1831-79) in his Treatise on Electricity and Magnetism. It was in this work that all electromagnetic phenomena and all optical phenomena were first accounted for, including waves, see section (Electromagnetic wave). It also included the first theoretical prediction for the speed of light.

At the end of the 19th century, when some erroneously considered physics to be very nearly complete, new physical phenomena began to be observed that could not be explained. These demanded a whole new set of theories that ultimately led to the discovery of general relativity and quantum mechanics; which, even now in the 21st century are still yielding exciting new discoveries. However, as this article is primarily concerned with classical wave phenomena, we will not pursue these topics further.

Historic data source: 'Dictionary of The History of Science [Byn-84].

Contents

Introduction

A wave is a time evolution phenomenon that we generally model mathematically using partial differential equations (PDEs) which have a dependent variable \(u(x,t)\) (representing the wave value), an independent variable time \(t\) and one or more independent spatial variables \(x\in\mathbb{R}^{n}\ ,\) where \(n\) is generally equal to \(1,2 \;\textrm{or}\; 3\ .\) The actual form that the wave takes is strongly dependent upon the system initial conditions, the boundary conditions on the solution domain and any system disturbances.

Waves occur in most scientific and engineering disciplines, for example: fluid mechanics, optics, electromagnetism, solid mechanics, structural mechanics, quantum mechanics, etc. The waves for all these applications are described by solutions to either linear or _nonlinear_PDEs. We do not focus here on methods of solution for each type ofwave equation, but rather we concentrate on a small selection of relevant topics. However, first, it is legitimate to ask: what actually is a wave? This is not a straight forward question to answer.

Now, whilst most people have a general notion of what a wave is, based on their everyday experience, it is not easy to formulate a definition that will satisfy everyone engaged in or interested in this wide ranging subject. In fact, many technical works related to waves eschew a formal definition altogether and introduce the concept by a series of examples; for example, Physics of waves [Elm-69] and Hydrodynamics [Lam-93]. Nevertheless, it is useful to at least make an attempt and a selection of various definitions from normally authoritative sources is given below:

The variety of definitions given above, and their clearly differing degrees of clarity, confirm that 'wave' is indeed not an easy concept to define!

Because this is an introductory article and the subject of linear and non-linear waves is so wide ranging, we can only include sufficient material here to provide an overview of the phenomena and related issues. Relativistic issues will not be addressed. To this end we will discuss, as proxies for the wide range of known wave phenomena, the linear wave equation and the nonlinear Korteweg-de Vries equation in some detail by way of examples. To supplement this discussion we provide brief details of other types of wave equation and their application; and, finally, we introduce a number of PDE wave solution methods and discuss some general properties of waves. Where appropriate, references are included to works that provide further detailed discussion.

Physical wave types

A non-exhaustive list is given below of physical wave types with examples of occurrence and references where more details may be found.

\[\tag{1} {\displaystyle \left[\begin{array}{c} {\displaystyle \frac{\partial h}{\partial t}}\\ \\\frac{ {\displaystyle \partial u}}{{\displaystyle \partial t} }\end{array}\right]+\left[\begin{array}{c} {\displaystyle \frac{\partial\left(hu\right)}{\partial x}}\\ \\{\displaystyle \frac{\partial\left({\textstyle \frac{1}{2}}u^{2}+gh\right)}{\partial x}}\end{array}\right]=}\left[\begin{array}{c} 0\\ \\-g{\displaystyle \frac{\partial b}{\partial x}}\end{array}\right].\ \]

Where,

\(b\left(x\right)\) = fluid bed topography

\(h\left(x,t\right)\) = fluid surface height above bed

\(u\left(x,t\right)\) = fluid velocity - horizontal

\(g\) = acceleration due to gravity.

For this situation, the celerity or speed of wave propagation can be approximated by \(c=\sqrt{gh}\ .\) For detailed discussion refer to [Joh-97].

\[\tag{2} c_{p} = \frac{\omega}{k}=\sqrt{\frac{g}{k}}, \]

\[\tag{3} c_{g} = \frac{d\omega\left(k\right)}{dk}=\frac{1}{2}c_{p}.\ \]

The result is that the ship's wake is a wedge-shaped envelope of waves having a semi-angle of \(\backsimeq19.5\) degrees and a feathered pattern with the ship at the vertex. The shape is a characteristic of such waves, regardless of the size of disturbance - from a small duckling paddling on a pond to large ocean liner cruising across an ocean. These patterns are referred to as Kelvin Ship Waves after Lord Kelvin (William Thomson) [Joh-97].

Linear waves

General

Linear waves are described by linear equations, i.e. those where in each term of the equation the dependent variable and its derivatives are at most first degree (raised to the first power).

This means that the superposition principle applies, and linear combinations of simple solutions can be used to form more complex solutions. Thus, all the linear system analysis tools are available to the analyst, with Fourier analysis: expressing general solutions in terms of sums or integrals of well known basic solutions, being one of the most useful. The classic linear wave is discussed in section (The linear wave equation) with some further examples given in section (Linear wave equation examples). Linear waves are modelled by PDEs that are linear in the dependent variable, \(u\ ,\) and its_first_ and higher derivatives, if they exist.

The linear wave equation

The following represents the classical wave equation in one dimension and describes undamped linear waves in an isotropic medium

\[\tag{4} {\displaystyle \frac{1}{c^{2}}\frac{\partial^{2}u}{\partial t^{2}}}={\displaystyle \frac{\partial^{2}u}{\partial x^{2}}.}\]

It is second order in \(t\) and \(x\ ,\) and therefore requires two_initial condition functions_ (ICs) and two boundary condition functions (BCs). For example, we could specify\[\tag{5} \begin{array}{lcl} \textrm{ICs:}\quad u\left(x,t=0\right)=f\left(x\right),\quad u_{t}\left(x,t=0\right)=g\left(x\right) , \end{array}\]

\[\tag{6} \begin{array}{lcl}\textrm{BCs:}\quad u\left(x=a,t\right)=u_{a},\quad u\left(x=b,t\right)=u_{b}. \end{array}\]

Consequently, equations (4), (5) and (6) constitute a complete description of the PDE problem.

We assume \(f\) to have a continuous second derivative (written \(f\in C^{2}\)) and \(g\) to have a continuous first derivative (\(g\in C^{1}\)). If this is the case, then \(u\) will have continuous second derivatives in \(x\) and \(t\ ,\) i.e. (\(u\in C^{2}\)), and will be a correct solution to equation (4) with any consistent set of appropriate ICs and BCs [Stra-92].

Extending equation (4) to three dimensions, the classical wave equation becomes,\[\tag{7} \frac{1}{c^{2}}\frac{\partial^{2}u}{\partial t^{2}}=\nabla^{2}u,\]

where \(\nabla^{2}=\nabla\cdot\nabla\) represents the _Laplacian_operator. Because the Laplacian is co-ordinate free, it can be applied within any co-ordinate system and for any number of dimensions. Given below are examples of wave equations in 3 dimensions for Cartesian,cylindrical and spherical co-ordinate systems \[ \begin{array}{lccl} \textrm{Cartesian co-ordinates:} & {\displaystyle \frac{1}{c^{2}}\frac{\partial^{2}u}{\partial t^{2}}} & = & {\displaystyle \frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}},}\\ \textrm{Cylindrical co-ordinates}: & {\displaystyle \frac{1}{c^{2}}\frac{\partial^{2}u}{\partial t^{2}}} & = & {\displaystyle \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u}{\partial r}\right)+\frac{1}{r^{2}}\frac{\partial^{2}u}{\partial\theta^{2}}+\frac{\partial^{2}u}{\partial z^{2}},}\\ \textrm{Spherical co-ordinates}: & {\displaystyle \frac{1}{c^{2}}\frac{\partial^{2}u}{\partial t^{2}}} & = & {\displaystyle \frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial u}{\partial r}\right)+\frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial u}{\partial\theta}\right)+\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}u}{\partial\phi^{2}}.}\end{array}\]

These equations occur in one form or another, in numerous applications in all areas of the physical sciences; see for example section (Linear wave equation examples ).

The d'Alembert solution

The solution to equations (4), (5) and (6) was first reported by the French mathematician Jean-le-Rond d'Alembert (1717-1783) in 1747 in a treatise on Vibrating Strings [Caj-61] [Far-93]. D'Alembert's remarkable solution, which used a method specific to the wave equation (based on the chain rule for differentiation), is given below

\[\tag{8} u(x,t)=\frac{1}{2}\left[f(x-ct)+f(x+ct)\right]+\frac{1}{2c}\int_{x-ct}^{x+ct}g(\xi)d\xi.\]

It can also be obtained by the Fourier Transform method or by the separation of variables (SOV) method, which are more general than the the method used by d'Alembert [Krey-93].

The d'Alembertian \(\square=\nabla^{2}-{\displaystyle \frac{1}{c^{2}}\frac{\partial^{2}}{\partial t^{2}}}\ ,\) also known as the d'Alembert operator or wave operator, allows a succinct notation for the wave equation, i.e. \(\square u=0\ .\) It first arose in d'Alembert's work on vibrating strings and plays a useful role in modern theoretical physics.

Linear wave equation examples

Acoustic (sound) wave

We will consider the acoustic or sound wave as a small amplitude disturbance of ambient conditions where second order effects can be ignored. We start with the Euler continuity and momentum equations

\[\tag{9} \frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho v\right) = 0,\]

\[\tag{10} \frac{\partial\left(\rho v\right)}{\partial t}+\nabla\cdot\left(\rho vv\right)-\rho g+\nabla p+\nabla\cdot T = 0,\]

where

\(T\) = stress tensor (Pa)

\(g\) = gravitational acceleration (m/s\(^{2}\))

\(p\) = pressure (Pa)

\(t\) = time (s)

\(v\) = fluid velocity (m/s)

\(\rho\) = fluid density (kg/m\(^{3}\))

We assume an inviscid dry gas situation where gravitational effects are negligible. This means that the third and fifth terms of equation (10) can be ignored. If we also assume that we can represent velocity by \(v=u_{0}+u\ ,\) where \(u_{o}\) is ambient velocity which we set to zero and \(u\) represents a small velocity disturbance, the second term in equation (10) can be ignored (because it becomes a second order effect). Thus, equations (9) and (10) reduce to

\[\tag{11} \frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho u\right) = 0,\]

\[\tag{12} \frac{\partial\left(\rho u\right)}{\partial t}+\nabla p = 0.\]

Now, taking the divergence of equation (12) and the time derivative of equation (11), we obtain\[ \frac{\partial^{2}\rho}{\partial t^{2}}-\nabla^{2}p=0.\]

To complete the analysis we need to apply an _equation-of-state_relating \(p\) and \(\rho\) when we obtain the linear acoustic wave equation

\[\tag{13} \frac{1}{c^{2}}\frac{\partial^{2}p}{\partial t^{2}}=\nabla^{2}p,\]

where \[\tag{14} c^{2}=\frac{\partial p}{\partial\rho}\ .\]

We now consider three cases:

For atmospheric air at standard conditions we have \(p=101325\)Pa, \(T_{0}=293.15\)K, \(R=8.3145\)J/mol/K, \(\gamma=1.4\) and \(MW=0.028965\)kg/mol, which gives

\[\tag{15} \textrm{isothermal:}\quad c = 290\textrm{m/s,}\]

\[\tag{16} \textrm{isentropic:}\quad \; c = 343\textrm{m/s.}\]

For liquid distilled water at \(20\)C we have \(\beta=2.18\times10^{9}\)Pa and \(\rho=1,000\)kg/m\(^{3},\) which gives

\[\tag{17} \textrm{liquid}:\quad c=1476\textrm{m/s.}\]

Waves in solids

Waves in solids are more complex than acoustic waves in fluids. Here we are dealing with displacement \(\varrho\ ,\) and the resulting waves can be either longitudinal, P-waves, or shear(transverse), S-waves. Starting with Newton's second Law we arrive at the vector wave equation [Elm-69, chapter 7]

\[\tag{18} \left(\lambda+\mu\right)\nabla\left(\nabla\cdot\varrho\right)+\mu\nabla^{2}\varrho=\rho\frac{\partial^{2}\varrho}{\partial t^{2}},\]

from which, using the fundamental identity from vector calculus, \(\nabla\times\left(\nabla\times\varrho\right)=\nabla\left(\nabla\cdot\varrho\right)-\nabla^{2}\varrho\ ,\) we obtain

\[\tag{19} \left(\lambda+2\mu\right)\nabla\left(\nabla\cdot\varrho\right)+\mu\nabla\times\left(\nabla\times\varrho\right)=\rho\frac{\partial^{2}\varrho}{\partial t^{2}}.\]

Now, for irrotational waves, which vibrate only in the direction of propagation \(x\ ,\) \(\nabla\times\varrho=0\Rightarrow\nabla\left(\nabla\cdot\varrho\right)=\nabla^{2}\varrho\) and equation (19) reduces to the familiar linear wave equation

\[\tag{20} \frac{1}{c^{2}}\frac{\partial^{2}\varrho}{\partial t^{2}}=\nabla^{2}\varrho,\]

where \(c=\sqrt{\left(\lambda+2\mu\right)/\rho}=\sqrt{\left(K+\frac{4}{3}\mu\right)/\rho}\) is the wave speed, \(\lambda=E\upsilon/\left(1+\upsilon\right)\left(1-2\upsilon\right)\) is the Lamé modulus, \(\mu={\displaystyle \frac{E}{2\left(1+\upsilon\right)}}\) is the shear modulus and \(K=E/3\left(1-2\upsilon\right)\) is the bulk modulus of the solid material. Here, \(E\) and \(\upsilon\) are Young's modulus and Poisson's ratio for the solid respectively. Irrotational waves are of the longitudinal type, or P-waves.

For solenoidal waves, which can vibrate independently in the \(y\) and \(z\) directions but not in the direction of propagation \(x\ ,\) we have \(\nabla\cdot\varrho=0\) and equation (18) reduces to the linear wave equation

\[\tag{21} \frac{1}{c^{2}}\frac{\partial^{2}\varrho}{\partial t^{2}}=\nabla^{2}\varrho,\]

where the wave speed is given by \(c=\sqrt{\mu/\rho}\) . Solenoidal waves are of the transverse type, or S-waves.

For a typical mild-steel at \(20\)C with \(\rho=7,860\)kg/m\(^{3}\ ,\) \(E=210\times10^{9}\)N/m\(^{2}\) and \(\upsilon=0.29\) we find that the P-wave speed is \(5917\)m/s and the S-wave speed is \(3,218\)m/s. For further discussion refer to [Cia-88].

Electromagnetic waves

The fundamental equations of electromagnetism are the Maxwell Equations, which in differential form and SI units, are usually written as:

\[\tag{22} \nabla\cdot E = \frac{1}{\epsilon_{0}}\rho,\]

\[\tag{23} \nabla\cdot B = 0,\]

\[\tag{24} \nabla\times E = -\frac{\partial B}{\partial t},\]

\[\tag{25} \nabla\times B = \mu_{0}J+\mu_{0}\epsilon_{0}\frac{\partial E}{\partial t},\]

where

\(B =\) magnetic field (T)

\(E =\) electric field (V/m)

\(J =\) current density (A/m\(^{2}\))

\(\; t =\) time (s)

\(\epsilon_{0} =\) permittivity of free space (\(8.8541878\times10^{-12}\simeq10^{-9}/36\pi\) F/m)

\(\mu_{0} =\) permeability of free space (\(4\pi\times10^{-7}\) H/m)

\(\; \rho =\) charge density (C/m\(^{3}\))

If we assume that \(J=0\) and \(\rho=0\ ,\) then on taking the curl of equation (24) and again using the fundamental identity from vector calculus, \(\nabla\times\left(\nabla\times E\right)=\nabla\left(\nabla\cdot E\right)-\nabla^{2}E\ ,\) we obtain

\[\tag{26} \frac{1}{c_{0}^{2}}\frac{\partial^{2}E}{\partial t^{2}}=\nabla^{2}E.\]

Similarly, taking the curl of equation (25) we obtain

\[\tag{27} \frac{1}{c_{0}^{2}}\frac{\partial^{2}B}{\partial t^{2}}=\nabla^{2}B.\]

Equations (26) and (27) are the linear electric and magnetic wave equations respectively, where \(c_{0}=1/\sqrt{\mu_{0}\epsilon_{0}}\simeq3\times10^{8}\) m/s, the speed of light in a vacuum. They take the familiar form of linear wave equation (4). For further discussion refer to [Sha-75].

Nonlinear waves

General

Nonlinear waves are described by nonlinear equations, and therefore the superposition principle does not generally apply. This means that nonlinear wave equations are more difficult to analyze mathematically and that no general analytical method for their solution exists. Thus, unfortunately, each particular wave equation has to be treated individually. An example of solving the Korteweg-de Vries equation by direct integration is given below.

Some advanced methods that have been used successfully to obtain closed-form solutions are listed in section (Closed form PDE solution methods), and example solutions to well known evolution equations are given in section (Nonlinear wave equation solutions).

Closed form PDE solution methods

There are no general methods guaranteed to find closed form solutions to non-linear PDEs. Nevertheless, some problems can yield to a _trial-and-error_approach. This hit-and-miss method seeks to deduce candidate solutions by looking for clues from the equation form, and then systematically investigating whether or not they satisfy the particular PDE. If the form is close to one with an already known solution, this approach may yield useful results. However, success is problematical and relies on the analyst having a keen insight into the problem.

We list below, in alphabetical order, a non-exhaustive selection of advanced solution methods that can assist in determining closed form solutions to nonlinear wave equations. We will not discuss further these methods and refer the reader to the references given for details. All these methods are greatly enhanced by use of a symbolic computer program such as: Maple V, Mathematica, Macysma, etc.

Some techniques for obtaining traveling wave solutions

The following are examples of techniques that transform PDEs into ODEs which are subsequently solved to obtain traveling wave solutions to the original equations.

Some example applications of these and other methods can be found in [Gri-11].

Nonlinear wave equation solutions

A non-exhaustive selection of well known 1D nonlinear wave equations and their closed-form solutions is given below. The closed form solutions are given by way of example only, as nonlinear wave equations often have many possible solutions.

- Applications: gas dynamics and traffic flow.

- Solution\[u=\varphi\left(\xi\right),\;\xi=x-\varphi\left(\xi\right)t.\]

where

\(u\left(x,t=0\right)=\varphi\left(x\right)\), arbitrary initial condition.

- Applications: acoustic and hydrodynamic waves.

- Solution\[u(x,t)=2ak\left[1-\tanh k\left(x-Vt\right)\right] .\]

where

\(k=\) wavenumber,

\(V=\)velocity,

\(a=\) arbitrary constant.

- Applications: heat and mass transfer, population dynamics, ecology.

- Solution\[u(x,t)=\frac{1}{4}\left\{ 1-\tanh k\left[x-Vt\right]\right\} ^{2}.\]

where

\(k={\displaystyle \frac{1}{2\sqrt{6}}}\) (wavenumber),

\(V={\displaystyle \frac{5}{\sqrt{6}}}\) (velocity).

Note: wavenumber and velocity are fixed values.

- Applications: various areas of physics

- Solution\[u\left(x,t\right)=\left\{ \begin{array}{l} {\displaystyle \frac{4}{\lambda}}\arctan\left[\exp\left(\pm{\displaystyle \frac{b\lambda\left(kx+\mu t+\theta_{0}\right)}{\sqrt{b\lambda\left(\mu^{2}-ak^{2}\right)}}}\right)\right],\quad b\lambda\left(\mu^{2}-ak^{2}\right)>0,\\ \\{\displaystyle \frac{4}{\lambda}}\arctan\left[\exp\left(\pm{\displaystyle \frac{b\lambda\left(kx+\mu t+\theta_{0}\right)}{\sqrt{b\lambda\left(ak^{2}-\mu^{2}\right)}}}\right)\right]-{\displaystyle \frac{\pi}{\lambda}},\quad b\lambda \left(\mu^{2}-ak^{2}\right)<0.\end{array}\right. \]

where

\(k=\) wavenumber,

\(\mu, \theta_{0}= \) arbitrary constants.

- Applications: various areas of physics, non-linear optics, superconductivity, plasma models.

- Solution\[u(x,t)=\sqrt{\frac{\alpha}{q}}\textrm{sech}\left(\sqrt{\alpha}\left(x-Vt\right)\right),\quad \alpha>0,q>0 .\]

where

\(V=\)velocity,

\(\alpha,q=\) arbitrary constants.

- Applications: various areas of physics, nonlinear mechanics, water waves.

- Solution\[u(x,t)=12bk^{2}\textrm{sech}^{2}k\left(x-Vt\right)\]

where

\(k=\)wavenumber,

\(V=\)velocity,

\(b=\) arbitrary constant.

- Applications: surface water waves

- Solution\[{\displaystyle \frac{1}{6}\left\{ 1+8k^{2}-V^{2}\right\} -2k^{2}\tanh^{2}k\left(x+Vt\right)}\]

where

\(k=\)wavenumber,

\(V=\)velocity.

This equation can be linearized in the general case. Some exact solutions are given in [Pol-04, pp252-255] and, by way of an example consider the following special case where \(f\left(u\right)=\alpha e^{\lambda u}\ :\)

Wave equation with exponential non-linearity: \(u_{tt}=\left(\alpha e^{\lambda u}u_{x}\right)_{x},\quad\alpha>0.\) [Pol-04, p223]

- Applications: traveling waves

- Solution\[u(x,t)={\displaystyle \frac{1}{\lambda}}\ln\left(\alpha ax^{2}+bx+c\right)-{\displaystyle \frac{2}{\lambda}}\ln\left(\alpha at+d\right)\ :\]

where

\(\alpha,\lambda,a,b,c,d=\) arbitrary constants.

Additional wide-ranging examples of traveling wave equations, with solutions, from the fields of mathematics, physics and engineering are given in Polyanin & Manzhirov [Pol-07] and Polyanin & Zaitsev [Pol-04]. Examples from the biological and medical fields can be found in Murray [Mur-02] and Murray [Mur-03]. A useful on-line resource is the DispersiveWiki [Dis-08].

The Korteweg-de Vries equation

The canonical form of the Korteweg-de Vries (KdV) equation is

\[\tag{28} \frac{\partial u}{\partial t}-6u\frac{\partial u}{\partial x}+\frac{\partial^{3}u}{\partial x^{3}}=0,\]

and is a non-dimensional version of the following equation originally derived by Korteweg and de Vries for a _moving (Lagrangian)_frame of reference [Jag-06], [Kor-95],

\[\tag{29} \frac{\partial\eta}{\partial\tau}=\frac{3}{2}\sqrt{\frac{g}{h_{o}}}\frac{\partial}{\partial\chi}\left[\frac{1}{2}\eta^{2}+\frac{2}{3}\alpha\eta+\frac{1}{3}\sigma\frac{\partial^{2}\eta}{\partial\chi^{2}}\right].\]

It is, historically, the most famous solitary wave equation and describes small amplitude, shallow water waves in a channel, where symbols have the following meaning:

\(g =\) gravitational acceleration (m/s\(^{2}\))

\(h_{o} =\) nominal water depth (m)

\(T =\) capillary surface tension of fluid (N/m)

\(\alpha =\) small arbitrary constant related to the uniform motion of the liquid (dimensionless)

\(\eta =\) wave height (m)

\(\rho =\) fluid density (kg/m\(^{3}\))

\(\tau =\) time (s)

\(\chi =\) distance (m)

After re-scaling and translating the dependent and independent variables to eliminate the physical constants using the transformations [Abl-91],

\[\tag{30} u=-\frac{1}{2}\eta-\frac{1}{3}\alpha;\quad x=-\frac{\chi}{\sqrt{\sigma}};\quad t=\frac{1}{2}\sqrt{\frac{g}{h_{o}\sigma}}\tau\]

where \(\sigma=h_{o}^{3}/3-Th_{o}/\left(\rho g\right)\ ,\) and \(Th_{o}/\left(\rho g\right)\) is called the Bond number (a measure of the relative strengths of surface tension and gravitational force), we arrive at the Korteweg-de Vries equation, i.e. equation (28).

The basic assumptions for the derivation of KdV waves in liquid, having wavelength \(\lambda\ ,\) are [Abl-91]:

The KdV equation was found to have solitary wave solutions [Lam-93], which confirmed John Scott-Russell's account of the solitary wave phenomena [Sco-44] discovered during his experimental investigations into water flow in channels to determine the most efficient design for canal boats [Jag-06]. Subsequently, the KdV equation has been shown to model various other nonlinear wave phenomena found in the physical sciences. John Scott-Russell, a Scottish engineer and naval architect, also described in poetic terms his first encounter with the solitary wave phenomena, thus:

"I was observing the motion of a boat which was rapidly drawn along a narrow channel by a pair of horses, when the boat suddenly stopped - not so the mass of water in the channel which it had put in motion; it accumulated round the prow of the vessel in a state of violent agitation, then suddenly leaving it behind, rolled forward with great velocity, assuming the form of a large solitary elevation, a rounded, smooth and well-defined heap of water, which continued its course along the channel apparently without change of form or diminution of speed. I followed it on horseback, and overtook it still rolling on at a rate of some eight or nine miles an hour, preserving its original figure some thirty feet long and a foot to a foot and a half in height. Its height gradually diminished, and after a chase of one or two miles I lost it in the windings of the channel. Such, in the month of August 1834, was my first chance interview with that singular and beautiful phenomenon which I have called the Wave of Translation" [Sco-44].

An experimental apparatus for re-creating the phenomena observed by Scott-Russell have been built at Herriot-Watt University. Scott-Russell also coined the term solitary wave and conducted some of the first experiments to investigate another nonlinear wave phenomena, the Doppler effect, publishing an independent explanation of the theory in 1848 [Sco-48].

It is interesting to note that, a KdV solitary wave in water that experiences a change in depth will retain its general shape. However, on encountering shallower water its velocity and height will increase and its width decrease; whereas, on encountering deeper water its velocity and height will decrease and its width increase [Joh-97, pp 268-277].

A closed form single soliton solution to the KdV equation (28) can be found using direct integration as follows.

Assume a travelling wave solution of the form\[\tag{31} u(x,t)=f(x-vt)=f(\xi).\]

Then on substituting into the canonical equation the PDE is transformed into the following ODE

\[\tag{32} -v\frac{df(\xi)}{d\xi}-6f\frac{df(\xi)}{d\xi}+\frac{d^{3}f(\xi)}{d\xi^{3}}=0.\]

Now integrate with respect to \(\xi\) and multiply by \({\displaystyle \frac{df(\xi)}{d\xi}}\) to obtain\[\tag{33} -vf(\xi)\frac{df(\xi)}{d\xi}-3f(\xi)^{2}\frac{df(\xi)}{d\xi}+\frac{df(\xi)}{d\xi}\left(\frac{d^{2}f(\xi)}{d\xi^{2}}\right)=A\frac{df(\xi)}{d\xi}.\]

Now integrate with respect to \(\xi\) once more, to obtain\[\tag{34} -\frac{1}{2}vf(\xi)^{2}-f(\xi)^{3}+\frac{1}{2}\left(\frac{df(\xi)}{d\xi}\right)^{2}=Af(\xi)+B.\]

Where \(A\) and \(B\) are arbitrary constants of integration which we set to zero. We justify this by assuming that we are modeling a physical system with properties such that \(f,f^{\prime}\) and \(f^{\prime\prime}\rightarrow0\) as \(\xi\rightarrow\pm\infty\ .\) After rearranging and evaluating the resulting integral, we find\[\tag{35} f\left(\xi\right)=\frac{v}{2}\textrm{sech}^{2}\left(\frac{\sqrt{v}}{2}\xi\right).\]

The solution is therefore\[\tag{36} u(x,t) = f(x-vt),\]

\[\tag{37} \quad= 2k^{2}\textrm{sech}^{2}\left(k\left[x-vt-x_{0}\right]\right),\]

where \(k={\displaystyle \frac{\sqrt{v}}{2}}\) represents _wavenumber_and the constant \(x_{0}\) has been included to locate the wave peak at \(t=0\ .\) Thus, we observe that the wave travels to the right with a speed that is equal to twice the peak amplitude. Hence, the taller a wave the faster it travels.

The KdV equation also admits many other solutions including multiple soliton solutions, see figure (15), and cnoidal (periodic) solutions.

Solutions of KdV equation can be systematically obtained from solutions \(\psi_{i}\) of of the free particle Schrödinger equation \[\tag{38} -\left(\frac{\partial^{2}}{\partial x^{2}}\psi_{i}\right)=E_{i}\psi_{i},\quad i=1,\cdots,n\]

using the the relationship\[\tag{39} u\left(x,t\right)=2\left(\frac{\partial^{2}}{\partial x^{2}}\ln\left(W_{n}\right)\right),\]

where we use the the Wronskian function \[\tag{40} W_{n}=W_{n}\left[\psi_{1},\psi_{2},\cdots,\psi_{n}\right].\]

The Wronskian is the determinant of a \(n\times n\) matrix [Dra-89]composed from the functions \(\psi_{i}(\xi_{i})\ ,\) where \(\xi_{i}\) for our purposes is given by\[\tag{41} \xi_{i} = k_{i}\left(x-v_{i}t\right),\quad E_{i}<0,\]

\[\tag{42} \xi_{i} = k_{i}\left(x+v_{i}t\right),\quad E_{i}>0.\]

For example, a two-soliton solution is given by\[\tag{43} u(x,t)=\frac{\left(k_{1}^{2}-k_{2}^{2}\right)\left\{ 2k_{2}^{2}\textrm{csch}\, k_{2}\left(x-v_{2}t\right)+2k_{1}^{2}\textrm{sech}\, k_{1}\left(x-v_{1}t\right)\right\} }{\left[k_{1}\tanh k_{1}\left(x-v_{1}t\right)+k_{2}\coth k_{2}\left(x-v_{2}t\right)\right]^{2}}\]

and a cnoidal wave solution is given by

\[\tag{44} u(x,t)=\frac{1}{6k}\left(4k^{2}(2m-1)-vk\right)-2k^{2}\textrm{cn}^{2}\left(kx-vkt+x_{0};m\right).\]

where 'cn' represents the Jacobi elliptic cosine function with modulus \( m, \left( 0<m<1 \right)\). Note: as \(m\rightarrow1\) the periodic solution tends to a single soliton solution.

Interestingly, the KdV equation is invariant under a Galilean transformation, i.e. its properties remain unchanged, see section (Galilean invariance).

Numerical solution methods

Linear and nonlinear evolutionary wave problems can very often be solved by application of general numerical techniques such as: finite difference, finite volume, finite element, spectral, least_squares, weighted residual (e.g. collocation and Galerkin) methods,_etc. These methods, which can all handle various boundary conditions,stiff problems and may involve explicit or _implicit_calculations, are well documented in the literature and will not be discussed further here. For general texts refer to [Bur-93],[Sch-94],[Sch-09], and for more detailed discussion refer to [Lev-02],[Mor-94],[Zie-77].

Some wave problems do, however, present significant problems when attempting to find a numerical solution. In particular we highlight problems that include shocks, sharp fronts or _large gradients_in their solutions. Because these problems often involve _inviscid_conditions (zero or vanishingly small viscosity), it is often only practical to obtain weak solutions. Some PDE problems do not have a mathematically rigorous solution, for example where _discontinuities_or jump conditions are present in the solution and/or characteristics intersect. Such problems are likely to occur when there is a hyperbolic (strongly convective) component present. In these situations weak solutions provide useful information. Detailed discussion of this approach is beyond the scope of this article and readers are referred to [Wes-01, chapters 9 and 10] for further discussion.

General methods are often not adequate for accurate resolution of steep gradient phenomena; they usually introduce non-physical effects such as smearing of the solution or spurious oscillations. Since publication of Godunov's order barrier theorem, which proved that linear methods cannot provide non-oscillatory solutions higher than first order [God-54],[God-59], these difficulties have attracted a lot of attention and a number of techniques have been developed that largely overcome these problems. To avoid spurious or non-physical oscillations where shocks are present, schemes that exhibit a total variation diminishing (TVD) characteristic are especially attractive. Two techniques that are proving to be particularly effective are MUSCL (Monotone Upstream-Centred Schemes for Conservation Laws) a flux/slope limiter method[van-79],[Hir-90],[Tan-97],[Lan-98],[Tor-99] and the WENO (Weighted Essentially Non-Oscillatory) method [Shu-98],[Shu-09]. MUSCL methods are usually referred to as _high resolution schemes_and are generally second-order accurate in smooth regions (although they can be formulated for higher orders) and provide good resolution, monotonic solutions around discontinuities. They are straight-forward to implement and are computationally efficient. For problems comprising both shocks and complex smooth solution structure, WENO schemes can provide higher accuracy than second-order schemes along with good resolution around discontinuities. Most applications tend to use a fifth order accurate WENO scheme, whilst higher order schemes can be used where the problem demands improved accuracy in smooth regions.

Initial conditions and boundary conditions

Consider the classic 1D linear wave equation

\[\tag{45} \dfrac{\partial^{2}u}{\partial t^{2}}=\frac{1}{c^{2}}\dfrac{\partial^{2}u}{\partial x^{2}}.\]

In order to obtain a solution we must first specify some auxiliary conditions to complete the statement of the PDE problem. The number of required auxiliary conditions is determined by the highest order derivative in each independent variable. Since equation (45) is second order in \(t\) and second order in \(x\ ,\) it requires two auxiliary conditions in \(t\) and two auxiliary conditions in \(x\ .\) To have a complete well posed problem, some additional conditions may have to be included - refer to section (Wellposedness).

The variable \(t\) is termed an initial value variable and therefore requires two initial conditions (ICs). It is an initial value variable since it starts at an initial value, \(t_{0}\ ,\) and moves forward over a finite interval \(t_{0}\leq t\leq t_{f}\) or a semi-infinite interval \(t_{0}\leq t\leq\infty\) without any additional conditions being imposed. Typically in a PDE application, the initial value variable is time, as in the case of equation (45).

The variable \(x\) is termed a boundary value variable and therefore requires two boundary conditions (BCs). It is a boundary value variable since it varies over a _finite interval_\(x_{0}\leq x\leq x_{f}\ ,\) a semi-infinite interval \(x_{0}\leq x\leq\infty\) or a fully infinite interval \(-\infty\leq x\leq\infty\ ,\) and at two different values of \(x\), conditions are imposed on \(u\) in equation (45). Typically, the two values of \(x\) correspond to boundaries of a physical system, and hence the name boundary conditions.

BCs can be of three types:

An important consideration is the possibility of discontinuities at the boundaries, produced for example by differences in initial and boundary conditions at the boundaries, which can cause computational difficulties, such as shocks - see section (Shock waves), particularly for hyperbolic PDEs such as equation (45) above.

Numerical dissipation and dispersion

General

Some dissipation and dispersion occur naturally in most physical systems described by PDEs. Errors in magnitude are termed _dissipation_and errors in phase are called dispersion. These terms are defined below. The term amplification factor is used to represent the change in the magnitude of a solution over time. It can be calculated in either the time domain, by considering solution harmonics, or in the complex frequency domain by taking Fourier transforms.

Dissipation and dispersion can also be introduced when PDEs are discretized in the process of seeking a numerical solution. This introduces numerical errors. The accuracy of a discretization scheme can be determined by comparing the numeric amplification factor \(G_{numeric},\) with the analytical or exact amplification factor \(G_{exact}\ ,\) over one time step.

For further reading refer to [Hir-88, chap. 8], [Lig-78, chap. 3], [Tan-97, chap. 4], [Wes-01, chap 8 and 9].

Dispersion relation

Physical waves that propagate in a particular medium will, in general, exhibit a specific group velocity as well as a specific phase velocity - see section (Group and phase velocity). This is because within a particular medium there is a fixed relationship between the wavenumber \(k\ ,\) and the frequency \(\omega\ ,\) of waves. Thus, frequency and wavenumber are not independent quantities and are related by a functional relationship, known as the dispersion relation , \(\omega(k)\).

We will demonstrate the process of obtaining the dispersion relation by example, using the advection equation

\[\tag{46} u_{t}+au_{x}=0.\]

Generally, each wavenumber \(k \,\ ,\) corresponds to \(s \,\) frequencies where \(s \,\) is the order of the PDE with respect to \(t\ .\) Now any linear PDE with constant coefficients admits a solution of the form\[\tag{47} u\left(x,t\right)=u_{0}e^{i\left(kx-\omega t\right)}.\]

Because we are considering a linear system, the principal of superposition applies and equation (47) can be considered to be a frequency component or harmonic of the Fourier series representation of a specific solution to the advection equation. On inserting this solution into a PDE we obtain the so called dispersion relation between \(\omega\) and \(k\) i.e.,\[\tag{48} \omega=\omega\left(k\right),\]

and each PDE will have its own distinct form. For example, we obtain the specific dispersion relation for the advection equation by substituting equation (47) into equation (46) to get \[ -i\omega u_{0}e^{i\left(kx-\omega t\right)} = -iaku_{0}e^{i\left(kx-\omega t\right)}\] \[\Downarrow \]\[\tag{49} \omega = ak.\]

This confirms that \(\omega\) and \(k\) cannot be determined independently for the advection equation, and therefore equation (47) becomes\[\tag{50} u\left(x,t\right)=u_{0}e^{ik\left(x-at\right)}.\]

Note: If the imaginary part of \(\omega\left(k\right)\) is zero, then the system is non-dissipative.

The physical meaning of equation (50) is that the initial value \(u\left(x,0\right)=u_{0}e^{ikx}\ ,\) is propagated from left to right, unchanged, at velocity \(a\ .\) Thus, there is no dissipation or attenuation and no dispersion.

A similar approach can be used to establish the dispersion relation for systems described by other forms of PDEs.

Amplification factor

As mentioned above, the accuracy of a numerical scheme can be determined by comparing the numeric amplification factor \(G_{numeric},\) with the exact amplification factor \(G_{exact}\ ,\) over one time step. The exact amplification factor can be determined by considering the change that takes place in the exact solution over a single time-step. For example, taking the advection equation (46) and assuming a solution of the form \(u\left(x,t\right)=u_{0}e^{ik\left(x-at\right)}\ ,\) we have \[ G_{exact} = \frac{u\left(x,t+\Delta t\right)}{u\left(x,t\right)}=\frac{u_{0}e^{ik\left(x-a\left(t+\Delta t\right)\right)}}{u_{0}e^{ik\left(x-at\right)}}.\]\[\tag{51} \therefore G_{exact} = e^{-iak\Delta t}.\]

We can also represent equation (51) in the form\[\tag{52} G_{exact}=\left|G_{exact}\right|e^{i\Phi_{exact}},\]

where \[\tag{53} \Phi_{exact}=\angle G=\tan^{-1}\left(\frac{\textrm{Im}\left\{ G\right\} }{\textrm{Re}\left\{ G\right\} }\right).\]

Thus, for this case \[\tag{54} \left|G_{exact}\right| = 1\]

and\[\tag{55} \Phi_{exact} = \tan^{-1}\left(\tan\left(-ak\Delta t\right)\right)=-ak\Delta t.\]

The amplification factor provides an indication of how the the solution will evolve because values of \(\left|\Phi\right|\rightarrow0\) are associated with low frequencies and values of \(\left|\Phi\right|\rightarrow\pi\) are associated with high frequencies. Also, because phase shift is associated with the imaginary part of \( G_{exact}\ ,\) if \(\Im\left\{ G_{exact}\right\} =0\ ,\) the system does not exhibit any phase shift and is purely dissipative. Conversely, if \(\Re\left\{ G_{exact}\right\} =1\ ,\) the system does not exhibit any amplitude attenuation and is purely dispersive

The numerical amplification factor \(G_{numeric}\) is calculated in the same way, except that the appropriate _numerical approximation_is used for \(u(x,t)\ .\) For stability of the numerical solution, \(\left|G_{numeric}\right|\leq1\) for all frequencies.

Numerical dissipation

Figure 1: Figure 1: Illustration of pure numeric dissipation effect on a single sinusoid, as it propagates along the spatial domain. Both exact and simulated dissipative waves begin with the same amplitude; however, the amplitude of the dissipative wave decreases over time, but stays in phase.

Figure 2: Figure 2: Effect of numerical dissipation on a step function applied to the advection equation \(u_{t}+u_{x}=0\ .\)

In a numerical scheme, a situation where waves of different frequencies are damped by different amounts, is called _numerical dissipation,see figure (1). Generally, this results in the higher frequency components being damped more than lower frequency components. The effect of dissipation therefore is that sharp gradients, discontinuities or shocks in the solution tend to be smeared out, thus losing resolution, see figure (2). Fortunately, in recent years, various high resolution schemes have been developed to obviate this effect to enable shocks to be captured with a high degree of accuracy, albeit at the expense ofcomplexity. Examples of particularly effective schemes are based upon_flux/slope limiters [Wes-01] and WENO methods [Shu-98]. Dissipation can be introduced by numerical discretization of a partial differential equation that models a non-dissipative process. Generally, dissipation improves stability and, in some numerical schemes it is introduced deliberately to aid stability of the resulting solution. Dissipation, whether real or numerically induced, tend to cause waves to lose energy.

The dissipation error as a result of discretization can be determined by comparing the magnitude of the _numeric amplification factor_\(\left|G_{numeric}\right|,\) with the magnitude of the exact amplification factor \(\left|G_{exact}\right|\ ,\) over one time step. The relative numerical diffusion error or relative numerical dissipation error compares real physical dissipation with the anomalous dissipation that results from numerical discretization. It can be defined as

\[\tag{56} \varepsilon_{D}=\frac{\left|G_{numeric}\right|}{\left|G_{exact}\right|},\]

and the total dissipation error resulting from \(n\) steps will be

\[\tag{57} \varepsilon_{Dtotal}=\left(\left|G_{numeric}\right|^{n}-\left|G_{exact}\right|^{n}\right)u_{0}.\]

If \(\varepsilon_{D}>1\) for a given value of \(\theta\) or Co, this discretization scheme will be unstable and a modification to the scheme will be necessary.

As mentioned above, if the imaginary part of \(\omega\left(k\right)\) is zero for a particular discretization, then the scheme is non-dissipative.

Numerical dispersion

Figure 3: Figure 3: Illustration of pure numeric dispersion effect on a single sinusoid, as it propagates along the spatial domain. Both exact and simulated dispersive waves start in phase; however, the phase of the dispersive wave lags the exact wave over time, but its amplitude is unaffected.

Figure 4: Figure 4: Effect of numerical dispersion on a step function applied to the advection equation \(u_{t}+u_{x}=0\ .\)

In a numerical scheme, a situation where waves of different frequencies move at different speeds without a change in amplitude, is called_numerical dispersion_ - see figure (3). Alternatively, the Fourier components of a wave can be considered to disperse relative to each other. It therefore follows that the effect of a dispersive scheme on a wave composed of different harmonics, will be to deform the wave as it propagates. However the energy contained within the wave is not lost and travels with the_group velocity_. Generally, this results in higher frequency components traveling at slower speeds than the lower frequency components. The effect of dispersion therefore is that often spurious oscillations or wiggles occur in solutions with sharp gradient, discontinuity or shock effects, usually with high frequency oscillations trailing the particular effect, see figure (4). The degree of dispersion can be determined by comparing the phase of the numeric amplification factor \(\left|G_{numeric}\right|,\) with the phase of the exact amplification factor \(\left|G_{exact}\right|\ ,\) over one time step. Dispersion represents phase shift and results from the imaginary part of the amplification factor. The relative numerical dispersion error compares real physical dispersion with the anomalous dispersion that results from numerical discretization. It can be defined as

\[\tag{58} \varepsilon_{P}=\frac{\Phi_{numeric}}{\Phi_{exact}},\]

where \(\Phi=\angle G=\tan^{-1}\left(\frac{\textrm{Im}\left\{ G\right\} }{\textrm{Re}\left\{ G\right\} }\right)\ .\) The total phase error resulting from \(n\) steps will be

\[\tag{59} \varepsilon_{Ptotal}=n\left(\Phi_{numeric}-\Phi_{exact}\right)\]

If \(\varepsilon_{P}>1\ ,\) this is termed a leading phase error. This means that the Fourier component of the solution has a wave speed_greater_ than the exact solution. Similarly, if \(\varepsilon_{P}<1\ ,\) this is termed a lagging phase error. This means that the Fourier component of the solution has a wave speed less than the exact solution.

Again, high resolution schemes can all but eliminate this effect, but at the expense of complexity. Although many physical processes are modeled by PDE's that are non-dispersive, when numerical discretization is applied to analyze them, some dispersion is usually introduced.

Group and phase velocity

The term group velocity refers to a wave packet consisting of a low frequency signal modulated (or multiplied) by a higher frequency wave. The result is a low frequency wave, consisting of a fundamental plus harmonics, that propagates with group velocity \(c_{g}\) along a continuum oscillating at a higher frequency. _Wave energy_and information signals propagate at this velocity, which is defined as being equal to the derivative of the real part of the frequency \(\omega\ ,\) with respect to wavenumber \(k\) (scalar or vector proportional to the number of wave lengths per unit distance), i.e. \[\tag{60} c_{g}=\frac{d\,\textrm{Re}\left\{ \omega\left(k\right)\right\} }{dk}.\]

If there are a number of spatial dimensions then the group velocity is equal to the gradient of frequency with respect to the wavenumber vector, i.e. \(c_{g}=\nabla\textrm{Re}\left\{ \omega\left(k\right)\right\} \ .\)

The complementary term to group velocity is phase velocity, \(c_{p}\ ,\) and this refers to the speed of propagation of an individual frequency component of the wave. It is defined as being equal to the real part of the ratio of frequency to wavenumber, i.e.\[\tag{61} c_{p}=\textrm{Re}\left\{ \frac{\omega}{k}\right\} .\]

It can also be viewed as the speed at which a particular phase of a wave propagates; for example, the speed of propagation of a wave crest. In one wave period \(T\) the crest advances one wave length \(\lambda\ ;\) therefore, the phase velocity is also given by \(c_{p}=\lambda/T\ .\) We see that this second form is equal to equation (61) due to the following relationships: wavenumber \(k=\frac{2\pi}{\lambda}\) and frequency \(\omega=2\pi f\) where \(f=\frac{1}{T}\ .\)

For a non-dispersive wave the phase error is zero and therefore \(c_{g}=c_{p}\ .\)

To calculate group and phase velocity for linear waves (or small amplitude waves) we assume a solution of the form \(u(x,t)=Ae^{i(kx-\omega t)}\ ,\) where \(A\) is a constant and \(x\) can be a scalar or vector, and substitute into the wave equation (or linearized wave equation) under consideration. For example, for \(u_{t}+u_{x}+u_{xxx}=0\) we obtain the dispersion relation \(\omega=k-k^{3}\ ,\) from which we calculate the group and phase velocities to be \(c_{g}=1-3k^{2}\) and \(c_{p}=1-k^{2}\) respectively. Thus, we observe that \(c_{g}\neq c_{p}\) and that therefore, this example is dispersive.

Wellposedness

For most practical situations our interest is primarily in solving partial differential equations numerically; and, before we embark on implementing a numerical procedure, we would usually like to have some idea as to the expected behaviour of the system being modeled, ideally from an analytical solution. However, an analytical solution is not usually available; otherwise we would not need a numerical solution. Nevertheless, we can usually carry out some basic analysis that may give some idea as to steady state, long term trend, bounds on key variables, and reduced order solution for ideal or special conditions, etc. One key estimate that we would like to know is whether the fundamental system is stable or well posed. This is particularly important because if our numerical solution produces seemingly unstable results we need to know if this is fundamental to the problem or whether it has been introduced by the solution method we have selected to implement. For most situations involving simulation this is not a concern as we would be dealing with a well analyzed and documented system. But there are situations where real physical systems can be unstable and we need to know these in advance. For a real system to become unstable there needs to be some form of energy source: kinetic, potential, reaction, etc., so this can provide a clue as to whether or not the system is likely to become unstable. If it is, then we may need to modify our computational approach so that we capture the essential behaviour correctly - although a complete solution may not be possible.

In general, solutions to PDE problems are sought to solve a particular problem or to provide insight into a class of problems. To this end_existence_, uniqueness and stability of the solution are of vital importance [Zwi-97, chapter 10]. Whilst at this introductory level we must restrict our discussion, it is desirable to emphasize that for a solution of an evolutionary PDE (together with appropriate ICs and BCs) to be useful we require that:

If these conditions are full-filled, then the problem is said to be_well posed_, in the sense of Hadamard [Had-23]. Numerical schemes for particular PDE systems can be analyzed mathematically to determine if the solutions remain bounded. By invoking Parseval's theorem of equality this analysis can be performed in the time domain or in the Fourier domain. A good introduction to this subject is given by LeVeque [Lev-07], and more advanced technical discussions can be found in the monographs by Tao [Tao-05] and Kreiss & Lorenz [Kre-04].

Characteristics

Characteristics are surfaces in the solution space of an evolutionary PDE problem that represent wave-fronts upon which information propagates. For example, consider the 1D advection equation problem

\[\tag{62} u_{t}=cu_{x},\quad u\left(x,t=0\right)=u_{0},\; t\geq0\]

where the characteristics are given by \(dx/dt=c\ .\) For this problem the characteristics are straight lines in the \(xt\)-plane with slope \(1/c\) and, along which, the dependent variable \(u\) is constant. The consequence of this is that the initial condition propagates from left to right at constant speed \(c\ .\) But, for other situations such as the inviscid Burgers equation problem,

\[\tag{63} u_{t}=uu_{x},\quad u\left(x,t=0\right)=u_{0},\; t\geq0,\]

the propagation speed is not constant and the shape of the characteristics depend upon the initial conditions. If the initial condition is monotonically increasing with \(x\ ,\) the characteristics will not overlap and the problem is well behaved. However, if the initial conditions are not monotonically increasing with \(x\ ,\) at some time \(t>0\) the characteristics will overlap and the solution will become multi-valued and a shock will develop. In this situation we can only find a weak solution(one where the problem is re-stated in integral form) by appealing to entropy considerations and the Rankine-Hugoniot jump condition. PDEs other than equations (62) and (63), such as those involving conservation laws, introduce additional complexity such as _rarefaction_or expansion waves. We will not discuss these aspects further here, and for additional discussion readers are referred to [Hir-90, chap. 16].

The method of characteristics

The method of characteristics (MOC) is a numerical method for solving evolutionary PDE problems by transforming them into a set of ODEs. The ODEs are solved along particular characteristics, using standard methods and the initial and boundary conditions of the problem. For more information refer to [Kno-00],[Ost-94],[Pol-07].

MOC is a quite general technique for solving PDE problems and has been particularly popular in the area of fluid dynamics for solving incompressible transient flow in pipelines. For an introduction refer to [Stre-97, chap. 12].

General topics

We conclude with a brief overview of some general aspects relating to linear and nonlinear waves.

Galilean invariance

Certain wave equations are Galilean invariant, i.e. the equation properties remain unchanged under a Galilean transformation. For example:

\[\tag{64} \tilde{u}=Au\left(\pm\lambda x+C_{1},\pm\lambda t+C_{2}\right),\]

where \(A\ ,\) \(C_{1}\ ,\) \(C_{2}\) and \(\lambda\) are arbitrary constants.

\[\tag{65} \tilde{u}=u\left(x-6\lambda t,t\right)-\lambda ,\]

where \(\lambda\) is an arbitrary constant.

Other invariant transformations are possible for many linear and nonlinear wave equations, for example the Lorentz transformation applied to Maxwell's equations, but these will not be discussed here.

Plane waves

Figure 5: Figure 5: Plane sinusoidal wave where its source is assumed to be at \( x = -\infty\ ,\) and its fronts are advancing from right to left.

A Plane wave is considered to exist far from its source and any physical boundaries so, effectively, it is located within an infinite domain. Its position vector remains perpendicular to a given plane and satisfies the 1D wave equation \[\tag{66} \frac{1}{c^2}\frac{\partial^{2}u}{\partial t^{2}}=\frac{\partial^{2}u}{\partial x^{2}}\]

with a solution of the form\[\tag{67} u=u_{0}\cos\left(\omega t-kx+\phi\right)\]

where \(c=\frac{\omega}{k}\) represents propagation velocity and \(\phi\) the phase of the wave. See figure (5).

Refraction and diffraction

Wave crests do not necessarily travel in a straight line as they proceed - this may be caused by refraction or diffraction.

Wave refraction is caused by segments of the wave moving at different speeds resulting from local changes in characteristic speed, usually due to a change in medium properties. Physically, the effect is that the overall direction of the wave changes, its wavelength either increases or decreases but its frequency remains unchanged. For example, in optics refraction is governed by Snell's law and in shallow water waves by the depth of water.

Wave diffraction is the effect whereby the direction of a wave changes as it interacts with objects in its path. The effect is greatest when the size of the object causing the wave to diffract is similar to the wavelength.

Reflection

Reflection results from a change of wave direction following a collision with a reflective surface or domain boundary. A _hard boundary_is one that is fixed which causes the wave to be reflected with opposite polarity, e.g. \(u(x-vt)\;\rightarrow\;-u(x+vt)\ .\) A _soft boundary_is one that changes on contact with the wave, which causes the wave to be reflected with the same polarity, e.g. \(u(x-vt)\;\rightarrow\; u(x+vt)\ .\) If the propagating medium is not isotropic, i.e. it is not spatially uniform, then a partial reflection can result, with an attenuated original wave continuing to propagate. The polarity of the partial reflection will depend upon the characteristics of the medium.

Consider a travelling wave situation where the domain has a soft boundary with incident wave \(\phi_{I}=I\exp\left(j\omega t-k_{1}x\right)\ ,\)reflected wave \(\phi_{R}=R\exp\left(j\omega t+k_{1}x\right)\) and transmitted wave \(\phi_{T}=T\exp\left(j\omega t-k_{2}x\right)\ .\) In addition, for simplicity, consider the medium on both sides of the boundary to be isotropic and non-dispersive, which implies that all three waves will have the same frequency. From the conservation of energy law we have \(\phi_{I}+\phi_{R}=\phi_{T}\) for all \(t\ ,\) which implies \(I+R=T.\) Also, on differentiating with respect to \(x\ ,\) we obtain \(-ik_{1}I+ik_{1}R=-ik_{2}T\ .\) Thus, on rearranging we have

\[\tag{68} \frac{T}{I} = \frac{2k_{1}}{k_{1}+k_{2}},\]

\[\tag{69} \frac{R}{I} = \frac{k_{1}-k_{2}}{k_{1}+k_{2}}.\]

Equations (68) and (69) indicate that:

Also, because \(c_{g}=c_{p}={\displaystyle \frac{\omega}{k}},\;\) if \(\;k_{1}>k_{2}\), then this implies that \(c_{g1}<c_{g2},\;\) see section (Group and phase velocity).

We mention two other quantities

\[\tag{70} \tau = \left|\frac{T}{I}\right| ,\]

\[\tag{71} \rho = \left|\frac{R}{I}\right| ,\]

the so-called coefficients of transmission and reflection respectively.

Resonance

Resonance describes a situation where a system oscillates at one of its natural frequencies, usually when the amplitude increases as a result of energy being supplied by a perturbing force. A striking example of this phenomena is the failure of the mile-long Tacoma Narrows Suspension Bridge. On 7 November 1940 the structure collapsed due to a nonlinear wave that grew in magnitude as a result of excitation by a 42 mph wind. A video of this disaster is available on line at:archive.org . Another less dramatic example of resonance that most people have experienced is the effect of sound feedback from loudspeaker to microphone.

A more complex form of resonance is autoresonance, a nonlinear phase-locking phenomenon which occurs when a resonantly driven nonlinear system becomes phase-locked (synchronized or in-step) with a driving perturbation or wave.

Doppler effect

The Doppler effect (or Doppler shift) relates to the change in frequency and wavelength of waves emitted from a source as perceived by an observer, where the source and observer are moving at a speed relative to each other. At each moment of time the source will radiate a wave and an observer will experience the following effects:

Perhaps the most famous discovery involving the Doppler effect, is that made in 1929 by Edwin Hubble in connection with the Earth's distance from receding galaxies: the redshift of light coming from distant galaxies is proportional to their distance. This is known as Hubble's law.

Transverse and longitudinal waves

Transverse waves oscillate in the plane perpendicular to the direction of wave propagation. They include: seismic S (secondary) waves, and electromagnetic waves, E (electric field) and H (magnetic field), both of which oscillate perpendicularly to each other as well to the direction of propagation of energy. Light, an electromagnetic wave, can be polarized (oriented in a specific direction) by use of a polarizing filter.

Longitudinal waves oscillate along the direction of wave propagation. They include sound waves (pressure, particle displacement, or particle velocity propagated in an elastic medium) and _seismic_P (earthquake or explosion) waves.

Surface water waves however, are an example of waves that involve a combination of both longitudinal and transverse motion.

Traveling waves

Traveling-wave solutions [Pol-08], [Gri-11], by definition, are of the form

\[\tag{72} u(x,t)=U(z),\quad z=kx-\lambda t\ ;\]

where \(\lambda/k\) plays the role of the wave propagation velocity (the value \(\lambda=0 \,\) corresponds to a stationary solution, and the value \(k=0 \,\) corresponds to a space-homogeneous solution).

Traveling-wave solutions are characterized by the fact that the profiles of these solutions at different time instants are obtained from one another by appropriate shifts (translations) along the \(\, x\)-axis. Consequently, a Cartesian coordinate system moving with a constant speed can be introduced in which the profile of the desired quantity is stationary. For \(\lambda>0 \,\) and \(k>0\ ,\) the wave described by equation (72) travels along the \(x\)-axis to the right (in the direction of increasing \(x \,\)).

The term traveling-wave solution is also used in situations where the variable \(t \,\) plays the role of a spatial coordinate, \(y \,\ .\)

Standing waves

Figure 6: Figure 6: A standing wave\[\Re \left( \phi\left(x,t\right) \right)\ .\]

_Standing waves' occur when two traveling waves of equal amplitude_and speed, but opposite direction, are superposed. The effect is that the wave amplitude varies with time but it does not move spatially. For example, consider two waves \(\phi_{1}\left(x,t\right)=\Phi_{1}\exp i\left(\omega t-kx\right)\) and \(\phi_{2}\left(x,t\right)=\Phi_{2}\exp i\left(\omega t+kx\right)\ ,\) where \(\phi_{1}\) moves to the right and \(\phi_{2}\) moves to the left. By definition we have \(\Phi_{2}=\Phi_{2}\ ,\) and by simple algebraic manipulation we obtain \[ \phi\left(x,t\right) = \phi_{1}\left(x,t\right)+\phi_{2}\left(x,t\right) ,\] \[ = \Phi_{1}\left[\exp i\left(\omega t-kx\right)+\exp i\left(\omega t+kx\right)\right] ,\] \[ \Downarrow \]\[\tag{73} = 2\Phi_{1}\exp i\omega t\;\cos kx.\]

A standing wave is illustrated in figures (6) and (7) by a plot of the real part of equation (73), i.e. \(\Re \left( \phi\left(x,t\right) \right)=2\Phi_{1}\cos\omega t\cos kx\) with \(k=1\ ,\) \(\omega=1\) and \(\Phi_{1}=\frac{1}{2}\ .\)

Figure 7: Figure 7: Animated standing wave\[\Re \left( \phi\left(x,t\right) \right)\ .\]

The points at which \(\phi=0\) are called nodes and the points at which \(\phi=2\left|\Phi\right|\) are called antinodes. These points are fixed and occur at \(kx=\left(2n+1\right)\frac{ {\displaystyle \pi}}{{\displaystyle 2} }\) and \(kx=\left(2n+1\right)\pi\) respectively \(\left(n=\pm1,\pm2,\cdots\right)\ .\)

Clearly, the existence of nonlinear standing waves can be demonstrated by application of Fourier analysis.

Waveguides

The idea of a waveguide is to constrain a wave such that its energy is directed along a specific path. The path may be fixed or capable of being varied to suit a particular application.

The operation of a waveguide is analyzed by solving the appropriate wave equation, subject to the prevailing boundary conditions. There will be multiple solutions, or modes, which are determined by the eigenfunctions associated with the particular wave equations, and the velocity of the wave as it propagates along the waveguide will be determined by the eigenvalues of the solution.

For detailed analysis and further discussion refer to [Lio-03],[Oka-06].

Wave-front

Figure 8: Figure 8: Circular wave-fronts emanating from a point source.

As a wave propagates through a medium, the wavefront represents the outward normal surface formed by points in space, at a particular instant, as the wave travels outwards from its origin.

One of the simplest form of wavefront to envisage is an expanding circle where its radius \(r\ ,\) expands with velocity \(v\ ,\) i.e. \(r=vt.\) Simple circular sinusoidal wave-fronts propagating from a point source are shown in figure (8). They can be described by\[\tag{74} u\left(r,t\right)=Re\left\{ \exp\left[i(kr-\omega t+\pi/2-\psi\right]\right\} ,\]

where \(k=\) wavenumber, \(r=\textrm{radius}\ ,\) \(t=\textrm{time}\ ,\) \(\omega=\textrm{frequency}\ ,\) \(\psi=\textrm{phase angle}\ .\)

Figure 9: Figure 9: A snapshot from a simulation of the Indian Ocean tsunami that occurred on 26th December 2004 resulting from an earthquake off the west coast of Sumatra. The non-circular wave-fronts are clearly visible, which indicates curved rays. See animation here

.

Depending upon the particular wave equation and medium in which the wave travels, the wavefront may not appear to be an expanding circle. The path upon which any point on the wave front has traveled is called a ray, and this can be a straight line or, more likely, a curve in space. In general, the wavefront is perpendicular to the ray path, and the ray curvature will depend on the circumstances of the particular physical situation. For example, its curvature will be influenced by: an anisotropic medium, refraction, diffraction, etc.

Consider a water wave where wave height is very much smaller than water depth \(h\ .\) Its speed of propagation \(c\ ,\) or _celerity,_is given by \(c=\sqrt{gh}\ ;\) thus, for an ocean with varying depth the velocity will vary at different locations (refraction). This can result in waves having non-circular wave-fronts and hence curved rays. This situation, which occurs in many different applications, is illustrated in figure (9) where the curved wave-fronts are due to a combination of effects due to refraction, diffraction, reflection and a non-point disturbance.

Huygens' principle

Figure 10: Figure 10: Advancing envelope of wave-fronts \(\Phi_{q_{0}}\left(t\right)\ .\)

We can consider all points of a wave-front of light in a vacuum or transparent medium to be new sources of wavelets that expand in every direction at a rate depending on their velocities.

This idea was originally proposed by the Dutch mathematician, physicist, and astronomer, Christiaan Huygens, in 1690, and is a powerful method for studying various optical phenomena[Enc-09]. Thus, the points on a wave can be viewed as each emitting a circular wave which combine to propagate the wave-front \(\Phi_{q_{0}}\left(t\right)\ .\) The wave-front can be thought of as an advancing line tangential to these circular waves - see figure (10). The points on a wave-front propagate from the wave source along so-called_rays_. The Huygens' principle applies generally to wave-fronts and the laws of reflection and refraction can both be derived from Huygens' principle. These results can also be obtained from Maxwell's equations.

For detailed analysis and proof of Huygens' principle, refer to [Arn-91].

Shock waves

There are an extremely large number of types and forms of shock wave phenomena, and the following are representative of some subject areas where shocks occur:

For further discussion relating to shock phenomena see ([Ben-00],[Whi-99]).

We briefly introduce two topics below by way of example.

Blast Wave - Sedov-Taylor Detonation

A blast wave can be analyzed from the following equations, \[\tag{75} \frac{\partial\rho}{\partial t}+v\frac{\partial\rho}{\partial r}+\rho\left(\frac{\partial v}{\partial r}+\frac{2v}{r}\right) = 0,\]

\[\tag{76} \quad \frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}+\frac{1}{\rho}\frac{\partial p}{\partial r} = 0,\]

\[\tag{77} \quad \frac{\partial\left(p/\rho^{\gamma}\right)}{\partial t}+v\frac{\partial\left(p/\rho^{\gamma}\right)}{\partial r} = 0,\]

where \(\rho\ ,\) \(v\ ,\) \(p\ ,\) \(r\ ,\) \(t\) and \(\gamma\) represent _density_of the medium in which the blast takes place (air), velocity of the blast front, blast pressure, blast radius, time and isentropic exponent (ratio of specific heats) of the medium respectively. Now, if we assume that,

Figure 11: Figure 11: Time-lapse photographs with distance scales (100 m) of the first atomic bomb explosion in the New Mexico desert - 5.29 A.M. on 16th July, 1945. Times from instant of detonation are indicated in bottom left corner of each photograph (Top first - left column: 0.006s, 0.016s; right column: 0.025s, 0.09s).

then, after some analysis, similarity considerations lead to the following equation [Tay-50b] \[\tag{78} E=c{\displaystyle \frac{R^{5}\rho}{t^{2}}},\]

where \(c\) is a similarity constant, \(R\) is the radius of the wave front and \(E\) is the total energy released by the explosion.

Back in 1945 Sir Geoffrey Ingram Taylor was asked by the British MAUD (Military Application of Uranium Detonation) Committee to deduce information regarding the power of the first atomic explosion in New Mexico. He derived this result, which was based on his earlier classified work[Tay-41], and was able to estimate, using only photographs of the blast (released into the public domain in 1947), that the yield of the bomb was equivalent to between \(16.8\) and \(22.9\) kilo-tons of TNT for values of \(\gamma\) equal to \(1.4\) and \(1.3\) respectively. Each of these photographs, crucially, contained a distance scale and precise time, see figure (11). Taylor used a value for the similarity constant of \(c=0.856\) that he obtained by a step-by-step method. However the correct analytical value for this constant was later shown to be \(0.8501\) [Sed-59].

This result was classified secret but, five years later he published the details [Tay-50a],[Tay-50b], much to the consternation of the British government. J. von Neumann and L. I. Sedov published similar independently derived results [Bet-47],[Sed-46]. For further discussion relating to the theory refer to [Kam-00],[Deb-58].

Sonic boom

Figure 12: Figure 12: The N-wave sonic boom.

As an aircraft proceeds in smooth flight at a speed greater than the speed of sound - the sound barrier - a shock wave is formed at its nose and finishes at its tail. The speed of sound is given by \(c=\sqrt{\gamma RT/MW}\ ,\) where \(\gamma\ ,\) \(R\ ,\) \(T\) and MW represent ratio of specific heats, universal gas constant, temperature and molecular weight respectively, and \(c\simeq330\)m/s at sea level for dry air at \(0^{o}\)C . The shock forms a high pressure, cone-shaped surface propagating with the aircraft. The half-angle (between direction of flight and the shock wave) \(\theta\) is given by \(\sin\left(\theta\right)=1/M\ ,\) where \(M=v_{aircraft}/c\) is known as the Mach number of the aircraft. Clearly, as \(v_{aircraft}\) increases, the cone becomes more pointed (\(\, \theta\) becomes smaller).

Figure 13: Figure 13: The U-wave sonic boom.

As the aircraft continues under steady flight conditions at high speed, there will be an abrupt rise in pressure at the aircraft's nose, which falls towards the tail when it then becomes negative. This is the so-called N-wave [Nak-08] - a pressure wave measured at sufficient distance such that it has lost its fine structure, see figure (12). A sonic boom occurs when the abrupt changes in pressure are of sufficient magnitude. Thus, steady supersonic flight results in two booms: one resulting from the rapid rise in pressure at the nose, and another when the pressure returns to normal as the tail passes the point vacated by the nose. This is the cause of the distinctive _double boom_from supersonic aircraft. At ground level typically, \(10<P_{max}<500\)Pa and \(\tau\simeq0.001-0.005\)s. The duration \(T\) varies from around 100 ms for a fighter plane to 500 ms for the Space Shuttle or Concorde.

Figure 14: Figure 14:A USAF B1B makes a high speed pass at the Pensacola Beach airshow - Florida, July 12, 2002. Copyright © Gregg Stansbery, Stansbery Photography - reproduced with permission.

Another form of sonic boom is the focused boom. These can result from high speed aircraft flight maneuvering operations. These result in the so-called U-waves which have positive shocks at the front and rear of the boom, see figure (13). Generally, U-waves result in higher peak over-pressures than N-waves - typically between 2 and 5 times. At ground level typically, \(20<P_{max}<2500\)Pa (although they can be much higher). The highest overpressure ever recorded was 6800 Pa [144 lbs/sq-ft] (source: USAF Fact Sheet 96-03). For further discussion related to sonic booms refer to [Kao-04].

As an aircraft passes through, or close to the sound barrier, water vapor in the air is compressed by the shock wave and becomes visible as a large cloud of condensation droplets formed as the air cools due to low pressure at the tail. A smaller shock wave can also form on top of the canopy. This phenomena is illustrated in figure (14).

Solitary waves and solitons

The correct term for a wave which is localized and retains its form over a long period of time is: solitary wave. However, a _soliton_is a solitary wave having the additional property that other solitons can pass through it without changing its shape. But, in the literature it is customary to refer to the solitary wave as a soliton, although this is strictly incorrect [Tao-08].

Figure 15: Figure 15: Evolution of a two-soliton solution of the KdV equation. Image illustrates the collision of two solitons that are both moving from left to right. The faster (taller) soliton overtakes the slower (shorter) soliton.

Solitons are stable, nonlinear pulses which exhibit a fine balance between non-linearity and dispersion. They often result from real physical phenomena that can be described by PDEs that are completely integrable, i.e. they can be solved exactly. Such PDEs describe: shallow water waves, nonlinear optics, electrical network pulses, and many other applications that arise in mathematical physics. Where multiple solitons moving at different velocities occur within the same domain, collisions can take place with the unexpected phenomenon that, first they combine, then the faster soliton emerges to proceed on its way. Both solitons then continue to proceed in the the same direction and eventually reach a situation where their speeds and shapes are unchanged. Thus, we have a situation where a faster soliton can overtake a slower soliton. There are two effects that distinguishes this phenomena from that which occurs in a linear wave system. The first is that the maximum height of the combined solitons is not equal to the sum of the individual soliton heights. The second is that, following the collision, there is a phase shift between the two solitons, i.e. the linear trajectory of each soliton before and after the collision is seen to be shifted horizontally - see figure (15).

Some additional discussion is given in section (The Korteweg-de Vries equation) and detailed technical overviews of the subject can be found in the works by Ablowitz & Clarkson [Abl-91], Drazin & Johnson [Dra-89] and Johnson [Joh-97]. Soliton theory is still an active area of research and a discussion on the various types of soliton solution that are known is given by Gerdjikov & Kaup [Ger-05].

Soliton types

Soliton types generally fall into thee types:

More details may be found in Drazin and Johnson [Dra-89].

Tsunami

The word tsunami is a Japanese term derived from the characters 津 (tsu) meaning harbor and 波 (nami) meaning wave. It is now generally accepted by the international scientific community to describe a series of traveling waves in water produced by the displacement of the sea floor associated with submarine earthquakes, volcanic eruptions, or landslides. They are also known as tidal waves.

Tsunami are usually preceded by a leading-depression N-wave(LDN), one in which the trough reaches the shoreline first. Eyewitnesses in Banda Aceh who observed the effects of the December 2004_Sumatra Tsunami_, see figure (9), resulting from a magnitude 9.3 seabed earthquake, described a series of three waves, beginning with a leading depression N wave [Bor-05]. Recent estimates indicate that this powerful tsunami resulted in excess of 275,000 deaths and extensive damage to property and infrastructure around the entire coast line of the Indian ocean [Kun-07].

Tsunami are long-wave phenomena and, because the wavelengths of tsunami in the ocean are long with respect to water depth, they can be considered shallow water waves. Thus, \(c_{p}=c_{g}=\sqrt{gh}\) and for a depth of 4km we see that the wave velocity is around 200 m/s. Hence, tsunami waves are often modelled using the shallow water equations, the Boussinesq equation, or other suitable equations that bring out in sufficient detail the required wave characteristics. However, one of the major challenges is to model shoreline _inundation_realistically, i.e. the effect of the wave when it encounters the shore - also known as run-up. As the wave approaches the shoreline, the water depth decreases sharply resulting in a greatly increased surge of water at the point where the wave strikes land. This requires special modeling techniques to be used, such as robust Riemann solvers [Tor-01],[Ran-06] or the _level-set_method [Set-99],[Osh-03], which can handle situations where dry regions become flooded and vice versa.

Acknowledgements

The authors would like to thank reviewers Prof. Andrei Polyanin and Dr. Alexei Zhurov for their positive and constructive comments.

References

Internal references