Week 5 - Finite Elements

Problem 9.1:

Consider the damped wave equation $ \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2} $. Take the solution domain to be interval $[0,1]$.

(a) Use the Galerkin method to find an approximating system of differential equations.

\begin{equation} \begin{split} R(x,t) = \frac{\partial^2 u}{\partial t^2} - v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2} \end{split} \end{equation}\begin{equation} \begin{split} u(x,t) = \sum_i R(x)\varphi_i(x) \end{split} \end{equation}\begin{equation} \begin{split} \int_0^1 R(x)\varphi_j(x) dx = 0 \end{split} \end{equation}\begin{equation} \begin{split} \int_o^1 R(x)\varphi_j(x) dx &= \int_0^1 (\frac{\partial^2 u}{\partial t^2} - v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2}) \varphi_j(x) dx \\ &= \sum_i \int_0^1 (\frac{\partial^2 a_i}{\partial t^2}\varphi_i - a_i v^2 \frac{\partial^2 \varphi_i}{\partial x^2} + \frac{\partial a_i}{\partial t}\gamma \frac{\partial^2 \varphi_i}{\partial x^2}) \varphi_j(x) dx \\ &= \sum_i [\frac{\partial^2 a_i}{\partial t^2} \int_0^1 \varphi_i \varphi_j(x) dx + (a_i v^2 + \frac{\partial a_i}{\partial t}\gamma)(\int_0^1 \frac{\partial \varphi_i}{\partial x} \frac{\partial \varphi_j}{\partial x}dx - \frac{\partial \varphi_i}{\partial x} \varphi_j \mid_0^1)] \end{split} \end{equation}

And we get \begin{equation} \mathbf{A} \frac{\partial^2 \overrightarrow{a}}{\partial t^2} + \gamma \mathbf{B} \frac{\partial \overrightarrow{a}}{\partial t} + v^2\mathbf{B}\overrightarrow{a} = 0 \end{equation}

where, $$ \mathbf{A}_{ij} = \int_0^1 \varphi_i \varphi_j(x) dx $$ and $$ \mathbf{B}_{ij} = \int_0^1 \frac{\partial \varphi_i}{\partial x} \frac{\partial \varphi_j}{\partial x}dx - \frac{\partial \varphi_i}{\partial x} \varphi_j \mid_0^1 .$$

(b) Evaluate the matrix coefficients for linear hat basis functions, using elements with a fixed size of $h$.

In [15]:
import sympy, math

h = 10**-1
N = math.ceil(1/h) + 1
In [16]:
x = sympy.symbols('x')
phi = [sympy.Piecewise((0, x <= (i-1)*h), (0, x >= (i+1)*h), 
                       ((x - (i-1)*h)/h, x < i*h),
                      (((i+1)*h - x)/h, True)) 
       for i in range(N)]
In [17]:
A = sympy.Matrix([['%.3f' %sympy.integrate(phi_i*phi_j, (x, 0, 1)) 
                   for phi_i in phi] for phi_j in phi])
In [18]:
B = sympy.Matrix([['%.3f' %(sympy.integrate(phi_i.diff(x)*phi_j.diff(x), (x, 0, 1)) - 
                    (phi_i.diff(x)*phi_j).subs(x,1) + (phi_i.diff(x)*phi_j). subs(x,0))
                   for phi_i in phi] for phi_j in phi])
In [19]:
A
Out[19]:
Matrix([
[0.033, 0.083,   0.0,   0.0,   0.0,   0.0,   0.0,    0.0,    0.0,    0.0,    0.0],
[0.083, 0.067, 0.083,   0.0,   0.0,   0.0,   0.0,    0.0,    0.0,    0.0,    0.0],
[  0.0, 0.083, 0.067,  1.35,   0.0,   0.0,   0.0,    0.0,    0.0,    0.0,    0.0],
[  0.0,   0.0,  1.35, 0.067,  1.35,   0.0,   0.0,    0.0,    0.0,    0.0,    0.0],
[  0.0,   0.0,   0.0,  1.35, 0.067, 2.933,   0.0,    0.0,    0.0,    0.0,    0.0],
[  0.0,   0.0,   0.0,   0.0, 2.933, 0.067,   9.0,    0.0,    0.0,    0.0,    0.0],
[  0.0,   0.0,   0.0,   0.0,   0.0,   9.0, 0.067,    9.0,    0.0,    0.0,    0.0],
[  0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   9.0,  0.067, 20.267,    0.0,    0.0],
[  0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0, 20.267,  0.067, 20.267,    0.0],
[  0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,    0.0, 20.267,  0.067, 38.333],
[  0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,    0.0,    0.0, 38.333,  0.033]])
In [20]:
B
Out[20]:
Matrix([
[ 0.0, 10.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,   0.0,   0.0],
[10.0, 20.0, 10.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,   0.0,   0.0],
[ 0.0, 10.0, 20.0, 30.0,  0.0,  0.0,  0.0,  0.0,  0.0,   0.0,   0.0],
[ 0.0,  0.0, 30.0, 20.0, 30.0,  0.0,  0.0,  0.0,  0.0,   0.0,   0.0],
[ 0.0,  0.0,  0.0, 30.0, 20.0, 40.0,  0.0,  0.0,  0.0,   0.0,   0.0],
[ 0.0,  0.0,  0.0,  0.0, 40.0, 20.0, 60.0,  0.0,  0.0,   0.0,   0.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0, 60.0, 20.0, 60.0,  0.0,   0.0,   0.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 60.0, 20.0, 80.0,   0.0,   0.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 80.0, 20.0,  80.0,   0.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 80.0,  20.0, 100.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 100.0,  20.0]])

(c) Now find the matrix coefficients for Hermite polynomial interpolation basis functions, once again using elements with a fixed size of $h$. A symbolic math environment is useful for this problem.

In [ ]:
hermite = [sympy.Piecewise((0, x <= (i-1)*h), (0, x >= (i+1)*h), 
                           ((x - (i-1)*h)/h, x < i*h),
                          (((i+1)*h - x)/h, True)) 
           for i in range(N)]

A_herm = sympy.Matrix([[sympy.integrate(herm_i*herm_j, (x, 0, 1)) 
                       for herm_i in hermite] for herm_j in hermite])

B_herm = sympy.Matrix([[(sympy.integrate(herm_i.diff(x)*herm_j.diff(x), (x, 0, 1)) - 
                        (herm.diff(x)*herm_j).subs(x,1) + (herm_i.diff(x)*herm_j). subs(x,0))
                       for herm_i in hermite] for herm_j in hermite])

Problem 9.2:

Model the bending of a beam (equation 9.29) under an applied load. Use Hermite polynomial interpolation, and boundary conditions fixing the displacement and slope at one end.

In [ ]: