WikiWaves - User contributions [en]

Traffic Waves

2010-11-14T00:13:38Z

Rua: /* Speed of the shock */

{{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) 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_fan1.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)} + \int_{s(t)}^{x_2}\right)
\rho(x,t)\mathrm{d}x = 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) 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]]
|}

[[Category:Simple Nonlinear Waves]]

Numerical Solution of the KdV

2010-11-09T03:58:40Z

Rua: /* Numerical Solution of the KdV */

{{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 \\
&=&-k\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>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.]]

|}