This page shows how the Korteweg-de Vries equation can be solved on a periodic domain using the method of lines, with the spatial derivatives computed using the pseudo-spectral method. In this method, the derivatives are computed in the frequency domain by first applying the FFT to the data, then multiplying by the appropriate values and converting back to the spatial domain with the inverse FFT. This method of differentiation is implemented by the diff function in the module scipy.fftpack.

We discretize the spatial domain, and compute the spatial derivatives using the diff function defined in the scipy.fftpack module. In the following code, this function is given the alias psdiff to avoid confusing it with the numpy function diff. By discretizing only the spatial dimension, we obtain a system of ordinary differential equations, which is implemented in the function kdv(u, t, L). The function kdv_solution(u0, t, L) uses scipy.integrate.odeint to solve this system.

   1 
   2 import numpy as np
   3 from scipy.integrate import odeint
   4 from scipy.fftpack import diff as psdiff
   5 
   6 
   7 def kdv_exact(x, c):
   8     """Profile of the exact solution to the KdV for a single soliton on the real line."""
   9     u = 0.5*c*np.cosh(0.5*np.sqrt(c)*x)**(-2)
  10     return u
  11 
  12 def kdv(u, t, L):
  13     """Differential equations for the KdV equation, discretized in x."""
  14     # Compute the x derivatives using the pseudo-spectral method.
  15     ux = psdiff(u, period=L)
  16     uxxx = psdiff(u, period=L, order=3)
  17 
  18     # Compute du/dt.    
  19     dudt = -6*u*ux - uxxx
  20 
  21     return dudt
  22 
  23 def kdv_solution(u0, t, L):
  24     """Use odeint to solve the KdV equation on a periodic domain.
  25     
  26     `u0` is initial condition, `t` is the array of time values at which
  27     the solution is to be computed, and `L` is the length of the periodic
  28     domain."""
  29 
  30     sol = odeint(kdv, u0, t, args=(L,), mxstep=5000)
  31     return sol
  32 
  33 
  34 if __name__ == "__main__":
  35     # Set the size of the domain, and create the discretized grid.
  36     L = 50.0
  37     N = 64
  38     dx = L / (N - 1.0)
  39     x = np.linspace(0, (1-1.0/N)*L, N)
  40 
  41     # Set the initial conditions.
  42     # Not exact for two solitons on a periodic domain, but close enough...
  43     u0 = kdv_exact(x-0.33*L, 0.75) + kdv_exact(x-0.65*L, 0.4)
  44 
  45     # Set the time sample grid.
  46     T = 200
  47     t = np.linspace(0, T, 501)
  48 
  49     print "Computing the solution."
  50     sol = kdv_solution(u0, t, L)
  51 
  52 
  53     print "Plotting."
  54 
  55     import matplotlib.pyplot as plt
  56 
  57     plt.figure(figsize=(6,5))
  58     plt.imshow(sol[::-1, :], extent=[0,L,0,T])
  59     plt.colorbar()
  60     plt.xlabel('x')
  61     plt.ylabel('t')
  62     plt.axis('normal')
  63     plt.title('Korteweg-de Vries on a Periodic Domain')
  64     plt.show()

The following plot is created by the above code:

Cookbook/KdV (last edited 2011-02-22 17:32:23 by Warren Weckesser)