%matplotlib inline
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

Galerkin expansion of damped wave equation

$$ \text{take the equation for a damped wave}\\ \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})\phi_j = 0 \\ \qquad \text{since:} \quad u(x,t) = \sum_{i} a_i(t) \phi_i(x) $$$$ \text{we can rearrange to get} \sum_{i} \frac{\partial^2 a_i(t)}{\partial t^2} \int_{0}^{1} \phi_i \phi_j dx - \frac{\partial a_i(t)}{\partial t} \gamma \int_{0}^{1} \frac{\partial ^2 \phi_i}{\partial x^2} \phi_j \; dx - a_i(t) v^2 \int_{0}^{1} \frac{\partial ^2 \phi_i}{\partial x^2} \phi_j \; dx $$$$ A\cdot\frac{\partial^2 a_i(t)}{\partial t^2}-\gamma B \cdot\frac{\partial a_i(t)}{\partial t} - v^2 B \cdot a_i(t) $$$$ \text{where} \quad A = \int_{0}^{1} \phi_i \phi_j dx \quad \text{and} \quad B = \int_{0}^{1} \frac{\partial^2 phi_i}{\partial x^2} \phi_j dx $$$$\text{using elements with spacing h we can evaluate these integrals for j=i}\\$$\begin{equation} \phi_i = \left \{ \begin{aligned} &\frac{x-x_{i-1}}{h} && x_{i-1}< x < x_i \\ \\ &\frac{x_{i+1}-x}{h} && x_i < x < x_{i+1} \\ \\ &0 && \text{otherwise} \end{aligned} \right. \end{equation}

$$\text{if} \quad i = j \\ A_{i=j} = \int_{0}^{1} \phi_i^2 dx \quad \ = \int_{x_{i-1}}^{x_i} [\frac{x-x_{i-1}}{h}]^2 dx + \int_{x_i}^{x_{i+1}} [\frac{x_{i+1}-x}{h}]^2 dx = 2*\frac{h^3}{3 h^2} = \frac{2h}{3}\\ \text{if} \quad i \pm j \\ A_{i, i\pm 1} = \int_{0}^{1} \phi_i \phi_{i-1} dx \quad \ = \int_{x_{i-1}}^{x_i} \frac{x-x_{i}}{h} \big(1-\frac{x-x_{i-1}}{h}\big) = \frac{h}{6}\\ \text{this is just 0.25% scaling of the first result} \; A_{i=j}\\ \text{all that is left are the boundaries which are just 0.5% of} \; A_{i=j} \implies A_{0,0} \: and \: A_{n-1,n-1} = \frac{h}{3} $$$$ \text{for B we first use integration by parts noting that} \quad \frac{\partial \phi_i}{\partial x} = \frac{1}{h}\\ \int_0^1 \frac{\partial^2 \phi_i}{\partial x^2} \phi_j = \frac{\partial \phi_i}{\partial x} \phi_j \Big|_0^1 - \int_0^1 \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_j }{\partial x} \\ B_{i=j} = 0 - \int_{0}^{1} \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_j }{\partial x} dx = - \int_{0}^{1} \frac{\partial \phi_i}{\partial x}^2 = - \int_{0}^{1} \frac{1}{h}^2 $$

Plot the hermite basis functions

# solve for h using the 
h = symbols('h')

u0 = Matrix([[1],[0],[0],[0]])
u1 = Matrix([[0],[1],[0],[0]])
u2 = Matrix([[0],[0],[1],[0]])
u3 = Matrix([[0],[0],[0],[1]])
A = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, h, h**2, h**3], [0, 1, 2*h, 3*h**2]])

a0 = A.inv()*u0
a1 = A.inv()*u1
a2 = A.inv()*u2
a3 = A.inv()*u3

A.inv()

Matrix([
[      1,       0,       0,       0],
[      0,       1,       0,       0],
[-3/h**2,    -2/h,  3/h**2,    -1/h],
[ 2/h**3, h**(-2), -2/h**3, h**(-2)]])

# make lines using h = 3
t0 = a0.subs(h,3)
t1 = a1.subs(h,3)
t2 = a2.subs(h,3)
t3 = a3.subs(h,3)

xs = np.arange(0.,606,)/100

u0 = t0[0]+(t0[1]*xs)+(t0[2]*xs**2)+(t0[3]*xs**3)
u1 = t1[0]+(t1[1]*xs)+(t1[2]*xs**2)+(t1[3]*xs**3)
u2 = t2[0]+(t2[1]*xs)+(t2[2]*xs**2)+(t2[3]*xs**3)
u3 = t3[0]+(t3[1]*xs)+(t3[2]*xs**2)+(t3[3]*xs**3)

fig, ((ax1, ax2)) = plt.subplots(nrows=1, ncols=2,)
ax1.set_ylim(-3, 3)
ax1.set_xlim(0, 3)
ax2.set_ylim(-3, 3)
ax2.set_xlim(0, 6)

ax1.set_title("Hermite basis functions")
ax1.plot(xs,u0,color = 'red')
ax1.plot(xs,u1,color = 'orange')
ax1.plot(xs,u2,color = 'blue')
ax1.plot(xs,u3,color = 'black')

ax2.set_title("expanded x axis")
ax2.plot(xs,u0,color = 'red')
ax2.plot(xs,u1,color = 'orange')
ax2.plot(xs,u2,color = 'blue')
ax2.plot(xs,u3,color = 'black')

plt.show

<function matplotlib.pyplot.show>

# get the equations for each basis function with our newly found coefficients
x = symbols('x')
ex0 = t0[0]+(t0[1]*x)+(t0[2]*x**2)+(t0[3]*x**3)
ex1 = t1[0]+(t1[1]*x)+(t1[2]*x**2)+(t1[3]*x**3)
ex2 = t2[0]+(t2[1]*x)+(t2[2]*x**2)+(t2[3]*x**3)
ex3 = t3[0]+(t3[1]*x)+(t3[2]*x**2)+(t3[3]*x**3)