{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Lecture 22: PDEs1 \n", "### Phys 503, J.R. Gladden, University of Mississippi" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Partial Differential Equations\n", "Partial differential equations mathematically describe a system which depends on multiple variables and their derivatives. Examples are:\n", "- Wave equation (acoustics, optics)\n", "- Laplace equation (electrostatics)\n", "- Schrodinger equation (quantum mechanics)\n", "- Navier-Stokes equation (fluid flow)\n", "\n", "\n", "An example of a typical 2nd order PDE would be the 2D Laplace equation for the electrostatic potential (more on this in Lecture 22):\n", "$$ \\mu \\left( \\frac{\\partial^2 \\Phi}{\\partial x^2} + \\frac{\\partial^2 \\Phi}{\\partial y^2} \\right)=0 $$\n", "\n", "A another classic example is the time evolution of the 1D wave equation:\n", "$$ \\frac{\\partial^2 \\Psi}{\\partial t^2} = c^2 \\frac{\\partial^2 \\Psi}{\\partial x^2} $$" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### The Solver\n", "This solver supplies both vector and scalar approaches to solving the 1D wave equation. Imports first..." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "%matplotlib wx\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import time\n", "import matplotlib.animation as animation" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def solver0(I,f,c,L,n,dx,tstop,dt=0.0,method='v'):\n", " '''\n", "\tWave equation solver.\n", "\tJRG, April 2014\n", "\tInput:\n", "\t\tI - callable initial condition function \n", "\t\tf - callable forcing function\n", "\t\tc - wavespeed\n", "\t\tL - length of domain\n", "\t\tn - number of space steps\n", "\t\tdx- space step\n", "\t\ttstop - stop time (assume start at t=0.0)\n", "\t\tdt - set to 0.0 to allow program to choose time step\n", "\t\tmethod - (v)ector or (s)calar.\n", "\tOutput:\n", "\t\t- computed dt [float]\n", "\t\t- fullU (solution at all time steps) [list of arrays]\n", "\t\t- method used [string]\n", " '''\n", "\t\n", "\t#Set up an empty list container for solutions at ALL time steps\n", " fullU=[]\n", " #if dt is not specified, make largest that is stable\n", " \n", " if dt <= 0: dt=dx/float(c)\n", "# print dt\n", " \n", " C2=(c*dt/dx)**2\n", " dt2=dt*dt\n", " \n", " up=np.zeros(n+1)\t\t#solution array at t=t+dt\n", " u=up.copy()\t\t\t# solution at t\n", " um=up.copy()\t\t#solution at t-dt\n", "\n", " t=0.0\n", " \n", " #Set up intitial conditions\n", " for i in range(n+1):\n", " u[i]=I(x[i])\n", " for i in range(1,n):\n", " um[i]=u[i] + 0.5*C2*(u[i-1] - 2*u[i] + u[i+1]) + dt2*f(x[i],t)\n", " \n", " #Set up Dirichlet boundary conditions\n", " um[0]=0.\n", " um[n]=0.\n", " \n", " while t <= tstop:\n", " t_old = t; t+=dt\n", " if method == 's':\n", " for i in range(1,n):\n", " up[i] = -um[i] +2*u[i] + C2*(u[i-1] - 2*u[i] + u[i+1]) + dt2*f(x[i],t_old)\n", " elif method == 'v':\n", " up[1:n] = -um[1:n] +2*u[1:n] +C2*(u[0:n-1] - 2*u[1:n] +u[2:n+1]) + dt2*f(x[1:n],t_old)\n", " else:\n", " print \"Method needs to be 's' for scalar or 'v' for vector. \\n Exiting ...\"\n", " break\n", " up[0] = 0.; up[n] = 0.\n", " fullU.append(u)\n", "\n", " um=u.copy(); u=up.copy()\n", " return dt, fullU, method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualization\n", "Here we provide a dynamic animation and static waterfall options to visualize the wave evolution." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def animate(fullU):\n", "\tanimFig=plt.figure()\n", "\tline,=plt.plot(x,fullU[0],'b-')\n", "\tspeed=1.0*dt\n", "\t\n", "\tdef init():\n", "\t\tline.set_ydata(np.zeros(len(x))+np.nan)\n", "\t\tplt.ylim([-1,1])\n", "\t\tplt.xlabel('X',fontsize=18)\n", "\t\tplt.ylabel('$\\psi(x)$',fontsize=20)\n", "\t\tplt.title('Solution to the 1D Wave Equation')\n", "\t\treturn line,\n", "\t\t\n", "\tdef update(i): \n", "\t\tline.set_ydata(fullU[i])\n", "\t\ttime=float(i)*dt\n", "\t\t#text(max(x)*0.9,max(fullU[0])*0.9,'Time=%2.2f s'%time)\n", "\t\tax=plt.axes()\n", "\t\t#txt=text(0.8,0.9,'Time=%2.2f s'%time,transform = ax.transAxes)\n", "\t\treturn line, #, txt\n", "\t\t\n", "\tani = animation.FuncAnimation(\n", "\t\tanimFig,\n", "\t\tupdate,\n", "\t\trange(len(fullU)),\n", "\t\tblit=True,\n", "\t\tinterval=20,\n", "\t\tinit_func=init\n", "\t\t)\n", "\tplt.show()\n", "\t\n", "def waterfallPlot(fullU,dstep=25,offsetFactor=0.5):\n", "\tnumPlots=0\n", "\toffset=0.0\n", "\tfor i in range(len(fullU)):\n", "\t\tif offset