---
title: 8-1. Fourier approximations - KdV equation
tags: Spectral method
---

**Reference**
1. Spectral methods by *J. Shen*, *T. Tang* and *L.-L. Wang*

---

## Korteweg–de Vries (KdV) Equation

Consider the KdV equation in the whole space:
$$
u_t + uu_{y}+u_{yyy} = 0, \quad y\in(-\infty, \infty), \quad t\ge 0,
$$
with initial condition
$$
u(y,0) = u_0(y), \quad y\in(-\infty, \infty),
$$
which has the exact solution
$$
u(y,t) = 12 \kappa^2 sech^2(\kappa(y-y_0)-4\kappa^3 t),
$$
where $y_0$ is the center of the initial profile $u(y,0)$, and $\kappa$ is a constant related to the traveling phase speed.

Since $u(y,t)$ decays exponentially to zero as $|y|\to\infty$, we can truncate the infinite interval to a finite one $[-\pi L, \pi L]$ with $L > 0$, and approximate the boundary conditions by the periodic boundary conditions on $[-\pi L, \pi L]$. It is expected that the initial-boundary valued problem with periodic boundary conditions can provide a good approximation to the original initial-valued problem as long as the soliton does not reach the boundaries.

For convenience, we map the interval $[-\pi L, \pi L]$ to $[0, 2\pi]$ through the coordinate transform:
$$
x = \frac{y}{L}+\pi, \quad y = L(x-\pi), \quad x\in[0, 2\pi], \quad y\in[-\pi L, \pi L],
$$
and denote
$$
v(x,t) = u(y,t), \quad v_0(x) = u_0(y).
$$
The transformed KdV equation reads
$$
v_t + \frac{1}{L}vv_{x}+\frac{1}{L^3}v_{xxx} = 0, \quad x\in [0, 2\pi], \quad t\ge 0,
$$
$$
v(x,0) = v_0(x).
$$
Writing $v(x,t) = \sum_k \hat{v}_k(t) e^{ikx}$, taking the inner product of the equation with $e^{ikx}$ and using the fact that $vv_{x} = \frac{1}{2}(v^2)_x$, we obtain
$$
\frac{d\hat{v}_k}{dt} - \frac{ik^3}{L^3}\hat{v}_k + \frac{ik}{2L}\widehat{(v^2)}_k = 0, \quad k=0,\pm 1, \cdots,
$$
with the initial condition
$$
\hat{v}_k(0) = \frac{1}{2\pi}\int^{2\pi}_0 v_0(x)e^{ikx}dx.
$$

As suggestd by Trefethen we solve the system using method of integrating factor. The equations are multiplied by the integrating factor $e^{-ik^3 t/L^3}$ to
$$
\frac{d}{dt}\left[e^{-ik^3t/L^3}\hat{v}_k\right] = -\frac{ik}{2L}e^{-ik^3t/L^3}\widehat{(v^2)}_k, \quad k=0,\pm 1, \cdots.
$$
This system can then be solved by standard ODE solver such as Runge-Kutta method, as shown in the following example.

% p27.m - Solve KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(-ik^3t)*u-hat.

% Set up grid and two-soliton initial data:
  N = 256; dt = .4/N^2; x = (2*pi/N)*(-N/2:N/2-1)';
  A = 25; B = 16; clf, drawnow, set(gcf,'renderer','zbuffer')
  u = 3*A^2*sech(.5*(A*(x+2))).^2 + 3*B^2*sech(.5*(B*(x+1))).^2; 
  v = fft(u); k = [0:N/2-1 0 -N/2+1:-1]'; ik3 = 1i*k.^3;

% Solve PDE and plot results:
  tmax = 0.006; nplt = floor((tmax/25)/dt); nmax = round(tmax/dt);
  udata = u; tdata = 0; h = waitbar(0,'please wait...');
  for n = 1:nmax
    t = n*dt; g = -.5i*dt*k;
    E = exp(dt*ik3/2); E2 = E.^2;
    a = g.*fft(real( ifft(     v    ) ).^2);
    b = g.*fft(real( ifft(E.*(v+a/2)) ).^2);     % 4th-order
    c = g.*fft(real( ifft(E.*v + b/2) ).^2);     % Runge-Kutta
    d = g.*fft(real( ifft(E2.*v+E.*c) ).^2);
    v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
    if mod(n,nplt) == 0 
      u = real(ifft(v)); waitbar(n/nmax)
      udata = [udata u]; tdata = [tdata t];
    end
  end
  waterfall(x,tdata,udata'), colormap(1e-6*[1 1 1]); view(-20,25)
  xlabel x, ylabel t, axis([-pi pi 0 tmax 0 2000]), grid off
  set(gca,'ztick',[0 2000]), close(h), pbaspect([1 1 .13])

One can also re-write the equations by make the following substitution
$$
g_k=e^{-ik^3t/L^3}, \quad \tilde{u}_k = g_k\hat{v}_k,
$$
the system is then written as
$$
\frac{d}{dt}\tilde{u}_k = -\frac{ik}{2L}g_k*fft\left(\left[ifft(g^{-1}*\tilde{u})\right]^2\right), \quad k=0,\pm 1, \cdots.
$$
The ODE solvers can be directly applied.

### Julia code(昆毅):

using LinearAlgebra
using Plots
using FFTW

# p27.m - Solve KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
#         FFT with integrating factor v = exp(-ik^3t)*u-hat.

# Set up grid and two-soliton initial data:
N = 256; dt = .4/N^2; x = (2*pi/N)*(-N/2:N/2-1)';
A = 25; B = 16;
u=zeros(length(x))
for i = 1:length(x)
    u[i] = 3*A^2*sech(.5*(A*(x[i]+2))).^2 + 3*B^2*sech(.5*(B*(x[i]+1))).^2
end
v = fft(u); k = [0:N/2-1;0;-N/2+1:-1]; ik3 = 1im*k.^3;

# Solve PDE and plot results:
tmax = 0.006; nplt = floor((tmax/25)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;
for n = 1:nmax
    t = n*dt; g = -.5im*dt*k;
    E = exp.(dt*ik3/2); E2 = E.^2;
    a = g.*fft(real( ifft(     v    ) ).^2);
    b = g.*fft(real( ifft(E.*(v+a/2)) ).^2);     # 4th-order
    c = g.*fft(real( ifft(E.*v + b/2) ).^2);     # Runge-Kutta
    d = g.*fft(real( ifft(E2.*v+E.*c) ).^2);
    v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
    
    if mod(n,nplt) == 0
        u = real(ifft(v));
        udata = [udata u]; 
        tdata = [tdata t];
    end
end
  
p = plot(x',udata[:,1])
anim = Animation()
for i = range(2, stop=26, length=250)
    T=round(tdata[Int(floor(i))],digits=6)
    plot(x', udata[:,Int(floor(i))],label="t=$T",ylim=(0,2000))
    frame(anim)
end
gif(anim, "KdV Equation.gif")

@Woodystick (1). 研究看看 label 時間那邊能不能顯示到小數點後幾位就好, 比如說6或7位, 顯示這麼長實在很醜. (2). 然後也要固定 y 軸範圍, 不然跳來跳去的. (3). 還有 x-t-u 三維的圖也研究看看怎麼做.



### Python code (沛聖)
The following code solves the boundary value problem
$$
\frac{d}{dt} e^{-ik^3t} v_k(t) = -\frac{ik}{2} e^{-ik^3t} (\widehat{u^2})_k,\quad k=0,\pm 1,\dots,\pm (N/2-1),N/2,
$$
where $v_k(t) = \widehat{u}(k)$, using Runge-Kutta 4th Order Method
#### Code

import numpy as np
from numpy.fft import fft,ifft
from numpy.linalg import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# Initialize variables
N = 256
dt = .4/N**2
x = 2*np.pi/N * np.arange(-N/2,N/2)
k = np.concatenate( (np.arange(N/2) , np.array([0]) , np.arange(-N/2+1,0)) )
ik3 = 1j*k**3

# Initialize initial condition
A = 25
B = 16
u = 3*(A**2)/np.cosh(0.5*(A*(x+2)))**2 + 3*(B**2)/np.cosh(0.5*(B*(x+1)))**2
v = fft(u)

# Initialize parameter for for-loop
tmax = 0.006
nplt = np.floor(tmax/25/dt)
nmax = round(tmax/dt)

# Create container to store u and t 
udata = np.array(u)
tdata = np.array([0])

# Iterate over  [0,tmax]
for n in range(1,nmax+1):
    t = n*dt
    g = -.5j*dt*k
    E = np.exp(dt*ik3/2)
    E2 = E**2
    a = g*fft( (ifft(v).real)**2 )
    b = g*fft( (ifft(E*(v+a/2)).real)**2 )
    c = g*fft( (ifft(E*v + b/2).real)**2 )
    d = g*fft( (ifft(E2*v + E*c).real)**2 )
    v = E2*v + (E2*a + 2*E*(b+c) +d)/6
    if n%nplt == 0:
        u = ifft(v).real
        udata = np.vstack( (udata,u) )
        tdata = np.hstack( (tdata,t) )

# Create the waterfall plot
x_grid , t_grid = np.meshgrid(x,tdata)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_wireframe(x_grid,t_grid,udata)
ax.set_xlabel('x')
ax.set_ylabel('t')
plt.show()

#### Results

### C++ code (沂君)

//compile with g++ EX8_1_graph.cpp -L.\. -lfftw3 -o E81G
#include <iostream>
#include <cmath>
#include <fstream>
#include <complex>
#include <vector>

#include "fftw3.h"

#define REAL 0
#define IMAG 1
#define PI 3.141592653589793

using namespace std;
void copy( fftw_complex c2[], fftw_complex c1[], int N);  //copy c1 into c2 
void norm_real( fftw_complex c[], int N);  //normalize and set imaginary part to 0
fftw_complex* mult( fftw_complex c1[], fftw_complex c2[], int N);  //multiply
fftw_complex* multc( fftw_complex c[], int num, int N);  //multiply a constant
fftw_complex* add( fftw_complex c1[], fftw_complex c2[], int N);  //add
fftw_complex* div( fftw_complex c[], int num, int N);  //divide by a constant
void u_to_file( fftw_complex v[], double t, double x[], int N, ofstream &outfile);  //output real u at t

int main(){

    ofstream outfile("./EX8_output.txt");
    outfile.precision(15);
    outfile << "# X T U" << endl;

    //% Set up grid and two-soliton initial data:
    int N = 256, A=25, B=16;
    double dt = 0.4/powl(N,2);
    double x[N], k[N];
    fftw_complex u[N], v[N], ik3[N];
    fftw_complex *in, *out;
    vector<int> tdata;

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    for(int i=0; i<N; i++)
    {
        x[i] = (-N/2+i)*2*PI/N;
        u[i][REAL] = 3*powl(A,2)*powl(1/coshl(0.5*(A*(x[i]+2))),2) + 3*powl(B,2)*powl(1/coshl(0.5*(B*(x[i]+1))),2);
        u[i][IMAG] = 0;
    }

    u_to_file(u, 0, x, N,outfile);
    outfile << endl << endl;

    fftw_plan plan1 = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_MEASURE);
    fftw_plan plan2 = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, FFTW_MEASURE);

    copy(in, u, N);
    fftw_execute(plan1);
    copy(v, out, N);

    for(int i=0; i<N; i++)
    {
        if(i == (N/2)) k[i] = 0;
	else if(i < (N/2)) k[i] = i;
	else if(i > (N/2)) k[i] = -N + i;
    
        ik3[i][REAL] = 0;
        ik3[i][IMAG] = powl(k[i],3);
    }

    //% Solve PDE and plot results:
    double tmax = 0.006, nmax = round(tmax/dt), t;
    int nplt = floor((tmax/25)/dt);
    fftw_complex g[N], E[N], E2[N];
    fftw_complex *a, *b, *c, *d;

    for(int n=1; n<=nmax; n++)
    {
        t = n*dt;

        for(int i=0; i<N; i++)
        {
            g[i][REAL] = 0;
            g[i][IMAG] = -0.5*dt*k[i];
            
            complex<double> c1(dt*ik3[i][REAL]*0.5,dt*ik3[i][IMAG]*0.5);
            E[i][REAL] = real(exp(c1));
            E[i][IMAG] = imag(exp(c1));

            complex<double> c2(E[i][REAL],E[i][IMAG]);
            E2[i][REAL] = real(pow(c2,2));
            E2[i][IMAG] = imag(pow(c2,2));
        }
    
        //a
        copy(in, v, N);
        fftw_execute(plan2);
        norm_real(out, N);
        copy(in, mult(out, out, N), N);
        fftw_execute(plan1);
        a = mult(g, out, N);

        //b
        copy(in, mult(E, add(v, div(a, 2, N), N), N), N);
        fftw_execute(plan2);
        norm_real(out, N);
        copy(in, mult(out, out, N), N);
        fftw_execute(plan1);
        b = mult(g, out, N);

        //c
        copy(in, add(mult(E, v, N), div(b, 2, N), N), N);
        fftw_execute(plan2);
        norm_real(out, N);
        copy(in, mult(out, out, N), N);
        fftw_execute(plan1);
        c = mult(g, out, N);

        //d
        copy(in, add(mult(E2, v, N), mult(E, c, N), N), N);
        fftw_execute(plan2);
        norm_real(out, N);
        copy(in, mult(out, out, N), N);
        fftw_execute(plan1);
        d = mult(g, out, N);

        //v
        copy(v,add(mult(E2, v, N) ,div(add(add(mult(E2, a, N), multc(mult(E, add(b, c, N), N), 2, N) , N), d, N), 6, N) ,N) ,N);

        if(n%nplt == 0)
        {
            copy(in, v, N);
            fftw_execute(plan2);
            norm_real(out, N);
            copy(u, out, N);
            tdata.push_back(t);
            u_to_file(u, t, x, N,outfile);
            outfile << endl << endl;
        }

    }

    //Plot results:
    FILE* gp = _popen("gnuplot","w");
    if (gp == NULL)
    {
        cout<<("Cannot open gnuplot!\n")<<endl;;
        return 0;
    }
    
    fprintf(gp,"set xlabel 'X'\nset ylabel 'T'\nset zlabel 'U'\nset zrange [0:2000]\nset view 30,351\n");
    fprintf(gp,"set term png\nset output 'EX8_1_g.png'\n");
    fprintf(gp,"splot 'EX8_output.txt' using 1:2:3 with lines title 'u'\n");
    
    fprintf(gp,"set term gif size 800,400 \nset term gif animate delay 20 \nset output 'EX8_1.gif'\n");
    fprintf(gp,"set ylabel 'U'\n");
    for(int i=0; i<=nmax/39; i++)
    {
        fprintf(gp,"plot 'EX8_output.txt' using 1:3 index %i with lines lw 1.5 title 't=%f'\n",i+1,39*(i+1)*dt);
    }
    fprintf(gp,"set output\n");
    _pclose(gp);
    
    fftw_destroy_plan(plan1);
    fftw_destroy_plan(plan2);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

void copy( fftw_complex c2[], fftw_complex c1[], int N)
{
    for(int i=0; i<N; i++)
    {
        c2[i][REAL] = c1[i][REAL];
        c2[i][IMAG] = c1[i][IMAG];
    }
}

void norm_real( fftw_complex c[], int N)
{
    for(int i=0; i<N; i++)
    {
        c[i][REAL] = c[i][REAL]/N;
        c[i][IMAG] = 0;
    }
}

fftw_complex* mult( fftw_complex c1[], fftw_complex c2[], int N)
{
    fftw_complex *c3;
    c3 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    for(int i=0; i<N; i++)
    {
        complex<double> c_1(c1[i][REAL],c1[i][IMAG]);
        complex<double> c_2(c2[i][REAL],c2[i][IMAG]);
        
        c3[i][REAL] = real(c_1*c_2);
        c3[i][IMAG] = imag(c_1*c_2);
    }

    return c3;
}

fftw_complex* multc( fftw_complex c[], int num, int N)
{
    fftw_complex *c3;
    c3 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    for(int i=0; i<N; i++)
    {
        c3[i][REAL] = c[i][REAL]*num;
        c3[i][IMAG] = c[i][IMAG]*num;
    }

    return c3;
}

fftw_complex* add( fftw_complex c1[], fftw_complex c2[], int N)
{
    fftw_complex *c3;
    c3 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    for(int i=0; i<N; i++)
    {
        c3[i][REAL] = c1[i][REAL] + c2[i][REAL];
        c3[i][IMAG] = c1[i][IMAG] + c2[i][IMAG];
    }

    return c3;
}

fftw_complex* div( fftw_complex c[], int num, int N)
{
    fftw_complex *c3;
    c3 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    for(int i=0; i<N; i++)
    {
        c3[i][REAL] = c[i][REAL]/num;
        c3[i][IMAG] = c[i][IMAG]/num;
    }

    return c3;
}

void u_to_file( fftw_complex v[], double t, double x[], int N, ofstream &outfile)
{
	for(int i=0; i<N-1; i++)
	{
	    outfile << x[i] << " " << t << " " << v[i][REAL] << endl; 
	}
}

#### Results: