9.1

a

The Galerkin method is as follows: $$ \int R(\vec{x})w_i(\vec{x})d\vec{x}=0 $$

Thus, using the Galerkin method, we have: $$ \int_0^1(\frac{\partial^2 u}{\partial t^2}-v^2\frac{\partial^2u}{\partial x^2}-\gamma\frac{\partial}{\partial t}\frac{\partial^2u}{\partial x^2})\phi_j dx=0) $$ Recalling (eqn. 9.13), we define

$$ u(x,t)=\sum_i a_i(t)\phi_i(x) $$

Plugging in...

$$ \sum_i\int_0^1\frac{d^2a_i}{dt^2}\phi_i\phi_j-v^2a_i\phi_j\frac{d^2\phi_i}{dx^2}-\gamma\frac{\partial}{\partial t}(a_i\phi_j\frac{d^2\phi_i}{dx^2})dx=0 $$$$ \sum_i\int_0^1\frac{d^2a_i}{dt^2}\phi_i\phi_j-v^2a_i\phi_j\frac{d^2\phi_i}{dx^2}-\gamma\frac{da_i}{dt}\phi_j\frac{d^2\phi_i}{dx^2})dx=0 $$

Using integration by parts (per eqns. 9.24 and 9.25...) brings us to:

$$ \sum_i\frac{d^2a_i}{dt^2}\int_0^1\phi_i\phi_jdx+\sum_ia_i[\int_0^1v^2\frac{d\phi_i}{dx}\frac{d\phi_j}{dx}dx-v^2\phi_j\frac{d\phi_i}{dx}|^1_0]+\sum_i\gamma\frac{da_i}{dt}[\int_0^1\frac{d\phi_i}{dx}\frac{d\phi_j}{dx}dx-\phi_j\frac{d\phi_i}{dx}|^1_0] $$

Where the integral in the first summation can be defined as a matrix $\mathbf{A_{ij}}$, the integral in the second summation can be defined as a matrix $\mathbf{B_{ij}}$ and the summation in the third summation as a matrix $\mathbf{C_{ij}}$ such that:

$$ \mathbf{A}\cdot\frac{d^2\vec{a}}{dt^2}+\mathbf{B}\cdot\vec{a}+\mathbf{C}\cdot\gamma\frac{d\vec{a}}{dt} $$

9.1b

import sympy as sym
%config IPCompleter.greedy=True

x, x_i, x_minus, x_plus, h, h_plus, h_minus = sym.symbols("x x_i x_minus x_plus h h_plus h_minus")

phi = sym.symbols('phi', cls=sym.Function)
phi = sym.Piecewise(((x-x_minus)/(x_i-x_minus), (x < x_i) & (x_minus<=x)), ((x_plus-x)/(x_plus-x_i), (x < x_plus) & (x_i<=x)), (0, True))
phi

$\displaystyle \begin{cases} \frac{x - x_{minus}}{x_{i} - x_{minus}} & \text{for}\: x \geq x_{minus} \wedge x < x_{i} \\\frac{- x + x_{plus}}{- x_{i} + x_{plus}} & \text{for}\: x \geq x_{i} \wedge x < x_{plus} \\0 & \text{otherwise} \end{cases}$

Let us evaluate $\mathbf{A_{ij}}=\int_0^1\phi_i\phi_jdx$. Let us consider the matrix diagonal, where $\phi_i=\phi_j$. We can thus break up the integral into four pieces:

$$ \int_0^{x_{i-1}}\phi_i\phi_idx+\int_{x_{i-1}}^{x_i}\phi_i\phi_idx+\int_{x_i}^{x_{i+1}}\phi_i\phi_idx+\int^1_{x_{i+1}}\phi_i\phi_idx $$

By definition, the first and fourth integrals are equal to 0. Using SymPy to evaluate

int1 = sym.integrate(phi*phi,(x,x_minus,x_i))
sym.simplify(int1.subs((x_i-x_minus),h,simultaneous=True))

$\displaystyle \begin{cases} \frac{x_{i}}{3} - \frac{x_{minus}}{3} & \text{for}\: x_{i} > x_{minus} \\\frac{\frac{x_{i}^{3}}{3} + x_{plus} \left(- x_{i}^{2} + x_{i} x_{plus} - x_{plus} \min\left(x_{minus}, \max\left(x_{i}, x_{plus}\right)\right) + \min\left(x_{minus}, \max\left(x_{i}, x_{plus}\right)\right)^{2}\right) - \frac{\min\left(x_{minus}, \max\left(x_{i}, x_{plus}\right)\right)^{3}}{3}}{x_{i}^{2} - 2 x_{i} x_{plus} + x_{plus}^{2}} & \text{otherwise} \end{cases}$

int2= sym.integrate(phi*phi,(x,x_i,x_plus))
sym.simplify(int2.subs((x_plus-x_i),h,simultaneous=True))

$\displaystyle \begin{cases} - \frac{x_{i}}{3} + \frac{x_{plus}}{3} & \text{for}\: x_{i} < x_{plus} \\\frac{- \frac{x_{i}^{3}}{3} + x_{minus} \left(x_{i}^{2} - x_{i} x_{minus} + x_{minus} \min\left(x_{i}, \max\left(x_{minus}, x_{plus}\right)\right) - \min\left(x_{i}, \max\left(x_{minus}, x_{plus}\right)\right)^{2}\right) + \frac{\min\left(x_{i}, \max\left(x_{minus}, x_{plus}\right)\right)^{3}}{3}}{x_{i}^{2} - 2 x_{i} x_{minus} + x_{minus}^{2}} & \text{otherwise} \end{cases}$

By definition, $(x_i-x_{i-1})$ and $(x_{i+1}-x_{i})$ both evaluate to $h$. Thus, we can see for the integration, we have the value of the diagonal along $A_{ii\forall i=j}=\frac{2h}{3}$.

Let us now consider when $\mathbf{A}$ when $i\neq j$. In general, we know $\phi_i(x_j)=0$. However, in the case where $j=i\pm1$, we can see (from figure 9.1) that one might expect some overlap in the equation.

As matrices $\mathbf{B}$ and $\mathbf{C}$ are nearly identical, we can consider the equation:

$$ \int_0^1\frac{d\phi_i}{dx}\frac{d\phi_j}{dx}dx-\phi_j\frac{d\phi_i}{dx}|^1_0 $$

recalling that both terms would be multiplied by $v^2$ (in the case of matrix $\mathbf{B}$) or $\gamma$ (in the case of matrix $\mathbf{C}$.

To begin, per definition: $$ \frac{d\phi_i}{dx}=\pm\frac{1}{h} $$ Under boundary conditions, we also know $\phi_j=0$ at $x=0,1$, so evaluating the second part of the matrix is simple, at it is equal to 0.

Again, considering the special case of $j=1$, we once again break up the integral::

$$ \int_0^{x_{i-1}}B+\int_{x_{i-1}}^{x_i}B+\int_{x_i}^{x_{i+1}}B+\int_{x_{i+1}s}^1B $$

. Once again, we know that, definitionally, the first and fourth terms are qual to 0, and so we only evaluate the middle two.

intb=sym.integrate(phi.diff(x)*phi.diff(x),(x,x_minus,x_i))
sym.simplify(intb.subs((x_i-x_minus),h,simultaneous=True))

$\displaystyle \begin{cases} \frac{x_{i} - x_{minus}}{h^{2}} & \text{for}\: x_{i} > x_{minus} \\\frac{x_{i} - \min\left(x_{minus}, \max\left(x_{i}, x_{plus}\right)\right)}{\left(x_{i} - x_{plus}\right)^{2}} & \text{otherwise} \end{cases}$

Though sympy still wont't quite simplify as much as I'd like, we can see that, this evaluates to $h/h^2$, or 1/h.

intb2=sym.integrate(phi.diff(x)*phi.diff(x),(x,x_i,x_plus))
sym.simplify(intb2.subs((x_plus-x_i),h,simultaneous=True))

$\displaystyle \begin{cases} \frac{- x_{i} + x_{plus}}{h^{2}} & \text{for}\: x_{i} < x_{plus} \\\frac{- x_{i} + \min\left(x_{i}, \max\left(x_{minus}, x_{plus}\right)\right)}{\left(x_{i} - x_{minus}\right)^{2}} & \text{otherwise} \end{cases}$

We once again have $1/h$. Adding these together we see that the diagonal entires of matrices $B_{ii}$=$\frac{2}{h}$.

9.1c

We once again define $u = a_0+a_1x+a_2x^2+a_3x^3$

u, udot, a, h = sym.symbols ("u udot a h")
Amatrix = sym.Matrix([[1, 0, 0, 0],[0,1,0,0],[1,h,h*h,h*h*h],[0,1,2*h,3*h*h]])
Uvector1=sym.Matrix([[1],[0],[0],[0]])
Uvector2=sym.Matrix([[0],[1],[0],[0]])
Uvector3=sym.Matrix([[0],[0],[1],[0]])
Uvector4=sym.Matrix([[0],[0],[0],[1]])

Amatrix.inv()*Uvector3

$\displaystyle \left[\begin{matrix}0\\0\\\frac{3}{h^{2}}\\- \frac{2}{h^{3}}\end{matrix}\right]$

For a wave equation with boundary conditions such that $u(0)=0$ and $u(1)=1$, with $\dot u = 0$ at $x=0,1$, we have the above solution, using the third Hermite polynomial.

9.2

a0, a1, a2, a3, E, I, force, L= sym.symbols("a0 a1 a2 a3 E I force L")

u = a0 + a1*x + a2*x*x + a3*x*x*x

We know from the boundary conditions that:

$$ u(0) = 0\\ \dot u(0) = 0 $$

We cannot make assumptions about $u(L)$ or $\dot u(L) $

We then want to integrate: $$ \int_0^L(\frac{1}{2}E*I(\frac{d^2u}{dx^2})-u(x)*f(x)))dx $$

udotdot = u.diff(x,x)

integrated = E*I/2*sym.integrate(udotdot-u*force,(x,0,L))
integrated

$\displaystyle \frac{E I \left(- \frac{L^{4} a_{3} force}{4} - \frac{L^{3} a_{2} force}{3} + L^{2} \left(- \frac{a_{1} force}{2} + 3 a_{3}\right) + L \left(- a_{0} force + 2 a_{2}\right)\right)}{2}$

eqn1 = u.subs(x,0)
eqn2 = u.diff(x).subs(x,0)
eqn3 = u.subs(x,L)
eqn4 = u.diff(x).subs(x,L)
solutions = sym.solve((eqn3,eqn4),(a2,a3))
solutions[a2]

$\displaystyle - \frac{2 L a_{1} + 3 a_{0}}{L^{2}}$

solutions[a3]

$\displaystyle \frac{L a_{1} + 2 a_{0}}{L^{3}}$

Now, because we know the initial conditions, we know that $a_0=0$ and $a_1=0$, so we seem to have come to an impass...