$$ \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}} $$ Exam INF5620, 2013

Exam INF5620, 2013

Hans Petter Langtangen

Dec 16, 2013

There is also a PDF version of this document.

About the exam. Six problems are given for this exam. For each problem, the candidate must prepare a 15 min oral presentation. Try to communicate a good overview and understanding of the topic, but compose the talk so that you can demonstrate knowledge about details too. Some of the problems require graphical illustrations - you can either sketch graphics on the whiteboard or bring your laptop or iPad to show or compute the graphics. Otherwise there are no aids besides a whiteboard. You will have a printout of this document with the exam problems available in the exam room. (Experience with this type of exam and various aids tells that learning the content by heart gives by far the best delivery and communication of understanding.)

We will throw a die and the number of eyes determines the problem to be presented. Thereafter, you will be given some questions, either about parts of your presention or facts from the other problems. After each presentation, the next candidate can throw the die and get about 10 min to collect the thoughts before presenting the assigned problem.

Usually some candidates decide not to show up on the exam so everybody must meet 0845 in the morning on the day they want to take the exam to get their specific time for the exam that day (0900, 0930, 1000, 1030, and so on). Candidates can choose their times in the sequence they appear in the list of candidates.

Problem 1: Numerical artifacts in time integration schemes

a) Consider the ODE problem

$$ T' = -kT,\quad T(0)=T_0.$$

Sketch graphically what kind of numerical problems (artifacts, non-physical features) that can arise from applying the Forward Euler, Backward Euler, and Crank-Nicolson schemes to solve this ODE problem. Explain how mathematical analysis can provide understanding of the numerical problems.

b) We consider a diffusion problem

$$ T_{t} = kT_{xx},$$ modeling the temperature in a solid body. Two pieces of the same material, with different temperatures, are brought together at \( t=0 \). The condition at \( t=0 \) can be formulated as \( T(x,0)=T_0 \) for \( x\in [0,L/2) \) and \( T(x,0)=T_1\neq T_0 \) for \( x\in [L/2,L] \). The PDE will then predict how the initially discontinuous temperature develops in time and space.

Illustrate what kind of numerical problems (artifacts, non-physical features) that may arise from the Forward Euler, Backward Euler, and Crank-Nicolson schemes applied to this PDE problem (with finite difference discretization in space). Explain how mathematical analysis can provide understanding of the numerical problems. Point out what is similar to the ODE problem in a) and what is new in this PDE problem.

Hint. A demo program for experimentation with an initial discontinuity is available: demo_osc.py.

Problem 2: 2D/3D wave equation with finite differences

a) Set up a wave equation problem in 2D with zero normal derivative as boundary condition. Assume a variable wave velocity.

Mention a physical problem where this mathematical model arises. Explain the physical interpretation of the unknown function.

b) Present a finite difference discretization. Explain in particular how the boundary conditions and the initial conditions are incorporated in the scheme. You can choose appropriate initial conditions.

c) Explain (in princple) how the 2D discretization can be extended to 3D.

d) Set up the stability condition in 3D. Also quote results about the accuracy of the method in 3D and define the accuracy measure(s) precisely.

Hint. The simplest accuracy measure in 3D is the truncation error. More information about the accuracy arises from the numerical dispersion relation: \( \tilde\omega(k,\Delta x, \Delta y, \Delta z, \Delta t) \). It is enough to define what \( \tilde\omega \) is and outline the structure of the formula for \( \tilde\omega \).

e) Explain how you can verify the implementation of the method.

f) The scheme for the wave equation is perfect for parallel computing. Why? What are the principal ideas behind a parallel version of the scheme?

Problem 3: Analysis of wave problems

a) Consider a vibration problem

$$ u''(t) + \omega^2 u(t) = 0,\quad u(0)=I,\ u'(0)=V\tp $$ Here, \( \omega \), \( I \), and \( V \) are given data. We discretize the ODE by centered diferences,

$$ [D_tD_t u + \omega^2 u =0 ]^n\tp $$ Explain how the stability and the accuracy of this scheme can be investigated via exact solutions of the discrete equations, and quote the main results. Illustrate the numerical problems that can arise from this scheme.

b) We now consider a 1D wave equation

$$ u_{tt} = c^2u_{xx},$$ with some appropriate boundary and initial conditions. Explain how the stability and accuracy of a centered difference scheme,

$$ [D_tD_t u = c^2 D_xD_x u]^n_i $$ can be investigated via exact solutions of the discrete equations. Quote the main results.

c) Explain how the analysis can help us to understand why a smooth initial condition gives relatively small numerical artifacts, and why a less smooth initial condition gives rise to significant numerical artifacts. The movies below show a wave propagating with unit Courant number in a medium and the wave enters another medium with 1/4 of the wave velocity (implying a Courant number of 1/4 in that medium). The propagation of waves in the left medium is exact, while the propagation in the other medium is subject to numerical errors.

Little noise.

Significant noise.

d) Explain how a truncation error analysis is carried out for the problem in a). Find correction terms such that the accuracy of the scheme becomes \( \Oof{\Delta t^4} \).

e) Explain how a truncation error analysis is carried out for the problem in b). Find correction terms such that the accuracy of the scheme becomes \( \Oof{\Delta t^4,\Delta x^4} \).

Problem 4: Finite elements for a 1D wave equation

We consider the 1D wave equation problem on \( \Omega = [0,L] \):

$$ u_{tt} = c^2u_{xx} + f,\quad u(x,0)=I(x),\ u_t(x,0)=0,\ u_x(0)=u_x(L)=0\tp $$

a) Explain how the initial condition can be approximated by the finite element method using the principles of least squares, projection (Galerkin), and interpolation (collocation).

b) Discretize in time by a centered difference: \( u_{tt}(x_i,t_n)\approx [D_tD_tu]^n_i \). Derive a variational formulation of the time-discrete wave equation problem using the Galerkin method. Derive formulas for the element matrix corresponding to the term with \( u_{xx} \) in the PDE.

c) Show how the element matrix associated with the \( u_{xx} \) term is computed for P1 elements. Explain the assembly principle and what the resulting global matrix look like when all cells have equal length.

d) Set up the discrete equations for this wave problem on operator form (assume P1 elements). Briefly explain the idea of an analysis of the scheme based on exact solution of the discrete equations. State the main results. Compare the main results with those of the finite difference method (Problem 3: Analysis of wave problems).

Problem 5: Nonlinear diffusion

We look at the PDE problem

$$ \begin{align*} u_t &= \nabla\cdot ((1+\dfc_0u^4)\nabla u),\quad &\x\in\Omega,\ t >0\\ u(\x,0) &= I(\x),\quad &\x\in\Omega\\ u(\x,t) &= g(\x),\quad &\x\in\partial\Omega_D,\ t>0\\ -(1+\dfc_0u^4)\frac{\partial}{\partial n} u(\x, t) &= h(u - T_s),\quad &\x\in\partial\Omega_R \end{align*} $$ Here, \( u(\x,t) \) is the temperature in some solid material, \( I \) is the initial temperature, \( g \) is the controlled temperature at a part \( \partial\Omega_D \) of the boundary, while at the rest of the boundary, \( \partial\Omega_R \), we apply Newton's cooling law with \( h \) as a heat transfer coefficient and \( T_s \) as the temperature in the surrounding air.

a) Perform a Crank-Nicolson time discretization and define a Picard iteration method for the resulting spatial problems. (Do not pay attention to the boundary conditions.)

b) Perform a Backward-Euler time discretization and derive the variational form for the resulting spatial problems. Use a Picard iteration method to linearize the variational form.

c) Apply Newton's method to the nonlinear variational form \( F=0 \) in b). Set up expressions for the right-hand side (\( -F \)) and the coefficient matrix (Jacobian of \( F \)) in the linear system that must be solved in each Newton iteration.

Problem 6: Finite element calculations with P2 elements

We address the problem

$$ -u''(x) = 2,\quad x\in (0,1),\quad u(0)=\beta,\ u(1)=\gamma,$$ and seek a numerical solution \( u \) in some vector space \( V \) (modulo boundary conditions) with basis \( \{\psi_i\}_{i=0}^N \).

a) Explain the principles of the least squares method, the Galerkin method, and the collocation method. Describe a method to incorporate the boundary conditions.

b) Let \( V=\hbox{span}\{\sin\pi x\} \), and compute the solution corresponding to the least squares method, the Galerkin method, and the collocation method (the latter with \( x=0.5 \) as collocation point). Set \( \beta = \gamma =2 \) for simplicity. What are the errors in each of the approximate method?

c) Now we want to use P2 elements on a uniform mesh. Explain how to calculate the element matrix and vector for cells in the interior of the mesh (those not affected by boundary conditions) and set up the results. Describe how the element matrix and vector are assembled into the global linear sytem.

d) Use only one element and explain how the boundary conditions affect the element matrix and vector. Why does this numerical solution coincide with the exact one? Two equal-sized P1 elements lead to exact values at the nodes. Sketch these P2 and P1 solutions.