$$ \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\vexd}[1]{{v_{\small\mbox{e}, #1}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\Std}[1]{\hbox{Std}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \newcommand{\X}{\boldsymbol{X}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\w}{\boldsymbol{w}} \newcommand{\V}{\boldsymbol{V}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\strain}{\boldsymbol{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\iz}{\boldsymbol{i}_z} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifd}{{I_d}} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\sequencei}[1]{\left\{ {#1}_i \right\}_{i\in\If}} \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\psib}{\boldsymbol{\psi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\boldsymbol{x}_{#1}} \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$ Third "compulsory" exercise: 1D wave equation with finite elements

Third "compulsory" exercise: 1D wave equation with finite elements

INF5620

2013

Project 1: 1D wave equation with finite elements

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

Time-dependent finite element discretization, which builds on finite element discretization in space and previous material on approximation of differential equations and approximation of functions.
Analysis of wave equations, which builds on analysis of vibration equations.

The finite element solution is written as

$$ \begin{equation} u^n(x) \equiv u(x,t_n) = \sum_{j=0}^N c_j^n\basphi_j(x), \label{fem:wave:u1expansion} \end{equation} $$ where $ N=N_n $, $ N_n $ being the number of nodes. Number the cells and nodes from left to right. The expansion \eqref{fem:wave:u1expansion} requires the Dirichlet boundary condition at $ x=0 $ to be incorporated by modifying the linear systems. Alternatively, we can seek

$$ \begin{equation} u^n(x) \equiv u(x,t_n) = U_0(t_n)\basphi_0(x) + \sum_{j=0}^N c_j^n\basphi_{j+1}(x)\tp \label{fem:wave:u2expansion} \end{equation} $$

a) Set up equations for the coefficients $ c_j^0 $ associated with the initial condition $ u(x,0)=I(x) $. Assume $ U_0(0)=I(0) $ and $ I_x(L)=0 $. Use two methods: Galerkin and collocation/interpolation.

b) Use a finite difference method in time to formulate a sequence of spatial problems for $ u^n(x) $, $ n=0,1,\ldots,N $. A special problem is needed for $ n=1 $ in order to incorporate the boundary condition with $ u_t(x,t) $.

c) Use a Galerkin method to derive a variational form of each of the spatial problems. Integrate the $ u_{xx} $ term by parts.

d) Derive the linear system to be solved at each time level and express it 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_0^n,c_1^n,\ldots,c^n_{N}) $, $ c_j^n $ being the coefficients in the expansion $ u^n = \sum_{j=0}^{N} c_j\varphi_j(x) $. 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) Assume 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 equation number $ i $ and group the $ c_j $ terms to mimic the differences above. Then substitute $ c_j $ by $ u_j $ (or $ u_{j+1} $, depending on how the Dirichlet is implemented) to 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.

Remarks

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.