WikiWaves - User contributions [en]

Category:Nonlinear PDE's Course

2026-04-13T06:21:12Z

Meylan:

These pages contain notes for the course in Nonlinear PDE's taught by Mike Meylan

# [[Method of Characteristics for Linear Equations]]
# [[Traffic Waves]]
# [[Nonlinear Shallow Water Waves]]
# [[Introduction to KdV]]
# [[Numerical Solution of the KdV]]
# [[Conservation Laws for the KdV]]
# [[Introduction to the Inverse Scattering Transform]]
# [[Properties of the Linear Schrodinger Equation]]
# [[Connection betwen KdV and the Schrodinger Equation]]
# [[Example Calculations for the KdV and IST]]
# [[Reaction-Diffusion Systems]]
# [[Burgers Equation]]

{| class="wikitable"
|+ 2026 course
! Week !! Date !! Topic !! Assessment
|-
| Week 1 || 2 Mar 2026 || [[Method of Characteristics for Linear Equations]] ||
|-
| Week 2 || 9 Mar 2026 || [[Traffic Waves]] ||
|-
| Week 3 || 16 Mar 2026 || [[Nonlinear Shallow Water Waves]] ||
|-
| Week 4 || 23 Mar 2026 || [[Introduction to KdV]] ||
|-
| Week 5 || 30 Mar 2026 || [[Numerical Solution of the KdV]] || Assignment 1
|-
| Week 6 || 6 Apr 2026 || [[Conservation Laws for the KdV]] ||
|-
| Week 7 || 13 Apr 2026 || [[Introduction to the Inverse Scattering Transform]] ||
|-
| Week 8 || 20 Apr 2026 || [[Properties of the Linear Schrodinger Equation]] ||
|-
| Week 9 || 27 Apr 2026 || No lecture ||
|-
| Week 9 || 4 May 2026 || [[Connection betwen KdV and the Schrodinger Equation]] ||
|-
| Week 10 || 11 May 2026 || [[Example Calculations for the KdV and IST]] ||
|-
| Week 11 || 18 May 2026 || [[Reaction-Diffusion Systems]] ||
|-
| Week 12 || 25 May 2026 || [[Burgers Equation]] || Assignment 2
|-
| Week 13 || 1 Jun 2026 || No topic ||
|-
| Week 14 || 8 Jun 2026 || No topic || Assignment 3
|-
| Week 15 || 15 Jun 2026 || No topic ||
|-
| Week 16 || 22 Jun 2026 || No topic || Final Exam
|}

File:Characteristic linear1.jpg

2026-03-23T04:32:37Z

Meylan:

Code to make can be found [http://www.math.auckland.ac.nz/~meylan/code/nonlinear/linear_characteristics.m linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-23T04:32:12Z

Meylan:

Code to make can be found http://www.math.auckland.ac.nz/~meylan/code/nonlinear/linear_characteristics.m linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-22T21:59:51Z

Meylan:

Code to make can be found [https://www.dropbox.com/scl/fi/dalyshvrswogg9mznth1j/linear_characteristics.m?rlkey=7v0wx8jk1hw4nfd7fap3wvsn1&st=bqfzmw0x&dl=0 linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-22T21:53:14Z

Meylan:

Code to make can be found [https://www.dropbox.com/scl/fi/dalyshvrswogg9mznth1j/linear_characteristics.m?rlkey=7v0wx8jk1hw4nfd7fap3wvsn1&st=rr1glcls&dl=0 linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-22T21:47:52Z

Meylan:

Code to make can be found [https://www.dropbox.com/scl/fi/dalyshvrswogg9mznth1j/linear_characteristics.m?rlkey=fgvnj0sv2o2fbcg9b72454djz&st=htmx6k9n&dl=0 linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-22T21:45:19Z

Meylan:

Code to make can be found [https://www.dropbox.com/scl/fi/cwyjexugz76puve3y98w3/phase_portrait_burgers.m?rlkey=tbye33u8q1ho78qhqj3h5zdkh&st=8x225n66&dl=0 linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-22T21:40:42Z

Meylan:

Code to make can be found [https://www.dropbox.com/scl/fi/dalyshvrswogg9mznth1j/linear_characteristics.m?rlkey=fgvnj0sv2o2fbcg9b72454djz&st=caom7qo6&dl=0 linear_characteristic.m]

File:Characteristic linear1.jpg

2026-03-22T21:39:43Z

Meylan:

Code to make can be found [https://www.dropbox.com/scl/fi/dalyshvrswogg9mznth1j/linear_characteristics.m?rlkey=fgvnj0sv2o2fbcg9b72454djz&st=7kfgors7&dl=0]

Category:Nonlinear PDE's Course

2026-03-02T06:04:58Z

Meylan:

These pages contain notes for the course in Nonlinear PDE's taught by Mike Meylan

# [[Method of Characteristics for Linear Equations]]
# [[Traffic Waves]]
# [[Nonlinear Shallow Water Waves]]
# [[Introduction to KdV]]
# [[Numerical Solution of the KdV]]
# [[Conservation Laws for the KdV]]
# [[Introduction to the Inverse Scattering Transform]]
# [[Properties of the Linear Schrodinger Equation]]
# [[Connection betwen KdV and the Schrodinger Equation]]
# [[Example Calculations for the KdV and IST]]
# [[Reaction-Diffusion Systems]]
# [[Burgers Equation]]

{| class="wikitable"
|+ 2026 course
! Week !! Date !! Topic !! Assessment
|-
| Week 1 || 2 Mar 2026 || [[Method of Characteristics for Linear Equations]] ||
|-
| Week 2 || 9 Mar 2026 || [[Traffic Waves]] ||
|-
| Week 3 || 16 Mar 2026 || [[Nonlinear Shallow Water Waves]] ||
|-
| Week 4 || 23 Mar 2026 || [[Introduction to KdV]] ||
|-
| Week 5 || 30 Mar 2026 || [[Numerical Solution of the KdV]] || Assignment 1
|-
| Week 6 || 6 Apr 2026 || [[Conservation Laws for the KdV]] ||
|-
| Week 7 || 13 Apr 2026 || [[Introduction to the Inverse Scattering Transform]] ||
|-
| Week 8 || 20 Apr 2026 || [[Properties of the Linear Schrodinger Equation]] ||
|-
| Week 9 || 27 Apr 2026 || [[Connection betwen KdV and the Schrodinger Equation]] ||
|-
| Week 10 || 4 May 2026 || [[Example Calculations for the KdV and IST]] ||
|-
| Week 11 || 11 May 2026 || [[Reaction-Diffusion Systems]] ||
|-
| Week 12 || 18 May 2026 || [[Burgers Equation]] || Assignment 2
|-
| Week 13 || 25 May 2026 || No topic ||
|-
| Week 14 || 1 Jun 2026 || No topic || Assignment 3
|-
| Week 15 || 8 Jun 2026 || No topic ||
|-
| Week 16 || 15 Jun 2026 || No topic || Final Exam
|}

Category:Nonlinear PDE's Course

2026-03-02T06:02:31Z

Meylan:

These pages contain notes for the course in Nonlinear PDE's taught by Mike Meylan

# [[Method of Characteristics for Linear Equations]]
# [[Traffic Waves]]
# [[Nonlinear Shallow Water Waves]]
# [[Introduction to KdV]]
# [[Numerical Solution of the KdV]]
# [[Conservation Laws for the KdV]]
# [[Introduction to the Inverse Scattering Transform]]
# [[Properties of the Linear Schrodinger Equation]]
# [[Connection betwen KdV and the Schrodinger Equation]]
# [[Example Calculations for the KdV and IST]]
# [[Reaction-Diffusion Systems]]
# [[Burgers Equation]]

{| class="wikitable"
|+ 2026 course
! Week !! Date !! Topic !! Assessment
|-
| Week 1 || 2/3/26 || [[Method of Characteristics for Linear Equations]] ||
|-
| Week 2 || 2/10/26 || [[Traffic Waves]] ||
|-
| Week 3 || 2/17/26 || [[Nonlinear Shallow Water Waves]] ||
|-
| Week 4 || 2/24/26 || [[Introduction to KdV]] ||
|-
| Week 5 || 3/3/26 || [[Numerical Solution of the KdV]] || Assignment 1
|-
| Week 6 || 3/10/26 || [[Conservation Laws for the KdV]] ||
|-
| Week 7 || 3/17/26 || [[Introduction to the Inverse Scattering Transform]] ||
|-
| Week 8 || 3/24/26 || [[Properties of the Linear Schrodinger Equation]] ||
|-
| Week 9 || 3/31/26 || [[Connection betwen KdV and the Schrodinger Equation]] ||
|-
| Week 10 || 4/7/26 || [[Example Calculations for the KdV and IST]] ||
|-
| Week 11 || 4/14/26 || [[Reaction-Diffusion Systems]] ||
|-
| Week 12 || 4/21/26 || [[Burgers Equation]] || Assignment 2
|-
| Week 13 || 4/28/26 || No topic ||
|-
| Week 14 || 5/5/26 || No topic || Assignment 3
|-
| Week 15 || 5/12/26 || No topic ||
|-
| Week 16 || 5/19/26 || No topic || Final Exam
|}

Category:Nonlinear PDE's Course

2026-03-02T06:01:37Z

Meylan:

These pages contain notes for the course in Nonlinear PDE's taught by Mike Meylan

# [[Method of Characteristics for Linear Equations]]
# [[Traffic Waves]]
# [[Nonlinear Shallow Water Waves]]
# [[Introduction to KdV]]
# [[Numerical Solution of the KdV]]
# [[Conservation Laws for the KdV]]
# [[Introduction to the Inverse Scattering Transform]]
# [[Properties of the Linear Schrodinger Equation]]
# [[Connection betwen KdV and the Schrodinger Equation]]
# [[Example Calculations for the KdV and IST]]
# [[Reaction-Diffusion Systems]]
# [[Burgers Equation]]

{| class="wikitable"
|+ 2026 course
! Week !! Date !! Topic !! Assessment
|-
| Week 1 || 2/3/26 || [[Method of Characteristics for Linear Equations]] ||
|-
| Week 2 || 2/10/26 || [[Traffic Waves]] ||
|-
| Week 3 || 2/17/26 || [[Nonlinear Shallow Water Waves]] ||
|-
| Week 4 || 2/24/26 || [[Introduction to KdV]] ||
|-
| Week 5 || 3/3/26 || [[Numerical Solution of the KdV]] || Assignment 1
|-
| Week 6 || 3/10/26 || [[Conservation Laws for the KdV]] ||
|-
| Week 7 || 3/17/26 || [[Introduction to the Inverse Scattering Transform]] ||
|-
| Week 8 || 3/24/26 || [[Properties of the Linear Schrodinger Equation]] ||
|-
| Week 9 || 3/31/26 || [[Connection betwen KdV and the Schrodinger Equation]] ||
|-
| Week 10 || 4/7/26 || [[Example Calculations for the KdV and IST]] ||
|-
| Week 11 || 4/14/26 || [[Reaction-Diffusion Systems]] ||
|-
| Week 12 || 4/21/26 || [[Burgers Equation]] || Assignment 2
|-
| Week 13 || 4/28/26 || ||
|-
| Week 14 || 5/5/26 || || Assignment 3
|-
| Week 15 || 5/12/26 || ||
|-
| Week 16 || 5/19/26 || || Final Exam
|}

Category:Nonlinear PDE's Course

2026-03-02T05:00:51Z

Meylan:

Template:Frequency domain equations for a floating plate

2025-11-12T23:01:52Z

Meylan:

If we make the assumption of [[Frequency Domain Problem]] that everything is proportional to
<math>\exp (-\mathrm{i}\omega t)\,</math> the equations become
<center><math>
\begin{align}
-\mathrm{i}\omega\zeta &= \partial_z\phi , &z=0 \\
-\rho g\zeta - \mathrm{i}\omega\rho \phi &= D \partial_x^4 \zeta -\omega^2 \rho_i h \zeta, &z=0 \\
\Delta \phi &= 0, &-h<z<0 \\
\partial_z \phi &= 0, &z=-h,
\end{align}
</math></center>
where
<math>\zeta</math> is the surface displacement and <math>\phi</math> is the velocity potential in the frequency domain.

These equations can be simplified by defining <math>\alpha = \omega^2/g</math>,
<math>\beta = D/\rho g</math> and <math>\gamma = \rho_i h/\rho</math> to obtain
<center><math>
\begin{align}
\Delta \phi &= 0, &-h < z \leq 0 \\
\partial_z \phi &= 0, &z = - h \\
\beta \partial_x^4 \zeta + \left( 1 - \gamma\alpha \right) \zeta &= -\mathrm{i} \sqrt{g \alpha}\phi, &z = 0 \\
-\mathrm{i}\omega\zeta &= \partial_z\phi , &z=0 .
\end{align}
</math></center>

Template:Frequency domain equations for a floating plate

2025-11-12T23:01:33Z

Meylan:

If we make the assumption of [[Frequency Domain Problem]] that everything is proportional to
<math>\exp (-\mathrm{i}\omega t)\,</math> the equations become
<center><math>
\begin{align}
-\mathrm{i}\omega\zeta &= \partial_z\phi , &z=0 \\
-\rho g\zeta - \mathrm{i}\omega\rho \phi &= D \partial_x^4 \zeta -\omega^2 \rho_i h \zeta, &z=0 \\
\Delta \phi &= 0, &-h<z<0 \\
\partial_z \phi &= 0, &z=-h,
\end{align}
</math></center>
where
<math>\zeta</math> is the surface displacement and <math>\phi</math> is the velocity potential in the frequency domain.

These equations can be simplified by defining <math>\alpha = \omega^2/g</math>,
<math>\beta = D/\rho g</math> and <math>\gamma = \rho_i h/\rho</math> to obtain
<center><math>
\begin{align}
\Delta \phi &= 0, &-h < z \leq 0 \\
\partial_z \phi &= 0, &z = - h \\
\beta \partial_x^4 \zeta + \left( 1 - \gamma\alpha \right) \zeta &= -\mathrm{i} \sqrt{\alpha/g}\phi, &z = 0 \\
-\mathrm{i}\omega\zeta &= \partial_z\phi , &z=0 .
\end{align}
</math></center>

Template:Frequency domain equations for a floating plate

2025-11-12T22:55:04Z

Meylan:

If we make the assumption of [[Frequency Domain Problem]] that everything is proportional to
<math>\exp (-\mathrm{i}\omega t)\,</math> the equations become
<center><math>
\begin{align}
-\mathrm{i}\omega\zeta &= \partial_z\phi , &z=0 \\
-\rho g\zeta - \mathrm{i}\omega\rho \phi &= D \partial_x^4 \zeta -\omega^2 \rho_i h \zeta, &z=0 \\
\Delta \phi &= 0, &-h<z<0 \\
\partial_z \phi &= 0, &z=-h,
\end{align}
</math></center>
where
<math>\zeta</math> is the surface displacement and <math>\phi</math> is the velocity potential in the frequency domain.

These equations can be simplified by defining <math>\alpha = \omega^2/g</math>,
<math>\beta = D/\rho g</math> and <math>\gamma = \rho_i h/\rho</math> to obtain
<center><math>
\begin{align}
\Delta \phi &= 0, &-h < z \leq 0 \\
\partial_z \phi &= 0, &z = - h \\
\beta \partial_x^4 \zeta + \left( 1 - \gamma\alpha \right) \zeta &= -\mathrm{i} \sqrt{\alpha}\phi, &z = 0 \\
-\mathrm{i}\omega\zeta &= \partial_z\phi , &z=0 .
\end{align}
</math></center>

Template:Linear elastic plate on water time domain

2025-11-12T22:53:51Z

Meylan:

We begin with the [[Linear and Second-Order Wave Theory| linear equations]] for a fluid.
The kinematic condition is the same
<center><math>
\frac{\partial\zeta}{\partial t} = \frac{\partial\Phi}{\partial z} , \ z=0; </math></center>
but the dynamic condition needs to be modified to include the effect of the the plate
<center><math> -\rho g\zeta + \rho \frac{\partial\Phi}{\partial t}
= D \frac{\partial^4 \eta}{\partial x^4} + \rho_i h \frac{\partial^2 \eta}{\partial t^2}
, \ z=0; </math></center>
We also have [[Laplace's Equation|Laplace's equation]]
<center><math>
\Delta \Phi = 0,\,\,-h<z<0
</math></center>
and the usual non-flow condition at the bottom surface
<center><math>
\partial_z \Phi = 0,\,\,z=-h,
</math></center>
where
<math>\zeta</math> is the surface displacement, <math>\Phi</math> is the velocity potential,
and <math>\rho</math> is the fluid density.

Introduction to the Inverse Scattering Transform

2025-09-10T03:23:03Z

Meylan:

Introduction to KdV

2025-08-14T05:32:41Z

Meylan: /* Travelling Wave Solution */

{{nonlinear waves course
| chapter title = Introduction to KdV
| next chapter = [[Numerical Solution of the KdV]]
| previous chapter = [[Nonlinear Shallow Water Waves]]
}}

{{complete pages}}

The KdV (Korteweg-De Vries) equation is one of the most important non-linear
pde's. It was originally derived to model shallow water waves with weak
nonlinearities, but it has a wide variety of applications. The derivation of
the KdV is given in [[KdV Equation Derivation]]. The KdV equation
is written as
<center><math>
\partial _{t}u+6u\partial _{x}u+\partial _{x}^{3}u=0.
</math></center>

More information about it can be found at [http://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation Korteweg de Vries equation]

==Travelling Wave Solution==

The KdV equation posesses travelling wave solutions. One particular
travelling wave solution is called a soltion and it was discovered
experimentally by [http://en.wikipedia.org/wiki/John_Scott_Russell John Scott Russell]
in 1834. However, it was not understood
theoretically until the work of [http://en.wikipedia.org/wiki/Diederik_Korteweg Korteweg] and
[http://en.wikipedia.org/wiki/Gustav_de_Vries de Vries] in 1895.

We begin with the assumption that the wave travels with contant form, i.e.
is of the form
<center><math>
u\left( x,t\right) =f\left( x-ct\right)
</math></center>
Note that in this equation the parameter <math>c</math> is an unknown as is the
function <math>f.</math>
Only very special values of <math>c</math> and <math>f</math> will give travelling
waves.
We introduce the coordinate <math>\zeta = x - ct</math>.
If we substitute this expression into the KdV equation we obtain
<center><math>
-c\frac{\mathrm{d}f}{\mathrm{d}\zeta}+6f\frac{\mathrm{d}f}{\mathrm{d}\zeta}+\frac{\mathrm{d}^{3}f}{\mathrm{d}\zeta^{3}}=0
</math></center>
We can integrate this with respect to <math>\zeta</math> to obtain
<center><math>
-cf+3f^{2}+\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}=A_{1}
</math></center>
where <math>A</math> is a constant of integration.

If we think about this equation as Newton's second law in a potential well <math>
V(f) </math> then the equation is
<center><math>
\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}= -\frac{\mathrm{d}V}{\mathrm{d}f}
</math></center>
then the potential well is given by
<center><math>
V\left( f\right) =-A_{0}-A_{1}f-c\frac{f^{2}}{2}+f^{3}
</math></center>
Therefore our equation for <math>f</math> may be thought of as the motion of a particle
in a cubic well.

The constant <math>A_0</math> has no effect on our solution so we can set it to be zero.
We can choose the constant <math>A_{1}=0</math> and then we have a
maximum at <math>f=0</math>. There is a solution which rolls from this at <math>t=-\infty </math>
and then runs up the other side and finally returns to the maximum at <math>
t=\infty .</math> This corresponds to a solitary wave solution.

We can also think about the equation as a first order system using <math>
f^{^{\prime }}=v.</math> This gives us
<center><math>\begin{matrix}
\frac{\mathrm{d}f}{\mathrm{d}\zeta} &=&v \\
\frac{\mathrm{d}v}{\mathrm{d}\zeta} &=&A_{1}+cf-3f^{2}
\end{matrix}</math></center>
If we chose <math>A_{1}=0</math> then we obtain two equilibria at <math>(f,v)=\left(
0,0\right) </math> and <math>(c/3,0).</math> If we analysis these equilibria we find the
first is a saddle and the second is a nonlinear center (it is neither repelling nor
attracting).
The Jacobian matrix for the saddle point is
<center><math>
J_{\left( 0,0\right) }=\left(
\begin{array}{cc}
0 & 1 \\
c & 0
\end{array}
\right)
</math></center>
which has eigenvalues at <math>\pm \sqrt{c}</math> and the incident directions are
<center><math>
\left(
\begin{array}{c}
1\\
\sqrt{c}
\end{array}
\right) \leftrightarrow \sqrt{c},\left(
\begin{array}{c}
-1 \\
\sqrt{c}
\end{array}
\right) \leftrightarrow -\sqrt{c}
</math></center>
There is a
homoclinic connection which connects the equilibrium point at the origin. This holoclinic
connection represents the solitary wave. Within this homoclinic connection
lie periodic orbits which represent the cnoidal waves.

{| class="wikitable"
|-
! Phase portrait
! Solitary Wave
! Cnoidal Wave
|-
| [[Image:Kdv phase portrait.jpg|thumb|350px|Phase portrait]]
| [[Image:Kdv wave solitary2.gif|thumb|350px|Solitary Wave]]
| [[Image:Kdv wave cn.gif|thumb|350px|Cnoidal Wave]]
|}

We can also integrate the equation
<center><math>
-cf+3f^{2}+\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}=A_{1}
</math></center>
by multiplying by <math>f^{\prime }</math>and integrating. This gives us
<center><math>
\frac{\left( f^{^{\prime }}\right) ^{2}}{2}=A_{0}+A_{1}f + c\frac{
f^{2}}{2}-f^{3} = -V(f)
</math></center>
It is no coincidence that the right hand side is the potential energy, because this
is nothing more that the equation for conservation of energy (or the first
integral of the Lagrangian system) which does not depend on <math>\zeta</math>.

This is a separable equation and the only challenge is to integrate
<center><math>
\int \frac{1}{\sqrt{-2V(f)}} \mathrm{d}f.
</math></center>

==Formula for the solitary wave==

We know that the solitary wave solution is found when <math>
A_{0}=A_{1}=0.</math> This gives us
<center><math>
\left( f^{^{\prime }}\right) ^{2}=f^{2}\left( c-2f\right)
</math></center>
This can be solved by separation of variables to give
<center><math>
\int \frac{\mathrm{d}f}{f\sqrt{c-2f}}=\int \mathrm{d}\zeta
</math></center>
We then substitute
<center><math>
f=\frac{1}{2}c\,\mathrm{sech}^{2}\left( s\right)
</math></center>
and note that
<center><math>
c-2f=c\left( 1-\mathrm{sech}^{2}\left( s\right) \right) =c\tanh ^{2}\left( f\right)
</math></center>
and that
<center><math>
\frac{\mathrm{d}f}{\mathrm{d}s}=-c\frac{\sinh \left( f\right) }{\cosh ^{3}\left( f\right) }
</math></center>
This means that
<center><math>\begin{matrix}
\int \frac{\mathrm{d}f}{f\sqrt{c-2f}} &=&-\frac{2}{\sqrt{c}}\int \frac{\sinh \left(
f\right) }{\mathrm{sech}^{2}\left( s\right) \tanh \left( f\right) \cosh ^{3}\left(
f\right) }\mathrm{d}s \\
&=&-\frac{2}{\sqrt{c}}s
\end{matrix}</math></center>
This gives us
<center><math>
-\frac{2}{\sqrt{c}}s=\zeta+a
</math></center>
Therefore
<center><math>
f=\frac{1}{2}c\,\mathrm{sech}^{2}\left( \frac{\sqrt{c}}{2}\left( \zeta+a\right) \right)
</math></center>
Of course we assumed that <math>x=x-ct</math> so the formula for the solitary wave is
given by
<center><math>
f\left( x-ct\right) =\frac{1}{2}c\,\mathrm{sech}^{2}\left[ \frac{\sqrt{c}}{2}\left(
x-ct+a\right) \right]
</math></center>
Note that a solution exists for each
<math>c</math>, and that the amplitude is proportional to <math>c.</math> All of this was
discovered experimentally by
[http://en.wikipedia.org/wiki/John_Scott_Russell John Scott Russell].

==Formula for the cnoidal wave==

If we consider the case when the solution oscillates between two values <math>F_2 < F_3</math>
(which we can assume are also roots of <math>V(f)</math> without loss of generality) then
we can integrate the equation to obtain
<center><math>
f\left( \zeta\right) =F_2+(F_3 - F_2) \mathrm{cn}^{2}\left( \gamma \zeta ;k\right)
</math></center>
where <math>cn</math> is a Jacobi Elliptic function and
<math>\gamma</math> and <math>k</math> are constants which depend on <math>c</math>.
Derivation of this equation is found [[KdV Cnoidal Wave Solutions]].
We can write this equation as
<center><math>
f\left( x-ct\right) =a+b \mathrm{cn}^{2}\left( \gamma (x-ct);k\right)
</math></center>
where <math>a=k^2\gamma^2</math> and <math>c = 6b + 4(2k^2 -1)\gamma^2</math>.
These waves are known as [http://en.wikipedia.org/wiki/Cnoidal_wave cnoidal waves].

In the limit the two solutions agree. We also obtain a sinusoidal solution in the limit of
small amplitude.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|SzZ-KhvvPio}}

=== Part 2 ===

{{#ev:youtube|QltlSQQBtrs}}

=== Part 3 ===

{{#ev:youtube|NJ7h3Z9QtvU}}

=== Part 4 ===

{{#ev:youtube|NictSlSgRbM}}

Introduction to the Inverse Scattering Transform

2025-08-07T00:56:57Z

Meylan:

{{nonlinear waves course
| chapter title = Introduction to the Inverse Scattering Transform
| next chapter = [[Properties of the Linear Schrodinger Equation]]
| previous chapter = [[Conservation Laws for the KdV]]
}}

= Miura transform =

The inverse scattering transformation gives a way to solve the KdV equation
exactly. You can think about is as being an analogous transformation to the
Fourier transformation, except it works for a non linear equation. We want to
be able to solve
<center><math>\begin{matrix}
\partial_{t}u+6u\partial_{x}u+\partial_{x}^{3}u & =0\\
u(x,0) & =f\left( x\right)
\end{matrix}</math></center>
with <math>\left\vert u\right\vert \rightarrow0</math> as <math>x\rightarrow\pm\infty.</math>

The Miura transformation is given by
<center><math>
u=-v^{2}-\partial_x v\,
</math></center>
and if <math>v</math> satisfies the mKdV
<center><math>
\partial_{t}v-6v^{2}\partial_{x}v+\partial_{x}^{3}v=0
</math></center>
then <math>u</math> satisfies the KdV (but not vice versa). We can think about the Miura
transformation as being a nonlinear ODE solving for <math>v</math> given <math>u.</math> This
nonlinear ODE is also known as the Riccati equation and there is a well known
transformation which linearises this equation. It we write
<center><math>
v=\frac{\left( \partial_{x}w\right) }{w}
</math></center>
then we obtain the equation
<center><math>
\partial_{x}^{2}w+uw=0
</math></center>
The KdV is invariant under the transformation <math>x\rightarrow x+6\lambda t,</math>
<math>u\rightarrow u+\lambda.</math> Therefore we consider the associated eigenvalue
problem
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
The eigenfunctions and eigenvalues of this scattering problem play a key role
in the inverse scattering transformation. Note that this is Schrodinger's equation.

= Lax Pair =

A powerful way to study the KdV equation is via the idea of a '''Lax pair'''.

A Lax pair consists of two linear operators <math> L </math> and <math> P </math> such that the KdV equation is equivalent to the so-called '''Lax equation''':

:<math> \frac{dL}{dt} = [P, L] = PL - LP </math>

For the KdV equation, a classical Lax pair is:

:<math> L = -\partial_x^2 + u(x,t) </math>

:<math> P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

These are differential operators acting on a function <math> \psi(x,t) </math>. The idea is that the evolution of <math> u(x,t) </math> is such that the operator <math> L </math> undergoes an isospectral deformation: its spectrum does not change in time.

If <math> L\psi = \lambda \psi </math> and <math> \psi_t = P\psi </math>, then consistency requires:

:<math> \frac{d}{dt}(L\psi) = L_t \psi + L\psi_t = PL\psi </math>

so

:<math> L_t = [P, L] </math>

This is the Lax equation.

=== Derivation of the KdV Equation from the Lax Pair ===

Let us now compute <math> [P, L] = PL - LP </math> explicitly, using:

:<math> L = -\partial_x^2 + u </math>

:<math> P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

We compute the commutator acting on a test function <math> \psi </math>:

'''Step 1:''' Compute <math> PL\psi </math>.

<math>
\begin{align}
PL\psi &= P(-\psi_{xx} + u\psi) \\
&= -4\partial_x^3(-\psi_{xx} + u\psi) + 3u\partial_x(-\psi_{xx} + u\psi) + 3\partial_x u (-\psi_{xx} + u\psi)
\end{align}
</math>

We compute each term:

* <math> -4\partial_x^3(-\psi_{xx}) = -4\partial_x^3(-\psi_{xx}) = 4\psi_{xxxxx} </math>
* <math> -4\partial_x^3(u\psi) = -4[(u\psi)_{xxx}] </math>
* <math> 3u\partial_x(-\psi_{xx}) = -3u\psi_{xxx} </math>
* <math> 3u\partial_x(u\psi) = 3u(u\psi)_x = 3u^2\psi_x + 3uu_x\psi </math>
* <math> 3\partial_x u (-\psi_{xx}) = -3u_x\psi_{xx} </math>
* <math> 3\partial_x u (u\psi) = 3u_x u \psi + 3u u_x \psi = 6u u_x \psi </math>

Add all the terms together:

<math>
\begin{align}
PL\psi &= 4\psi_{xxxxx} -4(u\psi)_{xxx} -3u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi -3u_x\psi_{xx} + 6uu_x\psi
\end{align}
</math>

'''Step 2:''' Compute <math> LP\psi </math>.

<math>
\begin{align}
LP\psi &= (-\partial_x^2 + u)P\psi = -\partial_x^2(P\psi) + u P\psi
\end{align}
</math>

We focus on the leading terms in <math> P\psi </math>:

<math>
P\psi = -4\psi_{xxx} + 3u\psi_x + 3u_x\psi
</math>

Then:

<math>
\begin{align}
LP\psi &= -\partial_x^2(-4\psi_{xxx} + 3u\psi_x + 3u_x\psi) + u(-4\psi_{xxx} + 3u\psi_x + 3u_x\psi) \\
&= 4\psi_{xxxxx} -3\partial_x^2(u\psi_x) -3\partial_x^2(u_x\psi) -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Compute:

* <math> \partial_x^2(u\psi_x) = u_{xx}\psi_x + 2u_x\psi_{xx} + u\psi_{xxx} </math>
* <math> \partial_x^2(u_x\psi) = u_{xxx}\psi + 2u_{xx}\psi_x + u_x\psi_{xx} </math>

So:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -3(u_{xx}\psi_x + 2u_x\psi_{xx} + u\psi_{xxx}) -3(u_{xxx}\psi + 2u_{xx}\psi_x + u_x\psi_{xx}) \\
&\quad -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Simplify:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -3u_{xx}\psi_x -6u_x\psi_{xx} -3u\psi_{xxx} -3u_{xxx}\psi -6u_{xx}\psi_x -3u_x\psi_{xx} \\
&\quad -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Group like terms:

* <math> -3u\psi_{xxx} -4u\psi_{xxx} = -7u\psi_{xxx} </math>
* <math> -6u_x\psi_{xx} -3u_x\psi_{xx} = -9u_x\psi_{xx} </math>
* <math> -3u_{xx}\psi_x -6u_{xx}\psi_x = -9u_{xx}\psi_x </math>
* <math> -3u_{xxx}\psi </math>

So:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -7u\psi_{xxx} -9u_x\psi_{xx} -9u_{xx}\psi_x -3u_{xxx}\psi + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

'''Step 3:''' Compute <math> [P, L]\psi = PL\psi - LP\psi </math>

Subtract term by term (all <math> \psi </math> terms):

* The <math> 4\psi_{xxxxx} </math> terms cancel.
* Remaining terms:

<math>
\begin{align}
[P, L]\psi &= (-4(u\psi)_{xxx} + 3u\psi_{xxx}) + 3u_x\psi_{xx} + 9u_{xx}\psi_x + 3u_{xxx}\psi + 6uu_x\psi
\end{align}
</math>

Use:

<math> (u\psi)_{xxx} = u_{xxx}\psi + 3u_{xx}\psi_x + 3u_x\psi_{xx} + u\psi_{xxx} </math>

So:

<math>
\begin{align}
-4(u\psi)_{xxx} &= -4(u_{xxx}\psi + 3u_{xx}\psi_x + 3u_x\psi_{xx} + u\psi_{xxx}) \\
&= -4u_{xxx}\psi -12u_{xx}\psi_x -12u_x\psi_{xx} -4u\psi_{xxx}
\end{align}
</math>

Now add everything:

<math>
\begin{align}
[P,L]\psi &= (-4u_{xxx}\psi -12u_{xx}\psi_x -12u_x\psi_{xx} -4u\psi_{xxx}) + 3u\psi_{xxx} + 3u_x\psi_{xx} + 9u_{xx}\psi_x + 3u_{xxx}\psi + 6uu_x\psi \\
&= (-u_{xxx}\psi -3u_{xx}\psi_x -9u_x\psi_{xx} -u\psi_{xxx}) + 6uu_x\psi
\end{align}
</math>

So the final result is:

<math>
[P,L]\psi = ( -u_{xxx} + 6uu_x ) \psi
</math>

Thus, we have:

:<math> \frac{dL}{dt} = [P, L] \quad \Rightarrow \quad \frac{du}{dt} = -u_{xxx} + 6uu_x </math>

Which is exactly the KdV equation.

=== Summary ===

The KdV equation:

:<math> u_t + 6uu_x + u_{xxx} = 0 </math>

can be written as the Lax equation:

:<math> \frac{dL}{dt} = [P, L] </math>

where

:<math> L = -\partial_x^2 + u, \qquad P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

This reformulation reveals the integrable structure of the KdV equation and allows powerful solution methods such as the inverse scattering transform.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|P3uMk9OS8p4}}

Introduction to the Inverse Scattering Transform

2025-08-07T00:56:24Z

Meylan:

{{nonlinear waves course
| chapter title = Introduction to the Inverse Scattering Transform
| next chapter = [[Properties of the Linear Schrodinger Equation]]
| previous chapter = [[Conservation Laws for the KdV]]
}}

= Miura transform =

The inverse scattering transformation gives a way to solve the KdV equation
exactly. You can think about is as being an analogous transformation to the
Fourier transformation, except it works for a non linear equation. We want to
be able to solve
<center><math>\begin{matrix}
\partial_{t}u+6u\partial_{x}u+\partial_{x}^{3}u & =0\\
u(x,0) & =f\left( x\right)
\end{matrix}</math></center>
with <math>\left\vert u\right\vert \rightarrow0</math> as <math>x\rightarrow\pm\infty.</math>

The Miura transformation is given by
<center><math>
u=-v^{2}-\partial_x v\,
</math></center>
and if <math>v</math> satisfies the mKdV
<center><math>
\partial_{t}v-6v^{2}\partial_{x}v+\partial_{x}^{3}v=0
</math></center>
then <math>u</math> satisfies the KdV (but not vice versa). We can think about the Miura
transformation as being a nonlinear ODE solving for <math>v</math> given <math>u.</math> This
nonlinear ODE is also known as the Riccati equation and there is a well known
transformation which linearises this equation. It we write
<center><math>
v=\frac{\left( \partial_{x}w\right) }{w}
</math></center>
then we obtain the equation
<center><math>
\partial_{x}^{2}w+uw=0
</math></center>
The KdV is invariant under the transformation <math>x\rightarrow x+6\lambda t,</math>
<math>u\rightarrow u+\lambda.</math> Therefore we consider the associated eigenvalue
problem
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
The eigenfunctions and eigenvalues of this scattering problem play a key role
in the inverse scattering transformation. Note that this is Schrodinger's equation.

== Lax Pair ==

A powerful way to study the KdV equation is via the idea of a '''Lax pair'''.

=== The Lax Pair ===

A Lax pair consists of two linear operators <math> L </math> and <math> P </math> such that the KdV equation is equivalent to the so-called '''Lax equation''':

:<math> \frac{dL}{dt} = [P, L] = PL - LP </math>

For the KdV equation, a classical Lax pair is:

:<math> L = -\partial_x^2 + u(x,t) </math>

:<math> P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

These are differential operators acting on a function <math> \psi(x,t) </math>. The idea is that the evolution of <math> u(x,t) </math> is such that the operator <math> L </math> undergoes an isospectral deformation: its spectrum does not change in time.

If <math> L\psi = \lambda \psi </math> and <math> \psi_t = P\psi </math>, then consistency requires:

:<math> \frac{d}{dt}(L\psi) = L_t \psi + L\psi_t = PL\psi </math>

so

:<math> L_t = [P, L] </math>

This is the Lax equation.

=== Derivation of the KdV Equation from the Lax Pair ===

Let us now compute <math> [P, L] = PL - LP </math> explicitly, using:

:<math> L = -\partial_x^2 + u </math>

:<math> P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

We compute the commutator acting on a test function <math> \psi </math>:

'''Step 1:''' Compute <math> PL\psi </math>.

<math>
\begin{align}
PL\psi &= P(-\psi_{xx} + u\psi) \\
&= -4\partial_x^3(-\psi_{xx} + u\psi) + 3u\partial_x(-\psi_{xx} + u\psi) + 3\partial_x u (-\psi_{xx} + u\psi)
\end{align}
</math>

We compute each term:

* <math> -4\partial_x^3(-\psi_{xx}) = -4\partial_x^3(-\psi_{xx}) = 4\psi_{xxxxx} </math>
* <math> -4\partial_x^3(u\psi) = -4[(u\psi)_{xxx}] </math>
* <math> 3u\partial_x(-\psi_{xx}) = -3u\psi_{xxx} </math>
* <math> 3u\partial_x(u\psi) = 3u(u\psi)_x = 3u^2\psi_x + 3uu_x\psi </math>
* <math> 3\partial_x u (-\psi_{xx}) = -3u_x\psi_{xx} </math>
* <math> 3\partial_x u (u\psi) = 3u_x u \psi + 3u u_x \psi = 6u u_x \psi </math>

Add all the terms together:

<math>
\begin{align}
PL\psi &= 4\psi_{xxxxx} -4(u\psi)_{xxx} -3u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi -3u_x\psi_{xx} + 6uu_x\psi
\end{align}
</math>

'''Step 2:''' Compute <math> LP\psi </math>.

<math>
\begin{align}
LP\psi &= (-\partial_x^2 + u)P\psi = -\partial_x^2(P\psi) + u P\psi
\end{align}
</math>

We focus on the leading terms in <math> P\psi </math>:

<math>
P\psi = -4\psi_{xxx} + 3u\psi_x + 3u_x\psi
</math>

Then:

<math>
\begin{align}
LP\psi &= -\partial_x^2(-4\psi_{xxx} + 3u\psi_x + 3u_x\psi) + u(-4\psi_{xxx} + 3u\psi_x + 3u_x\psi) \\
&= 4\psi_{xxxxx} -3\partial_x^2(u\psi_x) -3\partial_x^2(u_x\psi) -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Compute:

* <math> \partial_x^2(u\psi_x) = u_{xx}\psi_x + 2u_x\psi_{xx} + u\psi_{xxx} </math>
* <math> \partial_x^2(u_x\psi) = u_{xxx}\psi + 2u_{xx}\psi_x + u_x\psi_{xx} </math>

So:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -3(u_{xx}\psi_x + 2u_x\psi_{xx} + u\psi_{xxx}) -3(u_{xxx}\psi + 2u_{xx}\psi_x + u_x\psi_{xx}) \\
&\quad -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Simplify:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -3u_{xx}\psi_x -6u_x\psi_{xx} -3u\psi_{xxx} -3u_{xxx}\psi -6u_{xx}\psi_x -3u_x\psi_{xx} \\
&\quad -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Group like terms:

* <math> -3u\psi_{xxx} -4u\psi_{xxx} = -7u\psi_{xxx} </math>
* <math> -6u_x\psi_{xx} -3u_x\psi_{xx} = -9u_x\psi_{xx} </math>
* <math> -3u_{xx}\psi_x -6u_{xx}\psi_x = -9u_{xx}\psi_x </math>
* <math> -3u_{xxx}\psi </math>

So:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -7u\psi_{xxx} -9u_x\psi_{xx} -9u_{xx}\psi_x -3u_{xxx}\psi + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

'''Step 3:''' Compute <math> [P, L]\psi = PL\psi - LP\psi </math>

Subtract term by term (all <math> \psi </math> terms):

* The <math> 4\psi_{xxxxx} </math> terms cancel.
* Remaining terms:

<math>
\begin{align}
[P, L]\psi &= (-4(u\psi)_{xxx} + 3u\psi_{xxx}) + 3u_x\psi_{xx} + 9u_{xx}\psi_x + 3u_{xxx}\psi + 6uu_x\psi
\end{align}
</math>

Use:

<math> (u\psi)_{xxx} = u_{xxx}\psi + 3u_{xx}\psi_x + 3u_x\psi_{xx} + u\psi_{xxx} </math>

So:

<math>
\begin{align}
-4(u\psi)_{xxx} &= -4(u_{xxx}\psi + 3u_{xx}\psi_x + 3u_x\psi_{xx} + u\psi_{xxx}) \\
&= -4u_{xxx}\psi -12u_{xx}\psi_x -12u_x\psi_{xx} -4u\psi_{xxx}
\end{align}
</math>

Now add everything:

<math>
\begin{align}
[P,L]\psi &= (-4u_{xxx}\psi -12u_{xx}\psi_x -12u_x\psi_{xx} -4u\psi_{xxx}) + 3u\psi_{xxx} + 3u_x\psi_{xx} + 9u_{xx}\psi_x + 3u_{xxx}\psi + 6uu_x\psi \\
&= (-u_{xxx}\psi -3u_{xx}\psi_x -9u_x\psi_{xx} -u\psi_{xxx}) + 6uu_x\psi
\end{align}
</math>

So the final result is:

<math>
[P,L]\psi = ( -u_{xxx} + 6uu_x ) \psi
</math>

Thus, we have:

:<math> \frac{dL}{dt} = [P, L] \quad \Rightarrow \quad \frac{du}{dt} = -u_{xxx} + 6uu_x </math>

Which is exactly the KdV equation.

=== Summary ===

The KdV equation:

:<math> u_t + 6uu_x + u_{xxx} = 0 </math>

can be written as the Lax equation:

:<math> \frac{dL}{dt} = [P, L] </math>

where

:<math> L = -\partial_x^2 + u, \qquad P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

This reformulation reveals the integrable structure of the KdV equation and allows powerful solution methods such as the inverse scattering transform.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|P3uMk9OS8p4}}

Introduction to the Inverse Scattering Transform

2025-08-07T00:55:54Z

Meylan:

{{nonlinear waves course
| chapter title = Introduction to the Inverse Scattering Transform
| next chapter = [[Properties of the Linear Schrodinger Equation]]
| previous chapter = [[Conservation Laws for the KdV]]
}}

The inverse scattering transformation gives a way to solve the KdV equation
exactly. You can think about is as being an analogous transformation to the
Fourier transformation, except it works for a non linear equation. We want to
be able to solve
<center><math>\begin{matrix}
\partial_{t}u+6u\partial_{x}u+\partial_{x}^{3}u & =0\\
u(x,0) & =f\left( x\right)
\end{matrix}</math></center>
with <math>\left\vert u\right\vert \rightarrow0</math> as <math>x\rightarrow\pm\infty.</math>

The Miura transformation is given by
<center><math>
u=-v^{2}-\partial_x v\,
</math></center>
and if <math>v</math> satisfies the mKdV
<center><math>
\partial_{t}v-6v^{2}\partial_{x}v+\partial_{x}^{3}v=0
</math></center>
then <math>u</math> satisfies the KdV (but not vice versa). We can think about the Miura
transformation as being a nonlinear ODE solving for <math>v</math> given <math>u.</math> This
nonlinear ODE is also known as the Riccati equation and there is a well known
transformation which linearises this equation. It we write
<center><math>
v=\frac{\left( \partial_{x}w\right) }{w}
</math></center>
then we obtain the equation
<center><math>
\partial_{x}^{2}w+uw=0
</math></center>
The KdV is invariant under the transformation <math>x\rightarrow x+6\lambda t,</math>
<math>u\rightarrow u+\lambda.</math> Therefore we consider the associated eigenvalue
problem
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
The eigenfunctions and eigenvalues of this scattering problem play a key role
in the inverse scattering transformation. Note that this is Schrodinger's equation.

== Lax Pair ==

A powerful way to study the KdV equation is via the idea of a '''Lax pair'''.

=== The Lax Pair ===

A Lax pair consists of two linear operators <math> L </math> and <math> P </math> such that the KdV equation is equivalent to the so-called '''Lax equation''':

:<math> \frac{dL}{dt} = [P, L] = PL - LP </math>

For the KdV equation, a classical Lax pair is:

:<math> L = -\partial_x^2 + u(x,t) </math>

:<math> P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

These are differential operators acting on a function <math> \psi(x,t) </math>. The idea is that the evolution of <math> u(x,t) </math> is such that the operator <math> L </math> undergoes an isospectral deformation: its spectrum does not change in time.

If <math> L\psi = \lambda \psi </math> and <math> \psi_t = P\psi </math>, then consistency requires:

:<math> \frac{d}{dt}(L\psi) = L_t \psi + L\psi_t = PL\psi </math>

so

:<math> L_t = [P, L] </math>

This is the Lax equation.

=== Derivation of the KdV Equation from the Lax Pair ===

Let us now compute <math> [P, L] = PL - LP </math> explicitly, using:

:<math> L = -\partial_x^2 + u </math>

:<math> P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

We compute the commutator acting on a test function <math> \psi </math>:

'''Step 1:''' Compute <math> PL\psi </math>.

<math>
\begin{align}
PL\psi &= P(-\psi_{xx} + u\psi) \\
&= -4\partial_x^3(-\psi_{xx} + u\psi) + 3u\partial_x(-\psi_{xx} + u\psi) + 3\partial_x u (-\psi_{xx} + u\psi)
\end{align}
</math>

We compute each term:

* <math> -4\partial_x^3(-\psi_{xx}) = -4\partial_x^3(-\psi_{xx}) = 4\psi_{xxxxx} </math>
* <math> -4\partial_x^3(u\psi) = -4[(u\psi)_{xxx}] </math>
* <math> 3u\partial_x(-\psi_{xx}) = -3u\psi_{xxx} </math>
* <math> 3u\partial_x(u\psi) = 3u(u\psi)_x = 3u^2\psi_x + 3uu_x\psi </math>
* <math> 3\partial_x u (-\psi_{xx}) = -3u_x\psi_{xx} </math>
* <math> 3\partial_x u (u\psi) = 3u_x u \psi + 3u u_x \psi = 6u u_x \psi </math>

Add all the terms together:

<math>
\begin{align}
PL\psi &= 4\psi_{xxxxx} -4(u\psi)_{xxx} -3u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi -3u_x\psi_{xx} + 6uu_x\psi
\end{align}
</math>

'''Step 2:''' Compute <math> LP\psi </math>.

<math>
\begin{align}
LP\psi &= (-\partial_x^2 + u)P\psi = -\partial_x^2(P\psi) + u P\psi
\end{align}
</math>

We focus on the leading terms in <math> P\psi </math>:

<math>
P\psi = -4\psi_{xxx} + 3u\psi_x + 3u_x\psi
</math>

Then:

<math>
\begin{align}
LP\psi &= -\partial_x^2(-4\psi_{xxx} + 3u\psi_x + 3u_x\psi) + u(-4\psi_{xxx} + 3u\psi_x + 3u_x\psi) \\
&= 4\psi_{xxxxx} -3\partial_x^2(u\psi_x) -3\partial_x^2(u_x\psi) -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Compute:

* <math> \partial_x^2(u\psi_x) = u_{xx}\psi_x + 2u_x\psi_{xx} + u\psi_{xxx} </math>
* <math> \partial_x^2(u_x\psi) = u_{xxx}\psi + 2u_{xx}\psi_x + u_x\psi_{xx} </math>

So:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -3(u_{xx}\psi_x + 2u_x\psi_{xx} + u\psi_{xxx}) -3(u_{xxx}\psi + 2u_{xx}\psi_x + u_x\psi_{xx}) \\
&\quad -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Simplify:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -3u_{xx}\psi_x -6u_x\psi_{xx} -3u\psi_{xxx} -3u_{xxx}\psi -6u_{xx}\psi_x -3u_x\psi_{xx} \\
&\quad -4u\psi_{xxx} + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

Group like terms:

* <math> -3u\psi_{xxx} -4u\psi_{xxx} = -7u\psi_{xxx} </math>
* <math> -6u_x\psi_{xx} -3u_x\psi_{xx} = -9u_x\psi_{xx} </math>
* <math> -3u_{xx}\psi_x -6u_{xx}\psi_x = -9u_{xx}\psi_x </math>
* <math> -3u_{xxx}\psi </math>

So:

<math>
\begin{align}
LP\psi &= 4\psi_{xxxxx} -7u\psi_{xxx} -9u_x\psi_{xx} -9u_{xx}\psi_x -3u_{xxx}\psi + 3u^2\psi_x + 3uu_x\psi
\end{align}
</math>

'''Step 3:''' Compute <math> [P, L]\psi = PL\psi - LP\psi </math>

Subtract term by term (all <math> \psi </math> terms):

* The <math> 4\psi_{xxxxx} </math> terms cancel.
* Remaining terms:

<math>
\begin{align}
[P, L]\psi &= (-4(u\psi)_{xxx} + 3u\psi_{xxx}) + 3u_x\psi_{xx} + 9u_{xx}\psi_x + 3u_{xxx}\psi + 6uu_x\psi
\end{align}
</math>

Use:

<math> (u\psi)_{xxx} = u_{xxx}\psi + 3u_{xx}\psi_x + 3u_x\psi_{xx} + u\psi_{xxx} </math>

So:

<math>
\begin{align}
-4(u\psi)_{xxx} &= -4(u_{xxx}\psi + 3u_{xx}\psi_x + 3u_x\psi_{xx} + u\psi_{xxx}) \\
&= -4u_{xxx}\psi -12u_{xx}\psi_x -12u_x\psi_{xx} -4u\psi_{xxx}
\end{align}
</math>

Now add everything:

<math>
\begin{align}
[P,L]\psi &= (-4u_{xxx}\psi -12u_{xx}\psi_x -12u_x\psi_{xx} -4u\psi_{xxx}) + 3u\psi_{xxx} + 3u_x\psi_{xx} + 9u_{xx}\psi_x + 3u_{xxx}\psi + 6uu_x\psi \\
&= (-u_{xxx}\psi -3u_{xx}\psi_x -9u_x\psi_{xx} -u\psi_{xxx}) + 6uu_x\psi
\end{align}
</math>

So the final result is:

<math>
[P,L]\psi = ( -u_{xxx} + 6uu_x ) \psi
</math>

Thus, we have:

:<math> \frac{dL}{dt} = [P, L] \quad \Rightarrow \quad \frac{du}{dt} = -u_{xxx} + 6uu_x </math>

Which is exactly the KdV equation.

=== Summary ===

The KdV equation:

:<math> u_t + 6uu_x + u_{xxx} = 0 </math>

can be written as the Lax equation:

:<math> \frac{dL}{dt} = [P, L] </math>

where

:<math> L = -\partial_x^2 + u, \qquad P = -4\partial_x^3 + 3u\partial_x + 3\partial_x u </math>

This reformulation reveals the integrable structure of the KdV equation and allows powerful solution methods such as the inverse scattering transform.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|P3uMk9OS8p4}}

Main Page

2025-05-18T03:37:38Z

Meylan: Removed protection from "Main Page"

__NOTOC__
{| id="mp-topbanner" style="width:100%; background:#fcfcfc; margin-top:1.2em; border:1px solid #ccc;"
| style="width:56%; color:#000;" |
{| style="width:280px; border:none; background:none;"
| style="width:280px; text-align:center; white-space:nowrap; color:#000;" |
<div style="font-size:162%; border:none; margin:0; padding:.1em; color:#000;">Welcome to '''Wikiwaves''',</div>
<div style="top:+0.2em; font-size:95%;">the water waves wiki.</div>
<div id="articlecount" style="width:100%; text-align:center; font-size:85%;">[[Special:Statistics|{{NUMBEROFARTICLES}}]] articles</div>
|}
| style="width:17%; font-size:95%;" |
* [[:Category:Simple Linear Waves|Simple waves]]
* [[:Category:Linear Water-Wave Theory|Linear water waves]]
* [[:Category:Nonlinear Water-Wave Theory|Nonlinear water waves]]
| style="width:23%; font-size:95%;" |
* [[:Category:Geophysics|Geophysics]]
* [[:Category:Offshore Engineering|Offshore engineering]]
* [[:Category:Ocean Wave Interaction with Ships and Offshore Structures|Wave and wave body interactions]]
|}
{| id="mp-strapline" style="width:100%; background:none; margin:-.8em 0 -.7em 0;"
| style="font-size:95%; padding:10px 0; margin:0; text-align:left; white-space:nowrap; color:#000;" |[[FAQ]]
| style="font-size:95%; padding:10px 0; margin:0; text-align:right; white-space:nowrap; color:#000;" |
[[:Category:Reference|References]] '''·''' [[Standard Notation|Standard notation]] '''·''' [[test | Test page]]
|}

{| id="mp-upper" style="margin:0 0 0 0; background:none;"
| class="MainPageBG" style="width:55%; border:1px solid #cef2e0; background:#f5fffa; vertical-align:top; color:#000;" |
{| id="mp-left" cellpadding="2" cellspacing="5" style="width:100%; vertical-align:top; background:#f5fffa;"
! <h2 id="mp-i-h2" style="margin:0; background:#cef2e0; font-size:120%; font-weight:bold; border:1px solid #a3bfb1; text-align:left; color:#000; padding:0.2em 0.4em;">Introduction</h2>
|-
| style="color:#000;" |
<div>
Water waves are beautiful and their mathematics is fascinating, deep, and full of amazing results. Their understanding is also of great practical importance for the design of marine structures.

This wiki is devoted to the mathematics of water waves and to providing an online reference to explain all aspects of their theory. The wiki is a work in progress and strongly reflects the research interests of [[Michael Meylan]] who is doing most of work! However, you are strongly encouraged to contribute in any way you see fit (see [[Why you should contribute]]).

For a brief course in wave-body interactions begin with [[:Category:Ocean Wave Interaction with Ships and Offshore Structures|Wave and Wave Body Interactions]].
</div>
|}
| style="border:1px solid transparent;" |
| class="MainPageBG" style="width:45%; border:1px solid #cedff2; background:#f5faff; vertical-align:top;" |
{| id="mp-right" cellpadding="2" cellspacing="5" style="width:100%; vertical-align:top; background:#f5faff;"
! <h2 id="mp-fp-h2" style="margin:0; background:#cedff2; font-size:120%; font-weight:bold; border:1px solid #a3b0bf; text-align:left; color:#000; padding:0.2em 0.4em;">Featured pages</h2>
|-
| style="color:#000;" |
<div>
*[[:Category:Eigenfunction Matching Method|Eigenfunction Matching Method]]: The eigenfunction matching method.
*[[:Category:Wave Scattering in the Marginal Ice Zone|Wave Scattering in the Marginal Ice Zone]]: A description of the geophysical problem in water wave scattering.
*[[:Category:Interaction Theory|Interaction Theory]]: Presents the theory of multiple body interactions.
</div>
|-
! <h2 id="mp-s-h2" style="margin:0; background:#cedff2; font-size:120%; font-weight:bold; border:1px solid #a3b0bf; text-align:left; color:#000; padding:0.2em 0.4em;">Wikiwaves announcements</h2>
|-
| style="color:#000;" |
<div>
* At the moment we are trying to link this wiki to computer code.
* We appreciate all and every edit - even fixing a typing mistake.
</div>
|}
|}

===About Us===
[[Image:NZIMA.jpg|thumb|right]]

This website is organised by [[Michael Meylan]], Dave Smith, Ravi Pethiyagoda and [[WikiWaves:Administrators|others]]. It was initially supported by a grant from the
[http://www.nzima.auckland.ac.nz/ New Zealand Institute of Mathematics].
A significant amount of the initial content was derived from the [http://ocw.mit.edu/index.html MIT opencourseware].

Main Page

2025-05-05T01:40:23Z

Meylan: /* About Us */

Weakly Nonlinear Wave Theory for Periodic Waves (Stokes Expansion)

2025-05-04T23:02:16Z

Meylan: /* Introduction */

{{incomplete pages}}

== Introduction ==

The solution for Stokes waves is valid in deep or intermediate water depth. It is assumed that the wave steepness (given below) is much smaller than one.

<center><math> \varepsilon = a k \ll 1 \, </math></center>

<center><math> k h = O (1) \, </math></center>

where <math> k \, </math> is the wavenumber and <math> h \, </math> is the water depth which is assumed constant.

* Nondimensional Variables

<center><math> X = x k, \quad Z = z k, \quad Y = y k, \, </math></center>

<center><math> \bar{t} = t \sigma, \quad \bar{\eta} = \eta / a, \, </math></center>

<center><math> \Phi = \frac{\sigma\phi}{a g}, \quad \bar{C} = \frac{C}{a g}, \quad D = \frac{g k}{\sigma^2}, \, </math></center>

where <math> x, \ z, \ t, \ h, \ \phi \, </math> and <math> C \, </math> are dimensional variables and <math> X, \ Z, \ \bar{t}, \ \bar{\eta}\, \ \Phi \, </math> and <math> \bar{C} \, </math> are
corresponding nondimensional variables.

== Nondimensional Governing Equation & Boundary Conditions ==

{| border="0" align="center"
| <math> \frac{\partial^2\theta}{\partial X^2} + \frac{\partial^2\phi}{\partial Y^2} + \frac{\partial^2\theta}{\partial Z^2} = 0 \, </math> || width="100" | || width="150" | <math> - k h < Z < \varepsilon \bar{\eta} \, </math> || <math> (3.1.1) \, </math>
|-
| <math> \frac{\partial\Phi}{\partial Z} = 0 \, </math> || align="center" | at || width="150" | <math> Z = - k h \, </math> || <math> (3.1.2) \, </math>
|-
| <math> \frac{\partial\bar{\eta}}{\partial\bar{t}} + \varepsilon D \nabla_h \Phi \cdot \nabla_h \bar{\eta} = D \frac{\partial\Phi}{\partial Z} </math> || align="center" | at || <math> Z = \varepsilon \bar{\eta} \, </math> || <math> (3.1.3) \, </math>
|-
| <math> \frac{\partial\Phi}{\partial\bar{t}} + \frac{1}{2} \varepsilon D [ \nabla \Phi ]^2 = \bar{C} \, </math> || align="center" | at || <math> Z = \varepsilon \bar{\eta} \, </math> || <math> (3.1.4) \, </math>
|}

where <math> \nabla \, </math> and <math> \nabla_h \, </math> stand for gradient and horizontal gradient, respectively.

== Perturbation (Stokes Expansion) ==

Assuming the wave train is weakly nonlinear <math> ( \varepsilon = a k \ll 1 ) \, </math>, its potential and elevation can be perturbed in the order of <math> \varepsilon \, </math>.

<center><math> \Phi = \Phi^{(1)} + \varepsilon\Phi^{(2)} + \varepsilon^2\Phi^{(3)} + \ldots + \varepsilon^{j-1}\Phi^{(j)} + \ldots </math></center>

<center><math> \bar{\eta} = {\bar{\eta}}^{(1)} + \varepsilon {\bar{\eta}}^{(2)} + \varepsilon^2 {\bar{\eta}}^{(3)} + \ldots + \varepsilon^{j-1} {\bar{\eta}}^{(j)} + \dots </math></center>

<center><math> \bar{C} = {\bar{C}}^{(1)} + \varepsilon {\bar{C}}^{(2)} + \varepsilon^2 {\bar{C}}^{(3)} + \ldots + \varepsilon^{j-1} {\bar{C}}^{(j)} + \dots </math></center>

== Hierarchy Equations ==

Using the Taylor expansion, the free-surface boundary conditions (Equations(3.1.3) and (3.1.4) are expanded at the still water level <math> (Z=0) \, </math>. Then we substitute perturbation forms of potential and elevation into the Laplace Equation, bottom and free-surface boundary conditions. The equations are sorted and grouped according to the order in wave steepness <math> \varepsilon^{(j)} \, </math>. The governing equations for <math> j-th \, </math> order solutions is given by:

<center><math> \nabla^2\Phi^{(j)}=0 \quad -h k \leq Z \leq 0 \qquad (3.1.5) </math></center>

<center><math> \bar{\eta}^{(j)} + \frac{\partial{\Phi}^{(j)}}{\bar{\partial{t}}} = P^{(j)} \left[ \Phi^{(j-1)}, {\bar{\eta}}^{(j-1)} \right] \ \mbox{at} \quad Z = 0 \qquad (3.1.6) </math></center>

<center><math> \frac{{\bar{\partial{\eta}}}^{(j)}}{\bar{\partial{t}}} - D \frac{{\partial{\Phi}}^{(j)}}{\partial{Z}} = Q^{(j)} \left[ \Phi^{(j-1)}, {\bar{\eta}}^{(j-1)} \right] \quad \mbox{at} \quad Z = 0 \quad (3.1.7) </math></center>

<center><math> \frac{\partial\Phi^{(j)}}{\partial Z} = 0 \qquad \mbox{at} \quad Z = - k h \qquad (3.1.8) </math></center>

where the <math> P^{(j)} \, </math> and <math> Q^{(j)} \, </math> can be derived in terms of the solutions for the potential and elevation of order <math> (j-1) \, </math> or lower. Therefore, the above hierarchy equations must be solved sequentially from lower to higher order until the required accuracy is reached. To derive the third-order solution for a Stokes wave train, it is adequate to truncate the equations at <math> j = 3 \, </math>.

Up to <math> j = 3, \ P^{(j)} \, </math> and <math> Q^{(j)} \, </math> are given below.

<center><math> P^{(1)} = {\bar{C}}^{(1)} \quad \mbox{and} \quad Q^{(1)} = 0 </math></center>

<center><math> P^{(2)} = - {\bar{\eta}}^{(1)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}\partial{\bar{t}}} - \frac{D}{2} \left| \nabla\Phi^{(1)} \right|^2 + {\bar{C}}^{(2)} </math></center>

<center><math> Q^{(2)} = D {\bar{\eta}}^{(1)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}^2} - D \nabla_h \Phi^{(1)} \cdot \nabla_h {\bar{\eta}}^{(1)} </math></center>

<center><math> P^{(3)} = - D \nabla\Phi^{(1)} \cdot \nabla\Phi^{(2)} - {\bar{\eta}}^{(1)} \frac{\partial}{\partial{Z}} \left[ \frac{D}{2} \left| \nabla\Phi^{(1)} \right|^2 + \frac{\partial\Phi^{(2)}}{\partial\bar{t}} \right] </math></center>

<center><math> - {\bar{\eta}}^{(2)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}\partial{\bar{t}}} - \frac{1}{2} \left( {\bar{\eta}}^{(1)} \right)^2 \frac{\partial^3\Phi^{(1)}}{\partial{Z}^2
\partial{\bar{t}}} + {\bar{C}}^{(3)} </math></center>

<center><math> Q^{(3)} = D {\bar{\eta}}^{(2)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}^2} + \frac{D}{2} \left( {\bar{\eta}}^{(1)} \right)^2 \frac{\partial^3\Phi^{(1)}}{\partial{Z}^3} </math></center>

<center><math> +D {\bar{\eta}}^{(1)} \frac{\partial}{\partial{Z}} \left( \frac{\partial\Phi^2}{\partial{Z}} - \nabla_h \Phi^{(1)} \cdot \nabla_h {\bar{\eta}}^{(1)} \right) </math></center>

* Solving the non-dimensional Equations from lower order <math> (j=1) \, </math> to higher order <math> (j=3) \, </math> for the non-dimensional solutions ( '''wave advances in the x-direction''' ).

<center><math> \Phi^{(1)} = \frac{ \cosh ( k h + Z ) }{\cosh k h} \sin ( X - \bar{t} ), </math></center>

<center><math> {\bar{\eta}}^{(1)} = \cos ( X - \bar{t} ), \quad {\bar{C}}^{(1)} = 0 </math></center>

<center><math> \Phi^{(2)} = \frac{3}{8} \frac{\cosh (2kh+2Z)}{\sinh^3kh\cosh kh} \sin(2X-2\bar{t}) </math></center>

<center><math> {\bar{\eta}}^{(2)} = \frac{\alpha}{4} ( 3\alpha^2 - 1 ) \cos ( 2X - 2 \bar{t} ) \, </math></center>

<center><math> {\bar{C}}^{(2)} = \frac{1}{2\sinh 2 kh} \, </math></center>

<center><math> \Phi^{(3)} = \frac{1}{64} (\alpha^2-1) (\alpha^2+3) (9\alpha^2-13) \cdot \frac{\cosh(3kh+3Z)}{\cosh kh} \sin (3X - 3\bar{t}) </math></center>

<center><math> {\bar{\eta}}^{(3)} = \frac{3}{8} (\alpha^4 - 3\alpha^2 + 3) \cos (X-\bar{t}) + \frac{3}{64} (8\alpha^6 + (\alpha^2-1)^2) \cos (3X-3\bar{t}) </math></center>

<center><math> {\bar{C}}^{(3)} = 0 \, </math></center>

<center><math> D^{-1} = \alpha^{-1} \left\{ 1 + \varepsilon^2 \left[ \frac{9}{8} \left( \alpha^2 -1 \right)^2 + \alpha^2 \right] \right\} </math></center>

where <math> \alpha = \coth k h \, </math>

* The non-dimensional solutions are then transferred back to the dimensional form.

First-order:

<center><math> \phi^{(1)} = A \frac{\cosh [k(z+h)]}{\cosh (kh)} \sin\theta </math></center>

<center><math> \eta^{(1)} = a \cos\theta \, </math></center>

where

<center><math> \theta = kx - \sigma t + \beta \, </math> and <math> a = A \frac{\sigma}{g} \, </math></center>

Second-order:

<center><math> \phi^{(2)} = \frac{3akA}{8\alpha} (\alpha^2-1)^2 \cosh [ 2k (z+h)] \sin 2 \theta </math></center>

<center><math> \eta^2 = \frac{1}{4} \alpha (3\alpha^2-1) a^2 k \cos 2 \theta </math></center>

Bernoulli Constant: <math> C_o = \frac{1}{4\alpha} a^2 kg (\alpha^2-1) </math>

Third-order:

<center><math> \phi^{(3)} = \frac{1}{64} (\alpha^2-1)(\alpha^2+3)(9\alpha^2-13) \cdot \frac{\cosh[3k(z+h)]}{\cosh3kh} a^2 k^2 A \sin 3 \theta </math></center>

<center><math> \eta^{(3)} = - \frac{3}{8} ( \alpha^4 - 3 \alpha^2 + 3 ) a^3 k^2 \cos \theta + \frac{3}{64} ( 8 \alpha^6 + ( \alpha^2-1 )^2 ) a^3 k^2 \cos 3 \theta </math></center>

Nonlinear Dispersion Relation:

<center><math> \sigma^2 = g k \tanh (kh) \left\{ 1 + k^2 a^2 \left[ \frac{9}{8} \left( \alpha^2 - 1 \right)^2 + \alpha^2 \right] \right\} </math></center>

== Convergence ==

For the fast convergence of the perturbed coefficient, <math> \varepsilon \, </math>, must be much smaller than unity, which is consistent with weakly nonlinear assumption. However, when the ratio of depth to wave length is small, the Stokes perturbation may not be valid.

Convergence rate:

<center><math> R_\phi = \frac{\left| \phi^{(2)} \right|_{mag}}{\left| \phi^{(1)} \right|_{mag}}, \, </math></center>

<math> R_\phi \, </math> is the ratio of the potential magnitude of second-order to that of first order solution at <math> Z = 0 \, </math>.

<center><math> R_\phi = \frac{ 3 \varepsilon ( \alpha^2 - 1 )^2 }{ 8 a } \, </math></center>

For fast convergence, <math> R_\phi \, </math> should be <math> \ll 1 \, </math>. This is true when <math> k h \sim O (1) \, </math>. When <math> k h \ll 1 \, </math>, we have: <math> \alpha \sim (
kh )^{-1} \, </math>, hence <math> R_\phi \sim O \left( \frac{3}{8} \varepsilon ( kh )^{-3} \right) \, </math>

<math> R_\phi \, </math> may be much greater than unity <math> U_r \, </math> sell number <math> U_r = \frac{a}{h} \frac{1}{(kh)^2} = \frac{\varepsilon}{(kh)^3} \, </math>

For <math> R_\phi \ll 1 </math>, then <math> U_r \ll \frac{8}{3} \, </math>.

A few striking '''features''' of a nonlinear wave train can be described for the above equation:

* The crests are steeper and troughs are flatter; (see '''applet''' (''Nonlinear Wave Surface'')).

* Phase velocity increases with the increase in wave steepness.

* Non-closed trajectories of particles movement. (see '''applet''' (''N-Trajectory'')).

* Nonlinear wave characteristics (up to 2nd order).

Wave advancing in the x-direction

Particle velocity <math> \vec{V} = \vec{i} u + \vec{k} w \, </math>

<center><math> u = \frac{akg}{\sigma} \frac{\cosh[k(z+h)]}{\cosh kh} \cos\theta + \frac{3}{4} \frac{a^2 k^2 g}{\sigma} \frac{\cosh[2k(z+h)]}{\sinh^3 kh \cosh kh} \cos 2 \theta </math></center>

<center><math> w = \frac{akg}{\sigma} \frac{\sinh[k(z+h)]}{\cosh kh} \sin\theta + \frac{3}{4} \frac{a^2 k^2 g}{\sigma} \frac{\sinh[2k(z+h)]}{\sinh^3 kh \cosh kh} \sin 2 \theta </math></center>

Acceleration <math> \vec{a} = \vec{i} a_x + \vec{k} a_z \, </math>

<center><math> a_x = \frac{\partial u}{\partial t} + {\vec{V}}^{(1)} \cdot \nabla u^{(1)} = akg \frac{\cosh[k(z+h)]}{\cosh kh} \sin\theta + a^2 k^2 g \left[ \frac{3}{2} \frac{\cosh [2k(z+h)]}{\sinh^3kh \cosh kh} - \frac{1}{\sinh 2kh} \right] \sin 2 \theta </math></center>

<center><math> a_z = \frac{\partial w}{\partial t} + {\vec{V}}^{(1)} \cdot \nabla w^{(1)} = -akg \frac{\sinh[k(z+h)]}{\cosh kh} \cos\theta + a^2 k^2 g \left[ \frac{3}{2} \frac{\sinh [2k(z+h)]}{\sinh^3kh \cosh kh} - \frac{1}{\sinh 2kh} \right] \cos 2 \theta </math></center>

== Particle Trajectory ==

Denoting the mean position of a particle by <math> (x,z) \, </math>, and its instantaneous displacement from the mean position by <math> (\zeta, \zeta) \, </math>, the Lagrangian velocities of the particle are hence <math> u ( x + \zeta, z + \zeta ) \, </math> and <math> w ( x + \zeta, z + \zeta ) \, </math>, they are related to the Eulurian velocities through a Taylor Expansion:

<center><math> u (x+\zeta, z+\zeta) = u^{(1)}(x,z) + u^{(2)}(x,z) + \frac{\partial u^{(1)}}{\partial x} \zeta^{(1)} + \frac{\partial u^{(1)}}{\partial z} \zeta^{(1)} + O ( \varepsilon^2 ) u^{(1)} </math></center>

<center><math> w (x+\zeta, z+\zeta) = w^{(1)}(x,z) + w^{(2)}(x,z) + \frac{\partial w^{(1)}}{\partial x} \zeta^{(1)} + \frac{\partial w^{(1)}}{\partial z} \zeta^{(1)} + O (\varepsilon^2) w^{(1)} </math></center>

where <math> u^{(1)}(x,z),u^{(2)}(x,z),w^{(1)}(x,z) \, </math> and <math> w^{(2)}(x,z) \, </math> are first- and second-order horizontal and vertical velocities.

<center><math> ( \zeta, \xi) \, </math> are calculated by integrating the related Lagrangian velocities. <math> \zeta (t) = \int_0^t u (x+\zeta,z+\xi,\tau)\mathrm{d}t + \zeta_0; \ \xi(t) = \int_0^t w(x+\zeta,z+\xi,\tau)\mathrm{d}\tau + \xi_0 </math></center>

We intend to compute <math> (\zeta, \xi) \, </math> up to second order in wave steepness

<center><math> \zeta(t) = \zeta^{(1)} (x,z,t) + \zeta^{(2)} (x,z,t) + {\overline{\zeta}}^{(2)} (x,z) + O (\varepsilon^2) \zeta^{(1)} </math></center>

<center><math> \xi(t) = \xi^{(1)} (x,z,t) + \xi^{(2)} (x,z,t) + O ( \varepsilon^2) \xi^{(1)} </math></center>

where superscripts stand for orders and overbar denotes a secular term. At lending-order, the solution is the same as that in LWT,

<center><math> \zeta^{(1)} = \int_0^t u^{(1)} (x+\zeta, z+\xi, \tau) \mathrm{d}\tau + \zeta_0^{(1)} = - a \frac{\cosh [ k(z+h)]}{\sinh kh} \sin\theta </math></center>

<center><math> \xi^{(1)} = \int_0^t w^{(1)} (x+\zeta, z+\xi, \tau) \mathrm{d}\tau + \xi_0 = a \frac{\sinh [ k(z+h)]}{\sinh kh} \cos\theta </math></center>

<center><math> \frac{\partial{u}^{(1)}}{\partial{x}} \zeta^{(1)} + \frac{\partial{u}^{(1)}}{\partial{z}} \xi^{(1)} = \frac{a^2k^2g}{\sigma} \left[ \frac{\cosh^2 [k(z+h)]}{\sinh kh \cosh kh} \sin^2 \theta + \frac{\sinh^2 [k(z+h)]}{\sinh kh \cosh kh} \cos^2\theta \right] </math></center>
<center><math> = \frac{a^2k^2\sigma}{\sigma} \left\{ - \frac{1}{\sinh 2kh} \cos 2 \theta + \frac{ \cosh [ 2k(z+h)]}{\sinh 2kh} \right\} </math></center>

<center><math> \frac{\partial w^{(1)}}{\partial{x}} \zeta^{(1)} + \frac{\partial{w}^{(1)}}{\partial{z}} \xi^{(1)} = 0 \, </math></center>

<center><math> u^{(2)} (x,z) + \frac{\partial{u}^{(1)}}{\partial{x}} \zeta^{(1)} + \frac{\partial{u}^{(1)}}{\partial{z}} \xi^{(1)} = \frac{a^2k^2g}{\sigma} \left[ \frac{3}{4} \frac{\cosh [2k(z+h)]}{\sinh^3 kh \cosh kh} - \frac{1}{\sinh 2 kh} \right] \cos 2 \theta </math></center>
<center><math> + \frac{a^2 k^2 g}{\sigma} \frac{\cosh [2 k(z+h)]}{\sinh 2 kh}, </math></center>

<center><math> w^{(2)} (x,z) + \frac{ \partial{w}^{(1)}}{\partial x} \zeta^{(1)} + \frac{ \partial{w}^{(1)}}{\partial z} \xi^{(1)} = \frac{3}{4} \frac{a^2 k^2 g}{\sigma} \frac{ \sinh [ 2 k (z+h)]}{\sinh^3 kh \cosh kh} \sin 2 \theta </math></center>

The leading-order trajectory of a particle is an ellipse of the center at <math> (x,z) \, </math> and a major-axis <math> a \frac{\cosh [ k(z+h)]}{\sinh kh} </math> and minor-axis <math> a
\frac{\sinh[k(z+h)]}{\sinh kh} </math>.

<center><math> \frac{ \left(\zeta^{(1)} \right)^2}{\left[a\frac{\cosh[k(z+h)]}{\sinh kh} \right]^2} + \frac{\left(\xi^{(1)}\right)^2}{\left[a\frac{\sinh[k(z+h)]}{\sinh kh} \right]^2} = 1
</math></center>

The second-order solutions for the displacement are calculated by integrating the related second-order lagrangian velocities.

<center><math> \zeta^{(2)} = a^2 k \left[ - \frac{3}{8} \frac{\cosh [ 2k (z+h)]}{\sinh^4 kh} + \frac{1}{4\sinh^2 kh} \right] \sin 2\theta </math></center>

<center><math> {\bar{\zeta}}^{(2)} = a^2 k \sigma \frac{ \cosh [ 2k (z+h)]}{2\sinh^2 kh} t </math></center>

<center><math> \xi^{(2)} = \frac{3}{8} a^2 k \frac{ \sinh [ 2k (z+h) ] }{ \sinh^4 kh } \cos 2 \theta </math></center>

The secular term <math> \left( \bar{\zeta^{(2)}} \right) \, </math> in the horizontal displacement indicates the particles will continuously move in the wave direction. Hence, the trajectory of a
particle is no longer an ellipse. Because the horizontal mean position of a particle is not fixed at <math> x \, </math> but change with time, we re-define the horizontal mean position by

<center><math> x' = x + a^2 k \sigma \frac{ \cosh[2k(z+h)]}{2\sinh^2kh} t </math> and <math> \theta' = k x' - \sigma t \, </math>.

Correspondingly, the displacement with respect to the instantaneous mean position <math> (x', z) \, </math> is given by,

<center><math> \zeta^{(1)} = - a \frac{\cosh[k(z+h)]}{\sinh kh} \sin\theta', \quad \xi^{(1)} = a \frac{\sinh [k (z+h)]}{\sinh kh} \cos \theta', </math></center>

<center><math> \zeta^{(2)} = a^2 k \left[ - \frac{3}{8} \frac{\cosh [ 2k (z+h)]}{\sinh^4 kh} + \frac{1}{4\sinh^2 kh} \right] \sin 2 \theta', </math></center>

<center><math> \xi^{(2)} = \frac{3}{8} a^2 k \frac{\sinh[2k(z+h)]}{\sinh^4kh} \cos 2 \theta' </math></center>

The trajectory of a particle based on the solution is plotted in Applet (N-trajectory). The time average Lagrangian velocity of a particle is equal to the derivative of the secular term <math> \left(
\bar{\zeta}^{(2)}\right) \, </math> with respect to time.

<center><math> \bar{u}_l = a^2 k \sigma \frac{ \cosh [2k(z+h)]}{2\sinh^2 kh} </math></center>

The integral of the average Lagrangian velocity with respect to water depth renders the average mass flux induced by a periodic wave train over a unit width

<center> Mass flux <math> = \int_{-h}^0 \rho \bar{u}_l \mathrm{d}z = \frac{1}{2} \rho a^2 \sigma \alpha = \frac{1}{2} \rho a^2 kg/\sigma = E/C_p, </math></center>

which is consistent with the result derived using Eulurian approach.

== Dynamic Pressure ==

Using the Bernoulli equation, dynamic pressure head induced by a periodic wave train can be calculated up to second-order,

<center><math> \frac{p}{\rho g} = - \frac{1}{g} \left[ \frac{\partial\phi^{(1)}}{\partial t} + \frac{ \partial\phi^{(1)} }{\partial t} \right] - \frac{1}{2g} \left[ \left| \nabla \phi^{(1)} \right|^2
\right] + C_0, \, </math></center>

<center><math> \frac{p^{(1)}}{\rho g} = a \frac{\cosh [ k (z+h)]}{\cosh (kh)} \cos \theta, \, </math></center>

<center><math> \frac{p^{(2)}}{\rho g} = \frac{3a^2 k}{4\alpha} \left( \alpha^2 -1 \right)^2 \cosh [2k(z+h)] \cos 2 \theta - \frac{a^2k}{4} \alpha ( \alpha^2 -1 ) \cos 2 \theta, \, </math></center>

<center><math> \frac{ {\bar{p}}^{(2)}}{\rho g} = \frac{\alpha}{4} a^2 k ( \alpha^2 - 1 ) \left[ 1 - \cosh [ 2 k (z+h)] \right]. \, </math></center>

== Radiation Stress ==

'''Radiation stress''': defined as the time average of excess quasi momentum flux due to the presence of a periodic wave train.

Up to second order, a wave train advancing in the x-axis, <math> S = \begin{bmatrix} S_{xx} & S_{xy} \\ S_{yx} & S_{yy} \end{bmatrix} </math>

<center><math> S_{xx} = \overline{\int_{-h}^{\eta} \left( p + \rho u^2 \right) \mathrm{d}z} - \int_{-h}^0 p_0 \mathrm{d}z = \int_{-h}^0 \rho \overline{u^2} \mathrm{d}z + \int_{-h}^0 \overline{ \left( p - p_0 \right) } \mathrm{d}z + \overline{ \int_0^\eta p \mathrm{d}z } </math></center>

Noticing that <math> \overline{ \left( p + \rho w^2 \right) } = - \rho g z = p+0, \, </math>

<center><math> S_{xx} = \int_{-h}^0 \rho \overline{ \left( u^{(1)2} - w^{(1)2} \right) } \mathrm{d}z + \overline{ \int_0^\eta p^{(1)} \mathrm{d}z } = \frac{ \rho g a^2 k h }{ \sinh 2 k h } + \frac{1}{4} \rho g a^2 </math></center>

<center><math> S_{yy} = \overline{ \int_{-h}^\eta \left( p + \rho v^2 \right) \mathrm{d}z } - \int_{-h}^0 p_0 \mathrm{d}z = \int_{-h}^0 \overline{ (p - p_0) } \mathrm{d}z + \overline{ \int_0^\eta p \mathrm{d}z } = \frac{ \rho g a^2 kh}{ 2 \sinh 2 k h } </math></center>

<center><math> S_{xy} = S_{yx} = 0 \, </math></center>

<math> S = E \begin{bmatrix} \frac{2kh}{\sinh 2kh} + \frac{1}{2} & 0 \\ 0 \frac{kh}{\sinh 2kh} \end{bmatrix} </math>, where <math> E = \frac{1}{2} \rho g a^2 \, </math>, the energy density.

{| border="0" align="center"
| width="300" | In deep water || In shallow water
|-
| <math> S= E \begin{bmatrix} 1/2 & 0 \\ 0 & 0 \end{bmatrix} </math> || <math> S= E \begin{bmatrix} 3/2 & 0 \\ 0 & 1/2 \end{bmatrix} </math>
|}

In the case of a wave train having an angle, <math> \theta \, </math> with respect to the x-axis,

<center><math> S = E \begin{bmatrix} \left( \frac{kh}{\sinh 2kh} + \frac{1}{2} \right) \frac{(3+\cos 2\theta)}{2} - \frac{1}{2} & \left( 1 + \frac{2kh}{\sinh 2kh} \right) \frac{\sin 2\theta}{4} \\ \left(1 + \frac{2kh}{\sinh 2kh} \right) \frac{\sin 2 \theta}{4} & \left( \frac{kh}{\sinh 2kh} + \frac{1}{2} \right) \frac{(3-\cos 2 \theta)}{2} - \frac{1}{2} \end{bmatrix} </math></center>

-----

This article is based on the lecture notes of Dr. Jun Zhang and the original article can be found
[http://ceprofs.tamu.edu/jzhang/Ocen675/ch3-1.ppt here]

[[Category:Nonlinear Water-Wave Theory]]

Michael Meylan

2024-11-16T05:25:33Z

Meylan:

Michael Meylan is a Professor at the [http://www.newcastle.edu.au The University of Newcastle]. The wikiwaves site is largely his work. His home page can be found at [https://www.newcastle.edu.au/profile/mike-meylan]

[[Category:People|Meylan, Michael]]

Reaction-Diffusion Systems

2024-10-24T23:15:35Z

Meylan: /* The discrete Fourier transform */

{{nonlinear waves course
| chapter title = Reaction-Diffusion Systems
| next chapter = [[Burgers Equation]]
| previous chapter = [[Example Calculations for the KdV and IST]]
}}

{{complete pages}}

We present here a brief theory of reaction diffusion waves.

== Law of Mass Action ==

The law of mass action states that equation rates are proportional to the concentration
of reacting species and the ratio in which they combined. It is discussed in detail in
[[Billingham and King 2000]]. We will present here a few simple examples.

=== Example 1: Simple Decay ===

Suppose we have of chemical <math>P</math> which decays to <math>A</math>, i.e.
<center>
<math>P \to A</math>
</center>
with rate <math>k[P]</math> where <math>[P]</math> denotes concentration. Then if we
set <math>p=[P]</math> and <math>a = [A] </math> we obtain the equations
<center>
<math> \frac{\mathrm{d}p}{\mathrm{d}t} = -kp\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}a}{\mathrm{d}t} = kp</math>
</center>
which has solution
<center>
<math>
p = p_0 e^{-kt}\,\,\,\textrm{and}\,\,\, a = a_0 + p_0(1-e^{-kt})
</math>
</center>
where <math>a_0</math> and <math>p_0</math> are the values of <math>a</math>
<math>p</math> repectively at <math>t=0</math>.

=== Example 2: Quadratic Autocatalysis ===

This example will be important when we consider reaction diffusion problems.
We consider the reaction
<center>
<math>A + B \to 2B</math>
</center>
with rate proportional to <math>k[A][B]</math>. If we define <math>a = [A]</math>
and <math>b = [B]</math> we obtain the following equations
<center>
<math> \frac{\mathrm{d}a}{\mathrm{d}t} = -kab\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}b}{\mathrm{d}t} = kab</math>
</center>
We can solve these equations by observing that
<center>
<math>
\frac{\mathrm{d}(a+b)}{\mathrm{d}t} = 0
</math>
</center>
so that <math>a + b = a_0 + b_0</math>. We can then eliminate <math>a</math> to obtain
<center>
<math>
\frac{\mathrm{d}b}{\mathrm{d}t} = k(a_0 + b_0 - b)b
</math>
</center>
which is separable with solution
<center>
<math>
b = \frac{b_0(a_0 + b_0)e^{k(a_0 + b_0)t}}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
and
<center>
<math>
a = \frac{a_0(a_0 + b_0)}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
Note that <math>a\to 0</math> and <math>b\to a_0 + b_0</math>
as <math>t\to \infty.</math>

== Diffusion ==

The equation for spatially homogeneous diffusion of a chemical with concentration
<math>c</math> is
<center>
<math>
\partial_t c = D\nabla^2 c
</math>
</center>
which is the [http://en.wikipedia.org/wiki/Heat_equation heat equation]. We will consider this in
only one spatial dimension. Consider it on the boundary <math>-\infty < x < \infty</math>. In this case
we can solve by the [http://en.wikipedia.org/wiki/Fourier_transform Fourier transform] and obtain
<center>
<math>
\partial_t \hat{c} = -D k^2 \hat{c}
</math>
</center>
where <math>\hat{c}</math> is the Fourier transform of <math>c</math>. This has solution
<center>
<math>
\hat{c} = \hat{c}_0 e^{-D k^2 t}
</math>
</center>
We can find the inverse transform using convolution and obtain
<center>
<math>
c(x,t) = \frac{1}{\sqrt{4\pi D t}} \int_{-\infty}^{\infty} c_0(x) e^{(x-s)^2/4Dt}\mathrm{d}s
</math>
</center>

=== Solution of the dispersion equation using FFT ===

We can solve the dispersion equation using the discrete Fourier transform and
its closely related numerical implementation the FFT (Fast Fourier Transform).
We have already met the FFT [[Numerical Solution of the KdV]] but we consider it here in more detail.
We consider the concentration
on the finite domain <math>-L \leq x \leq L</math> and use a
[http://en.wikipedia.org/wiki/Fourier_series Fourier series] expansion
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(t) e^{\mathrm{i} k_n x}
</math>
</center>
where <math>k_n = \pi n /L </math>. If we substitute this into the diffusion equation we obtain
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0)e^{-k_n^2 D t} e^{\mathrm{i} k_n x}
</math>
</center>
Note that this is not the same solution as we obtained on the infinite domain because
of the boundary conditions on the finite domain. The coefficients <math>\hat{c}_n(0)</math>
are found using the initial conditions so that
<center>
<math>
\hat{c}_n(0) = \frac{1}{2L} \int_{-L}^{L} e^{-\mathrm{i} k_n x} c_0(x) \mathrm{d}x
</math>
</center>

The key to the numerical solution of this equation is the use of the FFT. We begin by discretising the
domain into a series of <math>N</math> points <math>x_m = -L + 2Lm/N </math>. We then use this to
approximate the integral above and obtain
<center>
<math>
\hat{c}_n(0) = \frac{1}{N} \sum_{m=0}^{N-1} e^{-\mathrm{i} k_n x_m} c_0(x_m)
</math>
</center>
<center>
<math>
= \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} e^{\mathrm{i} \pi n} c_0(x_m)
</math>
</center>

We also get
<center>
<math>
c(x_m,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{\mathrm{i} k_n x_m}
</math>
</center>
<center>
<math>
= \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{2\mathrm{i} \pi nm/N} e^{-\mathrm{i} \pi n}
</math>
</center>
but we know that
<center>
<math>
\hat{c}_n(0) e^{-\mathrm{i} \pi n} = \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} c_0(x_m)
</math>
</center>

=== The discrete Fourier transform ===
The
[http://en.wikipedia.org/wiki/Discrete_Fourier_transform discrete Fourier transform]
of a sequence of ''2N'' complex numbers ''c''<sub>0</sub>, ..., ''c''<sub>''2N''−1</sub> is transformed into the sequence of ''N'' complex numbers <math>\hat{c}</math><sub>0</sub>, ..., <math>\hat{c}</math><sub>''N''−1</sub> by the DFT according to the formula:
<center>
<math>\hat{c}_m = \sum_{n=0}^{N-1} c_n e^{-2\pi \mathrm{i}mn/N} \quad \quad m = 0, \dots, N-1</math>
</center>

We denote the transform by the symbol <math>\mathcal{F}</math>, as in <math>\mathbf{X} = \mathcal{F} \left \{ \mathbf{x} \right \} </math> or <math>\mathcal{F} \left ( \mathbf{x} \right )</math> or <math>\mathcal{F} \mathbf{x}</math>.

The '''inverse discrete Fourier transform (IDFT)''' is given by
<center>
<math>c_n = \frac{1}{N} \sum_{m=0}^{N-1} \hat{c}_m e^{2\pi \mathrm{i}mn/N} \quad \quad n = 0,\dots,N-1.</math>
</center>

Therefore we can write
<center>
<math>
c(x_m,t) = \mathcal{F} \left\{ e^{-k_n^2 D t} \mathcal{F}^{-1} \left\{ c_0(x_m) \right\}
\right\}
</math>
</center>
Note that the choice of where to put the <math>1/N</math> is arbitrary in the definition of FFT and IFFT and
does not exactly match here. Of course since it appears once the formula above is correct regardless.
The only difficulty is that we need to define carefully the values of
<math>k_n</math>

The real power of this method lies with the [http://en.wikipedia.org/wiki/FFT Fast Fourier Transform]
or '''FFT''' algorithm. A naive implementation of the discrete Fourier transform above (or its inverse)
will involve order <math>N^2</math> operations. Using FFT algorithms, this can be reduced to order <math>N \log(N)</math>. This is an incredible speed up, for example
if N = 1024, FFT algorithms are more efficient by a factor of 147. This is the reason FFT
algorithms are used so extensively.

== Reaction Diffusion Equations ==

We consider an auto catalytic reaction where the chemical species also diffuse. In this
case the equations are
<center>
<math>\partial_t a = D\partial_x^2 a - kab</math>
</center>
<center>
<math>\partial_t b = D\partial_x^2 b + kab</math>
</center>
We can non-dimensionalise these equations scaling the variables as
<center>
<math>
z = x/x^*\,\,\,\tau = t/t^*\,\,\,\alpha = a/a_0\,\,\,\beta = b/a_0
</math>
</center>
So that the equations become
<center>
<math>
\frac{1}{t^*}\partial_\tau \alpha = \frac{a_0}{(x^*)^2}D\partial_z^2 \alpha
- k a_0^2 \alpha\beta
</math>
</center>
<center>
<math>
\frac{1}{t^*}\partial_\tau \beta = \frac{a_0}{(x^*)^2}D\partial_z^2 \beta
+ k a_0^2 \alpha\beta
</math>
</center>
If we choose
<center>
<math>
x^* = \sqrt{\frac{D}{ka_0}}\,\,\,t^*=\frac{1}{ka_0}
</math>
</center>
then we obtain the system
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
-\alpha\beta
</math>
</center>
<center>
<math>
\partial_\tau \beta =\partial_z^2 \beta
+ \alpha\beta
</math>
</center>

=== Solution via split step method ===

[[Image:Reaction diffusion.gif|thumb|right|500px|Solutions for <math>\alpha(z,0) =1
\, \beta(z,0) = \exp(-10z^2)</math>]]

We can solve this equations numerically using a
[http://en.wikipedia.org/wiki/Split-step_method split step method]. We assume
that at time <math>\tau</math> we know <math>\alpha(z,\tau)</math>
and <math>\beta(z,\tau)</math>. We then solve first the following equation
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
</math>
</center>
from <math>\tau</math> to <math>\tau + \Delta\tau</math>
(which we can do exactly using the spectral methods just discussed for
the dispersion equation). We write this solution as
<math>\tilde{\alpha}(z,\tau + \Delta\tau)</math>
Then we solve
<center>
<math>
\partial_\tau \alpha = -\alpha\beta
</math>
</center>
by assuming that <math>\beta</math> is constant and subject to the boundary
condition that <math>\alpha(z,\tau) = \tilde{\alpha}(z,\tau + \Delta\tau)</math>.
This gives
<center>
<math>
\alpha(z,\tau + \Delta\tau) = e^{-\beta(z,\tau) \Delta\tau} \tilde{\alpha}(z,\tau+ \Delta\tau)\,
</math>
</center>
and we do likewise for the equation for <math>\beta</math>. Note that
while both steps are exact the result from the split step method is an
approximation with error which becomes smaller as the step size becomes
smaller.

We can easily implement this split step method in matlab and we obtain
a pair of travelling waves.

== Travelling Waves solution ==

When we solve the equations we found the solution formed travelling waves and
we now consider this phenomena in detail.

We define a new coordinate <math>y = z - v\tau</math> (so we will consider only
waves travelling to the right, although we could analyse waves travelling to
the left in a similar fashion). We seek stationary solutions in
<math>\alpha(y)</math> and <math>\beta(y)</math> which satisfy
<center>
<math>
\frac{\mathrm{d}^2 \alpha}{\mathrm{d}y^2} + v \frac{\mathrm{d} \alpha}{\mathrm{d}y} = \alpha\beta
</math>
</center>
and
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} = -\alpha\beta
</math>
</center>
If we add these equations we obtain
<center>
<math>
\frac{\mathrm{d}^2 (\alpha+\beta)}{\mathrm{d}y^2} + v \frac{\mathrm{d} (\alpha+\beta)}{\mathrm{d}y} = 0
</math>
</center>
so that <math>\alpha + \beta = c_0 + c_1 e^{-vy}</math>. Boundary conditions
are that as <math>y\to\infty </math> <math>\alpha = 1</math> and
<math>\beta = 0</math> and if <math>y\to-\infty</math> then <math>\alpha = 0</math> and
<math>\beta = 1</math>. Therefore <math>\alpha + \beta = 1</math>.
This means that, since <math>\alpha \geq 0</math>, we must have
<math>0\leq \beta \leq 1</math>.
We can then obtain the following equation
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} + \beta(1-\beta)= 0
</math>
</center>
which we can write as the system of first order equations.

[[Image:R_d_phase_portrait.jpg|thumb|right|500px|Phase portrait for out system showing
the equilibrium points and the heteroclinic connection]]

We define the variable <math>\gamma = \frac{\mathrm{d}\beta}{\mathrm{d}y}</math>
and we obtain
<center><math>
\begin{align}
\frac{\mathrm{d}\beta}{\mathrm{d} y} &= \gamma&\\
\frac{\mathrm{d}\gamma}{\mathrm{d} y} &= -v\gamma + \beta(\beta -1)& \\
\end{align}
</math></center>
This dynamical system has equilibrium points at <math>(0,0)</math>
and <math>(1,0)</math>. We can analyse these equilibrium points by
linearization. The Jacobian matrix is
<center>
<math>
J =\begin{pmatrix}
0 & 1 \\
-1 + 2\beta & -v
\end{pmatrix}
</math>
</center>

We can easily see that the Jacobian evaluated at our first equilibrium point is
<center>
<math>
J_{(0,0)} =\begin{pmatrix}
0 & 1 \\
-1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\mu_{\pm} = -1/2 (v \mp \sqrt{v^2-4})</math>. Therefore
this point is a nodal sink (possibly a spiral)

[[Image:R_d_wave.gif|right|Travelling wave solution for v=2 and the position on
the heteroclinic connection]]

Additionally,
<center>
<math>
J_{(1,0)} =\begin{pmatrix}
0 & 1 \\
1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\lambda_{\pm} = -1/2 (v \mp \sqrt{v^2+4})</math>.
This is a a saddle point. The unstable and stable
separatrices leave the equilibrium point at <math>(1,0)</math> in the directions
<math> \begin{pmatrix}\lambda_{\pm} \\ 1\end{pmatrix}</math>. The only path on which <math>\beta</math> is bounded
as <math>y\to-\infty</math> are the unstable separatrices. Also, only the
unstable separatrix which enters the region <math>\beta<1</math> is physically meaningful.

To find a travelling wave we need to find a heteroclinic connection
between the two equilibrium points which also has to satisfy the conditions
that <math>0\leq \beta \leq 1</math>.

We need to show that the heteroclinic connection does not cross the <math>\beta</math> axis.
Consider the region
<center>
<math>
R = \left\{(\beta,\gamma)\,|\, \beta<1,\,-k\beta<\gamma<0\right\}
</math>
</center>
On the line <math>\beta = 1,d\beta/dy<0</math> and hence all flow it into <math>R</math>.
On the line <math>\gamma = 0, d\gamma/dy < 0</math> for <math>0<\beta<1</math>. On the line
<math>\gamma = -k \beta</math> we know that <math>d\beta/dy < 0</math> so that integral paths
enter the region if and only if <math>d\gamma/d\beta < \gamma/\beta</math>. We know that
<center>
<math>
d\gamma/d\beta - \gamma/\beta = -v - \frac{\beta(1-\beta)}{\gamma} -\frac{\gamma}{\beta}
= \frac{1}{k} (k^2 - vk +1 -\beta)
</math>
</center>
when <math>\gamma = -k\beta</math>. Therefore we need to find a value of <math>k</math>
so that <math>k^2 - vk +1 < 0</math>, which is possible provided <math>v\geq 2</math>, for example
<math>k = \dfrac{v}{2}</math>.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|rF4X42jP0v8}}

=== Part 2 ===

{{#ev:youtube|D0NwYlM-uOg}}

=== Part 3 ===

{{#ev:youtube|5IEZJtJaDHk}}

=== Part 4 ===

{{#ev:youtube|t_OjTSwVgdo}}

=== Part 5 ===

{{#ev:youtube|kTHVZaYezLk}}

=== Part 6 ===

{{#ev:youtube|MVuSg5_sfYI}}

=== Part 7 ===

{{#ev:youtube|hxKMOHyy6Bw}}

=== Part 8 ===

{{#ev:youtube|yQ-O2KIqu44}}

[[Category:Simple Nonlinear Waves]]

Connection betwen KdV and the Schrodinger Equation

2023-09-14T05:30:01Z

Meylan: /* Scattering Data */

{{nonlinear waves course
| chapter title = Connection betwen KdV and the Schrodinger Equation
| next chapter = [[Example Calculations for the KdV and IST]]
| previous chapter = [[Properties of the Linear Schrodinger Equation]]
}}

If we substitute the relationship
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
into the KdV after some manipulation we obtain
<center><math>
\partial_{t}\lambda w^{2}+\partial_{x}\left( w\partial_{x}Q-\partial
_{x}wQ\right) =0
</math></center>
where <math>Q=\partial_{t}w+\partial_{x}^{3}w-3\left( \lambda-u\right)
\partial_{x}w.</math> If we integrate this equation then we obtain the result that
<center><math>
\partial_{t}\lambda=0
</math></center>
provided that the eigenfunction <math>w</math> is bounded (which is true for the bound
state eigenfunctions). This shows that the discrete eigenvalues are unchanged
and <math>u\left( x,t\right) </math> evolves according to the KdV. Many other
properties can be found

==Scattering Data==

For the discrete spectrum the eigenfunctions behave like
<center><math>
w_{n}\left( x,t\right) =c_{n}\left( t\right) e^{-k_{n}x}
</math></center>
as <math>x\rightarrow\infty</math> with
<center><math>
\int_{-\infty}^{\infty}\left( w_{n}\left( x,t\right) \right) ^{2}\mathrm{d}x=1
</math></center>
The continuous spectrum looks like
<center><math>
w\left( x,t\right) \approx \mathrm{e}^{-\mathrm{i}kx}+r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}
,\ \ \ x\rightarrow-\infty
</math></center>
<center><math>
w\left( x,t\right) \approx a\left( k,t\right) \mathrm{e}^{-\mathrm{i}kx},\ \ \ x\rightarrow
\infty
</math></center>
where <math>r</math> is the reflection coefficient and <math>a</math> is the transmission
coefficient. This gives us the scattering data at <math>t=0</math>
<center><math>
S\left( \lambda,0\right) =\left( \left\{ k_{n},c_{n}\left( 0\right)
\right\} _{n=1}^{N},r\left( k,0\right) ,a\left( k,0\right) \right)
</math></center>
The scattering data evolves as
<center><math>
k_{n}(t)=k_{n}(0) = k_{n}
</math></center>
<center><math>
c_{n}\left( t\right) =c_{n}\left( 0\right) e^{4k_{n}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
We can recover <math>u</math> from scattering data. We write
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}
x}+\frac{1}{2\pi}\int_{-\infty}^{\infty}r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>
Then solve
<center><math>
K\left( x,y;t\right) +F\left( x+y;t\right) +\int_{x}^{\infty}K\left(
x,z;t\right) F\left( z+y;t\right) \mathrm{d}z=0
</math></center>
This is a linear integral equation called the Gelfand-Levitan-Marchenko equation. We then find <math>u</math> from
<center><math>
u\left( x,t\right) =2\partial_{x}K\left( x,x,t\right)
</math></center>

==Reflectionless Potential==

In general the IST is difficult to solve. However, there is a simplification
we can make when we have a reflectionless potential (which we will see gives
rise to the soliton solutions). The reflectionless potential is the case when
<math>r\left( k,0\right) =0</math> for all values of <math>k</math> for some <math>u.</math> In this case
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}x}
</math></center>
then
<center><math>
K\left( x,y,t\right) +\sum_{n=1}^{N}c_{n}^{2}\left( t\right)
e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}K\left( x,z,t\right)
\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
From the equation we can see that
<center><math>
K\left( x,y,t\right) =-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}
</math></center>
If we substitute this into the equation,
<center><math>
-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}+\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}
-\sum_{m=1}^{N}v_{m}\left( x,t\right) e^{-k_{m}z}\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
which leads to
<center><math>
-\sum_{n=1}^{N} v_{n}\left( x,t\right) e^{-k_{n}y}
+\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }
-\sum_{n=1}^{N}\sum_{m=1}^{N}\frac{ c_{n}^{2}\left(
t\right) }{k_{n}+k_{m}}v_{m}\left( x\right) e^{-k_{m}x}e^{-k_{n}\left(
y+x\right) }=0
</math></center>
and we can eliminate the sum over <math>n</math> and the
<math>e^{-k_{n}y}</math> to obtain
<center><math>
-v_{n}\left( x,t\right) +c_{n}^2\left( t\right) e^{-k_{n}x}-\sum_{m=1}
^{N}\frac{c_{n}^2\left( t\right) }{k_{n}+k_{m}}
v_{m}\left( x,t\right) e^{-\left( k_{m}+k_{n}\right) x}=0
</math></center>
which is an algebraic (finite dimensional system) for the unknows <math>v_{n}.</math>.

We can write this as
<center><math>
\left( \mathbf{I}+\mathbf{C}\right) \vec{v}=\vec{f}
</math></center>
where <math>f_{m}=c_{m}^2\left( t\right) e^{-k_{m}x}</math> and
the elements of <math>\mathbf{C}</math> are given by
<center><math>
c_{mn} = \frac{c_{n}^2\left( t\right)}
{k_{n}+k_{m}}e^{-\left( k_{m}+k_{n}\right) x}
</math></center>
This gives us
<center><math>
K\left( x,y,t\right) =-\sum_{m=1}^{N} \left(
\mathbf{I}+\mathbf{C}\right) ^{-1}\vec{f}e^{-k_{m}y}
</math></center>
We then find <math>u(x,t)</math> from <math>K</math>.

=== Single Soliton Example ===

If <math>n=1</math> (a single soliton
solution) we get
<center><math>\begin{matrix}
K\left( x,x,t\right) & =-\frac{c_{1}^2\left( t\right)
e^{-2k_{1}x}}{1+\frac{c_{1}^2\left( t\right) }
{2k_{1}}e^{-2 k_{1} x}}\\
& =\frac{-1}{e^{2k_{1}x - 8k_{1}^{3}t-\alpha} + 1/2k_1}
\end{matrix}</math></center>
where <math>e^{-\alpha}=1/c_{0}^{2}\left( 0\right) .</math> Therefore
<center><math>\begin{align}
u\left( x,t\right) & =2\partial_{x}K\left( x,x,t\right) \\
& =\frac{4k_{1}e^{2k_{1}x-8k_{1}^{3}t-\alpha}}{\left( e^{2k_{1}
x-8k_{1}^{3}t-\alpha} + 1/2k_1\right) ^{2}}\\
& =\frac{8k_{1}^{2}}{\left( \sqrt{2k_{1}}e^{\theta}+e^{-\theta}/\sqrt
{2k_{1}}\right) ^{2}}\\
& =2k_{1}^{2}\mbox{sech}^{2}\left\{ k_{1}\left( x-x_{0}\right) -4k_{1}
^{3}t\right\}
\end{align}</math></center>
where <math>\theta=k_{1}x-4k^{3}t-\alpha/2</math> and <math>\sqrt{2k}e^{-\alpha/2}=e^{-kx_{0}
}</math>. This is of course the single soliton solution.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|KHGoPCoyP28}}

=== Part 2 ===

{{#ev:youtube|8mVq6MQWO3I}}

=== Part 3 ===

{{#ev:youtube|meTgaaGKsfQ}}

=== Part 4 ===

{{#ev:youtube|AfAQyjzvJUU}}

Category:Nonlinear PDE's Course

2023-08-31T03:38:21Z

Meylan:

Numerical Solution of the KdV

2023-08-17T02:58:20Z

Meylan: /* Fourier Transform */

{{nonlinear waves course
| chapter title = Numerical Solution of the KdV
| next chapter = [[Conservation Laws for the KdV]]
| previous chapter = [[Introduction to KdV]]
}}

== Introduction ==

We present here a method to solve the KdV equation numerically. There are
many different methods to solve the KdV and we use here a spectral method
which has been found to work well. Spectral methods work by using the
Fourier transform (or some varient of it) to calculate the derivative.

==Fourier Transform==

Recall that the Fourier transform is given by
<center><math>
\mathcal{F}\left[ f(x)\right] =\hat{f}\left( k\right) =\int_{-\infty
}^{\infty }f\left( x\right) \mathrm{e}^{-\mathrm{i}kx}\mathrm{d}x
</math></center>
and the inverse Fourier transform is given by
<center><math>
f\left( x\right) =\mathcal{F}^{-1}\left[ \hat{f}\left( k\right) \right] =
\frac{1}{2\pi }\int_{-\infty }^{\infty }\hat{f}\left( k\right) \mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>
(note that there are other ways of writing this transform). The most
important property of the Fourier transform is that
<center><math>
\begin{matrix}
\int_{-\infty }^{\infty }\left( \partial _{x}f\left( x\right) \right)
\mathrm{e}^{-\mathrm{i}kx}\mathrm{d}x &=&-\int_{-\infty }^{\infty }f\left( x\right)
\left( \partial_{x}\mathrm{e}^{-\mathrm{i}kx}\right) \mathrm{d}x \\
&=&ik\int_{-\infty }^{\infty }f\left( x\right) \mathrm{e}^{-\mathrm{i}kx}\mathrm{d}x
\end{matrix}
</math></center>
where we have assumed that the function <math>f\left( x\right) </math> vanishes at <math>\pm
\infty .</math> This means that the Fourier transform converts differentiation to
multiplication by <math>i k</math>.

==Solution for the Linearized KdV==

We begin with a simple example. Suppose we want to solve the linearized KdV
equation
<center><math>
\partial _{t}u+\partial _{x}^{3}u=0,\ \ -\infty < x <\infty
</math></center>
subject to solve initial conditions
<center><math>
u\left( x,0\right) =f\left( x\right)
</math></center>
We can solve this equation by taking the Fourier transform. and obtain
<center><math>
\partial _{t}\hat{u}-ik^{3}\hat{u}=0
</math></center>
So that
<center><math>
\hat{u}\left( k,t\right) =\hat{f}\left( k\right) \mathrm{e}^{\mathrm{i}k^{3}t}
</math></center>
Therefore
<center><math>
u\left( x,t\right) = \frac{1}{2\pi}
\int_{-\infty}^{\infty }\hat{f}\left( k\right)
\mathrm{e}^{\mathrm{i}k^{3}t}\mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>

==Numerical Implementation Using the FFT==

The Fast Fourier Transform (FFT) is a method to calculate the fourier
transform efficiently for discrete sets of points. These points need to be
evenly spaced in the <math>x</math> plane and are given by
<center><math>
x_{n}=x_{0}+n\Delta x,\ \ 0\leq n\leq N-1
</math></center>
(note they can start at any <math>x</math> value). For the FFT to be as efficient as
possible <math>N</math> should be a power of <math>2.</math> Corresponding to the discrete set of
points in the <math>x</math> domain is a discrete set of points in the <math>k</math> plane given
by
<center><math>
k_{n}=\left\{
\begin{matrix}
n\Delta k,\ \ 0\leq n\leq \frac{N}{2} \\
\left( n-N\right) \Delta k,\ \ \frac{N}{2}+1\leq n\leq N-1
\end{matrix}
\right.
</math></center>
where <math>\Delta k=2\pi /(N\Delta x).</math> Note that this numbering seems slighly
odd and is due to aliasing. We are not that interested in the frequency
domain solution but we need to make sure that we select the correct values
of <math>k</math> for our numerical code.

==Numerical Code for the Linear KdV==

Here is the code to solve the linear KdV using MATLAB

N = 1024;

t=0.1;

x = linspace(-10,10,N);

delta_x = x(2) - x(1);

delta_k = 2*pi/(N*delta_x);

k = [0:delta_k:N/2*delta_k,-(N/2-1)*delta_k:delta_k:-delta_k];

f = exp(-x.^2);

f_hat = fft(f);

u = real(ifft(f_hat.*exp(i*k.^3*t)));

==Numerical Solution of the KdV==

It turns out that a method to solve the KdV equation can be derived using
spectral methods. We begin with the KdV equation written as
<center><math>
\partial _{t}u+3\partial _{x}\left( u\right) ^{2}+\partial _{x}^{3}u=0.
</math></center>
The Fourier transform of the KdV is therefore
<center><math>
\partial _{t}\hat{u}+3ik\widehat{\left( u^{2}\right)} -ik^{3}\hat{u}=0
</math></center>
We solve this equation by a split step method. We write the equation as
<center><math>
\partial _{t}\hat{u}=-3ik\widehat{\left( u^{2}\right)} +ik^{3}\hat{u}
</math></center>
We can solve
<center><math>
\partial _{t}\hat{u}=ik^{3}\hat{u}
</math></center>
exactly while the equation
<center><math>
\partial _{t}\hat{u}=-3ik\widehat{\left( u^{2}\right)}
</math></center>
needs to be solved by time stepping. The idea of the split step method is to
solve alternatively each of these equations when stepping from <math>t</math> to <math>
t+\Delta t.</math> Therefore we solve
<center><math>\begin{matrix}
\hat{u}_{1}\left( k,t+\Delta t\right) &=&\hat{u}\left( k,t\right)
\mathrm{e}^{\mathrm{i}k^{3}\Delta t} \\
\hat{u}\left( k,t+\Delta t\right) &=&\hat{u}_{1}\left( k,t+\Delta t\right)
-3ik \Delta t\widehat{\left( u_{1}^{2}\right)}
\end{matrix}</math></center>
Note that we are using Euler's method to time step and that the solution
could be improved by using a better method, such as the Runge-Kutta 4
method.

The only slighly tricky thing is that we have both <math>u</math> and <math>\hat{u}</math> in the
equation, but we can simply use the Fourier transform to connect these. The
equation then becomes
<center><math>\begin{matrix}
\hat{u}_{1}\left( k,t+\Delta t\right) &=&\hat{u}\left( k,t\right)
\mathrm{e}^{\mathrm{i}k^{3}\Delta t} \\
\hat{u}\left( k,t+\Delta t\right) &=&\hat{u}_{1}\left( k,t+\Delta t\right)
- 3ik\Delta t\left( \mathcal{F}\left( \left( \mathcal{F}^{-1}\left[ \hat{u}_{1}\left( k,t+\Delta t\right)\right]
\right) ^{2}\right) \right)
\end{matrix}</math></center>

==Numerical Code to solve the KdV by the split step method==

Here is some code to solve the KdV using MATLAB

N = 256;

x = linspace(-10,10,N);

delta_x = x(2) - x(1);

delta_k = 2*pi/(N*delta_x);

k = [0:delta_k:N/2*delta_k,-(N/2-1)*delta_k:delta_k:-delta_k];

c=16;

u = 1/2*c*(sech(sqrt(c)/2*(x+8))).^2;

delta_t = 0.4/N^2;

tmax = 0.1; nmax = round(tmax/delta_t);

U = fft(u);

for n = 1:nmax

% first we solve the linear part

U = U.*exp(1i*k.^3*delta_t);

%then we solve the non linear part

U = U - delta_t*(3i*k.*fft(real(ifft(U)).^2));

end

== Example Calculations ==
We cosider the evolution of the KdV with two solitons as
initial condition as given below.
<center><math>
u(x,0) = 8\,\mathrm{sech}^{2}\left(2(x+8)\right)
+ \, 2 \mathrm{sech}^{2}\left(\left(
x +1\right) \right)</math></center>

{| class="wikitable"
|-
! Animation
! Three-dimensional plot.
|-
| [[Image:Two_soliton.gif|thumb|right|500px|Evolution of <math>u(x,t)</math> for two solitons.]]
| [[Image:Two_soliton.jpg|thumb|right|500px|Evolution of <math>u(x,t)</math> for two solitons.
Note the phase shift which occurs in the soliton interaction.]]

|}

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|ZjeuLNdxcRc}}

=== Part 2 ===

{{#ev:youtube|C7gI5PCKvFQ}}

=== Part 3 ===

{{#ev:youtube|X7-p8nBmKw4}}

=== Part 4 ===

{{#ev:youtube|xevnpuOoks0}}

Introduction to KdV

2023-08-10T03:57:26Z

Meylan: /* Travelling Wave Solution */

{{nonlinear waves course
| chapter title = Introduction to KdV
| next chapter = [[Numerical Solution of the KdV]]
| previous chapter = [[Nonlinear Shallow Water Waves]]
}}

{{complete pages}}

The KdV (Korteweg-De Vries) equation is one of the most important non-linear
pde's. It was originally derived to model shallow water waves with weak
nonlinearities, but it has a wide variety of applications. The derivation of
the KdV is given in [[KdV Equation Derivation]]. The KdV equation
is written as
<center><math>
\partial _{t}u+6u\partial _{x}u+\partial _{x}^{3}u=0.
</math></center>

More information about it can be found at [http://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation Korteweg de Vries equation]

==Travelling Wave Solution==

The KdV equation posesses travelling wave solutions. One particular
travelling wave solution is called a soltion and it was discovered
experimentally by [http://en.wikipedia.org/wiki/John_Scott_Russell John Scott Russell]
in 1834. However, it was not understood
theoretically until the work of [http://en.wikipedia.org/wiki/Diederik_Korteweg Korteweg] and
[http://en.wikipedia.org/wiki/Gustav_de_Vries de Vries] in 1895.

We begin with the assumption that the wave travels with contant form, i.e.
is of the form
<center><math>
u\left( x,t\right) =f\left( x-ct\right)
</math></center>
Note that in this equation the parameter <math>c</math> is an unknown as is the
function <math>f.</math>
Only very special values of <math>c</math> and <math>f</math> will give travelling
waves.
We introduce the coordinate <math>\zeta = x - ct</math>.
If we substitute this expression into the KdV equation we obtain
<center><math>
-c\frac{\mathrm{d}f}{\mathrm{d}\zeta}+6f\frac{\mathrm{d}f}{\mathrm{d}\zeta}+\frac{\mathrm{d}^{3}f}{\mathrm{d}\zeta^{3}}=0
</math></center>
We can integrate this with respect to <math>\zeta</math> to obtain
<center><math>
-cf+3f^{2}+\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}=A_{1}
</math></center>
where <math>A</math> is a constant of integration.

If we think about this equation as Newton's second law in a potential well <math>
V(f) </math> then the equation is
<center><math>
\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}= -\frac{\mathrm{d}V}{\mathrm{d}f}
</math></center>
then the potential well is given by
<center><math>
V\left( f\right) =-A_{0}-A_{1}f-c\frac{f^{2}}{2}+f^{3}
</math></center>
Therefore our equation for <math>f</math> may be thought of as the motion of a particle
in a cubic well.

The constant <math>A_0</math> has no effect on our solution so we can set it to be zero.
We can choose the constant <math>A_{1}=0</math> and then we have a
maximum at <math>f=0</math>. There is a solution which rolls from this at <math>t=-\infty </math>
and then runs up the other side and finally returns to the maximum at <math>
t=\infty .</math> This corresponds to a solitary wave solution.

We can also think about the equation as a first order system using <math>
f^{^{\prime }}=v.</math> This gives us
<center><math>\begin{matrix}
\frac{\mathrm{d}f}{\mathrm{d}\zeta} &=&v \\
\frac{\mathrm{d}v}{\mathrm{d}\zeta} &=&A_{1}+cf-3f^{2}
\end{matrix}</math></center>
If we chose <math>A_{1}=0</math> then we obtain two equilibria at <math>(f,v)=\left(
0,0\right) </math> and <math>(3/c,0).</math> If we analysis these equilibria we find the
first is a saddle and the second is a nonlinear center (it is neither repelling nor
attracting).
The Jacobian matrix for the saddle point is
<center><math>
J_{\left( 0,0\right) }=\left(
\begin{array}{cc}
0 & 1 \\
c & 0
\end{array}
\right)
</math></center>
which has eigenvalues at <math>\pm \sqrt{c}</math> and the incident directions are
<center><math>
\left(
\begin{array}{c}
1\\
\sqrt{c}
\end{array}
\right) \leftrightarrow \sqrt{c},\left(
\begin{array}{c}
-1 \\
\sqrt{c}
\end{array}
\right) \leftrightarrow -\sqrt{c}
</math></center>
There is a
homoclinic connection which connects the equilibrium point at the origin. This holoclinic
connection represents the solitary wave. Within this homoclinic connection
lie periodic orbits which represent the cnoidal waves.

{| class="wikitable"
|-
! Phase portrait
! Solitary Wave
! Cnoidal Wave
|-
| [[Image:Kdv phase portrait.jpg|thumb|350px|Phase portrait]]
| [[Image:Kdv wave solitary2.gif|thumb|350px|Solitary Wave]]
| [[Image:Kdv wave cn.gif|thumb|350px|Cnoidal Wave]]
|}

We can also integrate the equation
<center><math>
-cf+3f^{2}+\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}=A_{1}
</math></center>
by multiplying by <math>f^{\prime }</math>and integrating. This gives us
<center><math>
\frac{\left( f^{^{\prime }}\right) ^{2}}{2}=A_{0}+A_{1}f + c\frac{
f^{2}}{2}-f^{3} = -V(f)
</math></center>
It is no coincidence that the right hand side is the potential energy, because this
is nothing more that the equation for conservation of energy (or the first
integral of the Lagrangian system) which does not depend on <math>\zeta</math>.

This is a separable equation and the only challenge is to integrate
<center><math>
\int \frac{1}{\sqrt{-2V(f)}} \mathrm{d}f.
</math></center>

==Formula for the solitary wave==

We know that the solitary wave solution is found when <math>
A_{0}=A_{1}=0.</math> This gives us
<center><math>
\left( f^{^{\prime }}\right) ^{2}=f^{2}\left( c-2f\right)
</math></center>
This can be solved by separation of variables to give
<center><math>
\int \frac{\mathrm{d}f}{f\sqrt{c-2f}}=\int \mathrm{d}\zeta
</math></center>
We then substitute
<center><math>
f=\frac{1}{2}c\,\mathrm{sech}^{2}\left( s\right)
</math></center>
and note that
<center><math>
c-2f=c\left( 1-\mathrm{sech}^{2}\left( s\right) \right) =c\tanh ^{2}\left( f\right)
</math></center>
and that
<center><math>
\frac{\mathrm{d}f}{\mathrm{d}s}=-c\frac{\sinh \left( f\right) }{\cosh ^{3}\left( f\right) }
</math></center>
This means that
<center><math>\begin{matrix}
\int \frac{\mathrm{d}f}{f\sqrt{c-2f}} &=&-\frac{2}{\sqrt{c}}\int \frac{\sinh \left(
f\right) }{\mathrm{sech}^{2}\left( s\right) \tanh \left( f\right) \cosh ^{3}\left(
f\right) }\mathrm{d}s \\
&=&-\frac{2}{\sqrt{c}}s
\end{matrix}</math></center>
This gives us
<center><math>
-\frac{2}{\sqrt{c}}s=\zeta+a
</math></center>
Therefore
<center><math>
f=\frac{1}{2}c\,\mathrm{sech}^{2}\left( \frac{\sqrt{c}}{2}\left( \zeta+a\right) \right)
</math></center>
Of course we assumed that <math>x=x-ct</math> so the formula for the solitary wave is
given by
<center><math>
f\left( x-ct\right) =\frac{1}{2}c\,\mathrm{sech}^{2}\left[ \frac{\sqrt{c}}{2}\left(
x-ct+a\right) \right]
</math></center>
Note that a solution exists for each
<math>c</math>, and that the amplitude is proportional to <math>c.</math> All of this was
discovered experimentally by
[http://en.wikipedia.org/wiki/John_Scott_Russell John Scott Russell].

==Formula for the cnoidal wave==

If we consider the case when the solution oscillates between two values <math>F_2 < F_3</math>
(which we can assume are also roots of <math>V(f)</math> without loss of generality) then
we can integrate the equation to obtain
<center><math>
f\left( \zeta\right) =F_2+(F_3 - F_2) \mathrm{cn}^{2}\left( \gamma \zeta ;k\right)
</math></center>
where <math>cn</math> is a Jacobi Elliptic function and
<math>\gamma</math> and <math>k</math> are constants which depend on <math>c</math>.
Derivation of this equation is found [[KdV Cnoidal Wave Solutions]].
We can write this equation as
<center><math>
f\left( x-ct\right) =a+b \mathrm{cn}^{2}\left( \gamma (x-ct);k\right)
</math></center>
where <math>a=k^2\gamma^2</math> and <math>c = 6b + 4(2k^2 -1)\gamma^2</math>.
These waves are known as [http://en.wikipedia.org/wiki/Cnoidal_wave cnoidal waves].

In the limit the two solutions agree. We also obtain a sinusoidal solution in the limit of
small amplitude.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|SzZ-KhvvPio}}

=== Part 2 ===

{{#ev:youtube|QltlSQQBtrs}}

=== Part 3 ===

{{#ev:youtube|NJ7h3Z9QtvU}}

=== Part 4 ===

{{#ev:youtube|NictSlSgRbM}}

Traffic Waves

2023-07-27T05:04:42Z

Meylan: /* Shocks */

{{nonlinear waves course
| chapter title = Traffic Waves
| next chapter = [[Nonlinear Shallow Water Waves]]
| previous chapter = [[Method of Characteristics for Linear Equations]]
}}

{{complete pages}}

We consider here some simple equations which model traffic flow. This problem is discussed in
[[Billingham and King 2000]].

[[Category:Reference]]

== Equations ==

We consider a single lane of road, and we measure distance along the road with
the variable <math>x</math> and <math>t</math> is time.
We define the following variables
<center><math>
\begin{matrix}
&\rho(x,t) &: &\mbox{car density (cars/km)} \\
& v(\rho) &: &\mbox{car velocity (km/hour)} \\
& q(x,t) =\rho v &: &\mbox{car flow rate (cars/hour)} \\
\end{matrix}
</math></center>

If we consider a finite length of road <math>x_1\leq x \leq x_2</math> then the net flow of cars
in and out must be balanced by the change in density. This means that
<center><math>
\frac{\partial}{\partial t} \int_{x_1}^{x_2} \rho(x,t) \mathrm{d}x = -q(x_2,t) + q(x_1,t)
</math></center>
We now consider continuous densities (which is obviously an approximation) and
set <math>x_2 = x_1 + \Delta x</math> and we obtain
<center><math>
\frac{\partial}{\partial t} \rho(x_1,t) = -\frac{q(x_2,t) + q(x_1,t)}{\Delta x}
</math></center>
and if we take the limit as <math>\Delta x \to 0</math> we obtain the differential equation
<center><math>
\frac{\partial \rho}{\partial t} + \frac{\partial q}{\partial x} = 0
</math></center>
Note that this equation has been derived purely from the need to conserve cars (it is a conservation equation) and
is not possible to solve this equation until we have derived a connection between <math>\rho</math>
and <math>q</math>.

== Equation for <math>\rho</math> only ==

At the moment we assume that we have some expression for <math>v(\rho)</math>
If we substitute the expression for <math>q = v\rho</math> into our differential equation we obtain
<center><math>
\frac{\partial \rho}{\partial t} + \frac{\partial }{\partial x} \left(v(\rho)\rho\right) = 0
</math></center>
which gives us
<center><math>
\frac{\partial \rho}{\partial t} + \left(v^{\prime}(\rho)\rho + v(\rho)\right)
\frac{\partial \rho }{\partial x} = 0
</math></center>
or
<center><math>
\frac{\partial \rho}{\partial t} + c(\rho)\frac{\partial \rho }{\partial x} = 0
</math></center>
where <math>c(\rho) = \left(v^{\prime}(\rho)\rho + v(\rho)\right)</math>
is the '''kinematic wave speed'''. Note that this is not the speed of the cars, but
the speed at which disturbances in the density travel.

== A simple relationship between <math>\rho</math> and <math>q</math> ==

The relationship between <math>\rho</math> and <math>q</math> is an equation of state and
there is no ''exact'' equation since it depends on many unknowns. One of the
simplest relationship between <math>\rho</math> and <math>q</math> is derived from
the following assumptions

* When the density <math>\rho = 0</math> the speed is <math>v=v_0</math>
* When the density is <math>\rho = \rho_{\max} </math> the speed is <math>v=0</math>
* The speed is a linear function of <math>\rho</math> between these two values.

This also gives good fit with measured data. We will either consider the general case or use this simple
relationship. Using this we obtain
<center><math>
v(\rho) = v_0\frac{\rho_{\max} - \rho}{\rho_{\max}}
</math></center>
The flux of cars is given by
<center><math>
q = \rho v(\rho) = v_0\frac{\rho(\rho_{\max} - \rho)}{\rho_{\max}}
</math></center>
and the wave speed is
<center><math>
c(\rho) = \left(v^{\prime}(\rho)\rho + v(\rho)\right) = -\frac{v_0}{\rho_{\max}}\rho + v_0\frac{\rho_{\max} - \rho}{\rho_{\max}}
= v_0\frac{\rho_{\max} - 2\rho}{\rho_{\max}}
</math></center>

{| class="wikitable"
|-
! <math>v</math>
! <math>q</math>
! <math>c</math>
|-
| [[Image:Velocity.jpg|thumb|350px|<math>v(\rho)</math> versus <math>\rho</math>]]
| [[Image:Q_flux.jpg|thumb|350px|<math>q(\rho)</math> versus <math>\rho</math>]]
| [[Image:C_speed.jpg|thumb|350px|<math>\rho(\rho)</math> versus <math>\rho</math>]]
|}

== Small Amplitude Disturbances ==

We can linearise the model by assuming that the variation in density is small so
that we can write
<center><math>
\rho = \rho_0 + \tilde{\rho}
</math></center>
where we assume that <math>\tilde{\rho}</math> is small. This allows us to write the equations as
<center><math>
\frac{\partial \tilde{\rho}}{\partial t} + c(\rho_0) \frac{\partial \tilde{\rho}}{\partial x} = 0
</math></center>
where the main difference between this and the full equation is that the wave speed <math>c</math> is
a constant. This is the linearised equation. Note that this linearisation does not give a good model because
traffic density does not vary only a small amount about some mean (as is the case for accoustic waves where the
density of air is roughly constant).

Under these assumptions the solution to the equation is
<center><math>
\tilde{\rho} = f(x - c(\rho_0)t)
</math></center>
where <math>f</math> is determined by the initial condition. This represents
disturbances which travel with speed <math>c(\rho_0)</math> in the positive
<math>x</math> direction.

We now consider the '''characteristic curves''' which are curves along which the density
<math>\rho</math> is a constant. These are give by
<center><math>
x = X(t) = x_0 + c(\rho_0) t.\,
</math></center>
which are just straight lines of constant slope. We will see shortly that the full (nonlinear)
equations also possess characteristics.

== Nonlinear Initial Value Problem ==

We wish to solve
<center><math>
\frac{\partial \rho}{\partial t} + c(\rho) \frac{\partial \rho}{\partial x} = 0
</math></center>
subject to the initial conditions
<center><math>
\rho(x,0) = \rho_0(x) \,
</math></center>
It turns out that the concept of characteristic curves is very important for this problem.

If we want <math>\rho(X(t),t)</math> to be a constant then we require
<center><math>
\frac{\mathrm{d}}{\mathrm{d}t}\rho(X(t),t) = \frac{\mathrm{d} X}{\mathrm{d}t} \frac{\partial \rho}{\partial x} +
\frac{\partial \rho}{\partial t} = 0.
</math></center>
Comparing this to the governing partial differential equation we can see that we require
<center><math>
\frac{\mathrm{d} X}{\mathrm{d}t} = c(\rho)
</math></center>
This means that the characteristics are straight lines (since <math>\rho</math> is constant) with
slope given by <math> c(\rho_0(x_0))</math> so that the equation for the characteristics is
<center><math>
X(t) = x_0 + c(\rho_0(x_0))t \,
</math></center>

This does not allow us to write down a solution to the initial value problem,
all we can do is write
<center><math>
\rho(x_0 + c(\rho_0(x_0))t,t) = \rho_0(x_0)\,
</math></center>
which allows us to calculate the solution stepping forward in time, but not to determine the solution given
a value of <math>(x,t)</math> (because we have no way of knowing what <math>c(\rho_0(x_0))</math> is.

The characteristics are a family of straight lines which will all have different slopes. If two characteristics
meet, our solution method will break down because there will be two values of the density <math>\rho</math>.
This gives rise to a '''shock'''. It turns
out that this the formation of shocks is a product of the equations themselves and not with the solution method.
We will see shortly that special methods are required to treat these shocks.

=== Case when no shocks are formed ===

The characteristic curves will fill the space without meeting provided that the wave speed <math>c(\rho_0)</math>
is a monotonically increasing function of the distance <math>x</math>. If we work with our previous model we
have
<center><math>
v(\rho) = v_0\frac{\rho_{\max} - \rho}{\rho_{\max}}
</math></center>
and
<center><math>
c(\rho) = v_0\frac{\rho_{\max} - 2\rho}{\rho_{\max}}
</math></center>
so that <math>c</math> is a monotonically decreasing function of density <math>\rho</math>. This means
that the wave speed <math>c(\rho_0)</math> will be
a monotonically increasing function of the distance <math>x</math> if an only if the density is a
monotonically decreasing function. In this case the solution can be calculated straightforwardly
by expansion of the initial density.

==== No shock example ====

We consider the case when <math>\rho_{\max} = v_0 = 1</math> and where the initial density is given
by <math>\rho_0 = 1/2(1- \tanh(x))</math>. The figures below show the initial density, the initial speed,
the characteristics and <math>\rho(x,t)</math> for this case.

{| class="wikitable"
|-
! <math>\rho</math>
! <math>c</math>
|-
| [[Image:Traffic example1 rho.jpg|thumb|350px| <math>\rho_0 = 1/2(1- \tanh(x))</math> ]]
| [[Image:Traffic_example1_c.jpg|thumb|350px|<math>c(\rho_0)</math>]]
|}

{| class="wikitable"
|-
! characteristics
! <math>\rho(x,t)</math>
|-
| [[Image:Traffic_example1_characteristics.jpg|thumb|350px|Characterisitics for <math>\rho_0 = 1/2(1- \tanh(x))</math>]]
| [[Image:Traffic_example11.gif|thumb|350px|<math>\rho(x,t) </math><math> for </math><math>\rho_0 = 1/2(1- \tanh(x))</math>]]
|}

==== Riemann problem and the expansion fan ====

We can consider a simple problem in which there is a jump in the initial density
<center><math>
\rho(x,0) =
\left\{
\begin{matrix}
\rho_{L},& x < 0 \\
\rho_{R},& x > 0
\end{matrix}
\right.
</math></center>
where <math>\rho_{L} > \rho_{R}</math> so that we do not form a shock. In this case
the characteristics on each side of <math>x=0</math> have a different slope and the
question is what happens between. It is easiest to think about the following problem
<center><math>
\rho(x,0) =
\left\{
\begin{matrix}
\rho_{L},& x < -\epsilon \\
\frac{\rho_{R}-\rho_{L}}{2\epsilon}x + \frac{\rho_{R}+\rho_{L}}{2} & -\epsilon \leq x \leq \epsilon \\
\rho_{R},& x > \epsilon
\end{matrix}
\right.
</math></center>
We can then see that we have lines of uniformly varying slope for <math>-\epsilon < x <\epsilon</math>
with slope between <math>c(\rho_L)</math> and <math>c(\rho_R)</math>. If we then take the limit
as <math>\epsilon \to 0</math> we obtain an expansion fan emanating from <math>x=0</math>.

If we assume that
<center><math>
c(\rho) = v_0\frac{\rho_{\max} - 2\rho}{\rho_{\max}}
</math></center>
then we know that on the lines of the expansion fan (which all start at <math>x=0</math>) we have
<math>c(\rho) = x/t</math>. We can rearrange this and solve for <math>x</math> and obtain
<math>\rho =\frac{ 1}{2} \rho_{\max} (1-x/v_0 t)</math>

The solution is then given by
<center><math>
\rho(x,t) =
\left\{
\begin{matrix}
\rho_{L},& x < c(\rho_L) t\\
\frac{ \rho_{\max}}{2} (1-x/v_0 t),& c(\rho_L) t \leq x \leq c(\rho_R) t\\
\rho_{R},& x > c(\rho_R) t
\end{matrix}
\right.
</math></center>
This solution is known as an '''expansion fan'''.

We consider the case when <math>\rho_{\max} = v_0 = 1</math> and where the initial density is given
by
<center><math>
\rho(x,0) =
\left\{
\begin{matrix}
0.6,& x < 0 \\
0.3,& x > 0.
\end{matrix}
\right.
</math></center> The figures below show the initial density, the initial speed,
the characteristics and <math>\rho(x,t)</math> for this case.

{| class="wikitable"
|-
! <math>\rho</math>
! <math>c</math>
|-
| [[Image:Expansion_fan_rho.jpg|thumb|350px| <math>\rho_0</math> ]]
| [[Image:Expansion_fan_c.jpg|thumb|350px|<math>c(\rho_0)</math>]]
|}

{| class="wikitable"
|-
! characteristics
! <math>\rho(x,t)</math>
|-
| [[Image:Expansion_fan_characteristics.jpg|thumb|350px|Characterisitics]]
| [[Image:Expansion_fan2.gif|thumb|350px|<math>\rho(x,t)</math>]]
|}

=== Shocks ===

So far we have only considered the case when <math>c(x_0)</math> is monotonically increasing so that
two characteristics never cross. We now consider the case when characteristics can meet.
A movie of this case is shown below.

{{#ev:youtube|Suugn-p5C1M}}

We can easily see that
the first characteristics to meet will be neighbouring characteristics. Consider two characteristics
with
<center><math>
X_1(t) = x_0 + c(x_0)t\,
</math></center>
<center><math>
X_2(t) = x_0 + \delta x + c(x_0+\delta x)t\,
</math></center>
Then these curves will meet at time <math>T</math> where
<center><math>
x_0 + c(x_0)T = x_0 + \delta x + c(x_0+\delta x)T\,
</math></center>
which implies that
<center><math>
T = -\frac{1}{c^{\prime}(x_0)}
</math></center>
Note the following
* If <math>c^{\prime}(x) > 0 </math> then no shock will form.
* The shock first forms at the minimum positive value of
<math> - \frac{1}{c^{\prime}(x)} </math> for <math> -\infty < x <\infty </math>.

==== Shock Fitting ====

If we calculate the solution using our formula
<center><math>
\rho(x_0 + c(\rho_0(x_0))t,t) = \rho(x_0)\,
</math></center>
then we find that the solution becomes multivalued in the case when a shock forms.
We then have to fit a shock. One way to do this is by imposing the condition that equal
areas are removed and added when we chose the position of the shock.
This corresponds to the condition that
the number of cars must be conserved

==== Speed of the shock ====

If we consider the case when there is a shock at <math>s(t)</math> with
<math>\rho = \rho^{-}</math> at <math>x=s^{-}</math>
and
<math>\rho = \rho^{+}</math> at <math>x=s^{+}</math>
(where <math>s^{-}</math>
is just less than s(t) and <math>s^{+}</math>
is just greater than s(t) ). If we substitute
this into the governing integral equation we obtain
<center><math>
\frac{\partial}{\partial t} \left( \int_{x_1}^{s(t)} \rho(x,t)\mathrm{d}x + \int_{s(t)}^{x_2}
\rho(x,t)\mathrm{d}x \right) = q(x_1,t) - q(x_2,t)
</math></center>
and hence
<center><math>
\int_{x_1}^{x_2}
\frac{\partial \rho(x,t)}{\partial t} \mathrm{d}x + \frac{\mathrm{d}s}{\mathrm{d}t}\rho^{-}
- \frac{\mathrm{d}s}{\mathrm{d}t}\rho^{+} = q(x_1,t) - q(x_2,t)
</math></center>
If we now take the limit as <math>x_1\to x_2</math> we obtain
<center><math>
\frac{\mathrm{d}s}{\mathrm{d}t}\rho^{-}
- \frac{\mathrm{d}s}{\mathrm{d}t}\rho^{+} = q(\rho^{-}) - q(\rho^{+})
</math></center>
so that
<center><math>
\frac{\mathrm{d}s}{\mathrm{d}t} = \frac{q(\rho^{-}) - q(\rho^{+})}
{\rho^{-} - \rho^{+}}
</math></center>

==== Shock example ====

We consider the case when <math>\rho_{\max} = v_0 = 1</math> and where the initial density is given
by <math>\rho_0 = 1/2(1 + \tanh(x))</math>. The figures below show the initial density, the initial speed,
the characteristics and <math>\rho(x,t)</math> for this case.

{| class="wikitable"
|-
! <math>\rho</math>
! <math>c</math>
|-
| [[Image:Traffic example2 rho.jpg|thumb|350px| <math>\rho_0 = 1/2(1+ \tanh(x))</math> ]]
| [[Image:Traffic_example2_c.jpg|thumb|350px|<math>c(\rho_0)</math>]]
|}

{| class="wikitable"
|-
! characteristics
! <math>\rho(x,t)</math>
|-
| [[Image:Traffic_example2_characteristics.jpg|thumb|350px|Characterisitics for <math>\rho_0 = 1/2(1+ \tanh(x))</math>]]
| [[Image:Traffic_example2.gif|thumb|350px|<math>\rho(x,t) </math><math> for </math><math>\rho_0 = 1/2(1+ \tanh(x))</math> Dotted solution is without shock fitting.]]
|}

==== Riemann problem ====

We now consider the Riemann problem
<center><math>
\rho(x,0) =
\left\{
\begin{matrix}
\rho_{L},& x < 0 \\
\rho_{R},& x > 0
\end{matrix}
\right.
</math></center>
where <math>\rho_{L} < \rho_{R}</math>. In this case a shock forms immediately and
the characteristics terminate at the shock. The shock moves with constant speed given by
the equation for the motion of the shock (or can be found by the equal areas rule). We obtain
<center><math>
\rho(x,t) =
\left\{
\begin{matrix}
\rho_{L},& x < \frac{1}{2} \left(c(\rho_{L}) + c(\rho_{R}) \right) t \\
\rho_{R},& x > \frac{1}{2} \left(c(\rho_{L}) + c(\rho_{R}) \right) t
\end{matrix}
\right.
</math></center>
We consider the case when <math>\rho_{\max} = v_0 = 1</math> and where the initial density is given
by
<center><math>
\rho(x,0) =
\left\{
\begin{matrix}
0.3,& x < 0 \\
0.6,& x > 0
\end{matrix}
\right.
</math></center>

The figures below show the initial density, the initial speed,
the characteristics, and <math>\rho(x,t)</math> for this case.

{| class="wikitable"
|-
! <math>\rho</math>
! <math>c</math>
|-
| [[Image:Shock_rho.jpg|thumb|350px| <math>\rho_0</math> ]]
| [[Image:Shock_c.jpg|thumb|350px|<math>c(\rho_0)</math>]]
|}

{| class="wikitable"
|-
! characteristics
! <math>\rho(x,t)</math>
|-
| [[Image:Shock_characteristics.jpg|thumb|350px|Characterisitics with shock shown in green.]]
| [[Image:Shock3.gif|thumb|350px|<math>\rho(x,t)</math> with the red dotted line showing the solution without shock fitting]]
|}

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|SWIXL97H0fw}}

=== Part 2 ===

{{#ev:youtube|j652tvrgm4c}}

=== Part 3 ===

{{#ev:youtube|QVtYsUGtb7I}}

=== Part 4 ===

{{#ev:youtube|Ik2Q7H_aIV8}}

=== Part 5 ===

{{#ev:youtube|lNeFNHAqgEA}}

=== Part 6 ===

{{#ev:youtube|Yo0sW3Nh2FI}}

=== Part 7 ===

{{#ev:youtube|gOEMfOQZAAM}}

=== Part 8 ===

{{#ev:youtube|ZwvTZcWLZKM}}

[[Category:Simple Nonlinear Waves]]

Wavemaker Theory

2023-05-05T01:13:05Z

Meylan: /* Introduction */

{{Ocean Wave Interaction with Ships and Offshore Structures
| chapter title = Wavemaker Theory
| next chapter = [[Ship Kelvin Wake]]
| previous chapter = [[Wave Momentum Flux]]
}}

{{complete pages}}

== Introduction ==

[[Image:Wave_maker.png|600px|right|thumb|Wavemaker]]

We will derive the potential in a two-dimensional wavetank due the motion of the wavemaker. The method is based on the [[:Category:Eigenfunction Matching Method|Eigenfunction Matching Method]].
A paddle is undergoing small amplitude horizontal oscillations with displacement
<center><math> \zeta (z,t) = \mathrm{Re} \left \{\frac{1}{-\mathrm{i}\omega} f(z) e^{-i\omega t} \right \} </math></center>
where <math> f(z) \, </math> is assumed known. Since the time <math>t=0 \,</math> is arbitrary we can assume that
<math>f(z)\,</math> is real but this is not necessary.
Because the oscillations are small the [[Linear and Second-Order Wave Theory| linear equations]] apply (which will be given formally below).
This excitation creates plane progressive waves with amplitude <math> A \, </math> down the tank. The principal objective of wavemaker theory is to determine <math> A \, </math>
as a function of <math> \omega, f(z) \, </math> and <math> h \, </math>. Time-dependent wavemaker theories can also be developed.

== Expansion of the solution ==

{{frequency definition}}

{{velocity potential in frequency domain}}

The equations therefore become
{{standard linear wave scattering equations without body condition}}
The boundary condition at the wavemaker is
<center>
<math>
\left. \partial_x\phi \right|_{x=0} = \partial_t \xi = f(z).
</math>
</center>
We must also apply the [[Sommerfeld Radiation Condition]]
as <math>x\rightarrow\infty</math>. This essentially implies
that the only wave at infinity is propagating away.

{{separation of variables for a free surface}}

== Expansion in Eigenfunctions ==

The wavemaker velocity potential <math> \phi \,</math> can be expressed simply in terms of eigenfunctions

<center><math> \phi = \sum_{n=0}^{\infty} a_n \phi_n (z) e^{-k_n x} </math></center>

and we can solve for the coefficients by matching at <math>x=0 \,</math>

<center><math> \left. \phi_x \right|_{x=0} = \sum_{n=0}^{\infty} -k_n a_n \phi_n (z) = f(z)
</math></center>

It follows that

<center><math> a_n = -\frac{1}{k_n A_n} \int_{-h}^0 \phi_n(z) f(z)\mathrm{d}z </math></center>

=== Far Field Wave ===
One of the primary objecives of wavemaker theory is to determine the far-field wave amplitude <math> A \, </math> in terms of <math> f(z) \, </math>.
The far-field wave component representing progagating waves is given by:

<center><math> \lim_{x\to\infty} \phi = a_0 \phi_0(z) e^{-k_0 x} =
a_0 \frac{\cos k_0(z+h)}{\cos k_0 h } e^{-k_0 x} </math></center>

Note that <math> k_0 \, </math> is imaginery. We therefore obtain the complex amplitude of the propagating wave at infinity, namely modulus and phase, in terms of the wave maker displacement <math> f(z) \, </math>.

For what type of <math> f(z) \, </math> are the non-wavelike modes zero? It is easy to verify by virtue of orthogonality that

<center><math> f(z) \ \sim \ \phi_0 (z) </math></center>

Unfortunately this is not a "practical" displacement since <math> \phi_0 (z) \, </math> depends on <math> \omega\, </math>, so one would need to build a flexible paddle.

== Matlab Code ==

A program to calculate the coefficients for the wave maker problems can be found here
[http://www.math.auckland.ac.nz/~meylan/code/eigenfunction_matching/wavemaker.m wavemaker.m]

=== Additional code ===

This program requires
[http://www.math.auckland.ac.nz/~meylan/code/dispersion/dispersion_free_surface.m dispersion_free_surface.m]
to run

-----

This article is based in part on the MIT open course notes and the original article can be found
[http://ocw.mit.edu/NR/rdonlyres/Mechanical-Engineering/2-24Spring-2002/737E217E-0582-45D5-B1F9-B2ECF977C66E/0/lecture6.pdf here]

[[Ocean Wave Interaction with Ships and Offshore Energy Systems]]

[[Category:Eigenfunction Matching Method]]
[[Category:Pages with Matlab Code]]
[[Category:Complete Pages]]

Wavemaker Theory

2023-05-05T01:12:28Z

Meylan: /* Introduction */

{{Ocean Wave Interaction with Ships and Offshore Structures
| chapter title = Wavemaker Theory
| next chapter = [[Ship Kelvin Wake]]
| previous chapter = [[Wave Momentum Flux]]
}}

{{complete pages}}

== Introduction ==

[[Image:Wave_maker.png|600px|right|thumb|Wavemaker]]

We will derive the potential in a two-dimensional wavetank due the motion of the wavemaker. The method is based on the [[:Category:Eigenfunction Matching Method|Eigenfunction Matching Method]].
A paddle is undergoing small amplitude horizontal oscillations with displacement
<center><math> \zeta (t) = \mathrm{Re} \left \{\frac{1}{\mathrm{i}\omega} f(z) e^{-i\omega t} \right \} </math></center>
where <math> f(z) \, </math> is assumed known. Since the time <math>t=0 \,</math> is arbitrary we can assume that
<math>f(z)\,</math> is real but this is not necessary.
Because the oscillations are small the [[Linear and Second-Order Wave Theory| linear equations]] apply (which will be given formally below).
This excitation creates plane progressive waves with amplitude <math> A \, </math> down the tank. The principal objective of wavemaker theory is to determine <math> A \, </math>
as a function of <math> \omega, f(z) \, </math> and <math> h \, </math>. Time-dependent wavemaker theories can also be developed.

== Expansion of the solution ==

{{frequency definition}}

{{velocity potential in frequency domain}}

The equations therefore become
{{standard linear wave scattering equations without body condition}}
The boundary condition at the wavemaker is
<center>
<math>
\left. \partial_x\phi \right|_{x=0} = \partial_t \xi = f(z).
</math>
</center>
We must also apply the [[Sommerfeld Radiation Condition]]
as <math>x\rightarrow\infty</math>. This essentially implies
that the only wave at infinity is propagating away.

{{separation of variables for a free surface}}

== Expansion in Eigenfunctions ==

The wavemaker velocity potential <math> \phi \,</math> can be expressed simply in terms of eigenfunctions

<center><math> \phi = \sum_{n=0}^{\infty} a_n \phi_n (z) e^{-k_n x} </math></center>

and we can solve for the coefficients by matching at <math>x=0 \,</math>

<center><math> \left. \phi_x \right|_{x=0} = \sum_{n=0}^{\infty} -k_n a_n \phi_n (z) = f(z)
</math></center>

It follows that

<center><math> a_n = -\frac{1}{k_n A_n} \int_{-h}^0 \phi_n(z) f(z)\mathrm{d}z </math></center>

=== Far Field Wave ===
One of the primary objecives of wavemaker theory is to determine the far-field wave amplitude <math> A \, </math> in terms of <math> f(z) \, </math>.
The far-field wave component representing progagating waves is given by:

<center><math> \lim_{x\to\infty} \phi = a_0 \phi_0(z) e^{-k_0 x} =
a_0 \frac{\cos k_0(z+h)}{\cos k_0 h } e^{-k_0 x} </math></center>

Note that <math> k_0 \, </math> is imaginery. We therefore obtain the complex amplitude of the propagating wave at infinity, namely modulus and phase, in terms of the wave maker displacement <math> f(z) \, </math>.

For what type of <math> f(z) \, </math> are the non-wavelike modes zero? It is easy to verify by virtue of orthogonality that

<center><math> f(z) \ \sim \ \phi_0 (z) </math></center>

Unfortunately this is not a "practical" displacement since <math> \phi_0 (z) \, </math> depends on <math> \omega\, </math>, so one would need to build a flexible paddle.

== Matlab Code ==

A program to calculate the coefficients for the wave maker problems can be found here
[http://www.math.auckland.ac.nz/~meylan/code/eigenfunction_matching/wavemaker.m wavemaker.m]

=== Additional code ===

This program requires
[http://www.math.auckland.ac.nz/~meylan/code/dispersion/dispersion_free_surface.m dispersion_free_surface.m]
to run

-----

This article is based in part on the MIT open course notes and the original article can be found
[http://ocw.mit.edu/NR/rdonlyres/Mechanical-Engineering/2-24Spring-2002/737E217E-0582-45D5-B1F9-B2ECF977C66E/0/lecture6.pdf here]

[[Ocean Wave Interaction with Ships and Offshore Energy Systems]]

[[Category:Eigenfunction Matching Method]]
[[Category:Pages with Matlab Code]]
[[Category:Complete Pages]]

Introduction to KdV

2022-08-12T02:31:03Z

Meylan: /* Formula for the solitary wave */

{{nonlinear waves course
| chapter title = Introduction to KdV
| next chapter = [[Numerical Solution of the KdV]]
| previous chapter = [[Nonlinear Shallow Water Waves]]
}}

{{complete pages}}

The KdV (Korteweg-De Vries) equation is one of the most important non-linear
pde's. It was originally derived to model shallow water waves with weak
nonlinearities, but it has a wide variety of applications. The derivation of
the KdV is given in [[KdV Equation Derivation]]. The KdV equation
is written as
<center><math>
\partial _{t}u+6u\partial _{x}u+\partial _{x}^{3}u=0.
</math></center>

More information about it can be found at [http://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation Korteweg de Vries equation]

==Travelling Wave Solution==

The KdV equation posesses travelling wave solutions. One particular
travelling wave solution is called a soltion and it was discovered
experimentally by [http://en.wikipedia.org/wiki/John_Scott_Russell John Scott Russell]
in 1834. However, it was not understood
theoretically until the work of [http://en.wikipedia.org/wiki/Diederik_Korteweg Korteweg] and
[http://en.wikipedia.org/wiki/Gustav_de_Vries de Vries] in 1895.

We begin with the assumption that the wave travels with contant form, i.e.
is of the form
<center><math>
u\left( x,t\right) =f\left( x-ct\right)
</math></center>
Note that in this equation the parameter <math>c</math> is an unknown as is the
function <math>f.</math>
Only very special values of <math>c</math> and <math>f</math> will give travelling
waves.
We introduce the coordinate <math>\zeta = x - ct</math>.
If we substitute this expression into the KdV equation we obtain
<center><math>
-c\frac{\mathrm{d}f}{\mathrm{d}\zeta}+6f\frac{\mathrm{d}f}{\mathrm{d}\zeta}+\frac{\mathrm{d}^{3}f}{\mathrm{d}\zeta^{3}}=0
</math></center>
We can integrate this with respect to <math>\zeta</math> to obtain
<center><math>
-cf+3f^{2}+\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}=A_{1}
</math></center>
where <math>A</math> is a constant of integration.

If think about this equation as Newton's second law in a potential well <math>
V(f) </math> for which the equation is
<center><math>
\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}= -\frac{\mathrm{d}V}{\mathrm{d}f}
</math></center>
then the potential well is given by
<center><math>
V\left( f\right) =-A_{0}-A_{1}f-c\frac{f^{2}}{2}+f^{3}
</math></center>
Therefore our equation for <math>f</math> may be thought of as the motion of a particle
in a cubic well.

The constant <math>A_0</math> has no effect on our solution so we can set it to be zero.
We can choose the constant <math>A_{1}=0</math> and then we have a
maximum at <math>f=0</math>. There is a solution which rolls from this at <math>t=-\infty </math>
and then runs up the other side and finally returns to the maximum at <math>
t=\infty .</math> This corresponds to a solitary wave solution.

We can also think about the equation as a first order system using <math>
f^{^{\prime }}=v.</math> This gives us
<center><math>\begin{matrix}
\frac{\mathrm{d}f}{\mathrm{d}\zeta} &=&v \\
\frac{\mathrm{d}v}{\mathrm{d}\zeta} &=&A_{1}+cf-3f^{2}
\end{matrix}</math></center>
If we chose <math>A_{1}=0</math> then we obtain two equilibria at <math>(f,v)=\left(
0,0\right) </math> and <math>(3/c,0).</math> If we analysis these equilibria we find the
first is a saddle and the second is a nonlinear center (it is neither repelling nor
attracting).
The Jacobian matrix for the saddle point is
<center><math>
J_{\left( 0,0\right) }=\left(
\begin{array}{cc}
0 & 1 \\
c & 0
\end{array}
\right)
</math></center>
which has eigenvalues at <math>\pm \sqrt{c}</math> and the incident directions are
<center><math>
\left(
\begin{array}{c}
1\\
\sqrt{c}
\end{array}
\right) \leftrightarrow \sqrt{c},\left(
\begin{array}{c}
-1 \\
\sqrt{c}
\end{array}
\right) \leftrightarrow -\sqrt{c}
</math></center>
There is a
homoclinic connection which connects the equilibrium point at the origin. This holoclinic
connection represents the solitary wave. Within this homoclinic connection
lie periodic orbits which represent the cnoidal waves.

{| class="wikitable"
|-
! Phase portrait
! Solitary Wave
! Cnoidal Wave
|-
| [[Image:Kdv phase portrait.jpg|thumb|350px|Phase portrait]]
| [[Image:Kdv wave solitary2.gif|thumb|350px|Solitary Wave]]
| [[Image:Kdv wave cn.gif|thumb|350px|Cnoidal Wave]]
|}

We can also integrate the equation
<center><math>
-cf+3f^{2}+\frac{\mathrm{d}^{2}f}{\mathrm{d}\zeta^{2}}=A_{1}
</math></center>
by multiplying by <math>f^{\prime }</math>and integrating. This gives us
<center><math>
\frac{\left( f^{^{\prime }}\right) ^{2}}{2}=A_{0}+A_{1}f + c\frac{
f^{2}}{2}-f^{3} = -V(f)
</math></center>
It is no coincidence that the right hand side is the potential energy, because this
is nothing more that the equation for conservation of energy (or the first
integral of the Lagrangian system) which does not depend on <math>\zeta</math>.

This is a separable equation and the only challenge is to integrate
<center><math>
\int \frac{1}{\sqrt{-2V(f)}} \mathrm{d}f.
</math></center>

==Formula for the solitary wave==

We know that the solitary wave solution is found when <math>
A_{0}=A_{1}=0.</math> This gives us
<center><math>
\left( f^{^{\prime }}\right) ^{2}=f^{2}\left( c-2f\right)
</math></center>
This can be solved by separation of variables to give
<center><math>
\int \frac{\mathrm{d}f}{f\sqrt{c-2f}}=\int \mathrm{d}\zeta
</math></center>
We then substitute
<center><math>
f=\frac{1}{2}c\,\mathrm{sech}^{2}\left( s\right)
</math></center>
and note that
<center><math>
c-2f=c\left( 1-\mathrm{sech}^{2}\left( s\right) \right) =c\tanh ^{2}\left( f\right)
</math></center>
and that
<center><math>
\frac{\mathrm{d}f}{\mathrm{d}s}=-c\frac{\sinh \left( f\right) }{\cosh ^{3}\left( f\right) }
</math></center>
This means that
<center><math>\begin{matrix}
\int \frac{\mathrm{d}f}{f\sqrt{c-2f}} &=&-\frac{2}{\sqrt{c}}\int \frac{\sinh \left(
f\right) }{\mathrm{sech}^{2}\left( s\right) \tanh \left( f\right) \cosh ^{3}\left(
f\right) }\mathrm{d}s \\
&=&-\frac{2}{\sqrt{c}}s
\end{matrix}</math></center>
This gives us
<center><math>
-\frac{2}{\sqrt{c}}s=\zeta+a
</math></center>
Therefore
<center><math>
f=\frac{1}{2}c\,\mathrm{sech}^{2}\left( \frac{\sqrt{c}}{2}\left( \zeta+a\right) \right)
</math></center>
Of course we assumed that <math>x=x-ct</math> so the formula for the solitary wave is
given by
<center><math>
f\left( x-ct\right) =\frac{1}{2}c\,\mathrm{sech}^{2}\left[ \frac{\sqrt{c}}{2}\left(
x-ct+a\right) \right]
</math></center>
Note that a solution exists for each
<math>c</math>, and that the amplitude is proportional to <math>c.</math> All of this was
discovered experimentally by
[http://en.wikipedia.org/wiki/John_Scott_Russell John Scott Russell].

==Formula for the cnoidal wave==

If we consider the case when the solution oscillates between two values <math>F_2 < F_3</math>
(which we can assume are also roots of <math>V(f)</math> without loss of generality) then
we can integrate the equation to obtain
<center><math>
f\left( \zeta\right) =F_2+(F_3 - F_2) \mathrm{cn}^{2}\left( \gamma \zeta ;k\right)
</math></center>
where <math>cn</math> is a Jacobi Elliptic function and
<math>\gamma</math> and <math>k</math> are constants which depend on <math>c</math>.
Derivation of this equation is found [[KdV Cnoidal Wave Solutions]].
We can write this equation as
<center><math>
f\left( x-ct\right) =a+b \mathrm{cn}^{2}\left( \gamma (x-ct);k\right)
</math></center>
where <math>a=k^2\gamma^2</math> and <math>c = 6b + 4(2k^2 -1)\gamma^2</math>.
These waves are known as [http://en.wikipedia.org/wiki/Cnoidal_wave cnoidal waves].

In the limit the two solutions agree. We also obtain a sinusoidal solution in the limit of
small amplitude.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|SzZ-KhvvPio}}

=== Part 2 ===

{{#ev:youtube|QltlSQQBtrs}}

=== Part 3 ===

{{#ev:youtube|NJ7h3Z9QtvU}}

=== Part 4 ===

{{#ev:youtube|NictSlSgRbM}}

File:Burgers phase.jpg

2021-10-13T03:42:05Z

Meylan:

Code to make can be found [http://www.math.auckland.ac.nz/~meylan/code/nonlinear/phase_portrait_burgers.m phase_portrait_burgers.m]

Also needs
* {{gif_add_frame}}

Reaction-Diffusion Systems

2021-10-08T05:06:59Z

Meylan:

{{nonlinear waves course
| chapter title = Reaction-Diffusion Systems
| next chapter = [[Burgers Equation]]
| previous chapter = [[Example Calculations for the KdV and IST]]
}}

{{complete pages}}

We present here a brief theory of reaction diffusion waves.

== Law of Mass Action ==

The law of mass action states that equation rates are proportional to the concentration
of reacting species and the ratio in which they combined. It is discussed in detail in
[[Billingham and King 2000]]. We will present here a few simple examples.

=== Example 1: Simple Decay ===

Suppose we have of chemical <math>P</math> which decays to <math>A</math>, i.e.
<center>
<math>P \to A</math>
</center>
with rate <math>k[P]</math> where <math>[P]</math> denotes concentration. Then if we
set <math>p=[P]</math> and <math>a = [A] </math> we obtain the equations
<center>
<math> \frac{\mathrm{d}p}{\mathrm{d}t} = -kp\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}a}{\mathrm{d}t} = kp</math>
</center>
which has solution
<center>
<math>
p = p_0 e^{-kt}\,\,\,\textrm{and}\,\,\, a = a_0 + p_0(1-e^{-kt})
</math>
</center>
where <math>a_0</math> and <math>p_0</math> are the values of <math>a</math>
<math>p</math> repectively at <math>t=0</math>.

=== Example 2: Quadratic Autocatalysis ===

This example will be important when we consider reaction diffusion problems.
We consider the reaction
<center>
<math>A + B \to 2B</math>
</center>
with rate proportional to <math>k[A][B]</math>. If we define <math>a = [A]</math>
and <math>b = [B]</math> we obtain the following equations
<center>
<math> \frac{\mathrm{d}a}{\mathrm{d}t} = -kab\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}b}{\mathrm{d}t} = kab</math>
</center>
We can solve these equations by observing that
<center>
<math>
\frac{\mathrm{d}(a+b)}{\mathrm{d}t} = 0
</math>
</center>
so that <math>a + b = a_0 + b_0</math>. We can then eliminate <math>a</math> to obtain
<center>
<math>
\frac{\mathrm{d}b}{\mathrm{d}t} = k(a_0 + b_0 - b)b
</math>
</center>
which is separable with solution
<center>
<math>
b = \frac{b_0(a_0 + b_0)e^{k(a_0 + b_0)t}}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
and
<center>
<math>
a = \frac{a_0(a_0 + b_0)}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
Note that <math>a\to 0</math> and <math>b\to a_0 + b_0</math>
as <math>t\to \infty.</math>

== Diffusion ==

The equation for spatially homogeneous diffusion of a chemical with concentration
<math>c</math> is
<center>
<math>
\partial_t c = D\nabla^2 c
</math>
</center>
which is the [http://en.wikipedia.org/wiki/Heat_equation heat equation]. We will consider this in
only one spatial dimension. Consider it on the boundary <math>-\infty < x < \infty</math>. In this case
we can solve by the [http://en.wikipedia.org/wiki/Fourier_transform Fourier transform] and obtain
<center>
<math>
\partial_t \hat{c} = -D k^2 \hat{c}
</math>
</center>
where <math>\hat{c}</math> is the Fourier transform of <math>c</math>. This has solution
<center>
<math>
\hat{c} = \hat{c}_0 e^{-D k^2 t}
</math>
</center>
We can find the inverse transform using convolution and obtain
<center>
<math>
c(x,t) = \frac{1}{\sqrt{4\pi D t}} \int_{-\infty}^{\infty} c_0(x) e^{(x-s)^2/4Dt}\mathrm{d}s
</math>
</center>

=== Solution of the dispersion equation using FFT ===

We can solve the dispersion equation using the discrete Fourier transform and
its closely related numerical implementation the FFT (Fast Fourier Transform).
We have already met the FFT [[Numerical Solution of the KdV]] but we consider it here in more detail.
We consider the concentration
on the finite domain <math>-L \leq x \leq L</math> and use a
[http://en.wikipedia.org/wiki/Fourier_series Fourier series] expansion
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(t) e^{\mathrm{i} k_n x}
</math>
</center>
where <math>k_n = \pi n /L </math>. If we substitute this into the diffusion equation we obtain
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0)e^{-k_n^2 D t} e^{\mathrm{i} k_n x}
</math>
</center>
Note that this is not the same solution as we obtained on the infinite domain because
of the boundary conditions on the finite domain. The coefficients <math>\hat{c}_n(0)</math>
are found using the initial conditions so that
<center>
<math>
\hat{c}_n(0) = \frac{1}{2L} \int_{-L}^{L} e^{-\mathrm{i} k_n x} c_0(x) \mathrm{d}x
</math>
</center>

The key to the numerical solution of this equation is the use of the FFT. We begin by discretising the
domain into a series of <math>N</math> points <math>x_m = -L + 2Lm/N </math>. We then use this to
approximate the integral above and obtain
<center>
<math>
\hat{c}_n(0) = \frac{1}{N} \sum_{m=0}^{N-1} e^{-\mathrm{i} k_n x_m} c_0(x_m)
</math>
</center>
<center>
<math>
= \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} e^{\mathrm{i} \pi n} c_0(x_m)
</math>
</center>

We also get
<center>
<math>
c(x_m,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{\mathrm{i} k_n x_m}
</math>
</center>
<center>
<math>
= \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{2\mathrm{i} \pi nm/N} e^{-\mathrm{i} \pi n}
</math>
</center>
but we know that
<center>
<math>
\hat{c}_n(0) e^{-\mathrm{i} \pi n} = \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} c_0(x_m)
</math>
</center>

=== The discrete Fourier transform ===
The
[http://en.wikipedia.org/wiki/Discrete_Fourier_transform discrete Fourier transform]
of a sequence of ''2N'' complex numbers ''c''<sub>0</sub>, ..., ''c''<sub>''2N''−1</sub> is transformed into the sequence of ''N'' complex numbers <math>\hat{c}</math><sub>0</sub>, ..., <math>\hat{c}</math><sub>''N''−1</sub> by the DFT according to the formula:
<center>
<math>\hat{c}_m = \sum_{n=0}^{N-1} c_n e^{-2\pi \mathrm{i}mn/N} \quad \quad m = 0, \dots, N-1</math>
</center>

We denote the transform by the symbol <math>\mathcal{F}</math>, as in <math>\mathbf{X} = \mathcal{F} \left \{ \mathbf{x} \right \} </math> or <math>\mathcal{F} \left ( \mathbf{x} \right )</math> or <math>\mathcal{F} \mathbf{x}</math>.

The '''inverse discrete Fourier transform (IDFT)''' is given by
<center>
<math>c_n = \frac{1}{2N} \sum_{m=0}^{N-1} \hat{c}_m e^{2\pi \mathrm{i}mn/N} \quad \quad n = 0,\dots,N-1.</math>
</center>

Therefore we can write
<center>
<math>
c(x_m,t) = \mathcal{F} \left\{ e^{-k_n^2 D t} \mathcal{F}^{-1} \left\{ c_0(x_m) \right\}
\right\}
</math>
</center>

The only difficulty is that we need to define carefully the values of
<math>k_n</math>

The real power of this method lies with the [http://en.wikipedia.org/wiki/FFT Fast Fourier Transform]
or '''FFT''' algorithm. A naive implementation of the discrete Fourier transform above (or its inverse)
will involve order <math>N^2</math> operations. Using FFT algorithms, this can be reduced to order <math>N \log(N)</math>. This is an incredible speed up, for example
if N = 1024, FFT algorithms are more efficient by a factor of 147. This is the reason FFT
algorithms are used so extensively.

== Reaction Diffusion Equations ==

We consider an auto catalytic reaction where the chemical species also diffuse. In this
case the equations are
<center>
<math>\partial_t a = D\partial_x^2 a - kab</math>
</center>
<center>
<math>\partial_t b = D\partial_x^2 b + kab</math>
</center>
We can non-dimensionalise these equations scaling the variables as
<center>
<math>
z = x/x^*\,\,\,\tau = t/t^*\,\,\,\alpha = a/a_0\,\,\,\beta = b/a_0
</math>
</center>
So that the equations become
<center>
<math>
\frac{1}{t^*}\partial_\tau \alpha = \frac{a_0}{(x^*)^2}D\partial_z^2 \alpha
- k a_0^2 \alpha\beta
</math>
</center>
<center>
<math>
\frac{1}{t^*}\partial_\tau \beta = \frac{a_0}{(x^*)^2}D\partial_z^2 \beta
+ k a_0^2 \alpha\beta
</math>
</center>
If we choose
<center>
<math>
x^* = \sqrt{\frac{D}{ka_0}}\,\,\,t^*=\frac{1}{ka_0}
</math>
</center>
then we obtain the system
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
-\alpha\beta
</math>
</center>
<center>
<math>
\partial_\tau \beta =\partial_z^2 \beta
+ \alpha\beta
</math>
</center>

=== Solution via split step method ===

[[Image:Reaction diffusion.gif|thumb|right|500px|Solutions for <math>\alpha(z,0) =1
\, \beta(z,0) = \exp(-10z^2)</math>]]

We can solve this equations numerically using a
[http://en.wikipedia.org/wiki/Split-step_method split step method]. We assume
that at time <math>\tau</math> we know <math>\alpha(z,\tau)</math>
and <math>\beta(z,\tau)</math>. We then solve first the following equation
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
</math>
</center>
from <math>\tau</math> to <math>\tau + \Delta\tau</math>
(which we can do exactly using the spectral methods just discussed for
the dispersion equation). We write this solution as
<math>\tilde{\alpha}(z,\tau + \Delta\tau)</math>
Then we solve
<center>
<math>
\partial_\tau \alpha = -\alpha\beta
</math>
</center>
by assuming that <math>\beta</math> is constant and subject to the boundary
condition that <math>\alpha(z,\tau) = \tilde{\alpha}(z,\tau + \Delta\tau)</math>.
This gives
<center>
<math>
\alpha(z,\tau + \Delta\tau) = e^{-\beta(z,\tau) \Delta\tau} \tilde{\alpha}(z,\tau+ \Delta\tau)\,
</math>
</center>
and we do likewise for the equation for <math>\beta</math>. Note that
while both steps are exact the result from the split step method is an
approximation with error which becomes smaller as the step size becomes
smaller.

We can easily implement this split step method in matlab and we obtain
a pair of travelling waves.

== Travelling Waves solution ==

When we solve the equations we found the solution formed travelling waves and
we now consider this phenomena in detail.

We define a new coordinate <math>y = z - v\tau</math> (so we will consider only
waves travelling to the right, although we could analyse waves travelling to
the left in a similar fashion). We seek stationary solutions in
<math>\alpha(y)</math> and <math>\beta(y)</math> which satisfy
<center>
<math>
\frac{\mathrm{d}^2 \alpha}{\mathrm{d}y^2} + v \frac{\mathrm{d} \alpha}{\mathrm{d}y} = \alpha\beta
</math>
</center>
and
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} = -\alpha\beta
</math>
</center>
If we add these equations we obtain
<center>
<math>
\frac{\mathrm{d}^2 (\alpha+\beta)}{\mathrm{d}y^2} + v \frac{\mathrm{d} (\alpha+\beta)}{\mathrm{d}y} = 0
</math>
</center>
so that <math>\alpha + \beta = c_0 + c_1 e^{-vy}</math>. Boundary conditions
are that as <math>y\to\infty </math> <math>\alpha = 1</math> and
<math>\beta = 0</math> and if <math>y\to-\infty</math> then <math>\alpha = 0</math> and
<math>\beta = 1</math>. Therefore <math>\alpha + \beta = 1</math>.
This means that, since <math>\alpha \geq 0</math>, we must have
<math>0\leq \beta \leq 1</math>.
We can then obtain the following equation
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} + \beta(1-\beta)= 0
</math>
</center>
which we can write as the system of first order equations.

[[Image:R_d_phase_portrait.jpg|thumb|right|500px|Phase portrait for out system showing
the equilibrium points and the heteroclinic connection]]

We define the variable <math>\gamma = \frac{\mathrm{d}\beta}{\mathrm{d}y}</math>
and we obtain
<center><math>
\begin{align}
\frac{\mathrm{d}\beta}{\mathrm{d} y} &= \gamma&\\
\frac{\mathrm{d}\gamma}{\mathrm{d} y} &= -v\gamma + \beta(\beta -1)& \\
\end{align}
</math></center>
This dynamical system has equilibrium points at <math>(0,0)</math>
and <math>(1,0)</math>. We can analyse these equilibrium points by
linearization. The Jacobian matrix is
<center>
<math>
J =\begin{pmatrix}
0 & 1 \\
-1 + 2\beta & -v
\end{pmatrix}
</math>
</center>

We can easily see that the Jacobian evaluated at our first equilibrium point is
<center>
<math>
J_{(0,0)} =\begin{pmatrix}
0 & 1 \\
-1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\mu_{\pm} = -1/2 (v \mp \sqrt{v^2-4})</math>. Therefore
this point is a nodal sink (possibly a spiral)

[[Image:R_d_wave.gif|right|Travelling wave solution for v=2 and the position on
the heteroclinic connection]]

Additionally,
<center>
<math>
J_{(1,0)} =\begin{pmatrix}
0 & 1 \\
1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\lambda_{\pm} = -1/2 (v \mp \sqrt{v^2+4})</math>.
This is a a saddle point. The unstable and stable
separatrices leave the equilibrium point at <math>(1,0)</math> in the directions
<math> \begin{pmatrix}\lambda_{\pm} \\ 1\end{pmatrix}</math>. The only path on which <math>\beta</math> is bounded
as <math>y\to-\infty</math> are the unstable separatrices. Also, only the
unstable separatrix which enters the region <math>\beta<1</math> is physically meaningful.

To find a travelling wave we need to find a heteroclinic connection
between the two equilibrium points which also has to satisfy the conditions
that <math>0\leq \beta \leq 1</math>.

We need to show that the heteroclinic connection does not cross the <math>\beta</math> axis.
Consider the region
<center>
<math>
R = \left\{(\beta,\gamma)\,|\, \beta<1,\,-k\beta<\gamma<0\right\}
</math>
</center>
On the line <math>\beta = 1,d\beta/dy<0</math> and hence all flow it into <math>R</math>.
On the line <math>\gamma = 0, d\gamma/dy < 0</math> for <math>0<\beta<1</math>. On the line
<math>\gamma = -k \beta</math> we know that <math>d\beta/dy < 0</math> so that integral paths
enter the region if and only if <math>d\gamma/d\beta < \gamma/\beta</math>. We know that
<center>
<math>
d\gamma/d\beta - \gamma/\beta = -v - \frac{\beta(1-\beta)}{\gamma} -\frac{\gamma}{\beta}
= \frac{1}{k} (k^2 - vk +1 -\beta)
</math>
</center>
when <math>\gamma = -k\beta</math>. Therefore we need to find a value of <math>k</math>
so that <math>k^2 - vk +1 < 0</math>, which is possible provided <math>v\geq 2</math>, for example
<math>k = \dfrac{v}{2}</math>.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|rF4X42jP0v8}}

=== Part 2 ===

{{#ev:youtube|D0NwYlM-uOg}}

=== Part 3 ===

{{#ev:youtube|5IEZJtJaDHk}}

=== Part 4 ===

{{#ev:youtube|t_OjTSwVgdo}}

=== Part 5 ===

{{#ev:youtube|kTHVZaYezLk}}

=== Part 6 ===

{{#ev:youtube|MVuSg5_sfYI}}

=== Part 7 ===

{{#ev:youtube|hxKMOHyy6Bw}}

=== Part 8 ===

{{#ev:youtube|yQ-O2KIqu44}}

[[Category:Simple Nonlinear Waves]]

Burgers Equation

2020-10-24T02:14:18Z

Meylan: /* Exact Solution of Burgers equations */

{{nonlinear waves course
| chapter title = Burgers Equation
| next chapter =
| previous chapter = [[Reaction-Diffusion Systems]]
}}

==Introduction==

We have already met the conservation law for the traffic equations
<center><math>
\partial _{t}\rho +c\left( \rho \right) \partial _{x}\rho =0
</math></center>
and seen how this leads to shocks. We can smooth this equation by adding
dispersion to the equation to give us
<center><math>
\partial _{t}\rho +c\left( \rho \right) \partial _{x}\rho =\nu \partial
_{x}^{2}\rho
</math></center>
where <math>\nu >0.</math>

The simplest equation of this type is to write
<center><math>
\partial _{t}u+u\partial _{x}u=\nu \partial _{x}^{2}u
</math></center>
(changing variables to <math>u</math> and this equation is known as Burgers equation.

==Travelling Wave Solution==

We can find a travelling wave solution by assuming that
<center><math>
u\left( x,t\right) =u\left( x-ct\right) =u\left( \xi \right)
</math></center>
This leads to the equations
<center><math>
-cu^{\prime }+u^{\prime }u-\nu u^{\prime \prime }=0
</math></center>
We begin by looking at the phase plane for this system, writing <math>w=u^{\prime
}</math> so that
<center><math>\begin{matrix}
\dfrac{\mathrm{d}u}{\mathrm{d}\xi } &=&w \\
\dfrac{\mathrm{d}w}{\mathrm{d}\xi } &=&\frac{1}{\nu }\left( w\left( u-c\right) \right)
\end{matrix}</math></center>
This is a degenerate system with the entire <math>u</math> axis being equilibria.

We can also solve this equation exactly as follows.
<center><math>
-cu^{\prime }+u^{\prime }u-\nu u^{\prime \prime }=0
</math></center>
can be integrated to give
<center><math>
-cu+\frac{1}{2}\left( u\right) ^{2}-\nu u^{\prime} =c_{1}
</math></center>
which can be rearranged to give
<center><math>
u^{\prime }=\frac{1}{2\nu }\left( \left( u\right) ^{2}-2cu-2c_{1}\right)
</math></center>
We define the two roots of the quadratic <math>\left( u\right) ^{2}-2\nu
u-2c_{1}=0</math> by <math>u_{1}</math> and <math>u_{2}</math>
and we assume that <math>u_{2} < u_{1}</math>. Note that there is only a bounded
solution if we have two real roots and for the bounded solution
<math>u_{2} < u < u_{1}</math>.
We note that the wave speed
is
<center><math>
c=\frac{1}{2}\left( u_{1}+u_{2}\right)
</math></center>
The equation can therefore be written as
<center><math>
2\nu u^{\prime }=\left( u-u_{1}\right) \left( u-u_{2}\right)
</math></center>
which has solution
<center><math>
u\left( \xi \right) =\frac{1}{2}\left( u_{1}+u_{2}\right) -\frac{1}{2}\left(
u_{1}-u_{2}\right) \tanh \left[ \left( \frac{\xi }{4\nu }\right) \left(
u_{1}-u_{2}\right) \right]
</math></center>

==Numerical Solution of Burgers equation==

We can solve the equation using our split step spectral method. The equation
can be written as
<center><math>
\partial _{t}u=-\frac{1}{2}\partial _{x}\left( u^{2}\right) +\nu \partial
_{x}^{2}u
</math></center>
We solve this by solving in Fourier space to give
<center><math>
\partial _{t}\hat{u}=-\frac{1}{2}ik \widehat{\left( u^{2}\right)} -\nu k^{2}\hat{u}
</math></center>
Then we solve each of the steps in turn
for a small time interval to give
<center><math>\begin{matrix}
\tilde{u}\left( k,t+\Delta t\right) &=&\hat{u}\left( k,t\right) -\frac{
\Delta t}{2}ik\mathcal{F}\left( \left[ \mathcal{F}^{-1}\hat{u}\left(
k,t\right) \right] ^{2}\right) \\
\hat{u}\left( k,t+\Delta t\right) &=&\tilde{u}\left( k,t+\Delta t\right)
\exp \left( -\nu k^{2}\Delta t\right)
\end{matrix}</math></center>

{| class="wikitable"
|-
! Phase plane for a travelling wave solution
! Numerical solution of Burgers equation
|-
| [[Image:Burgers_phase.jpg|thumb|right|500px|Phase plane for a travelling wave solution of Burgers equation]]
| [[Image:File-Burgers2.gif|thumb|right|500px| Numerical solution of Burgers equation]]

|}

==Exact Solution of Burgers equations==

We can find an exact solution to Burgers equation. We want to solve
<center><math>\begin{matrix}
\partial _{t}u+u\partial _{x}u &=&\nu \partial _{x}^{2}u \\
u\left( x,0\right) &=&F\left( x\right)
\end{matrix}</math></center>
Frist we write the equation as
<center><math>
\partial _{t}u+\partial _{x}\left( \frac{u^{2}}{2}-\nu \partial _{x}u\right)
=0
</math></center>
We want to find a function <math>\psi \left( x,t\right) </math> such that
<center><math>
\partial _{x}\psi =u,\ \ \partial _{t}\psi =\nu \partial _{x}u-\frac{u^{2}}{2
}
</math></center>
Note that because <math>\partial _{x}\partial _{t}\psi =\partial _{t}\partial
_{x}\psi </math> we will satisfy Burgers equation. This gives us the following
equation for <math>\psi </math>
<center><math>
\partial _{t}\psi =\nu \partial _{x}^{2}\psi -\frac{1}{2}\left( \partial
_{x}\psi \right) ^{2}
</math></center>
We introduce the ''Cole-Hopf '' transformation
<center><math>
\psi =-2\nu \log \left( \phi \right)
</math></center>
From this we can obtain the three results:
<center><math>
\begin{align}
\partial _{x}\psi &=-2\nu \frac{\partial _{x}\phi }{\phi } \\
\partial _{x}^{2}\psi &=2\nu \left( \frac{\partial _{x}\phi }{\phi }\right)
^{2}-\frac{2\nu }{\phi }\partial _{x}^{2}\phi \\
\partial _{t}\psi &=-2\nu \frac{\partial _{t}\phi }{\phi }
\end{align}
</math></center>
Therefore
<center><math>
\partial _{t}\psi =\nu \partial _{x}^{2}\psi -\frac{1}{2}\left( \partial
_{x}\psi \right) ^{2}
</math></center>
becomes
<center><math>
-2\nu \frac{\partial _{t}\phi }{\phi }=2\nu ^{2}\left( \frac{\partial
_{x}\phi }{\phi }\right) ^{2}
-2\nu^2 \frac{\partial_x^2\phi}{\phi}
-\frac{1}{2}\left( 2\nu \frac{\partial _{x}\phi
}{\phi }\right) ^{2}
</math></center>
or
<center><math>
\partial _{t}\phi =\nu \partial _{x}^{2}\phi
</math></center>
which is just the diffusion equation. Note that we also have to transform the
boundary conditions. We have
<center><math>
F\left( x\right) =u\left( x,0\right) =-2\nu \frac{\partial _{x}\phi \left(
x,0\right) }{\phi \left( x,0\right) }
</math></center>
We can write this as
<center><math>
\frac{\mathrm{d}}{\mathrm{d}x}\left( \log \left( \phi \right) \right) =-\frac{1}{2\nu }F\left(
x\right)
</math></center>
which has solution
<center><math>
\phi \left( x,0\right) =\Phi \left( x\right) =\exp \left( -\frac{1}{2\nu }
\int_{0}^{x}F\left( s\right) \mathrm{d}s\right)
</math></center>
We need to solve
<center><math>\begin{matrix}
\partial _{t}\phi &=&\nu \partial _{x}^{2}\phi \\
\phi \left( x,0\right) &=&\Phi \left( x\right)
\end{matrix}</math></center>
We take the Fourier transform and obtain
<center><math>\begin{matrix}
\partial _{t}\hat{\phi} &=&-k^{2}\nu \hat{\phi} \\
\hat{\phi}\left( k,0\right) &=&\hat{\Phi}\left( k\right)
\end{matrix}</math></center>
which has solution
<center><math>
\hat{\phi}\left( k,t\right) =\hat{\Phi}\left( k\right) e^{-k^{2}\nu t}
</math></center>
We can then use the convolution theorem to write
<center><math>\begin{matrix}
\phi \left( x,t\right) &=&\Phi \left( x\right) * \mathcal{F}^{-1}\left[
e^{-k^{2}\nu t}\right] \\
&=&\frac{1}{2\sqrt{\pi \nu t}}\int_{-\infty }^{\infty }\Phi \left( y\right)
\exp \left[ -\frac{\left( x-y\right) ^{2}}{4\nu t}\right] \mathrm{d}y
\end{matrix}</math></center>
Which can be expressed as
<center><math>
\phi \left( x,t\right) =\frac{1}{2\sqrt{\pi \nu t}}\int_{-\infty }^{\infty
}\exp \left[ -\frac{f}{2\nu }\right] \mathrm{d}y
</math></center>
where
<center><math>
f\left( x,y,t\right) =\frac{1}{2\nu }\int_{0}^{y}F\left( s\right) \mathrm{d}s+\frac{
\left( x-y\right) ^{2}}{2t}
</math></center>
To find <math>u</math> we recall that
<center><math>\begin{matrix}
u\left( x,t\right) &=&-2\nu \dfrac{\partial _{x}\phi \left( x,t\right) }{\phi
\left( x,t\right) } \\
&=&\dfrac{\int_{-\infty }^{\infty }\left( \frac{x-y}{t}\right) \exp \left[ -
\dfrac{f}{2\nu }\right] \mathrm{d}y}{\int_{-\infty }^{\infty }\exp \left[ -\frac{f}{
2\nu }\right] \mathrm{d}y}
\end{matrix}</math></center>

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|tVXQmxOG_6Y}}

=== Part 2 ===

{{#ev:youtube|hzgpMM_wWts}}

=== Part 3 ===

{{#ev:youtube|uH4B1XsGB-0}}

=== Part 4 ===

{{#ev:youtube|h6aDmCtJygM}}

=== Part 5 ===

{{#ev:youtube|CsnUKrLjtyQ}}

Reaction-Diffusion Systems

2020-10-21T01:28:05Z

Meylan: /* Lecture Videos */

{{nonlinear waves course
| chapter title = Reaction-Diffusion Systems
| next chapter = [[Burgers Equation]]
| previous chapter = [[Example Calculations for the KdV and IST]]
}}

{{complete pages}}

We present here a brief theory of reaction diffusion waves.

== Law of Mass Action ==

The law of mass action states that equation rates are proportional to the concentration
of reacting species and the ratio in which they combined. It is discussed in detail in
[[Billingham and King 2000]]. We will present here a few simple examples.

=== Example 1: Simple Decay ===

Suppose we have of chemical <math>P</math> which decays to <math>A</math>, i.e.
<center>
<math>P \to A</math>
</center>
with rate <math>k[P]</math> where <math>[P]</math> denotes concentration. Then if we
set <math>p=[P]</math> and <math>a = [A] </math> we obtain the equations
<center>
<math> \frac{\mathrm{d}p}{\mathrm{d}t} = -kp\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}a}{\mathrm{d}t} = kp</math>
</center>
which has solution
<center>
<math>
p = p_0 e^{-kt}\,\,\,\textrm{and}\,\,\, a = a_0 + p_0(1-e^{-kt})
</math>
</center>
where <math>a_0</math> and <math>p_0</math> are the values of <math>a</math>
<math>p</math> repectively at <math>t=0</math>.

=== Example 2: Quadratic Autocatalysis ===

This example will be important when we consider reaction diffusion problems.
We consider the reaction
<center>
<math>A + B \to 2B</math>
</center>
with rate proportional to <math>k[A][B]</math>. If we define <math>a = [A]</math>
and <math>b = [B]</math> we obtain the following equations
<center>
<math> \frac{\mathrm{d}a}{\mathrm{d}t} = -kab\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}b}{\mathrm{d}t} = kab</math>
</center>
We can solve these equations by observing that
<center>
<math>
\frac{\mathrm{d}(a+b)}{\mathrm{d}t} = 0
</math>
</center>
so that <math>a + b = a_0 + b_0</math>. We can then eliminate <math>a</math> to obtain
<center>
<math>
\frac{\mathrm{d}b}{\mathrm{d}t} = k(a_0 + b_0 - b)b
</math>
</center>
which is separable with solution
<center>
<math>
b = \frac{b_0(a_0 + b_0)e^{k(a_0 + b_0)t}}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
and
<center>
<math>
a = \frac{a_0(a_0 + b_0)}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
Note that <math>a\to 0</math> and <math>b\to a_0 + b_0</math>
as <math>t\to \infty.</math>

== Diffusion ==

The equation for spatially homogeneous diffusion of a chemical with concentration
<math>c</math> is
<center>
<math>
\partial_t c = D\nabla^2 c
</math>
</center>
which is the [http://en.wikipedia.org/wiki/Heat_equation heat equation]. We will consider this in
only one spatial dimension. Consider it on the boundary <math>-\infty < x < \infty</math>. In this case
we can solve by the [http://en.wikipedia.org/wiki/Fourier_transform Fourier transform] and obtain
<center>
<math>
\partial_t \hat{c} = -D k^2 \hat{c}
</math>
</center>
where <math>\hat{c}</math> is the Fourier transform of <math>c</math>. This has solution
<center>
<math>
\hat{c} = \hat{c}_0 e^{-D k^2 t}
</math>
</center>
We can find the inverse transform using convolution and obtain
<center>
<math>
c(x,t) = \frac{1}{\sqrt{4\pi D t}} \int_{-\infty}^{\infty} c_0(x) e^{(x-s)^2/4Dt}\mathrm{d}s
</math>
</center>

=== Solution of the dispersion equation using FFT ===

We can solve the dispersion equation using the discrete Fourier transform and
its closely related numerical implementation the FFT (Fast Fourier Transform).
We have already met the FFT [[Numerical Solution of the KdV]] but we consider it here in more detail.
We consider the concentration
on the finite domain <math>-L \leq x \leq L</math> and use a
[http://en.wikipedia.org/wiki/Fourier_series Fourier series] expansion
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(t) e^{\mathrm{i} k_n x}
</math>
</center>
where <math>k_n = \pi n /L </math>. If we substitute this into the diffusion equation we obtain
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0)e^{-k_n^2 D t} e^{\mathrm{i} k_n x}
</math>
</center>
Note that this is not the same solution as we obtained on the infinite domain because
of the boundary conditions on the finite domain. The coefficients <math>\hat{c}_n(0)</math>
are found using the initial conditions so that
<center>
<math>
\hat{c}_n(0) = \frac{1}{2L} \int_{-L}^{L} e^{-\mathrm{i} k_n x} c_0(x) \mathrm{d}x
</math>
</center>

The key to the numerical solution of this equation is the use of the FFT. We begin by discretising the
domain into a series of <math>N</math> points <math>x_m = -L + 2Lm/N </math>. We then use this to
approximate the integral above and obtain
<center>
<math>
\hat{c}_n(0) = \frac{1}{N} \sum_{m=0}^{N-1} e^{-\mathrm{i} k_n x_m} c_0(x_m)
</math>
</center>
<center>
<math>
= \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} e^{\mathrm{i} \pi n} c_0(x_m)
</math>
</center>

We also get
<center>
<math>
c(x_m,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{\mathrm{i} k_n x_m}
</math>
</center>
<center>
<math>
= \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{2\mathrm{i} \pi nm/N} e^{-\mathrm{i} \pi n}
</math>
</center>
but we know that
<center>
<math>
\hat{c}_n(0) e^{-\mathrm{i} \pi n} = \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} c_0(x_m)
</math>
</center>

=== The discrete Fourier transform ===
The
[http://en.wikipedia.org/wiki/Discrete_Fourier_transform discrete Fourier transform]
of a sequence of ''2N'' complex numbers ''c''<sub>0</sub>, ..., ''c''<sub>''2N''−1</sub> is transformed into the sequence of ''N'' complex numbers <math>\hat{c}</math><sub>0</sub>, ..., <math>\hat{c}</math><sub>''N''−1</sub> by the DFT according to the formula:
<center>
<math>\hat{c}_m = \sum_{n=0}^{N-1} c_n e^{-2\pi \mathrm{i}mn/N} \quad \quad m = 0, \dots, N-1</math>
</center>

We denote the transform by the symbol <math>\mathcal{F}</math>, as in <math>\mathbf{X} = \mathcal{F} \left \{ \mathbf{x} \right \} </math> or <math>\mathcal{F} \left ( \mathbf{x} \right )</math> or <math>\mathcal{F} \mathbf{x}</math>.

The '''inverse discrete Fourier transform (IDFT)''' is given by
<center>
<math>c_n = \frac{1}{2N} \sum_{m=0}^{N-1} \hat{c}_m e^{2\pi \mathrm{i}mn/N} \quad \quad n = 0,\dots,N-1.</math>
</center>

Therefore we can write
<center>
<math>
c(x_m,t) = \mathcal{F} \left\{ e^{-k_n^2 D t} \mathcal{F}^{-1} \left\{ c_0(x_m) \right\}
\right\}
</math>
</center>

The only difficulty is that we need to define carefully the values of
<math>k_n</math>

The real power of this method lies with the [http://en.wikipedia.org/wiki/FFT Fast Fourier Transform]
or '''FFT''' algorithm. A naive implementation of the discrete Fourier transform above (or its inverse)
will involve order <math>N^2</math> operations. Using FFT algorithms, this can be reduced to order <math>N \log(N)</math>. This is an incredible speed up, for example
if N = 1024, FFT algorithms are more efficient by a factor of 147. This is the reason FFT
algorithms are used so extensively.

== Reaction Diffusion Equations ==

We consider an auto catalytic reaction where the chemical species also diffuse. In this
case the equations are
<center>
<math>\partial_t a = D\partial_x^2 a - kab</math>
</center>
<center>
<math>\partial_t b = D\partial_x^2 b + kab</math>
</center>
We can non-dimensionalise these equations scaling the variables as
<center>
<math>
z = x/x^*\,\,\,\tau = t/t^*\,\,\,\alpha = a/a_0\,\,\,\beta = b/a_0
</math>
</center>
So that the equations become
<center>
<math>
\frac{1}{t^*}\partial_\tau \alpha = \frac{a_0}{(x^*)^2}D\partial_z^2 \alpha
- k a_0^2 \alpha\beta
</math>
</center>
<center>
<math>
\frac{1}{t^*}\partial_\tau \beta = \frac{a_0}{(x^*)^2}D\partial_z^2 \beta
+ k a_0^2 \alpha\beta
</math>
</center>
If we choose
<center>
<math>
x^* = \sqrt{\frac{D}{ka_0}}\,\,\,t^*=\frac{1}{ka_0}
</math>
</center>
then we obtain the system
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
-\alpha\beta
</math>
</center>
<center>
<math>
\partial_\tau \beta =\partial_z^2 \beta
+ \alpha\beta
</math>
</center>

=== Solution via split step method ===

[[Image:Reaction diffusion.gif|thumb|right|500px|Solutions for <math>\alpha(z,0) =1
\, \beta(z,0) = \exp(-10z^2)</math>]]

We can solve this equations numerically using a
[http://en.wikipedia.org/wiki/Split-step_method split step method]. We assume
that at time <math>\tau</math> we know <math>\alpha(z,\tau)</math>
and <math>\beta(z,\tau)</math>. We then solve first the following equation
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
</math>
</center>
from <math>\tau</math> to <math>\tau + \Delta\tau</math>
(which we can do exactly using the spectral methods just discussed for
the dispersion equation). We write this solution as
<math>\tilde{\alpha}(z,\tau + \Delta\tau)</math>
Then we solve
<center>
<math>
\partial_\tau \alpha = -\alpha\beta
</math>
</center>
by assuming that <math>\beta</math> is constant and subject to the boundary
condition that <math>\alpha(z,\tau) = \tilde{\alpha}(z,\tau + \Delta\tau)</math>.
This gives
<center>
<math>
\alpha(z,\tau + \Delta\tau) = e^{-\beta(z,\tau) \Delta\tau} \tilde{\alpha}(z,\tau+ \Delta\tau)\,
</math>
</center>
and we do likewise for the equation for <math>\beta</math>. Note that
while both steps are exact the result from the split step method is an
approximation with error which becomes smaller as the step size becomes
smaller.

We can easily implement this split step method in matlab and we obtain
a pair of travelling waves.

== Travelling Waves solution ==

When we solve the equations we found the solution formed travelling waves and
we now consider this phenomena in detail.

We define a new coordinate <math>y = z - v\tau</math> (so we will consider only
waves travelling to the right, although we could analyse waves travelling to
the left in a similar fashion). We seek stationary solutions in
<math>\alpha(y)</math> and <math>\beta(y)</math> which satisfy
<center>
<math>
\frac{\mathrm{d}^2 \alpha}{\mathrm{d}y^2} + v \frac{\mathrm{d} \alpha}{\mathrm{d}y} = \alpha\beta
</math>
</center>
and
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} = -\alpha\beta
</math>
</center>
If we add these equations we obtain
<center>
<math>
\frac{\mathrm{d}^2 (\alpha+\beta)}{\mathrm{d}y^2} + v \frac{\mathrm{d} (\alpha+\beta)}{\mathrm{d}y} = 0
</math>
</center>
so that <math>\alpha + \beta = c_0 + c_1 e^{-vy}</math>. Boundary conditions
are that as <math>y\to\infty </math> <math>\alpha = 1</math> and
<math>\beta = 0</math> and if <math>y\to-\infty</math> then <math>\alpha = 0</math> and
<math>\beta = 1</math>. Therefore <math>\alpha + \beta = 1</math>.
This means that, since <math>\alpha \geq 0</math>, we must have
<math>0\leq \beta \leq 1</math>.
We can then obtain the following equation
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} + \beta(1-\beta)= 0
</math>
</center>
which we can write as the system of first order equations.

[[Image:R_d_phase_portrait.jpg|thumb|right|500px|Phase portrait for out system showing
the equilibrium points and the heteroclinic connection]]

We define the variable <math>\gamma = \frac{\mathrm{d}\beta}{\mathrm{d}y}</math>
and we obtain
<center><math>
\begin{align}
\frac{\mathrm{d}\beta}{\mathrm{d} y} &= \gamma&\\
\frac{\mathrm{d}\gamma}{\mathrm{d} y} &= -v\gamma + \beta(\beta -1)& \\
\end{align}
</math></center>
This dynamical system has equilibrium points at <math>(0,0)</math>
and <math>(1,0)</math>. We can analyse these equilibrium points by
linearization. The Jacobian matrix is
<center>
<math>
J =\begin{pmatrix}
0 & 1 \\
-1 + 2\beta & -v
\end{pmatrix}
</math>
</center>

We can easily see that the Jacobian evaluated at our first equilibrium point is
<center>
<math>
J_{(0,0)} =\begin{pmatrix}
0 & 1 \\
-1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\mu_{\pm} = -1/2 (v \mp \sqrt{v^2-4})</math>. Therefore
this point is a nodal sink (possibly a spiral)

[[Image:R_d_wave.gif|right|Travelling wave solution for v=2 and the position on
the heteroclinic connection]]

Additionally,
<center>
<math>
J_{(1,0)} =\begin{pmatrix}
0 & 1 \\
1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\lambda_{\pm} = -1/2 (v \mp \sqrt{v^2+4})</math>.
This is a a saddle point. The unstable and stable
separatrices leave the equilibrium point at <math>(1,0)</math> in the directions
<math> \begin{pmatrix}\lambda_{\pm} \\ 1\end{pmatrix}</math>. The only path on which <math>\beta</math> is bounded
as <math>y\to-\infty</math> are the unstable separatrices. Also, only the
unstable separatrix which enters the region <math>\beta<1</math> is physically meaningful.

To find a travelling wave we need to find a heteroclinic connection
between the two equilibrium points which also has to satisfy the conditions
that <math>0\leq \beta \leq 1</math>.

We need to show that the heteroclinic connection does not cross the <math>\beta</math> axis.
Consider the region
<center>
<math>
R = \left\{(\beta,\gamma)\,|\, \beta<1,\,-k\beta<\gamma<0\right\}
</math>
</center>
On the line <math>\beta = 1,d\beta/dy<0</math> and hence all flow it into <math>R</math>.
On the line <math>\gamma = 0, d\gamma/dy < 0</math> for <math>0<\beta<1</math>. On the line
<math>\gamma = -k \beta</math> we know that <math>d\beta/dy < 0</math> so that integral paths
enter the region if and only if <math>d\gamma/d\beta < \gamma/\beta</math>. We know that
<center>
<math>
d\gamma/d\beta - \gamma/\beta = -v - \frac{\beta(1-\beta)}{\gamma} -\frac{\gamma}{\beta}
= \frac{1}{k} (k^2 - vk +1 -\beta)
</math>
</center>
when <math>\gamma = -k\beta</math>. Therefore we need to find a value of <math>k</math>
so that <math>k^2 - vk +1 < 0</math>, which is possible provided <math>v\geq 2</math>, for example
<math>k = \dfrac{v}{2}</math>.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|rF4X42jP0v8}}

=== Part 2 ===

{{#ev:youtube|MR7xSsJ_zLs}}

=== Part 3 ===

{{#ev:youtube|l9tXlpnlqn8}}

=== Part 4 ===

{{#ev:youtube|Ip353jAdBkk}}

=== Part 5 ===

{{#ev:youtube|kTHVZaYezLk}}

=== Part 6 ===

{{#ev:youtube|MVuSg5_sfYI}}

=== Part 7 ===

{{#ev:youtube|hxKMOHyy6Bw}}

=== Part 8 ===

{{#ev:youtube|t3TiemLMxDg}}

[[Category:Simple Nonlinear Waves]]

Reaction-Diffusion Systems

2020-10-21T01:08:01Z

Meylan: /* Travelling Waves solution */

{{nonlinear waves course
| chapter title = Reaction-Diffusion Systems
| next chapter = [[Burgers Equation]]
| previous chapter = [[Example Calculations for the KdV and IST]]
}}

{{complete pages}}

We present here a brief theory of reaction diffusion waves.

== Law of Mass Action ==

The law of mass action states that equation rates are proportional to the concentration
of reacting species and the ratio in which they combined. It is discussed in detail in
[[Billingham and King 2000]]. We will present here a few simple examples.

=== Example 1: Simple Decay ===

Suppose we have of chemical <math>P</math> which decays to <math>A</math>, i.e.
<center>
<math>P \to A</math>
</center>
with rate <math>k[P]</math> where <math>[P]</math> denotes concentration. Then if we
set <math>p=[P]</math> and <math>a = [A] </math> we obtain the equations
<center>
<math> \frac{\mathrm{d}p}{\mathrm{d}t} = -kp\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}a}{\mathrm{d}t} = kp</math>
</center>
which has solution
<center>
<math>
p = p_0 e^{-kt}\,\,\,\textrm{and}\,\,\, a = a_0 + p_0(1-e^{-kt})
</math>
</center>
where <math>a_0</math> and <math>p_0</math> are the values of <math>a</math>
<math>p</math> repectively at <math>t=0</math>.

=== Example 2: Quadratic Autocatalysis ===

This example will be important when we consider reaction diffusion problems.
We consider the reaction
<center>
<math>A + B \to 2B</math>
</center>
with rate proportional to <math>k[A][B]</math>. If we define <math>a = [A]</math>
and <math>b = [B]</math> we obtain the following equations
<center>
<math> \frac{\mathrm{d}a}{\mathrm{d}t} = -kab\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}b}{\mathrm{d}t} = kab</math>
</center>
We can solve these equations by observing that
<center>
<math>
\frac{\mathrm{d}(a+b)}{\mathrm{d}t} = 0
</math>
</center>
so that <math>a + b = a_0 + b_0</math>. We can then eliminate <math>a</math> to obtain
<center>
<math>
\frac{\mathrm{d}b}{\mathrm{d}t} = k(a_0 + b_0 - b)b
</math>
</center>
which is separable with solution
<center>
<math>
b = \frac{b_0(a_0 + b_0)e^{k(a_0 + b_0)t}}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
and
<center>
<math>
a = \frac{a_0(a_0 + b_0)}{a_0 + b_0e^{k(a_0 + b_0)t}}
</math>
</center>
Note that <math>a\to 0</math> and <math>b\to a_0 + b_0</math>
as <math>t\to \infty.</math>

== Diffusion ==

The equation for spatially homogeneous diffusion of a chemical with concentration
<math>c</math> is
<center>
<math>
\partial_t c = D\nabla^2 c
</math>
</center>
which is the [http://en.wikipedia.org/wiki/Heat_equation heat equation]. We will consider this in
only one spatial dimension. Consider it on the boundary <math>-\infty < x < \infty</math>. In this case
we can solve by the [http://en.wikipedia.org/wiki/Fourier_transform Fourier transform] and obtain
<center>
<math>
\partial_t \hat{c} = -D k^2 \hat{c}
</math>
</center>
where <math>\hat{c}</math> is the Fourier transform of <math>c</math>. This has solution
<center>
<math>
\hat{c} = \hat{c}_0 e^{-D k^2 t}
</math>
</center>
We can find the inverse transform using convolution and obtain
<center>
<math>
c(x,t) = \frac{1}{\sqrt{4\pi D t}} \int_{-\infty}^{\infty} c_0(x) e^{(x-s)^2/4Dt}\mathrm{d}s
</math>
</center>

=== Solution of the dispersion equation using FFT ===

We can solve the dispersion equation using the discrete Fourier transform and
its closely related numerical implementation the FFT (Fast Fourier Transform).
We have already met the FFT [[Numerical Solution of the KdV]] but we consider it here in more detail.
We consider the concentration
on the finite domain <math>-L \leq x \leq L</math> and use a
[http://en.wikipedia.org/wiki/Fourier_series Fourier series] expansion
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(t) e^{\mathrm{i} k_n x}
</math>
</center>
where <math>k_n = \pi n /L </math>. If we substitute this into the diffusion equation we obtain
<center>
<math>
c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0)e^{-k_n^2 D t} e^{\mathrm{i} k_n x}
</math>
</center>
Note that this is not the same solution as we obtained on the infinite domain because
of the boundary conditions on the finite domain. The coefficients <math>\hat{c}_n(0)</math>
are found using the initial conditions so that
<center>
<math>
\hat{c}_n(0) = \frac{1}{2L} \int_{-L}^{L} e^{-\mathrm{i} k_n x} c_0(x) \mathrm{d}x
</math>
</center>

The key to the numerical solution of this equation is the use of the FFT. We begin by discretising the
domain into a series of <math>N</math> points <math>x_m = -L + 2Lm/N </math>. We then use this to
approximate the integral above and obtain
<center>
<math>
\hat{c}_n(0) = \frac{1}{N} \sum_{m=0}^{N-1} e^{-\mathrm{i} k_n x_m} c_0(x_m)
</math>
</center>
<center>
<math>
= \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} e^{\mathrm{i} \pi n} c_0(x_m)
</math>
</center>

We also get
<center>
<math>
c(x_m,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{\mathrm{i} k_n x_m}
</math>
</center>
<center>
<math>
= \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{2\mathrm{i} \pi nm/N} e^{-\mathrm{i} \pi n}
</math>
</center>
but we know that
<center>
<math>
\hat{c}_n(0) e^{-\mathrm{i} \pi n} = \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} c_0(x_m)
</math>
</center>

=== The discrete Fourier transform ===
The
[http://en.wikipedia.org/wiki/Discrete_Fourier_transform discrete Fourier transform]
of a sequence of ''2N'' complex numbers ''c''<sub>0</sub>, ..., ''c''<sub>''2N''−1</sub> is transformed into the sequence of ''N'' complex numbers <math>\hat{c}</math><sub>0</sub>, ..., <math>\hat{c}</math><sub>''N''−1</sub> by the DFT according to the formula:
<center>
<math>\hat{c}_m = \sum_{n=0}^{N-1} c_n e^{-2\pi \mathrm{i}mn/N} \quad \quad m = 0, \dots, N-1</math>
</center>

We denote the transform by the symbol <math>\mathcal{F}</math>, as in <math>\mathbf{X} = \mathcal{F} \left \{ \mathbf{x} \right \} </math> or <math>\mathcal{F} \left ( \mathbf{x} \right )</math> or <math>\mathcal{F} \mathbf{x}</math>.

The '''inverse discrete Fourier transform (IDFT)''' is given by
<center>
<math>c_n = \frac{1}{2N} \sum_{m=0}^{N-1} \hat{c}_m e^{2\pi \mathrm{i}mn/N} \quad \quad n = 0,\dots,N-1.</math>
</center>

Therefore we can write
<center>
<math>
c(x_m,t) = \mathcal{F} \left\{ e^{-k_n^2 D t} \mathcal{F}^{-1} \left\{ c_0(x_m) \right\}
\right\}
</math>
</center>

The only difficulty is that we need to define carefully the values of
<math>k_n</math>

The real power of this method lies with the [http://en.wikipedia.org/wiki/FFT Fast Fourier Transform]
or '''FFT''' algorithm. A naive implementation of the discrete Fourier transform above (or its inverse)
will involve order <math>N^2</math> operations. Using FFT algorithms, this can be reduced to order <math>N \log(N)</math>. This is an incredible speed up, for example
if N = 1024, FFT algorithms are more efficient by a factor of 147. This is the reason FFT
algorithms are used so extensively.

== Reaction Diffusion Equations ==

We consider an auto catalytic reaction where the chemical species also diffuse. In this
case the equations are
<center>
<math>\partial_t a = D\partial_x^2 a - kab</math>
</center>
<center>
<math>\partial_t b = D\partial_x^2 b + kab</math>
</center>
We can non-dimensionalise these equations scaling the variables as
<center>
<math>
z = x/x^*\,\,\,\tau = t/t^*\,\,\,\alpha = a/a_0\,\,\,\beta = b/a_0
</math>
</center>
So that the equations become
<center>
<math>
\frac{1}{t^*}\partial_\tau \alpha = \frac{a_0}{(x^*)^2}D\partial_z^2 \alpha
- k a_0^2 \alpha\beta
</math>
</center>
<center>
<math>
\frac{1}{t^*}\partial_\tau \beta = \frac{a_0}{(x^*)^2}D\partial_z^2 \beta
+ k a_0^2 \alpha\beta
</math>
</center>
If we choose
<center>
<math>
x^* = \sqrt{\frac{D}{ka_0}}\,\,\,t^*=\frac{1}{ka_0}
</math>
</center>
then we obtain the system
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
-\alpha\beta
</math>
</center>
<center>
<math>
\partial_\tau \beta =\partial_z^2 \beta
+ \alpha\beta
</math>
</center>

=== Solution via split step method ===

[[Image:Reaction diffusion.gif|thumb|right|500px|Solutions for <math>\alpha(z,0) =1
\, \beta(z,0) = \exp(-10z^2)</math>]]

We can solve this equations numerically using a
[http://en.wikipedia.org/wiki/Split-step_method split step method]. We assume
that at time <math>\tau</math> we know <math>\alpha(z,\tau)</math>
and <math>\beta(z,\tau)</math>. We then solve first the following equation
<center>
<math>
\partial_\tau \alpha = \partial_z^2 \alpha
</math>
</center>
from <math>\tau</math> to <math>\tau + \Delta\tau</math>
(which we can do exactly using the spectral methods just discussed for
the dispersion equation). We write this solution as
<math>\tilde{\alpha}(z,\tau + \Delta\tau)</math>
Then we solve
<center>
<math>
\partial_\tau \alpha = -\alpha\beta
</math>
</center>
by assuming that <math>\beta</math> is constant and subject to the boundary
condition that <math>\alpha(z,\tau) = \tilde{\alpha}(z,\tau + \Delta\tau)</math>.
This gives
<center>
<math>
\alpha(z,\tau + \Delta\tau) = e^{-\beta(z,\tau) \Delta\tau} \tilde{\alpha}(z,\tau+ \Delta\tau)\,
</math>
</center>
and we do likewise for the equation for <math>\beta</math>. Note that
while both steps are exact the result from the split step method is an
approximation with error which becomes smaller as the step size becomes
smaller.

We can easily implement this split step method in matlab and we obtain
a pair of travelling waves.

== Travelling Waves solution ==

When we solve the equations we found the solution formed travelling waves and
we now consider this phenomena in detail.

We define a new coordinate <math>y = z - v\tau</math> (so we will consider only
waves travelling to the right, although we could analyse waves travelling to
the left in a similar fashion). We seek stationary solutions in
<math>\alpha(y)</math> and <math>\beta(y)</math> which satisfy
<center>
<math>
\frac{\mathrm{d}^2 \alpha}{\mathrm{d}y^2} + v \frac{\mathrm{d} \alpha}{\mathrm{d}y} = \alpha\beta
</math>
</center>
and
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} = -\alpha\beta
</math>
</center>
If we add these equations we obtain
<center>
<math>
\frac{\mathrm{d}^2 (\alpha+\beta)}{\mathrm{d}y^2} + v \frac{\mathrm{d} (\alpha+\beta)}{\mathrm{d}y} = 0
</math>
</center>
so that <math>\alpha + \beta = c_0 + c_1 e^{-vy}</math>. Boundary conditions
are that as <math>y\to\infty </math> <math>\alpha = 1</math> and
<math>\beta = 0</math> and if <math>y\to-\infty</math> then <math>\alpha = 0</math> and
<math>\beta = 1</math>. Therefore <math>\alpha + \beta = 1</math>.
This means that, since <math>\alpha \geq 0</math>, we must have
<math>0\leq \beta \leq 1</math>.
We can then obtain the following equation
<center>
<math>
\frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} + \beta(1-\beta)= 0
</math>
</center>
which we can write as the system of first order equations.

[[Image:R_d_phase_portrait.jpg|thumb|right|500px|Phase portrait for out system showing
the equilibrium points and the heteroclinic connection]]

We define the variable <math>\gamma = \frac{\mathrm{d}\beta}{\mathrm{d}y}</math>
and we obtain
<center><math>
\begin{align}
\frac{\mathrm{d}\beta}{\mathrm{d} y} &= \gamma&\\
\frac{\mathrm{d}\gamma}{\mathrm{d} y} &= -v\gamma + \beta(\beta -1)& \\
\end{align}
</math></center>
This dynamical system has equilibrium points at <math>(0,0)</math>
and <math>(1,0)</math>. We can analyse these equilibrium points by
linearization. The Jacobian matrix is
<center>
<math>
J =\begin{pmatrix}
0 & 1 \\
-1 + 2\beta & -v
\end{pmatrix}
</math>
</center>

We can easily see that the Jacobian evaluated at our first equilibrium point is
<center>
<math>
J_{(0,0)} =\begin{pmatrix}
0 & 1 \\
-1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\mu_{\pm} = -1/2 (v \mp \sqrt{v^2-4})</math>. Therefore
this point is a nodal sink (possibly a spiral)

[[Image:R_d_wave.gif|right|Travelling wave solution for v=2 and the position on
the heteroclinic connection]]

Additionally,
<center>
<math>
J_{(1,0)} =\begin{pmatrix}
0 & 1 \\
1 & -v
\end{pmatrix}
</math>
</center>
which has eigenvalues <math>\lambda_{\pm} = -1/2 (v \mp \sqrt{v^2+4})</math>.
This is a a saddle point. The unstable and stable
separatrices leave the equilibrium point at <math>(1,0)</math> in the directions
<math> \begin{pmatrix}\lambda_{\pm} \\ 1\end{pmatrix}</math>. The only path on which <math>\beta</math> is bounded
as <math>y\to-\infty</math> are the unstable separatrices. Also, only the
unstable separatrix which enters the region <math>\beta<1</math> is physically meaningful.

To find a travelling wave we need to find a heteroclinic connection
between the two equilibrium points which also has to satisfy the conditions
that <math>0\leq \beta \leq 1</math>.

We need to show that the heteroclinic connection does not cross the <math>\beta</math> axis.
Consider the region
<center>
<math>
R = \left\{(\beta,\gamma)\,|\, \beta<1,\,-k\beta<\gamma<0\right\}
</math>
</center>
On the line <math>\beta = 1,d\beta/dy<0</math> and hence all flow it into <math>R</math>.
On the line <math>\gamma = 0, d\gamma/dy < 0</math> for <math>0<\beta<1</math>. On the line
<math>\gamma = -k \beta</math> we know that <math>d\beta/dy < 0</math> so that integral paths
enter the region if and only if <math>d\gamma/d\beta < \gamma/\beta</math>. We know that
<center>
<math>
d\gamma/d\beta - \gamma/\beta = -v - \frac{\beta(1-\beta)}{\gamma} -\frac{\gamma}{\beta}
= \frac{1}{k} (k^2 - vk +1 -\beta)
</math>
</center>
when <math>\gamma = -k\beta</math>. Therefore we need to find a value of <math>k</math>
so that <math>k^2 - vk +1 < 0</math>, which is possible provided <math>v\geq 2</math>, for example
<math>k = \dfrac{v}{2}</math>.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|Yi8fytUszrQ}}

=== Part 2 ===

{{#ev:youtube|_668xmJ9DAQ}}

=== Part 3 ===

{{#ev:youtube|yTat11ERJMg}}

=== Part 4 ===

{{#ev:youtube|qG8bPJwX96k}}

=== Part 5 ===

{{#ev:youtube|38ZWW2dM1Qs}}

=== Part 6 ===

{{#ev:youtube|AZHlhA5pQBg}}

=== Part 7 ===

{{#ev:youtube|VhlY5lQcTlM}}

=== Part 8 ===

{{#ev:youtube|VhlY5lQcTlM}}

[[Category:Simple Nonlinear Waves]]

Example Calculations for the KdV and IST

2020-10-01T03:38:19Z

Meylan:

{{nonlinear waves course
| chapter title = Example Calculations for the KdV and IST
| next chapter = [[Reaction-Diffusion Systems]]
| previous chapter = [[Connection betwen KdV and the Schrodinger Equation]]
}}

We consider here the two examples we treated in [[Properties of the Linear Schrodinger Equation]].

==Example1: <math>\delta</math> function potential==

We have already calculated the scattering data for the delta function
potential in [[Properties of the Linear Schrodinger Equation]]. The scattering data is
<center><math>
S\left( \lambda,0\right) =\left( k_{1},\sqrt{k_{1}},\frac{u_{0}}{2ik-u_{0}
},\frac{2ik}{2ik-u_{0}}\right)
</math></center>
The spectral data evolves as
<center><math>
k_{1}=k_{1}
</math></center>
<center><math>
c_{1}\left( t\right) =c_{1}\left( 0\right) e^{4k_{1}^{3}t}=\sqrt{k_{1}
}e^{4k_{1}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
so that
<center><math>
S\left( \lambda,t\right) =\left( k_{1},\sqrt{k_{1}}e^{4k_{1}^{3}t}
,\frac{u_{0}}{2ik-u_{0}}e^{8ik^{3}t},\frac{2ik}{2ik-u_{0}}\right)
</math></center>

==Example 2: Hat Function Potential==

We solve for the case when
<center><math>
u\left( x\right) =\left\{
\begin{matrix}
0, & x\notin\left[ -1,1\right] \\
20, & x\in\left[ -1,1\right]
\end{matrix}
\right.
</math></center>
We have already solved this case in [[Properties of the Linear Schrodinger Equation]].
For the even solutions we need to solve
<center><math>
\tan\kappa=\frac{k}{\kappa}
</math></center>
where <math>\kappa=\sqrt{b-k^{2}}</math>.

For the odd solutions we need to solve and
<center><math>
\tan\kappa=-\frac{\kappa}{k}
</math></center>

Recall that the solitons have amplitude <math>2k_{n}^{2}</math> or <math>-2\lambda_{n}</math>. This
can be seen in the height of the solitary waves.

We cannot work with a hat function numerically, because the jump in <math>u</math> leads
to high frequencies which dominate the response.. We can smooth our function
by a number of methods. We use here the function <math>\tanh\left( x\right) </math> so
we write
<center><math>
u\left( x\right) =\frac{20}{2}\left( \tanh\left( \nu\left( x+1\right)
\right) -\tanh\left( \nu\left( x-1\right) \right) \right)
</math></center>
where <math>\nu</math> is an appropriate constant to make the function increase in value
sufficiently rapidly but not too rapidly.

{| class="wikitable"
|-
! Animation
! Three-dimensional plot.
|-
| [[Image:Wide_function.gif|thumb|right|500px|Evolution of <math>u(x,t)</math>.]]
| [[Image:Widefunction.jpg|thumb|right|500px|Evolution of <math>u(x,t)</math> ]]

|}

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|K0zeheguRKo}}

Example Calculations for the KdV and IST

2020-10-01T03:32:02Z

Meylan: /* Example 2: Hat Function Potential */

{{nonlinear waves course
| chapter title = Example Calculations for the KdV and IST
| next chapter = [[Reaction-Diffusion Systems]]
| previous chapter = [[Connection betwen KdV and the Schrodinger Equation]]
}}

We consider here the two examples we treated in [[Properties of the Linear Schrodinger Equation]].

==Example1: <math>\delta</math> function potential==

We have already calculated the scattering data for the delta function
potential in [[Properties of the Linear Schrodinger Equation]]. The scattering data is
<center><math>
S\left( \lambda,0\right) =\left( k_{1},\sqrt{k_{1}},\frac{u_{0}}{2ik-u_{0}
},\frac{2ik}{2ik-u_{0}}\right)
</math></center>
The spectral data evolves as
<center><math>
k_{1}=k_{1}
</math></center>
<center><math>
c_{1}\left( t\right) =c_{1}\left( 0\right) e^{4k_{1}^{3}t}=\sqrt{k_{1}
}e^{4k_{1}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
so that
<center><math>
S\left( \lambda,t\right) =\left( k_{1},\sqrt{k_{1}}e^{4k_{1}^{3}t}
,\frac{u_{0}}{2ik-u_{0}}e^{8ik^{3}t},\frac{2ik}{2ik-u_{0}}\right)
</math></center>

==Example 2: Hat Function Potential==

We solve for the case when
<center><math>
u\left( x\right) =\left\{
\begin{matrix}
0, & x\notin\left[ -1,1\right] \\
20, & x\in\left[ -1,1\right]
\end{matrix}
\right.
</math></center>
We have already solved this case in [[Properties of the Linear Schrodinger Equation]].
For the even solutions we need to solve
<center><math>
\tan\kappa=\frac{k}{\kappa}
</math></center>
where <math>\kappa=\sqrt{b-k^{2}}</math>.

For the odd solutions we need to solve and
<center><math>
\tan\kappa=-\frac{\kappa}{k}
</math></center>

Recall that the solitons have amplitude <math>2k_{n}^{2}</math> or <math>-2\lambda_{n}</math>. This
can be seen in the height of the solitary waves.

We cannot work with a hat function numerically, because the jump in <math>u</math> leads
to high frequencies which dominate the response.. We can smooth our function
by a number of methods. We use here the function <math>\tanh\left( x\right) </math> so
we write
<center><math>
u\left( x\right) =\frac{20}{2}\left( \tanh\left( \nu\left( x+1\right)
\right) -\tanh\left( \nu\left( x-1\right) \right) \right)
</math></center>
where <math>\nu</math> is an appropriate constant to make the function increase in value
sufficiently rapidly but not too rapidly.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|K0zeheguRKo}}

{| class="wikitable"
|-
! Animation
! Three-dimensional plot.
|-
| [[Image:Wide_function.gif|thumb|right|500px|Evolution of <math>u(x,t)</math>.]]
| [[Image:Widefunction.jpg|thumb|right|500px|Evolution of <math>u(x,t)</math> ]]

|}

Connection betwen KdV and the Schrodinger Equation

2020-10-01T03:28:30Z

Meylan: /* Single Soliton Example */

{{nonlinear waves course
| chapter title = Connection betwen KdV and the Schrodinger Equation
| next chapter = [[Example Calculations for the KdV and IST]]
| previous chapter = [[Properties of the Linear Schrodinger Equation]]
}}

If we substitute the relationship
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
into the KdV after some manipulation we obtain
<center><math>
\partial_{t}\lambda w^{2}+\partial_{x}\left( w\partial_{x}Q-\partial
_{x}wQ\right) =0
</math></center>
where <math>Q=\partial_{t}w+\partial_{x}^{3}w-3\left( \lambda-u\right)
\partial_{x}w.</math> If we integrate this equation then we obtain the result that
<center><math>
\partial_{t}\lambda=0
</math></center>
provided that the eigenfunction <math>w</math> is bounded (which is true for the bound
state eigenfunctions). This shows that the discrete eigenvalues are unchanged
and <math>u\left( x,t\right) </math> evolves according to the KdV. Many other
properties can be found

==Scattering Data==

For the discrete spectrum the eigenfunctions behave like
<center><math>
w_{n}\left( x,t\right) =c_{n}\left( t\right) e^{-k_{n}x}
</math></center>
as <math>x\rightarrow\infty</math> with
<center><math>
\int_{-\infty}^{\infty}\left( w_{n}\left( x,t\right) \right) ^{2}\mathrm{d}x=1
</math></center>
The continuous spectrum looks like
<center><math>
w\left( x,t\right) \approx \mathrm{e}^{-\mathrm{i}kx}+r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}
,\ \ \ x\rightarrow-\infty
</math></center>
<center><math>
w\left( x,t\right) \approx a\left( k,t\right) \mathrm{e}^{-\mathrm{i}kx},\ \ \ x\rightarrow
\infty
</math></center>
where <math>r</math> is the reflection coefficient and <math>a</math> is the transmission
coefficient. This gives us the scattering data at <math>t=0</math>
<center><math>
S\left( \lambda,0\right) =\left( \left\{ k_{n},c_{n}\left( 0\right)
\right\} _{n=1}^{N},r\left( k,0\right) ,a\left( k,0\right) \right)
</math></center>
The scattering data evolves as
<center><math>
k_{n}=k_{n}
</math></center>
<center><math>
c_{n}\left( t\right) =c_{n}\left( 0\right) e^{4k_{n}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
We can recover <math>u</math> from scattering data. We write
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}
x}+\frac{1}{2\pi}\int_{-\infty}^{\infty}r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>
Then solve
<center><math>
K\left( x,y;t\right) +F\left( x+y;t\right) +\int_{x}^{\infty}K\left(
x,z;t\right) F\left( z+y;t\right) \mathrm{d}z=0
</math></center>
This is a linear integral equation called the Gelfand-Levitan-Marchenko equation. We then find <math>u</math> from
<center><math>
u\left( x,t\right) =2\partial_{x}K\left( x,x,t\right)
</math></center>

==Reflectionless Potential==

In general the IST is difficult to solve. However, there is a simplification
we can make when we have a reflectionless potential (which we will see gives
rise to the soliton solutions). The reflectionless potential is the case when
<math>r\left( k,0\right) =0</math> for all values of <math>k</math> for some <math>u.</math> In this case
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}x}
</math></center>
then
<center><math>
K\left( x,y,t\right) +\sum_{n=1}^{N}c_{n}^{2}\left( t\right)
e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}K\left( x,z,t\right)
\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
From the equation we can see that
<center><math>
K\left( x,y,t\right) =-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}
</math></center>
If we substitute this into the equation,
<center><math>
-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}+\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}
-\sum_{m=1}^{N}v_{m}\left( x,t\right) e^{-k_{m}z}\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
which leads to
<center><math>
-\sum_{n=1}^{N} v_{n}\left( x,t\right) e^{-k_{n}y}
+\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }
-\sum_{n=1}^{N}\sum_{m=1}^{N}\frac{ c_{n}^{2}\left(
t\right) }{k_{n}+k_{m}}v_{m}\left( x\right) e^{-k_{m}x}e^{-k_{n}\left(
y+x\right) }=0
</math></center>
and we can eliminate the sum over <math>n</math> and the
<math>e^{-k_{n}y}</math> to obtain
<center><math>
-v_{n}\left( x,t\right) +c_{n}^2\left( t\right) e^{-k_{n}x}-\sum_{m=1}
^{N}\frac{c_{n}^2\left( t\right) }{k_{n}+k_{m}}
v_{m}\left( x,t\right) e^{-\left( k_{m}+k_{n}\right) x}=0
</math></center>
which is an algebraic (finite dimensional system) for the unknows <math>v_{n}.</math>.

We can write this as
<center><math>
\left( \mathbf{I}+\mathbf{C}\right) \vec{v}=\vec{f}
</math></center>
where <math>f_{m}=c_{m}^2\left( t\right) e^{-k_{m}x}</math> and
the elements of <math>\mathbf{C}</math> are given by
<center><math>
c_{mn} = \frac{c_{n}^2\left( t\right)}
{k_{n}+k_{m}}e^{-\left( k_{m}+k_{n}\right) x}
</math></center>
This gives us
<center><math>
K\left( x,y,t\right) =-\sum_{m=1}^{N} \left(
\mathbf{I}+\mathbf{C}\right) ^{-1}\vec{f}e^{-k_{m}y}
</math></center>
We then find <math>u(x,t)</math> from <math>K</math>.

=== Single Soliton Example ===

If <math>n=1</math> (a single soliton
solution) we get
<center><math>\begin{matrix}
K\left( x,x,t\right) & =-\frac{c_{1}^2\left( t\right)
e^{-2k_{1}x}}{1+\frac{c_{1}^2\left( t\right) }
{2k_{1}}e^{-2 k_{1} x}}\\
& =\frac{-1}{e^{2k_{1}x - 8k_{1}^{3}t-\alpha} + 1/2k_1}
\end{matrix}</math></center>
where <math>e^{-\alpha}=1/c_{0}^{2}\left( 0\right) .</math> Therefore
<center><math>\begin{align}
u\left( x,t\right) & =2\partial_{x}K\left( x,x,t\right) \\
& =\frac{4k_{1}e^{2k_{1}x-8k_{1}^{3}t-\alpha}}{\left( e^{2k_{1}
x-8k_{1}^{3}t-\alpha} + 1/2k_1\right) ^{2}}\\
& =\frac{8k_{1}^{2}}{\left( \sqrt{2k_{1}}e^{\theta}+e^{-\theta}/\sqrt
{2k_{1}}\right) ^{2}}\\
& =2k_{1}^{2}\mbox{sech}^{2}\left\{ k_{1}\left( x-x_{0}\right) -4k_{1}
^{3}t\right\}
\end{align}</math></center>
where <math>\theta=k_{1}x-4k^{3}t-\alpha/2</math> and <math>\sqrt{2k}e^{-\alpha/2}=e^{-kx_{0}
}</math>. This is of course the single soliton solution.

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|KHGoPCoyP28}}

=== Part 2 ===

{{#ev:youtube|8mVq6MQWO3I}}

=== Part 3 ===

{{#ev:youtube|meTgaaGKsfQ}}

=== Part 4 ===

{{#ev:youtube|AfAQyjzvJUU}}

Connection betwen KdV and the Schrodinger Equation

2020-10-01T03:02:50Z

Meylan: /* Single Soliton Example */

{{nonlinear waves course
| chapter title = Connection betwen KdV and the Schrodinger Equation
| next chapter = [[Example Calculations for the KdV and IST]]
| previous chapter = [[Properties of the Linear Schrodinger Equation]]
}}

If we substitute the relationship
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
into the KdV after some manipulation we obtain
<center><math>
\partial_{t}\lambda w^{2}+\partial_{x}\left( w\partial_{x}Q-\partial
_{x}wQ\right) =0
</math></center>
where <math>Q=\partial_{t}w+\partial_{x}^{3}w-3\left( \lambda-u\right)
\partial_{x}w.</math> If we integrate this equation then we obtain the result that
<center><math>
\partial_{t}\lambda=0
</math></center>
provided that the eigenfunction <math>w</math> is bounded (which is true for the bound
state eigenfunctions). This shows that the discrete eigenvalues are unchanged
and <math>u\left( x,t\right) </math> evolves according to the KdV. Many other
properties can be found

==Scattering Data==

For the discrete spectrum the eigenfunctions behave like
<center><math>
w_{n}\left( x,t\right) =c_{n}\left( t\right) e^{-k_{n}x}
</math></center>
as <math>x\rightarrow\infty</math> with
<center><math>
\int_{-\infty}^{\infty}\left( w_{n}\left( x,t\right) \right) ^{2}\mathrm{d}x=1
</math></center>
The continuous spectrum looks like
<center><math>
w\left( x,t\right) \approx \mathrm{e}^{-\mathrm{i}kx}+r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}
,\ \ \ x\rightarrow-\infty
</math></center>
<center><math>
w\left( x,t\right) \approx a\left( k,t\right) \mathrm{e}^{-\mathrm{i}kx},\ \ \ x\rightarrow
\infty
</math></center>
where <math>r</math> is the reflection coefficient and <math>a</math> is the transmission
coefficient. This gives us the scattering data at <math>t=0</math>
<center><math>
S\left( \lambda,0\right) =\left( \left\{ k_{n},c_{n}\left( 0\right)
\right\} _{n=1}^{N},r\left( k,0\right) ,a\left( k,0\right) \right)
</math></center>
The scattering data evolves as
<center><math>
k_{n}=k_{n}
</math></center>
<center><math>
c_{n}\left( t\right) =c_{n}\left( 0\right) e^{4k_{n}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
We can recover <math>u</math> from scattering data. We write
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}
x}+\frac{1}{2\pi}\int_{-\infty}^{\infty}r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>
Then solve
<center><math>
K\left( x,y;t\right) +F\left( x+y;t\right) +\int_{x}^{\infty}K\left(
x,z;t\right) F\left( z+y;t\right) \mathrm{d}z=0
</math></center>
This is a linear integral equation called the Gelfand-Levitan-Marchenko equation. We then find <math>u</math> from
<center><math>
u\left( x,t\right) =2\partial_{x}K\left( x,x,t\right)
</math></center>

==Reflectionless Potential==

In general the IST is difficult to solve. However, there is a simplification
we can make when we have a reflectionless potential (which we will see gives
rise to the soliton solutions). The reflectionless potential is the case when
<math>r\left( k,0\right) =0</math> for all values of <math>k</math> for some <math>u.</math> In this case
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}x}
</math></center>
then
<center><math>
K\left( x,y,t\right) +\sum_{n=1}^{N}c_{n}^{2}\left( t\right)
e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}K\left( x,z,t\right)
\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
From the equation we can see that
<center><math>
K\left( x,y,t\right) =-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}
</math></center>
If we substitute this into the equation,
<center><math>
-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}+\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}
-\sum_{m=1}^{N}v_{m}\left( x,t\right) e^{-k_{m}z}\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
which leads to
<center><math>
-\sum_{n=1}^{N} v_{n}\left( x,t\right) e^{-k_{n}y}
+\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }
-\sum_{n=1}^{N}\sum_{m=1}^{N}\frac{ c_{n}^{2}\left(
t\right) }{k_{n}+k_{m}}v_{m}\left( x\right) e^{-k_{m}x}e^{-k_{n}\left(
y+x\right) }=0
</math></center>
and we can eliminate the sum over <math>n</math> and the
<math>e^{-k_{n}y}</math> to obtain
<center><math>
-v_{n}\left( x,t\right) +c_{n}^2\left( t\right) e^{-k_{n}x}-\sum_{m=1}
^{N}\frac{c_{n}^2\left( t\right) }{k_{n}+k_{m}}
v_{m}\left( x,t\right) e^{-\left( k_{m}+k_{n}\right) x}=0
</math></center>
which is an algebraic (finite dimensional system) for the unknows <math>v_{n}.</math>.

We can write this as
<center><math>
\left( \mathbf{I}+\mathbf{C}\right) \vec{v}=\vec{f}
</math></center>
where <math>f_{m}=c_{m}^2\left( t\right) e^{-k_{m}x}</math> and
the elements of <math>\mathbf{C}</math> are given by
<center><math>
c_{mn} = \frac{c_{n}^2\left( t\right)}
{k_{n}+k_{m}}e^{-\left( k_{m}+k_{n}\right) x}
</math></center>
This gives us
<center><math>
K\left( x,y,t\right) =-\sum_{m=1}^{N} \left(
\mathbf{I}+\mathbf{C}\right) ^{-1}\vec{f}e^{-k_{m}y}
</math></center>
We then find <math>u(x,t)</math> from <math>K</math>.

=== Single Soliton Example ===

If <math>n=1</math> (a single soliton
solution) we get
<center><math>\begin{matrix}
K\left( x,x,t\right) & =-\frac{c_{1}^2\left( t\right)
e^{-2k_{1}x}}{1+\frac{c_{1}^2\left( t\right) }
{2k_{1}}e^{-2 k_{1} x}}\\
& =\frac{-1}{e^{2k_{1}x - 8k_{1}^{3}t-\alpha} + 1/2k_1}
\end{matrix}</math></center>
where <math>e^{-\alpha}=1/c_{0}^{2}\left( 0\right) .</math> Therefore
<center><math>\begin{align}
u\left( x,t\right) & =2\partial_{x}K\left( x,x,t\right) \\
& =\frac{4k_{1}e^{2k_{1}x-8k_{1}^{3}t-\alpha}}{\left( e^{2k_{1}
x-8k_{1}^{3}t-\alpha} + 1/2k_1\right) ^{2}}\\
& =\frac{8k_{1}^{2}}{\left( \sqrt{2k_{1}}e^{\theta}+e^{-\theta}/\sqrt
{2k_{1}}\right) ^{2}}\\
& =2k_{1}^{2}\mbox{sech}^{2}\left\{ k_{1}\left( x-x_{0}\right) -4k_{1}
^{3}t\right\}
\end{align}</math></center>
where <math>\theta=k_{1}x-4k^{3}t-\alpha/2</math> and <math>\sqrt{2k}e^{-\alpha/2}=e^{-kx_{0}
}</math>. This is of course the single soliton solution.

Connection betwen KdV and the Schrodinger Equation

2020-10-01T01:46:16Z

Meylan: /* Scattering Data */

{{nonlinear waves course
| chapter title = Connection betwen KdV and the Schrodinger Equation
| next chapter = [[Example Calculations for the KdV and IST]]
| previous chapter = [[Properties of the Linear Schrodinger Equation]]
}}

If we substitute the relationship
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
into the KdV after some manipulation we obtain
<center><math>
\partial_{t}\lambda w^{2}+\partial_{x}\left( w\partial_{x}Q-\partial
_{x}wQ\right) =0
</math></center>
where <math>Q=\partial_{t}w+\partial_{x}^{3}w-3\left( \lambda-u\right)
\partial_{x}w.</math> If we integrate this equation then we obtain the result that
<center><math>
\partial_{t}\lambda=0
</math></center>
provided that the eigenfunction <math>w</math> is bounded (which is true for the bound
state eigenfunctions). This shows that the discrete eigenvalues are unchanged
and <math>u\left( x,t\right) </math> evolves according to the KdV. Many other
properties can be found

==Scattering Data==

For the discrete spectrum the eigenfunctions behave like
<center><math>
w_{n}\left( x,t\right) =c_{n}\left( t\right) e^{-k_{n}x}
</math></center>
as <math>x\rightarrow\infty</math> with
<center><math>
\int_{-\infty}^{\infty}\left( w_{n}\left( x,t\right) \right) ^{2}\mathrm{d}x=1
</math></center>
The continuous spectrum looks like
<center><math>
w\left( x,t\right) \approx \mathrm{e}^{-\mathrm{i}kx}+r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}
,\ \ \ x\rightarrow-\infty
</math></center>
<center><math>
w\left( x,t\right) \approx a\left( k,t\right) \mathrm{e}^{-\mathrm{i}kx},\ \ \ x\rightarrow
\infty
</math></center>
where <math>r</math> is the reflection coefficient and <math>a</math> is the transmission
coefficient. This gives us the scattering data at <math>t=0</math>
<center><math>
S\left( \lambda,0\right) =\left( \left\{ k_{n},c_{n}\left( 0\right)
\right\} _{n=1}^{N},r\left( k,0\right) ,a\left( k,0\right) \right)
</math></center>
The scattering data evolves as
<center><math>
k_{n}=k_{n}
</math></center>
<center><math>
c_{n}\left( t\right) =c_{n}\left( 0\right) e^{4k_{n}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
We can recover <math>u</math> from scattering data. We write
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}
x}+\frac{1}{2\pi}\int_{-\infty}^{\infty}r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>
Then solve
<center><math>
K\left( x,y;t\right) +F\left( x+y;t\right) +\int_{x}^{\infty}K\left(
x,z;t\right) F\left( z+y;t\right) \mathrm{d}z=0
</math></center>
This is a linear integral equation called the Gelfand-Levitan-Marchenko equation. We then find <math>u</math> from
<center><math>
u\left( x,t\right) =2\partial_{x}K\left( x,x,t\right)
</math></center>

==Reflectionless Potential==

In general the IST is difficult to solve. However, there is a simplification
we can make when we have a reflectionless potential (which we will see gives
rise to the soliton solutions). The reflectionless potential is the case when
<math>r\left( k,0\right) =0</math> for all values of <math>k</math> for some <math>u.</math> In this case
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}x}
</math></center>
then
<center><math>
K\left( x,y,t\right) +\sum_{n=1}^{N}c_{n}^{2}\left( t\right)
e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}K\left( x,z,t\right)
\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
From the equation we can see that
<center><math>
K\left( x,y,t\right) =-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}
</math></center>
If we substitute this into the equation,
<center><math>
-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}+\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}
-\sum_{m=1}^{N}v_{m}\left( x,t\right) e^{-k_{m}z}\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
which leads to
<center><math>
-\sum_{n=1}^{N} v_{n}\left( x,t\right) e^{-k_{n}y}
+\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }
-\sum_{n=1}^{N}\sum_{m=1}^{N}\frac{ c_{n}^{2}\left(
t\right) }{k_{n}+k_{m}}v_{m}\left( x\right) e^{-k_{m}x}e^{-k_{n}\left(
y+x\right) }=0
</math></center>
and we can eliminate the sum over <math>n</math> and the
<math>e^{-k_{n}y}</math> to obtain
<center><math>
-v_{n}\left( x,t\right) +c_{n}^2\left( t\right) e^{-k_{n}x}-\sum_{m=1}
^{N}\frac{c_{n}^2\left( t\right) }{k_{n}+k_{m}}
v_{m}\left( x,t\right) e^{-\left( k_{m}+k_{n}\right) x}=0
</math></center>
which is an algebraic (finite dimensional system) for the unknows <math>v_{n}.</math>.

We can write this as
<center><math>
\left( \mathbf{I}+\mathbf{C}\right) \vec{v}=\vec{f}
</math></center>
where <math>f_{m}=c_{m}^2\left( t\right) e^{-k_{m}x}</math> and
the elements of <math>\mathbf{C}</math> are given by
<center><math>
c_{mn} = \frac{c_{n}^2\left( t\right)}
{k_{n}+k_{m}}e^{-\left( k_{m}+k_{n}\right) x}
</math></center>
This gives us
<center><math>
K\left( x,y,t\right) =-\sum_{m=1}^{N} \left(
\mathbf{I}+\mathbf{C}\right) ^{-1}\vec{f}e^{-k_{m}y}
</math></center>
We then find <math>u(x,t)</math> from <math>K</math>.

=== Single Soliton Example ===

If <math>n=1</math> (a single soliton
solution) we get
<center><math>\begin{matrix}
K\left( x,x,t\right) & =-\frac{c_{1}^2\left( t\right)
e^{-2k_{1}x}}{1+\frac{c_{1}^2\left( t\right) }
{2k_{1}}e^{-2 k_{1} x}}\\
& =\frac{-1}{e^{-2k_{1}x + 8k_{1}^{3}t-\alpha} + 1/2k_1}
\end{matrix}</math></center>
where <math>e^{-\alpha}=1/c_{0}^{2}\left( 0\right) .</math> Therefore
<center><math>\begin{align}
u\left( x,t\right) & =2\partial_{x}K\left( x,x,t\right) \\
& =\frac{4k_{1}e^{2k_{1}x-8k_{1}^{3}t-\alpha}}{\left( e^{2k_{1}
x-8k_{1}^{3}t-\alpha} + 1/2k_1\right) ^{2}}\\
& =\frac{8k_{1}^{2}}{\left( \sqrt{2k_{1}}e^{\theta}+e^{-\theta}/\sqrt
{2k_{1}}\right) ^{2}}\\
& =2k_{1}^{2}\mbox{sech}^{2}\left\{ k_{1}\left( x-x_{0}\right) -4k_{1}
^{3}t\right\}
\end{align}</math></center>
where <math>\theta=k_{1}x-4k^{3}t-\alpha/2</math> and <math>\sqrt{2k}e^{-\alpha/2}=e^{-kx_{0}
}</math>. This is of course the single soliton solution.

Properties of the Linear Schrodinger Equation

2020-09-24T01:03:35Z

Meylan: /* Case when \lambda>0 */

{{nonlinear waves course
| chapter title = Properties of the Linear Schrodinger Equation
| next chapter = [[Connection betwen KdV and the Schrodinger Equation]]
| previous chapter = [[Introduction to the Inverse Scattering Transform]]
}}

The linear Schrodinger equation
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
has two kinds of solutions for <math>u\rightarrow0</math> as <math>x\rightarrow\pm\infty.</math> The
first are waves and the second are bound solutions. It is well known that
there are at most a finite number of bound solutions (provided <math>u\rightarrow0</math>
as <math>x\pm\infty</math> sufficiently rapidly) and a continum of solutions for the
incident waves. This is easiest seen through the following examples

==Example 1: <math>\delta</math> function potential==

We consider here the case when <math>u\left( x,0\right) = u_0 \delta\left( x\right)
.</math> Note that this function can be thought of as the limit as of the potential
<center><math>
u\left( x\right) =\left\{
\begin{matrix}

0 & x\notin\left[ -\varepsilon,\varepsilon\right] \\
\frac{u_{0}}{2\varepsilon} & x\in\left[ -\varepsilon,\varepsilon\right]
\end{matrix}
\right.
</math></center>
In this case we need to solve
<center><math>
\partial_{x}^{2}w+ u_0\delta(x) w=-\lambda w
</math></center>
We consider the case of <math>\lambda<0</math> and <math>\lambda>0</math> separately. For the first
case we write <math>\lambda=-k^{2}</math> and we obtain
<center><math>
w\left( x\right) =\left\{
\begin{matrix}

ae^{kx}, & x<0\\
be^{-kx}, & x>0
\end{matrix}
\right.
</math></center>
We have two conditions at <math>x=0,</math> <math>w</math> must be continuous at <math>0</math> and
<math>\partial_{x}w\left( 0^{+}\right) -\partial_{x}w\left( 0^{-}\right)
+u_0 w\left( 0\right) =0.</math> This final condition is obtained by integrating `across' zero as follows
<center><math>\begin{align}
\int_{0^{-}}^{0^{+}} \left(\partial_x^2 w +u_0\delta(x) w + \lambda w \right) \ \mathrm{d}x = 0.
\end{align}
</math></center>

This gives the condition that <math>a=b</math> and
<math>k=u_{0}/2.</math> We need to normalise the eigenfunctions so that
<center><math>
\int_{-\infty}^{\infty}\left( w\left( x\right) \right) ^{2}\mathrm{d}x=1.
</math></center>
Therefore
<center><math>
2\int_{0}^{\infty}\left( ae^{-u_{0}x/2}\right) ^{2}\mathrm{d}x=1
</math></center>
which means that <math>a=\sqrt{u_{0}/2}.</math> Therefore, there is only one discrete
spectral point which we denote by <math>k_{1}=u_{0}/2</math>
<center><math>
w_{1}\left( x\right) =\left\{
\begin{matrix}

\sqrt{k_{1}}e^{k_{1}x}, & x<0\\
\sqrt{k_{1}}e^{-k_{1}x}, & x>0
\end{matrix}
\right.
</math></center>
The continuous eigenfunctions correspond to <math>\lambda=k^{2}>0</math> are of the form
<center><math>
w\left( x\right) =\left\{
\begin{matrix}

\mathrm{e}^{-\mathrm{i}kx}+r\mathrm{e}^{\mathrm{i}kx}, & x<0\\
t\mathrm{e}^{-\mathrm{i}kx}, & x>0
\end{matrix}
\right.
</math></center>
Again we have the conditions that <math>w</math> must be continuous at <math>0</math> and
<math>\partial_{x}w\left( 0^{+}\right) -\partial_{x}w\left( 0^{-}\right)
+u_{0}w\left( 0\right) =0.</math> This gives us
<center><math>\begin{matrix}
1+r & =t\\
-ikt+ik-ikr & =-tu_{0}
\end{matrix}</math></center>
which has solution
<center><math>\begin{matrix}
r & =\frac{u_{0}}{2ik-u_{0}}\\
t & =\frac{2ik}{2ik-u_{0}}
\end{matrix}</math></center>

==Example 2: Hat Function Potential==

The properties of the eigenfunction is perhaps seem most easily through the
following example
<center><math>
u\left( x\right) =\left\{
\begin{matrix}

0 & x\notin\left[ -\varsigma,\varsigma\right] \\
b & x\in\left[ -\varsigma,\varsigma\right]
\end{matrix}
\right.
</math></center>
where <math>b>0.</math>

===Case when <math>\lambda<0</math>===

If we solve this equation for the case when <math>\lambda<0,</math> <math>\lambda=-k^{2}</math> we
get
<center><math>
w\left( x\right) =\left\{
\begin{matrix}
a_{1}e^{kx}, & x<-\varsigma,\\
b_{1}\cos\kappa x+b_{2}\sin\kappa x, & -\varsigma< x <\varsigma,\\
a_{2}e^{-kx}, & x>\varsigma,
\end{matrix}
\right.
</math></center>
where <math>\kappa=\sqrt{b-k^{2}}</math> which means that <math>0\leq k\leq\sqrt{b}</math> (there is
no solution for <math>k>\sqrt{b}).</math> We then match <math>w</math> and its derivative at
<math>x=\pm\varsigma</math> to solve for <math>a</math> and <math>b</math>. This leads to two system of
equations, one for the even (<math>a_{1}=a_{2}</math> and <math>b_{2}=0</math> ) and one for the odd
solutions (<math>a_{1}=-a_{2}</math> and <math>b_{1}=0)</math>. The solution for the even solutions
is
<center><math>
w\left( x\right) =\left\{
\begin{matrix}
a_{1}e^{kx}, & x<-\varsigma,\\
b_{1}\cos\kappa x, & -\varsigma< x <\varsigma,\\
a_{1}e^{-kx}, & x>\varsigma,
\end{matrix}
\right.
</math></center>
If we impose the condition that the function and its derivative are continuous at
<math>x=\pm\varsigma</math> we obtain the following equation
<center><math>
\left(
\begin{matrix}

e^{-k\varsigma} & -\cos\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\sin\kappa\varsigma
\end{matrix}
\right) \left(
\begin{matrix}

a_{1}\\
b_{1}
\end{matrix}
\right) =\left(
\begin{matrix}

0\\
0
\end{matrix}
\right)
</math></center>
This has non trivial solutions when
<center><math>
\det\left(
\begin{matrix}

e^{-k\varsigma} & -\cos\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\sin\kappa\varsigma
\end{matrix}
\right) =0
</math></center>
which gives us the equation
<center><math>
-\kappa e^{-k\varsigma}\sin\kappa\varsigma+k\cos\kappa\varsigma
e^{-k\varsigma}=0
</math></center>
or
<center><math>
\tan\kappa\varsigma=\frac{k}{\kappa}
</math></center>
We know that <math>0<\kappa<\sqrt{b}</math> and if we plot this we see that we obtain a
finite number of solutions.

The solution for the odd solutions is
<center><math>
w\left( x\right) =\left\{
\begin{matrix}
a_{1}e^{kx}, & x <-\varsigma,\\
b_{2}\sin\kappa x, & -\varsigma< x <\varsigma,\\
-a_{1}e^{-kx} & x > \varsigma,
\end{matrix}
\right.
</math></center>
and again imposing the condition that the solution and its derivative is continuous
at <math>x=\pm\varsigma</math> gives
<center><math>
\left(
\begin{matrix}

e^{-k\varsigma} & \sin\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\cos\kappa\varsigma
\end{matrix}
\right) \left(
\begin{matrix}

a_{1}\\
b_{1}
\end{matrix}
\right) =\left(
\begin{matrix}

0\\
0
\end{matrix}
\right)
</math></center>
This can non trivial solutions when
<center><math>
\det\left(
\begin{matrix}

e^{-k\varsigma} & \sin\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\cos\kappa\varsigma
\end{matrix}
\right) =0
</math></center>
which gives us the equation
<center><math>
\kappa e^{-k\varsigma}a\cos\kappa\varsigma+k\sin\kappa\varsigma e^{-k\varsigma}=0
</math></center>
or
<center><math>
\tan\varsigma\kappa=-\frac{\kappa}{k}
</math></center>

=== Case when <math>\lambda>0</math> ===

When <math>\lambda>0</math> we write <math>\lambda=k^{2}</math> and we obtain solution
<center><math>
w\left( x\right) =\left\{
\begin{matrix}

\mathrm{e}^{-\mathrm{i}kx}+r\mathrm{e}^{\mathrm{i}kx}, & x <-\varsigma\\
b_{1}\cos\kappa x+b_{2}\sin\kappa x & -\varsigma< x <\varsigma\\
t\mathrm{e}^{-\mathrm{i}kx} & x>\varsigma
\end{matrix}
\right.
</math></center>
where <math>\kappa=\sqrt{b+k^{2}}.</math> Matching <math>w</math> and its derivaties at <math>x=\pm1</math> we
obtain
<center><math>
\left(
\begin{matrix}

-\mathrm{e}^{-\mathrm{i}k\varsigma} & \cos\kappa\varsigma & -\sin\kappa\varsigma & 0\\
ik\mathrm{e}^{-\mathrm{i}k\varsigma} & \kappa\sin\kappa\varsigma & \kappa\cos\kappa
\varsigma & 0\\
0 & \cos\kappa\varsigma & \sin\kappa\varsigma & -\mathrm{e}^{-\mathrm{i}k\varsigma}\\
0 & -\kappa\sin\kappa\varsigma & \kappa\cos\kappa\varsigma &
ik\mathrm{e}^{-\mathrm{i}k\varsigma}
\end{matrix}
\right) \left(
\begin{matrix}

r\\
b_{1}\\
b_{2}\\
t
\end{matrix}
\right) =\left(
\begin{matrix}

\mathrm{e}^{\mathrm{i}k}\\
ik\mathrm{e}^{-\mathrm{i}k}\\
0\\
0
\end{matrix}
\right)
</math></center>

== Lecture Videos ==

=== Part 1 ===

{{#ev:youtube|anAThvCcpNw}}

=== Part 2 ===

{{#ev:youtube|SDPIx42VjLQ}}

=== Part 3 ===

{{#ev:youtube|OUmjeLZWr3M}}

=== Part 4 ===

{{#ev:youtube|hIfcO3a8_XU}}

=== Part 5 ===

{{#ev:youtube|z13lKSTficA}}

=== Part 6 ===

{{#ev:youtube|2XlQpEscxE4}}

=== Part 7 ===

{{#ev:youtube|iMMQ4NUdXNc}}

=== Part 8 ===

{{#ev:youtube|0F_dINNxMlw}}

Properties of the Linear Schrodinger Equation

2020-09-21T02:25:36Z

Meylan: /* Case when \lambda>0 */

Properties of the Linear Schrodinger Equation

2020-09-21T02:24:08Z

Meylan: /* Example 1: \delta function potential */

{{nonlinear waves course
| chapter title = Properties of the Linear Schrodinger Equation
| next chapter = [[Connection betwen KdV and the Schrodinger Equation]]
| previous chapter = [[Introduction to the Inverse Scattering Transform]]
}}

The linear Schrodinger equation
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
has two kinds of solutions for <math>u\rightarrow0</math> as <math>x\rightarrow\pm\infty.</math> The
first are waves and the second are bound solutions. It is well known that
there are at most a finite number of bound solutions (provided <math>u\rightarrow0</math>
as <math>x\pm\infty</math> sufficiently rapidly) and a continum of solutions for the
incident waves. This is easiest seen through the following examples

==Example 1: <math>\delta</math> function potential==

We consider here the case when <math>u\left( x,0\right) = u_0 \delta\left( x\right)
.</math> Note that this function can be thought of as the limit as of the potential
<center><math>
u\left( x\right) =\left\{
\begin{matrix}

0 & x\notin\left[ -\varepsilon,\varepsilon\right] \\
\frac{u_{0}}{2\varepsilon} & x\in\left[ -\varepsilon,\varepsilon\right]
\end{matrix}
\right.
</math></center>
In this case we need to solve
<center><math>
\partial_{x}^{2}w+ u_0\delta(x) w=-\lambda w
</math></center>
We consider the case of <math>\lambda<0</math> and <math>\lambda>0</math> separately. For the first
case we write <math>\lambda=-k^{2}</math> and we obtain
<center><math>
w\left( x\right) =\left\{
\begin{matrix}

ae^{kx}, & x<0\\
be^{-kx}, & x>0
\end{matrix}
\right.
</math></center>
We have two conditions at <math>x=0,</math> <math>w</math> must be continuous at <math>0</math> and
<math>\partial_{x}w\left( 0^{+}\right) -\partial_{x}w\left( 0^{-}\right)
+u_0 w\left( 0\right) =0.</math> This final condition is obtained by integrating `across' zero as follows
<center><math>\begin{align}
\int_{0^{-}}^{0^{+}} \left(\partial_x^2 w +u_0\delta(x) w + \lambda w \right) \ \mathrm{d}x = 0.
\end{align}
</math></center>

This gives the condition that <math>a=b</math> and
<math>k=u_{0}/2.</math> We need to normalise the eigenfunctions so that
<center><math>
\int_{-\infty}^{\infty}\left( w\left( x\right) \right) ^{2}\mathrm{d}x=1.
</math></center>
Therefore
<center><math>
2\int_{0}^{\infty}\left( ae^{-u_{0}x/2}\right) ^{2}\mathrm{d}x=1
</math></center>
which means that <math>a=\sqrt{u_{0}/2}.</math> Therefore, there is only one discrete
spectral point which we denote by <math>k_{1}=u_{0}/2</math>
<center><math>
w_{1}\left( x\right) =\left\{
\begin{matrix}

\sqrt{k_{1}}e^{k_{1}x}, & x<0\\
\sqrt{k_{1}}e^{-k_{1}x}, & x>0
\end{matrix}
\right.
</math></center>
The continuous eigenfunctions correspond to <math>\lambda=k^{2}>0</math> are of the form
<center><math>
w\left( x\right) =\left\{
\begin{matrix}

\mathrm{e}^{-\mathrm{i}kx}+r\mathrm{e}^{\mathrm{i}kx}, & x<0\\
t\mathrm{e}^{-\mathrm{i}kx}, & x>0
\end{matrix}
\right.
</math></center>
Again we have the conditions that <math>w</math> must be continuous at <math>0</math> and
<math>\partial_{x}w\left( 0^{+}\right) -\partial_{x}w\left( 0^{-}\right)
+u_{0}w\left( 0\right) =0.</math> This gives us
<center><math>\begin{matrix}
1+r & =t\\
-ikt+ik-ikr & =-tu_{0}
\end{matrix}</math></center>
which has solution
<center><math>\begin{matrix}
r & =\frac{u_{0}}{2ik-u_{0}}\\
t & =\frac{2ik}{2ik-u_{0}}
\end{matrix}</math></center>

==Example 2: Hat Function Potential==

The properties of the eigenfunction is perhaps seem most easily through the
following example
<center><math>
u\left( x\right) =\left\{
\begin{matrix}

0 & x\notin\left[ -\varsigma,\varsigma\right] \\
b & x\in\left[ -\varsigma,\varsigma\right]
\end{matrix}
\right.
</math></center>
where <math>b>0.</math>

===Case when <math>\lambda<0</math>===

If we solve this equation for the case when <math>\lambda<0,</math> <math>\lambda=-k^{2}</math> we
get
<center><math>
w\left( x\right) =\left\{
\begin{matrix}
a_{1}e^{kx}, & x<-\varsigma,\\
b_{1}\cos\kappa x+b_{2}\sin\kappa x, & -\varsigma< x <\varsigma,\\
a_{2}e^{-kx}, & x>\varsigma,
\end{matrix}
\right.
</math></center>
where <math>\kappa=\sqrt{b-k^{2}}</math> which means that <math>0\leq k\leq\sqrt{b}</math> (there is
no solution for <math>k>\sqrt{b}).</math> We then match <math>w</math> and its derivative at
<math>x=\pm\varsigma</math> to solve for <math>a</math> and <math>b</math>. This leads to two system of
equations, one for the even (<math>a_{1}=a_{2}</math> and <math>b_{2}=0</math> ) and one for the odd
solutions (<math>a_{1}=-a_{2}</math> and <math>b_{1}=0)</math>. The solution for the even solutions
is
<center><math>
w\left( x\right) =\left\{
\begin{matrix}
a_{1}e^{kx}, & x<-\varsigma,\\
b_{1}\cos\kappa x, & -\varsigma< x <\varsigma,\\
a_{1}e^{-kx}, & x>\varsigma,
\end{matrix}
\right.
</math></center>
If we impose the condition that the function and its derivative are continuous at
<math>x=\pm\varsigma</math> we obtain the following equation
<center><math>
\left(
\begin{matrix}

e^{-k\varsigma} & -\cos\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\sin\kappa\varsigma
\end{matrix}
\right) \left(
\begin{matrix}

a_{1}\\
b_{1}
\end{matrix}
\right) =\left(
\begin{matrix}

0\\
0
\end{matrix}
\right)
</math></center>
This has non trivial solutions when
<center><math>
\det\left(
\begin{matrix}

e^{-k\varsigma} & -\cos\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\sin\kappa\varsigma
\end{matrix}
\right) =0
</math></center>
which gives us the equation
<center><math>
-\kappa e^{-k\varsigma}\sin\kappa\varsigma+k\cos\kappa\varsigma
e^{-k\varsigma}=0
</math></center>
or
<center><math>
\tan\kappa\varsigma=\frac{k}{\kappa}
</math></center>
We know that <math>0<\kappa<\sqrt{b}</math> and if we plot this we see that we obtain a
finite number of solutions.

The solution for the odd solutions is
<center><math>
w\left( x\right) =\left\{
\begin{matrix}
a_{1}e^{kx}, & x <-\varsigma,\\
b_{2}\sin\kappa x, & -\varsigma< x <\varsigma,\\
-a_{1}e^{-kx} & x > \varsigma,
\end{matrix}
\right.
</math></center>
and again imposing the condition that the solution and its derivative is continuous
at <math>x=\pm\varsigma</math> gives
<center><math>
\left(
\begin{matrix}

e^{-k\varsigma} & \sin\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\cos\kappa\varsigma
\end{matrix}
\right) \left(
\begin{matrix}

a_{1}\\
b_{1}
\end{matrix}
\right) =\left(
\begin{matrix}

0\\
0
\end{matrix}
\right)
</math></center>
This can non trivial solutions when
<center><math>
\det\left(
\begin{matrix}

e^{-k\varsigma} & \sin\kappa\varsigma\\
ke^{-k\varsigma} & -\kappa\cos\kappa\varsigma
\end{matrix}
\right) =0
</math></center>
which gives us the equation
<center><math>
\kappa e^{-k\varsigma}a\cos\kappa\varsigma+k\sin\kappa\varsigma e^{-k\varsigma}=0
</math></center>
or
<center><math>
\tan\varsigma\kappa=-\frac{\kappa}{k}
</math></center>

=== Case when <math>\lambda>0</math> ===

When <math>\lambda>0</math> we write <math>\lambda=k^{2}</math> and we obtain solution
<center><math>
w\left( x\right) =\left\{
\begin{matrix}

\mathrm{e}^{-\mathrm{i}kx}+r\mathrm{e}^{\mathrm{i}kx}, & x <-\varsigma\\
b_{1}\cos\kappa x+b_{2}\sin\kappa x & -\varsigma< x <\varsigma\\
a\mathrm{e}^{-\mathrm{i}kx} & x>\varsigma
\end{matrix}
\right.
</math></center>
where <math>\kappa=\sqrt{b+k^{2}}.</math> Matching <math>w</math> and its derivaties at <math>x=\pm1</math> we
obtain
<center><math>
\left(
\begin{matrix}

-\mathrm{e}^{-\mathrm{i}k\varsigma} & \cos\kappa\varsigma & -\sin\kappa\varsigma & 0\\
ik\mathrm{e}^{-\mathrm{i}k\varsigma} & \kappa\sin\kappa\varsigma & \kappa\cos\kappa
\varsigma & 0\\
0 & \cos\kappa\varsigma & \sin\kappa\varsigma & -\mathrm{e}^{-\mathrm{i}k\varsigma}\\
0 & -\kappa\sin\kappa\varsigma & \kappa\cos\kappa\varsigma &
ik\mathrm{e}^{-\mathrm{i}k\varsigma}
\end{matrix}
\right) \left(
\begin{matrix}

r\\
b_{1}\\
b_{2}\\
a
\end{matrix}
\right) =\left(
\begin{matrix}

\mathrm{e}^{\mathrm{i}k}\\
ik\mathrm{e}^{-\mathrm{i}k}\\
0\\
0
\end{matrix}
\right)
</math></center>

Introduction to the Inverse Scattering Transform

2020-09-15T03:26:15Z

Meylan:

Connection betwen KdV and the Schrodinger Equation

2020-09-09T07:22:01Z

Meylan: /* Scattering Data */

{{nonlinear waves course
| chapter title = Connection betwen KdV and the Schrodinger Equation
| next chapter = [[Example Calculations for the KdV and IST]]
| previous chapter = [[Properties of the Linear Schrodinger Equation]]
}}

If we substitute the relationship
<center><math>
\partial_{x}^{2}w+uw=-\lambda w
</math></center>
into the KdV after some manipulation we obtain
<center><math>
\partial_{t}\lambda w^{2}+\partial_{x}\left( w\partial_{x}Q-\partial
_{x}wQ\right) =0
</math></center>
where <math>Q=\partial_{t}w+\partial_{x}^{3}w-3\left( \lambda-u\right)
\partial_{x}w.</math> If we integrate this equation then we obtain the result that
<center><math>
\partial_{t}\lambda=0
</math></center>
provided that the eigenfunction <math>w</math> is bounded (which is true for the bound
state eigenfunctions). This shows that the discrete eigenvalues are unchanged
and <math>u\left( x,t\right) </math> evolves according to the KdV. Many other
properties can be found

==Scattering Data==

For the discrete spectrum the eigenfunctions behave like
<center><math>
w_{n}\left( x\right) =c_{n}\left( t\right) e^{-k_{n}x}
</math></center>
as <math>x\rightarrow\infty</math> with
<center><math>
\int_{-\infty}^{\infty}\left( w_{n}\left( x\right) \right) ^{2}\mathrm{d}x=1
</math></center>
The continuous spectrum looks like
<center><math>
w\left( x,t\right) \approx \mathrm{e}^{-\mathrm{i}kx}+r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}
,\ \ \ x\rightarrow-\infty
</math></center>
<center><math>
w\left( x,t\right) \approx a\left( k,t\right) \mathrm{e}^{-\mathrm{i}kx},\ \ \ x\rightarrow
\infty
</math></center>
where <math>r</math> is the reflection coefficient and <math>a</math> is the transmission
coefficient. This gives us the scattering data at <math>t=0</math>
<center><math>
S\left( \lambda,0\right) =\left( \left\{ k_{n},c_{n}\left( 0\right)
\right\} _{n=1}^{N},r\left( k,0\right) ,a\left( k,0\right) \right)
</math></center>
The scattering data evolves as
<center><math>
k_{n}=k_{n}
</math></center>
<center><math>
c_{n}\left( t\right) =c_{n}\left( 0\right) e^{4k_{n}^{3}t}
</math></center>
<center><math>
r\left( k,t\right) =r\left( k,0\right) e^{8ik^{3}t}
</math></center>
<center><math>
a\left( k,t\right) =a\left( k,0\right)
</math></center>
We can recover <math>u</math> from scattering data. We write
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}
x}+\frac{1}{2\pi}\int_{-\infty}^{\infty}r\left( k,t\right) \mathrm{e}^{\mathrm{i}kx}\mathrm{d}k
</math></center>
Then solve
<center><math>
K\left( x,y;t\right) +F\left( x+y;t\right) +\int_{x}^{\infty}K\left(
x,z;t\right) F\left( z+y;t\right) \mathrm{d}z=0
</math></center>
This is a linear integral equation called the Gelfand-Levitan-Marchenko equation. We then find <math>u</math> from
<center><math>
u\left( x,t\right) =2\partial_{x}K\left( x,x,t\right)
</math></center>

==Reflectionless Potential==

In general the IST is difficult to solve. However, there is a simplification
we can make when we have a reflectionless potential (which we will see gives
rise to the soliton solutions). The reflectionless potential is the case when
<math>r\left( k,0\right) =0</math> for all values of <math>k</math> for some <math>u.</math> In this case
<center><math>
F\left( x,t\right) =\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}x}
</math></center>
then
<center><math>
K\left( x,y,t\right) +\sum_{n=1}^{N}c_{n}^{2}\left( t\right)
e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}K\left( x,z,t\right)
\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
From the equation we can see that
<center><math>
K\left( x,y,t\right) =-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}
</math></center>
If we substitute this into the equation,
<center><math>
-\sum_{n=1}^{N}v_{n}\left( x,t\right) e^{-k_{n}y}+\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }+\int_{x}^{\infty}
-\sum_{m=1}^{N}v_{m}\left( x,t\right) e^{-k_{m}z}\sum_{n=1}^{N}c_{n}
^{2}\left( t\right) e^{-k_{n}\left( y+z\right) }\mathrm{d}z=0
</math></center>
which leads to
<center><math>
-\sum_{n=1}^{N} v_{n}\left( x,t\right) e^{-k_{n}y}
+\sum_{n=1}^{N}c_{n}^{2}\left( t\right) e^{-k_{n}\left( x+y\right) }
-\sum_{n=1}^{N}\sum_{m=1}^{N}\frac{ c_{n}^{2}\left(
t\right) }{k_{n}+k_{m}}v_{m}\left( x\right) e^{-k_{m}x}e^{-k_{n}\left(
y+x\right) }=0
</math></center>
and we can eliminate the sum over <math>n</math> and the
<math>e^{-k_{n}y}</math> to obtain
<center><math>
-v_{n}\left( x,t\right) +c_{n}^2\left( t\right) e^{-k_{n}x}-\sum_{m=1}
^{N}\frac{c_{n}^2\left( t\right) }{k_{n}+k_{m}}
v_{m}\left( x,t\right) e^{-\left( k_{m}+k_{n}\right) x}=0
</math></center>
which is an algebraic (finite dimensional system) for the unknows <math>v_{n}.</math>.

We can write this as
<center><math>
\left( \mathbf{I}+\mathbf{C}\right) \vec{v}=\vec{f}
</math></center>
where <math>f_{m}=c_{m}^2\left( t\right) e^{-k_{m}x}</math> and
the elements of <math>\mathbf{C}</math> are given by
<center><math>
c_{mn} = \frac{c_{n}^2\left( t\right)}
{k_{n}+k_{m}}e^{-\left( k_{m}+k_{n}\right) x}
</math></center>
This gives us
<center><math>
K\left( x,y,t\right) =-\sum_{m=1}^{N} \left(
\mathbf{I}+\mathbf{C}\right) ^{-1}\vec{f}e^{-k_{m}y}
</math></center>
We then find <math>u(x,t)</math> from <math>K</math>.

=== Single Soliton Example ===

If <math>n=1</math> (a single soliton
solution) we get
<center><math>\begin{matrix}
K\left( x,x,t\right) & =-\frac{c_{1}^2\left( t\right)
e^{-2k_{1}x}}{1+\frac{c_{1}^2\left( t\right) }
{2k_{1}}e^{-2 k_{1} x}}\\
& =\frac{-1}{e^{-2k_{1}x + 8k_{1}^{3}t-\alpha} + 1/2k_1}
\end{matrix}</math></center>
where <math>e^{-\alpha}=1/c_{0}^{2}\left( 0\right) .</math> Therefore
<center><math>\begin{align}
u\left( x,t\right) & =2\partial_{x}K\left( x,x,t\right) \\
& =\frac{4k_{1}e^{2k_{1}x-8k_{1}^{3}t-\alpha}}{\left( e^{2k_{1}
x-8k_{1}^{3}t-\alpha} + 1/2k_1\right) ^{2}}\\
& =\frac{8k_{1}^{2}}{\left( \sqrt{2k_{1}}e^{\theta}+e^{-\theta}/\sqrt
{2k_{1}}\right) ^{2}}\\
& =2k_{1}^{2}\mbox{sech}^{2}\left\{ k_{1}\left( x-x_{0}\right) -4k_{1}
^{3}t\right\}
\end{align}</math></center>
where <math>\theta=k_{1}x-4k^{3}t-\alpha/2</math> and <math>\sqrt{2k}e^{-\alpha/2}=e^{-kx_{0}
}</math>. This is of course the single soliton solution.