(9.1)

Consider the damped save 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 solutoin domain to be the interval [0, 1].

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

We will substitute $u = \sum\limits_i a_i(t) \varphi_i(x)$ into the residual weighting equation, which gives us:

$0 = \int_0^1 \left(\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}\right) \varphi_j(x)dx$

$ = \sum\limits_i \frac{d^2a_i}{dt^2} \int_0^1 \varphi_i(x) \varphi_j(x) dx - \sum\limits_i \left(v^2a_i + \gamma \frac{da_i}{dt}\right) \int_0^1 \frac{d^2\varphi_i}{dx^2} \varphi_j(x)dx$

We will define the integral $\int_0^1 \varphi_i(x) \varphi_j(x) dx$ as $A_{ij}$.

$\int_0^1 \frac{d^2\varphi_i}{dx^2} \varphi_j(x)dx = \frac{d\varphi_i}{dx}\varphi_j|_0^1 - \int_0^1 \frac{d\varphi_i}{dx} \frac{d\varphi_j}{dx}dx$

We will define the entire above term as $-B_{ij}$.

With these two terms defined, we get:

$0 = \matrix{A} \cdot \frac{d^2\vec{a}}{dt^2} + \gamma \matrix{B} \cdot \frac{d\vec{a}}{dt} + v^2 \matrix{B} \cdot \vec{a}$

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

We need to solve for each element of the $A$ and $B$ matrices. Luckily, many elements of each matrix will be the same.

$A_{i,i} = \int_{x_i-h}^{x_i} \left[\frac{x-(x_i-h)}{h}\right]^2dx + \int_{x_i}^{x_i+h} \left[\frac{x_i+h-x}{h}\right]^2dx = \frac{2h}{3}$

$A_{i-1,i} = A_{i+1,i} = A_{i,i-1} = A_{i,i+1} = \int_{x_i-h}^{x_i} \left[\frac{x_i-x}{h}\right] \left[\frac{x-(x_i-h)}{h}\right]dx = \frac{h}{6}$

$B_{i,i} = \int_{x_i-h}^{x_i} \frac{1}{h^2}dx + \int_{x_i}^{x_i+h} \frac{1}{h^2}dx = \frac{2}{h}$

$B_{i-1,i} = B_{i+1,i} = B_{i,i-1} = B_{i,i+1} = \int_{x_i-h}^{x_i} -\frac{1}{h^2}dx = -\frac{1}{h}$

(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.

Since there are two types of basis functions, with with value $u$ and one with slope $u$ as its expansion coefficient at the element boundary, we will use the notation $A_{u_{i-1},du_i}$ to represent the term for the overlap intergral between the basis function with coefficient $u_{i-1}$ at $x_{i-1}$ and the basis function with coefficient $\dot{u_i}$ at $x_i$. The nontrivial matrix elements are found from the following code:

In [1]:
# From Neil's femat.py code

from sympy import *

x, h = symbols('x h')

M = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, h, h**2, h**3], [0, 1, 2*h, 3*h**2]])
M_inv = M.inv()
v = [1, x, x**2, x**3]
psi = []
dpsi = []

for i in range (4):
    psi.append((M_inv[:,i].dot(v)))
    dpsi.append(diff(psi[i], x))
    
print('A_(u_i,u_i) = ', end='')
print(integrate(psi[2]*psi[2], (x, 0, h)) + integrate(psi[0]*psi[0], (x, 0, h)))
print('A_(du_i,du_i) = ', end='') 
print(integrate(psi[3]*psi[3], (x, 0, h)) + integrate(psi[1]*psi[1], (x, 0, h))) 
print('A_(u_i,du_i) = ', end='') 
print(integrate(psi[2]*psi[3], (x, 0, h)) + integrate(psi[0]*psi[1], (x, 0, h))) 
print('A_(u_(i-1),u_i) = ', end='') 
print(integrate(psi[0]*psi[2], (x, 0, h)))
print('A_(du_(i-1),du_i) = ', end='') 
print(integrate(psi[1]*psi[3], (x, 0, h)))
print('A_(du_(i-1),u_i) = ', end='') 
print(integrate(psi[1]*psi[2], (x, 0, h))) 
print('B_(u_i,u_i) = ', end='') 
print(integrate(dpsi[2]*dpsi[2], (x, 0, h)) + integrate(dpsi[0]*dpsi[0], (x, 0, h))) 
print('B_(du_i,du_i) = ', end='') 
print(integrate(dpsi[3]*dpsi[3], (x, 0, h)) + integrate(dpsi[1]*dpsi[1], (x, 0, h))) 
print('B_(u_i,du_i) = ', end='') 
print(integrate(dpsi[2]*dpsi[3], (x, 0, h)) + integrate(dpsi[0]*dpsi[1], (x, 0, h))) 
print('B_(u_(i-1),u_i) = ', end='') 
print(integrate(dpsi[0]*dpsi[2], (x, 0, h))) 
print('B_(du_(i-1),du_i) = ', end='') 
print(integrate(dpsi[1]*dpsi[3], (x, 0, h))) 
print('B_(du_(i-1),u_i) = ', end='') 
print(integrate(dpsi[1]*dpsi[2], (x, 0, h))) 

A_(u_i,u_i) = 26*h/35
A_(du_i,du_i) = 2*h**3/105
A_(u_i,du_i) = 0
A_(u_(i-1),u_i) = 9*h/70
A_(du_(i-1),du_i) = -h**3/140
A_(du_(i-1),u_i) = 13*h**2/420
B_(u_i,u_i) = 12/(5*h)
B_(du_i,du_i) = 4*h/15
B_(u_i,du_i) = 0
B_(u_(i-1),u_i) = -6/(5*h)
B_(du_(i-1),du_i) = -h/30
B_(du_(i-1),u_i) = -1/10

Cleaning the notation up a bit to be easier to read, we get:

$A_{u_i,u_i} = \frac{26h}{35}, \; A_{du_i,du_i} = \frac{2h^3}{105}, \; A_{u_i,du_i} = 0, \; A_{u_{i-1},u_i} = \frac{9h}{70}, \; A_{du_{i-1},du_i} = -\frac{h^3}{140}, \; A_{du_{i-1},u_i} = \frac{13h^2}{420}$

$B_{u_i,u_i} = \frac{12}{5h}, \; B_{du_i,du_i} = \frac{4h}{15}, \; B_{u_i,du_i} = 0, \; B_{u_{i-1},u_i} = -\frac{6}{5h}, \; B_{du_{i-1},du_i} = -\frac{h}{30}, \; B_{du_{i-1},u_i} = -\frac{1}{10}$

(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, and applying a force at the other end.

In [39]:
# From Neil's beam.py code

import matplotlib.pyplot as plt
from numpy import *
from sympy import *

x, h = symbols('x h')

M = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, h, h**2, h**3], [0, 1, 2*h, 3*h**2]])
M_inv = M.inv()

v = [1, x, x**2, x**3]
psi = zeros(1,4)
dpsi = zeros(1,4)
ddpsi = zeros(1,4)

for i in range(4):
    psi[i] = M_inv[:,i].dot(v)
    dpsi[i] = diff(psi[i], x)
    ddpsi[i] = diff(dpsi[i], x)
    
I = zeros(4, 4)

for i in range(4):
    for j in range(4):
        I[i,j] = integrate(ddpsi[i]*ddpsi[j], (x, 0, h)).subs(h, 1)
        
N = 100
x = linspace(0, 1, N)

y0 = eval(str(psi[0].subs(h, 1)))
y1 = eval(str(psi[1].subs(h, 1)))
y2 = eval(str(psi[2].subs(h, 1)))
y3 = eval(str(psi[3].subs(h, 1)))

plt.title("Functions")
plt.plot(x, y0, 'r', x, y1, 'g', x, y2, 'b', x, y3, 'k')
plt.draw()
In [36]:
y0 = eval(str(dpsi[0].subs(h, 1)))
y1 = eval(str(dpsi[1].subs(h, 1)))
y2 = eval(str(dpsi[2].subs(h, 1)))
y3 = eval(str(dpsi[3].subs(h, 1)))

plt.title("First Derivatives")
plt.plot(x, y0, 'r', x, y1, 'g', x, y2, 'b', x, y3, 'k')
plt.draw()
In [37]:
y0 = eval(str(ddpsi[0].subs(h, 1)))
y1 = eval(str(ddpsi[1].subs(h, 1)))
y2 = eval(str(ddpsi[2].subs(h, 1)))
y3 = eval(str(ddpsi[3].subs(h, 1)))

plt.title("Second Derivatives")
plt.plot(x, y0, 'r', x, y1, 'g', x, y2, 'b', x, y3, 'k')
plt.draw()
In [45]:
A = matrix(zeros(N, N), dtype='float')
A[0, 0] = I[0, 0]
A[0, 1] = I[0, 1]
A[0, 2] = I[0, 2]
A[0, 3] = I[0, 3]
A[1, 0] = I[1, 0]
A[1, 1] = I[1, 1]
A[1, 2] = I[1, 2]
A[1, 3] = I[1, 3]

for i in range (2, N-3, 2):
    A[i, i-2] = I[2, 0]
    A[i, i-1] = I[2, 1]
    A[i, i] = I[2, 2] + I[0, 0]
    A[i, i+1] = I[2, 3] + I[0, 1]
    A[i, i+2] = I[0, 2]
    A[i, i+3] = I[0, 3]
    A[i+1, i-2] = I[3, 0]
    A[i+1, i-1] = I[3, 1]
    A[i+1, i] = I[3, 2] + I[1, 0]
    A[i+1, i+1] = I[3, 3] + I[1, 1]
    A[i+1, i+2] = I[1, 2]
    A[i+1, i+3] = I[1, 3]

A[N-2, N-4] = I[2, 0]
A[N-2, N-3] = I[2, 1]
A[N-2, N-2] = I[2, 2]
A[N-2, N-1] = I[2, 3]
A[N-1, N-4] = I[3, 0]
A[N-1, N-3] = I[3, 1]
A[N-1, N-2] = I[3, 2]
A[N-1, N-1] = I[3, 3]

for i in range(N):
    A[i, 0] = 0
    A[0, i] = 0
    A[i, 1] = 0
    A[1, i] = 0
    
A[0,0] = 1
A[1,1] = 1

B = linalg.inv(A)

b = zeros(N)

for i in linspace(0, 1, 10):
    b[-2] = i
    c = dot(B, b)
    u = c[::2]
    plt.plot(u, 'k')
    
plt.show()