`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
}