{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# p20: Second order wave equation in 2-D via FFT"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the wave equation in 2-d\n",
"\n",
"$$\n",
"u_{tt} = u_{xx}+u_{yy}, \\qquad -1 < x,y < 1, \\qquad t >0\n",
"$$\n",
"\n",
"with $u = 0$ on the boundary and initial condition\n",
"\n",
"$$\n",
"u(x,y,0) = e^{-40((x-0.4)^2 + y^2)}, \\qquad u_{t}(x,y,0) = 0\n",
"$$\n",
"\n",
"In the following code, `i` index runs along $x$-axis and `j` index runs along $y$-axis; see how the `meshgrid` function is used with `indexing='ij'`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"from numpy import meshgrid,cos,pi,round,exp,real,remainder,zeros,fliplr,flipud,array,arange\n",
"from numpy.fft import fft, ifft\n",
"from matplotlib.pyplot import subplot, figure ,title,axis\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from matplotlib.pyplot import figure,subplot,plot,title,axis,xlabel,ylabel\n",
"from matplotlib import cm\n",
"from scipy.interpolate import RegularGridInterpolator"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Grid and inital Data:\n",
"N = 24; x = cos(pi*arange(0,N+1)/N); y = x;\n",
"t = 0.0; dt = (6.0)/(N**2)\n",
"xx, yy = meshgrid(x,y,indexing='ij')\n",
"plotgap = int (round( (1.0/3.0) / (dt))); dt = (1.0/3.0)/(plotgap); \n",
"vv = exp(-40*((xx-0.4)**2 + yy**2)); # Solution at n\n",
"vvold = vv; # Solution at n-1\n",
"\n",
"w1_hat = 1j*zeros(2*N);\n",
"w1_hat[0:N] = 1j*arange(0,N)\n",
"w1_hat[N+1:] = 1j*arange(-N+1,0)\n",
"w2_hat = 1j*zeros(2*N);\n",
"w2_hat[0:N+1] = arange(0,N+1)\n",
"w2_hat[N+1:] = arange(-N+1,0)\n",
"\n",
"#Time stepping Leapfrog Formula:\n",
"fig = figure(figsize=(12,12))\n",
"k = 1;\n",
"for n in range(0,(3*plotgap)+1):\n",
" t = n*dt;\n",
" if (remainder(n+0.5,plotgap) < 1):\n",
" ax = fig.add_subplot(2,2,k,projection ='3d')\n",
" f = RegularGridInterpolator((x,y),vv,method='cubic');\n",
" xxx = arange(-1.,1.+1./16,1./16);\n",
" X,Y = meshgrid(xxx,xxx,indexing='ij');\n",
" vvv = f((X,Y))\n",
" ax.plot_surface(X,Y,vvv,rstride=1,cstride=1,cmap=cm.jet,edgecolor='black')\n",
" ax.set_zlim3d([-0.15,1])\n",
" ax.set_xlim3d([-1,1])\n",
" ax.set_ylim3d([-1,1])\n",
" ax.view_init(elev=40., azim=250.)\n",
" title(\"$ t $= \" +str(t))\n",
" xlabel(\"x\"); ylabel(\"y\");\n",
" k = k+1;\n",
" \n",
" uxx = zeros((N+1,N+1)); uyy = zeros((N+1,N+1));\n",
" ii = arange(1,N);\n",
" \n",
" for i in range(1,N): # Compute uyy\n",
" v = vv[i,:]; \n",
" V = list(v) + list(flipud(v[ii]));\n",
" U = real(fft(V));\n",
" W1 = real(ifft(w1_hat * U))\n",
" W2 = real(ifft((-w2_hat**2) * U))\n",
" uyy[i,ii] = W2[ii]/(1-x[ii]**2) - (x[ii]*W1[ii])/(1-x[ii]**2)**(3.0/2);\n",
" for j in range(1,N): # Compute uxx\n",
" v = vv[:,j]; \n",
" V = list(v) + list(flipud(v[ii]));\n",
" U = real(fft(V))\n",
" W1 = real(ifft(w1_hat * U))\n",
" W2 = real(ifft(-(w2_hat**2) * U))\n",
" uxx[ii,j] = W2[ii]/(1-y[ii]**2) - y[ii]*W1[ii]/(1-y[ii]**2)**(3.0/2.0);\n",
" vvnew = 2*vv - vvold + dt**2 *(uxx+uyy) # Solution at n+1\n",
" vvold = vv; vv = vvnew;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"Make an animation of the solution."
]
}
],
"metadata": {
"anaconda-cloud": {},
"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.12.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}