{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# p27: Solve KdV equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the KdV equation using FFT\n", "$$\n", "u_t + u u_x + u_{xxx} = 0, \\qquad x \\in [-\\pi,\\pi]\n", "$$" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format='svg'\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from matplotlib.collections import LineCollection\n", "from numpy import pi,cosh,exp,round,zeros,arange,real\n", "from numpy.fft import fft,ifft\n", "from matplotlib.pyplot import figure" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-10-29T10:43:49.023670\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.4.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Set up grid and differentiation matrix:\n", "N = 256; dt = 0.4/N**2; x = (2*pi/N)*arange(-N/2,N/2);\n", "A, B = 25.0, 16.0\n", "u = 3*A**2/cosh(0.5*A*(x+2))**2 + 3*B**2/cosh(0.5*B*(x+1))**2\n", "v = fft(u); \n", "k = zeros(N); k[0:N//2] = arange(0,N/2); k[N//2+1:] = arange(-N/2+1,0,1)\n", "ik3 = 1j*k**3\n", "\n", "# Time-stepping by Runge-Kutta\n", "tmax = 0.006; nplt = int(round((tmax/25)/dt))\n", "nmax = int(round(tmax/dt))\n", "udata = []; udata.append(list(zip(x, u)))\n", "tdata = [0.0]\n", "for n in range(1,nmax+1):\n", " t = n*dt; g = -0.5j*dt*k\n", " E = exp(dt*ik3/2); E2 = E**2\n", " a = g * fft(real(ifft( v ))**2)\n", " b = g * fft(real(ifft( E*(v+a/2) ))**2)\n", " c = g * fft(real(ifft( E*v+b/2 ))**2)\n", " d = g * fft(real(ifft( E2*v+E*c ))**2)\n", " v = E2*v + (E2*a + 2*E*(b+c) + d)/6\n", " if n%nplt == 0:\n", " u = real(ifft(v))\n", " udata.append(list(zip(x, u)))\n", " tdata.append(t);\n", "\n", "fig = figure(figsize=(12,10))\n", "ax = fig.add_subplot(111,projection='3d')\n", "poly = LineCollection(udata)\n", "poly.set_alpha(0.5)\n", "ax.add_collection3d(poly, zs=tdata, zdir='y')\n", "ax.set_xlabel('x')\n", "ax.set_xlim3d(-pi, pi)\n", "ax.set_ylabel('t')\n", "ax.set_ylim3d(0, tmax)\n", "ax.set_zlabel('u')\n", "ax.set_zlim3d(0, 2000)\n", "ax.view_init(80,-115);" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 1 }