This is an archival dump of old wiki content --- see scipy.org for current material.
Please see http://scipy-cookbook.readthedocs.org/

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

The following plot is created by the above code:

kdv2.png

SciPy: Cookbook/KdV (last edited 2015-10-24 17:48:25 by anonymous)