I'm solving the 1D wave equation
\begin{equation}
\frac{\partial^2 \eta}{\partial t ^2} – \frac{\partial^2 \eta}{\partial x ^2} = 0
\end{equation}
with boundary conditions
\begin{equation}
\frac{\partial \eta}{\partial x} = 0 \qquad \qquad \text{on} \qquad \qquad x = 0, 1
\end{equation}
using Finite Element Method.
Weak formulation, where $w(x)$ is the weight function, yields:
\begin{equation}
\int_{0}^{1} w(x) \frac{\partial^2 \eta}{\partial t ^2} \mathrm{d}x – \int_{0}^{1} w(x) \frac{\partial^2 \eta}{\partial x ^2} \mathrm{d}x = 0
\end{equation}
\begin{equation}
\int_{0}^{1} w \frac{\partial^2 \eta}{\partial t ^2} \mathrm{d}x
– \int_{0}^{1} \frac{\partial }{\partial x }\left( w \frac{\partial \eta}{\partial x }\right) \mathrm{d}x
+ \int_{0}^{1} \frac{\partial w}{\partial x } \frac{\partial \eta}{\partial x } \mathrm{d}x = 0
\end{equation}
the middle integral vanishes through integration and use of B.C.'s
Expressing $\eta$ and $w$ in terms of basis functions:
\begin{align}
\eta = \eta_j (t) \phi_j(x) \\
w = \phi_i(x)
\end{align}
($\phi_i(x)$ are hat functions in my case) gives
\begin{equation}
\left(\int_{0}^{1} \phi_i \phi_j \mathrm{d}x \right) \frac{d^2 \eta_j}{d t ^2}
+ \left(\int_{0}^{1} \frac{d \phi_i}{d x } \frac{d \phi_j}{d x } \mathrm{d}x \right)\eta_j= 0
\end{equation}
which can be written in matrix form as follows
\begin{equation}
M_{ij} \frac{d^2 \eta_j}{d t ^2} = – S_{ij} \eta_j
\end{equation}
TASK: Introduce an auxiliary variable $p_j = \frac{d \eta_j}{d t }$, using Crank-Nicolson time discretization formulate the final system $A \mathbf{x} = \mathbf{b}$.
Introducing the auxiliary variable gives rise to a system of first order equations in time:
\begin{align*}
\frac{d \eta_j}{d t } &= p_j \\
M_{ij} \frac{d p_j}{d t} &= – S_{ij} \eta_j
\end{align*}
Using C-N time discretization:
\begin{align*}
\frac{\eta_j^{n+1} – \eta_j^{n}}{\Delta t } &= \frac{1}{2} (p_j^{n+1}+p_j^{n} )\\
M_{ij} \frac{p_j^{n+1} – p_j^{n}}{\Delta t } &= -\frac{1}{2} S_{ij} \left(\eta_j^{n+1} + \eta_j^{n} \right)
\end{align*}
Here's where I got stuck. I tried eliminating $\eta_j^{n+1}$ and $p_j^{n+1}$ from right hand sides of above equations and obtained:
\begin{align}
\left[M_{ij} + \left(\frac{\Delta t}{2} \right)^2 S_{ij} \right] p_j^{n+1} &= \left[M_{ij} – \left(\frac{\Delta t}{2} \right)^2 S_{ij} \right] p_j^{n} – \Delta t S_{ij} \eta_j^{n} \\
\left[M_{ij} + \left(\frac{\Delta t}{2} \right)^2 S_{ij} \right] \eta_j^{n+1} &= \left[M_{ij} – \left(\frac{\Delta t}{2} \right)^2 S_{ij} \right] \eta_j^{n} + \Delta t M_{ij} p_j^{n}
\end{align}
Out of these I can form a single matrix equation, or just implement as they are. However, neither approach has worked – the numerical solution does not make much sense.
Questions: Have missed something? Is my method valid? Does anybody have experience/advice on solving these kind of problems numerically?
Thanks.
Best Answer
I think I found a way...
Rearranging the Crank-Nicolson formulation as follows: \begin{align*} \eta_j^{n+1} - \frac{\Delta t }{2} p_j^{n+1} &= \eta_j^{n} + \frac{\Delta t }{2} p_j^{n} \\ M_{ij} p_j^{n+1} +\frac{\Delta t}{2} S_{ij} \eta_j^{n+1} &= M_{ij} p_j^{n} -\frac{\Delta t}{2} S_{ij} \eta_j^{n} \end{align*}
and solving the matrix form: \begin{equation} \begin{bmatrix} I_{(N+1)\times(N+1)} & - \frac{\Delta t }{2} I_{(N+1)\times(N+1)} \\ \frac{\Delta t }{2} S & M \end{bmatrix} \begin{bmatrix} \mathbf{\eta}^{n+1} \\ \mathbf{p}^{n+1} \end{bmatrix} = \begin{bmatrix} I_{(N+1)\times(N+1)} & \frac{\Delta t }{2} I_{(N+1)\times(N+1)} \\ -\frac{\Delta t }{2} S & M \end{bmatrix} \begin{bmatrix} \mathbf{\eta}^{n} \\ \mathbf{p}^{n} \end{bmatrix} \end{equation} where $I_{(N+1)\times(N+1)}$ is an $(N+1)\times(N+1)$ identity matrix, produces sensible results!