Eigenvalues of beams and plates using Finite Differences.

beam_eigs_FD.py

# -*- coding: utf-8 -*-
"""
Eigenmodes of an Euler-Bernoulli beam with fixed ends [1]_. The stencil
used is a central finite difference of second order [2]_.


References
----------
.. [1] Euler–Bernoulli beam theory. (2015, June 2). In Wikipedia, 
    The Free Encyclopedia. Retrieved 20:11, June 3, 2015, from
    http://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory
 
.. [2] Finite difference coefficient. (2015, April 20). In Wikipedia,
    The Free Encyclopedia. Retrieved 20:13, June 3, 2015, from 
    http://en.wikipedia.org/wiki/Finite_difference_coefficient

@author: Nicolas Guarin-Zapata
"""

from __future__ import division
import numpy as np
from scipy.sparse import diags
from scipy.linalg import eigh
import matplotlib.pyplot as plt
from matplotlib import rcParams
 
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16 
 
 
#%% Computation phase
L = 1
N = 101
x = np.linspace(0, L, N)
dx = x[1] - x[0]
 
# Just use the inner nodes in the calculation
stiff = diags([6., -4., -4., 1, 1], [0,-1, 1, -2, 2], shape=(N-2, N-2))/dx**4
stiff = stiff.toarray()
 
vals, vecs = eigh(stiff)

#%% Plot phase
nvals = 20
nvecs = 4
num = np.linspace(1, nvals, nvals)
plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
plt.plot(num, np.sqrt(vals[0:nvals]), 'o')
plt.xlabel(r"$N$", size=18)
plt.ylabel(r"$\omega\sqrt{\frac{\lambda}{EI}}$", size=18)

plt.subplot(1, 2 ,2)
for k in range(nvecs):
    b = (2*k + 3 )**2 * np.pi**2/4
    vec = vecs[:,k].copy()
    vec = np.append(vec, 0)
    vec = np.insert(vec, 0, 0)
    plt.plot(x, vec, label=r'$n=%i$'%(k+1))

plt.xlabel(r"$x$", size=18)
plt.ylabel(r"$w$", size=18)
plt.legend(ncol=2, framealpha=0.4)
plt.show()

plate_eigs_FD.py

# -*- coding: utf-8 -*-
"""
Eigenmodes of an plate with fixed ends [1]_. The stencil
used is a central finite difference of second order [2]_.


References
----------
.. [1] Vibration of plates. (2014, September 25). In Wikipedia,
    The Free Encyclopedia. Retrieved 20:30, June 3, 2015, from
    http://en.wikipedia.org/wiki/Vibration_of_plates
 
.. [2] Finite difference coefficient. (2015, April 20). In Wikipedia,
    The Free Encyclopedia. Retrieved 20:13, June 3, 2015, from 
    http://en.wikipedia.org/wiki/Finite_difference_coefficient

@author: Nicolas Guarin-Zapata
"""

from __future__ import division
import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
from matplotlib import rcParams
 
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16 
 
 
#%% Setup phase
Lx = 1.
Ly = 1.
Nx = 51
Ny = 51
x = np.linspace(0., Lx, Nx)
y = np.linspace(0., Ly, Ny)
dx = x[1] - x[0]
dy = y[1] - y[0]
dx = 1

# Just use the inner nodes in the calculation
diag_vals = [20.,
             -8., -8., -8., -8.,
             2., 2., 2., 2.,
             1., 1., 1., 1.]
diag_num = [0,
            -1, 1, -Nx, Nx,
            Nx + 1, Nx - 1, -Nx + 1, -Nx - 1,
            -2, 2, -2*Nx, 2*Nx]
stiff = diags(diag_vals, diag_num, shape=(Nx*Ny, Nx*Ny))/dx**4

nvals = 20

#%%Computation phase
vals, vecs = eigsh(stiff, k=nvals, which='SA')

#%% Plot phase
plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
num = np.linspace(1, nvals, nvals)
plt.plot(num, np.sqrt(vals[0:nvals]), 'o')
plt.xlabel(r"$N$", size=18)
plt.ylabel(r"$\omega L^2\sqrt{\frac{\rho\, h}{D}}$", size=18)

#plt.subplot(1, 2, 2)
plt.subplot(2, 4, 3)
position = [3, 4, 7, 8]
for k in range(4):
    plt.subplot(2, 4, position[k])
    vec = vecs[:,k].copy()
    vec.shape = (Nx, Ny)
    plt.contourf(x, y, vec, 11, cmap="RdYlBu")
    plt.xticks([])
    plt.yticks([])

#plt.colorbar()
plt.show()