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:
