There is also a PDF version of this document.
The purpose of this project is to derive and analyze a finite element method for the 1D wave equation $$ u_{tt} = c^2 u_{xx},\quad x\in [0,L],\ t\in (0,T],$$ with boundary and initial conditions $$ u(0,t) = U_0(t),\quad u_x(L,t)=0,\quad u(x,0)=I(x),\quad u_t(x,0)=V(x) \thinspace .$$
The analysis will give insight into how to adjust the default behavior of the finite element method such that its properties are as good as those for the finite difference method for this particular equation. With the necessary adjustments discovered in detailed 1D calculations, one gets a recipe for constructing an accurate finite element method for simulating 2D and 3D waves in complex geometries.
Relevant background material consists of
Let \( \uex(x,t) \) be the exact solution of the initial-boundary value problem. After discretization in time by finite differences, we get a set of time-discrete problems for \( \uex^n(x) \), where \( \uex^n(x)\approx \uex(x,t_n) \). These functions of space are then approximated by finite element expansions $$ \begin{equation} \uex^n(x) \approx u^n(x) = \sum_{j\in\If} c_j^n\basphi_j(x), \label{fem:wave:u1expansion} \end{equation} $$ in case the linear system for the coefficients \( \{c^n_j\}_{j\in\If} \) is modified such that \( c_0^j=U_0(t_n) \). The unknown coefficients now involve \( u \) at all the nodes and consequently \( \If = \{0,\ldots,N\} \).
Alternatively, one can use a boundary function to take care of the Dirichlet condition at \( x=0 \): $$ \begin{equation} u^n(x) \equiv u(x,t_n) = U_0(t_n)\basphi_0(x) + \sum_{j\in\If} c_j^n\basphi_{j+1}(x)\tp \label{fem:wave:u2expansion} \end{equation} $$ Now, \( \If = \{0,\ldots,N-1\} \) and the unknowns \( c_0^n,\ldots,c_{N-1}^n \) are the values of \( u \) at the \( N \) nodes \( x_1,\ldots,x_N \), i.e., all the nodes where we do not have Dirichlet conditions.
a) Set up equations for the coefficients \( c_j^0 \), \( j\in\If \), associated with the initial condition \( u(x,0)=I(x) \). Use two methods: Galerkin and collocation/interpolation.
b) Use a finite difference method in time to formulate a sequence of spatial problems for \( \uex^n(x) \), \( n=0,1,\ldots,N_t \). A special problem is needed for \( n=1 \) in order to incorporate the boundary condition with \( u_t(x,t) \).
c) Apply Galerkin's method to derive a variational form of each of the spatial problems. Integrate the term with the second-order derivative by parts. Express the variational form in terms of \( u^{n+1} \), \( u^n \), and \( u^{n-1} \).
d) Insert \eqref{fem:wave:u1expansion} or \eqref{fem:wave:u2expansion} in the variational form and derive the linear system to be solved at each time level. Express the system in the form $$ Mc^{n+1} = 2Mc^{n} - Mc^{n-1} - \Delta t^2 c^2 Kc^{n},$$ where \( M \) and \( K \) are matrices, \( c^n=\{c_j^n\}_{j\in\If} \). Set up the formulas for the matrix entries \( M_{i,j} \) and \( K_{i,j} \).
Remark. Unless \( M \) is diagonal, a (tridiagonal) linear system must be solved at each time step, whereas the finite difference method leads to an explicit formula for \( u^{n+1}_i \) at each spatial point at a new time level.
e) Use P1 elements and compute element matrices and vectors corresponding to the \( M \) and \( K \) matrices. Assemble the element contributions to a global matrices.
f) Interpret equation number \( i \) in the linear system as a finite difference approximation of \( u_{tt}=c^2u_{xx} \) using the following scheme: $$ \begin{equation} [D_tD_t( u + \frac{1}{6}\Delta x^2 D_xD_x u ) = c^2 D_xD_x u]_{i}^{n} \label{oblig2:fdm:fem:discrete} \thinspace . \end{equation} $$
Hint. Write out an arbitrary equation in the linear system and group the unknown coefficients to mimic the differences above. Then substitute the coefficients by their corresponding \( u \) values, using a notation \( u^n_i \) as the finite element approximation of \( \uex(x_i,t_n) \), and write the finite element equation in the same form as a finite difference scheme.
g) Perform an analysis of the scheme \eqref{oblig2:fdm:fem:discrete} in the same way as is done for the corresponding finite difference scheme in the course notes. That is, investigate a Fourier component \( u^n_p = \exp{(i(kp\Delta x - \tilde \omega n\Delta t))} \). Show that the stability criterion is $$ C \leq \frac{1}{\sqrt{3}},$$ and therefore that any hope for recovering the exact solution for \( C=1 \) becomes impossible.
h) Find the numerical dispersion relation (\( \tilde\omega \) expressed by other parameters) and plot the error in wave velocity \( \tilde c/c \) (\( \tilde c = \tilde \omega/k \), \( c=\omega/c \)) as a function of \( k\Delta x \) for various Courant numbers. Compare with the corresponding plot for the finite difference method for \( u_{tt}=c^2u_{xx} \) (computed with the aid of the program dispersion_relation_1D.py).
i) Use the Trapezoidal rule to compute the \( M \) matrix. Show that the finite element method with P1 elements now produces the same scheme at the interior mesh points as the standard finite difference method for \( u_{tt}=c^2u_{xx} \). Also show that the first and last equations, which are affected by the boundary conditions, also are identical in the two methods.
j) Instead of using the Trapezoidal rule to produce a diagonal \( M \) matrix, one can replace \( M \) by \( \mbox{diag}(Me) \), where \( e=(1,1,\ldots,1) \) and \( \mbox{diag}(w) \) means a diagonal matrix with the vector \( w \) on the diagonal. The interpretation of this approach is that we sum each row in \( M \) and place the sum on the diagonal. Show that this method produces the same results as Trapezoidal integration (in a 1D problem). Also show that if you replace each row in the element matrices associated with \( M \) by the row sum on the diagonal, the same result arises.
Filenames: wave1D_fem.py, wave1D_fem.pdf.
Say we want to solve the 3D wave equation \( u_{tt}=c^2\nabla^2u \) with finite elements and get the same stability as in the finite difference method. We can then compute \( M \) and \( K \) in the usual way and thereafter just replace \( M \) by \( \mbox{diag}(Me) \). When the mass matrix is lumped this way and becomes diagonal, we get the same representation of of the terms \( u^{n-1} \), \( u^n \), and \( u^{n+1} \) as in the finite difference method. The \( K \) matrix also resembles finite differences for P1 elements. With Q1 elements the \( K \) matrix has more nonzero entries per row, but this does not change the numerical properties of the scheme significantly. The key issue is to use a lumped mass matrix.