# 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{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}$$$$\begin{split} u(x,t) = \sum_i R(x)\varphi_i(x) \end{split}$$$$\begin{split} \int_0^1 R(x)\varphi_j(x) dx = 0 \end{split}$$$$\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}$$

And we get $$\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$$

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 [ ]: