{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Can you Hear the Shape of a Drum?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have seen that the PageRank Problem is defined in the form $A \\vec{v} = \\lambda \\vec{v}$, where the transition of users from web page to web page reaches a steady state: even though the matrix $A$ re-distributes users to some new sites, the number of users on each web page doesn't change. In general, this represents a class of problems that are important in disciplines that require modeling. \n",
"\n",
"In the PageRank problem, the state $\\vec{v}$ tells you how many users there are on each site at a particular time, and $\\lambda$ tells you the score for each page. When you use the ($A \\vec{v} = \\lambda \\vec{v}$) format for vibrational modes of a string or a membrane, the state $\\vec{v}$ tells you how much displacement there is at a particular location on the object, and $\\lambda$ tells you how much energy there is in that particular vibrational mode described by $\\vec{v}$. \n",
"\n",
"This notebook will help you construct the matrix $A$ given some geometry, and then you will write a small amount of code to solve the problem $A \\vec{v} = \\lambda \\vec{v}$ for $\\lambda$ and $\\vec{v}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%pylab inline\n",
"import numpy as np\n",
"from matplotlib import path"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define Some Helper Functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You will need to make edits to two functions below: **construct_1D_FDE** and **construct_2DSquare_FDE**. \n",
"\n",
"**construct_1D_FDE(l, N)**: This function should take in two variables (l, the length of a string; N, the number of points on the string to model, including the anchor points) and output a matrix, $A$, which describes the 3-point finite difference model of the vibration of the string. **$A$ should be ~~NXN~~ $(N-2) \\times (N-2)$.**\n",
"\n",
"Reminder: the 3-point difference formula is $$\\frac{d^2 u}{dx^2} \\approx \\frac{u(x+h)-2u(x)+u(x-h)}{h^2}$$\n",
"\n",
"**construct_2DSquare_FDE(l, N)**: This function should take in two variables (l, the side-length of a square membrane; N, the number of points on one side of a membrane to model, including the anchor points) and output a matrix, $A$, which describes the 5-point finite difference model of the vibration of the membrane. **$A$ should be ~~N^2XN^2~~ $(N-2)^2\\times (N-2)^2$.**\n",
"\n",
"Reminder: the 5-point difference formula is $$ \\nabla^2 u(x,y) = \\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2} \\approx \\frac{u(x+h,y) + u(x,y+h)-4u(x,y) + u(x,y-h)+ u(x-h,y)}{h^2}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def construct_1D_FDE(l, N):\n",
" # l = length of a string\n",
" # N = number of points on a string\n",
" ######## STUDENT: write code to generate matrix, A\n",
"\n",
" ######## END STUDENT EDITS\n",
" return A;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def construct_2DSquare_FDE(l,N):\n",
" # l = sidelength of a square\n",
" # N = number of points on a side\n",
" ######## STUDENT: write code to generate matrix, A\n",
"\n",
" ######## END STUDENT EDITS\n",
" \n",
" ######## Do not edit the section below\n",
" G = arange((N-2)*(N-2))+1;\n",
" G = np.reshape(G,(N-2,N-2)).T;\n",
" G = np.c_[zeros((N-2,1)),G,zeros((N-2,1))]\n",
" G = np.r_[zeros((1,N)),G,zeros((1,N))]\n",
" ######## Do not edit the section above\n",
"\n",
" return [A,G]"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"**The helper functions numgrid and delsq do not need to be edited.** They will be used to automatically generate the $A$ matrix for more arbitrary geometries than strings or squares. They are adapted from MATLAB developer Cleve Moler."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def delsq(G):\n",
" # Do not edit.\n",
" \"\"\"\n",
" DELSQ Construct five-point finite difference Laplacian.\n",
" delsq(G) is the sparse form of the two-dimensional,\n",
" 5-point discrete negative Laplacian on the grid G.\n",
" adapted from C. Moler, 7-16-91.\n",
" Copyright (c) 1984-94 by The MathWorks, Inc.\n",
" \"\"\"\n",
" [m,n] = G.shape\n",
" # Indices of interior points\n",
" G1 = G.flatten()\n",
" p = np.where(G1)[0]\n",
" N = len(p)\n",
" # Connect interior points to themselves with 4's.\n",
" i = G1[p]-1\n",
" j = G1[p]-1\n",
" s = 4*np.ones(p.shape)\n",
"\n",
" # for k = north, east, south, west\n",
" for k in [-1, m, 1, -m]:\n",
" # Possible neighbors in k-th direction\n",
" Q = G1[p+k]\n",
" # Index of points with interior neighbors\n",
" q = np.where(Q)[0]\n",
" # Connect interior points to neighbors with -1's.\n",
" i = np.concatenate([i, G1[p[q]]-1])\n",
" j = np.concatenate([j,Q[q]-1])\n",
" s = np.concatenate([s,-np.ones(q.shape)])\n",
" # sparse matrix with 5 diagonals\n",
" A = zeros((N,N));\n",
" for ind in range(0,i.shape[0]-1):\n",
" A[i[ind],j[ind]] = s[ind];\n",
" return A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**The helper functions plotDrumMode and points_in_drum do not need to be edited.** They will be used to visualize the vibrational modes of a membrane once you've solved the eigenvalue problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def plotDrumMode(V,modeNum,G,xx,yy):\n",
" # Do not edit.\n",
" numberOfPoints_x = xx.shape[0];\n",
" numberOfPoints_y = yy.shape[0];\n",
" V_n = V[:,modeNum];\n",
" a_n = zeros_like(xx);\n",
" for i in range(0,numberOfPoints_x-1):\n",
" for j in range(0,numberOfPoints_y-1):\n",
" V_ind = G[i,j]-1;\n",
" if (V_ind >= 0)&(V_ind < V_n.shape[0]):\n",
" a_n[i,j] = V_n[V_ind]\n",
" else:\n",
" a_n[i,j] = 0;\n",
" plt.figure(figsize=(5,5))\n",
" CS = plt.contour(xx, yy, a_n)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def points_in_drum(xx,yy,drumPath):\n",
" # Do not edit.\n",
" h = xx[0,1]-xx[0,0];\n",
" positions = np.vstack([xx.ravel(), yy.ravel()])\n",
" positionBooleanIn = drumPath.contains_points(positions.T,transform=None,radius=-0.000000000000001)\n",
" positionBooleanOnIn = drumPath.contains_points(positions.T,transform=None,radius=0.000000000000001)\n",
" pointsInPolygon = positions.T[positionBooleanIn]/h;\n",
" pointsOnPolygon = positions.T[positionBooleanOnIn^positionBooleanIn]/h;\n",
" G = np.zeros(xx.shape,dtype=np.int)\n",
" for i in range(pointsInPolygon.shape[0]):\n",
" G[pointsInPolygon[i,0],pointsInPolygon[i,1]] = i+1;\n",
" \n",
" return [pointsInPolygon,pointsOnPolygon,G]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def construct_2DPolygon_FDE(gridDensity,gridLength,drum_path):\n",
" # Do not edit.\n",
" N = gridDensity*gridLength;\n",
" h = 1.0/gridDensity;\n",
" x = linspace(0,gridLength,N+1);\n",
" xx,yy = meshgrid(x,x);\n",
" [pointsInPolygon,pointsOnPolygon,G] = points_in_drum(xx,yy,drum_path);\n",
" A_drum = delsq(G)/(h**2)\n",
" return [A_drum,G]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part a-d)\n",
"Use the construct_1D_FDE helper function to generate the matrix $A$ for a string length of 1 and 50 model points. Then use an eigenvalue solver to find the eigenvalues and eigenvectors for $A$. (You can use functions built into the linalg library to do this. I suggest the eigh function.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"stringLength = 1.0; # play with this value\n",
"numberOfPoints = 50; # play with this value\n",
"h = stringLength/(numberOfPoints-1);\n",
"x = arange(numberOfPoints)*h;\n",
"\n",
"A = construct_1D_FDE(stringLength,numberOfPoints);\n",
"# hint: if you implemented this code correctly, when stringLength=1.0 and numberOfPoints=5,\n",
"# you should get the 3x3 matrix that part a) asks for."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Solution to the eigenvalue problem:\n",
"##### Student utilize solver here.\n",
"[evals,evecs] = ;\n",
"# evecs = matrix whose columns are the eigenvectors of A\n",
"# evals = vector whose columns are the eigenvalues of A corresponding to the columns of evecs"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot the first and last eigenvectors\n",
"first_evec = evecs[:,0]\n",
"last_evec = evecs[:,-1]\n",
"first_eval = evals[0]\n",
"last_eval = evals[-1]\n",
"\n",
"x = arange(numberOfPoints)*h;\n",
"\n",
"plt.figure(figsize=(7,7))\n",
"plt.plot(x,np.r_[0,first_evec,0],'r-o');\n",
"plt.plot(x,np.r_[0,last_evec,0],'b-o');"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false,
"scrolled": false
},
"source": [
"## Part g)\n",
"Use the construct_2DSquare_FDE helper function to generate the matrix $A$ for a square membrane with side-length of 1 and 50 points along a side. Then use an eigenvalue solver to find the eigenvalues and eigenvectors for $A$. (Use the same eigenvalue solver you used above.) There is a little extra code to generate a matrix, G, which will be used to plot the results. You don't need to modify this code to get your solution working. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [],
"source": [
"sidelength = 1.0; # play with this value\n",
"numberOfPoints = 50; # play with this value\n",
"\n",
"x = linspace(0,sidelength,numberOfPoints) # Do not edit\n",
"[xx,yy] = meshgrid(x,x); # Do not edit\n",
"\n",
"[A_squareDrum,G] = construct_2DSquare_FDE(sidelength,numberOfPoints); # calls the helper function you defined above.\n",
"\n",
"######## Student: implement eigen-solution here to find eigen values of A_squareDrum\n",
"[D,V] = ;\n",
"# V = matrix whose columns are the eigenvectors of A_squareDrum\n",
"# D = vector whose columns are the eigenvalues of A_squareDrum corresponding to the columns of V"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plotDrumMode function takes your eigenvectors (formatted as column vectors; if you use [D,V] = linalg.eigh(A_squareDrum), you can pass V), a number corresponding to the mode you want to plot, and the variables defined in the \"do not edit\" section (G, xx, and yy). Plot the zero-th and first modes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plotDrumMode(V,0,G,xx,yy)\n",
"plotDrumMode(V,1,G,xx,yy)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Parts h-i)\n",
"Here are two polygon shapes that we will study, drum1 and drum1. The variables gridDensity and gridLength describe the density of model points and the side-length of the square model grid. You can modify these values to get higher spatial resolution results, but remember that this trades off with the amount of memory and time the code needs to run!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"drum1_path = path.Path([(0,0), (1,0), (3,2), (2,2),\n",
" (2,3),(1,2),(1,1),(0,1),(0,0)])\n",
"drum2_path = path.Path([(0,1), (1,0), (2,0), (2,2),\n",
" (3,2),(2,3),(1,2),(1,1),(0,1)])\n",
"\n",
"gridDensity = 5; # increase this to change the number of points in the model. More points = more time and memory, so be careful.\n",
"gridLength = 3.0;\n",
"\n",
"\n",
"[A_weirdDrum1,G1] = construct_2DPolygon_FDE(gridDensity,gridLength,drum1_path);\n",
"[A_weirdDrum2,G2] = construct_2DPolygon_FDE(gridDensity,gridLength,drum2_path);\n",
"\n",
"[D1,V1] = linalg.eigh(A_weirdDrum1);\n",
"[D2,V2] = linalg.eigh(A_weirdDrum2);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# defining drum1 and drum2 for easy plotting of the drum shape.\n",
"drum1 = np.array([[0, 0, 2, 2, 3, 2, 1, 1, 0],\n",
" [0, 1, 3, 2, 2, 1, 1, 0, 0]]);\n",
"drum2 = np.array([[1, 0, 0, 2, 2, 3, 2, 1, 1],\n",
" [0, 1, 2, 2, 3, 2, 1, 1, 0]]);\n",
"N = gridDensity*gridLength;\n",
"x = linspace(0,gridLength,N+1);\n",
"xx,yy = meshgrid(x,x);\n",
"\n",
"# plot a drum mode\n",
"modeNum = 0; # play with this value to see different vibrational modes.\n",
"plotDrumMode(V1,modeNum,G1,xx,yy)\n",
"\n",
"# plot the outline of the drum\n",
"plt.plot(drum1[0,:],drum1[1,:],'b')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# plot the drum mode\n",
"plotDrumMode(V2,modenum,G2,xx,yy)\n",
"\n",
"# plot the outline of the drum\n",
"plt.plot(drum2[0,:],drum2[1,:],'b')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare the eigenvalues for the modes of the two drum shapes. These correspond to the drum pitches, or frequencies. Do the drums sound the same according to your simulation? Why or why not?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"D1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"D2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Final student answer: "
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"For fun, you can go back and edit the drum shape paths to create differently-shaped membranes. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [Root]",
"language": "python",
"name": "Python [Root]"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}