Past exam of the finite element methods for partial differential equations course of the University of Oxford

Original solutions to past examinations of the finite element methods for partial differential equations course of the University of Oxford.

Solutions to all three questions of the finite element methods for partial differential equations course of the University of Oxford Part C examination, sat on 6 June 2023. Source: original released PDF, bundle pages 325–328; paper pages 1–4. The original numbering is retained; question statements are not reproduced.

All function spaces and bilinear forms below are real. Write $(p,q{)}_0={\int }_{\Omega }pq$ for the L2 inner product, and $\Vert p{\Vert }_1^2=\Vert p{\Vert }_0^2+\Vert \nabla p{\Vert }_0^2$ for the first-order Sobolev norm; the notation applies componentwise to vector fields. We use the Poincaré inequality explicitly allowed by the paper. On an arbitrary unbounded Lipschitz domain that inequality would need an additional assumption. The convergence assertions use the usual full conforming finite element spaces, or vector subspaces containing the indicated interpolants; merely being a vector subspace of the stated piecewise polynomial functions is insufficient. Part 1(d) explains the intended meaning of positivity for its nonsymmetric matrix.

For a sufficiently regular solution, put $w=-\Delta u$, with $\Delta$ the Laplacian. Then $-\Delta u-w=0$ and $u-\Delta w=u+{\Delta }^{2}u=f$. The zero Dirichlet boundary conditions on $u$ and $\Delta u$ become zero Dirichlet boundary conditions on $u$ and $w$. Conversely, eliminating $w$ from these two equations recovers the fourth-order equation with the biharmonic operator. For the weak solution we will construct, these interior identities hold as distributional identities; zero Sobolev traces are imposed on $u$ and $w$, without assuming classical fourth derivatives of $u$.

Let $V=H_0^1(\Omega )\times H_0^1(\Omega )$, the Cartesian product of two zero-boundary Sobolev spaces, with normThis is a Hilbert space. Applying integration by parts separately to the two equations gives the weak formulationEquivalently, for $U=(u,w)$ and $T=(v,z)$, define the bilinear form and linear functionalWe seek $U\in V$ with $b(U,T)=F(T)$ for every $T\in V$.

The skew-coupled elliptic system has a useful bounded bilinear form estimate with constant one. Indeed, write the integrand as the dot product of $(\nabla p,\nabla q,-q,p)$ and $(\nabla v,\nabla z,v,z)$. The Cauchy-Schwarz inequality givesFor equal arguments, the two zeroth-order coupling terms cancel. If $K$ is the constant in the supplied Poincaré inequality, thenThus $b$ is a coercive bilinear form, although it is not a symmetric bilinear form. The Lax-Milgram theorem applies and yieldsIn particular, $\Vert (u,w){\Vert }_V\le (1+K)\Vert f{\Vert }_0$. Testing with compactly supported smooth test functions identifies $-\Delta u=w$ and $-\Delta w=f-u$ as distributional identities. This gives the asserted relation to the original boundary value problem at the available regularity.

Take $V_h={\mathcal{U}}_h\times {\mathcal{W}}_h\subset V$. The Galerkin method asks for $U_h=(u_h,w_h)\in V_h$ such thatSince a finite-dimensional subspace is closed, $V_h$ is itself a Hilbert space. The bounded bilinear form and coercive bilinear form estimates in part (a) restrict to $V_h$ with the same constants. The Lax-Milgram theorem therefore gives exactly one $U_h$ for every $h$.

Put $E=U-U_h=(e_u,e_w)$. Subtracting the continuous and discrete weak formulations gives Galerkin orthogonality:For any $Y_h=(v_h,z_h)\in V_h$, the difference $Y_h-U_h$ also belongs to $V_h$, soThe last step uses the constant-one Cauchy-Schwarz inequality estimate from part (a). Bounding the four terms separately by a sum would obscure the precise constant requested here.

The two minimization problems over ${\mathcal{U}}_h$ and ${\mathcal{W}}_h$ are independent. Each has a best approximation because its vector subspace is finite-dimensional. Consequently, ifthen the preceding inequality gives exactlyBoth errors belong to the zero-boundary Sobolev space. Applying the Poincaré inequality componentwise now yieldsWhen $U=U_h$ the desired estimate is automatic; otherwise divide by the nonzero norm. ThusThe constant $C=1+K$ is independent of $h$. This is the Céa lemma estimate for this particular skew-coupled elliptic system.

Use the continuous piecewise-linear nodal interpolants $I_{h}u$ and $I_{h}w$ on the shape-regular mesh. Since $u,w\in H^2(\Omega )$ in two dimensions, the Sobolev embedding theorem supplies continuous representatives, so the nodal values are well-defined. Their zero Sobolev traces imply zero boundary nodal values. Hence $I_{h}u\in {\mathcal{U}}_h$ and $I_{h}w\in {\mathcal{W}}_h$ for the usual full conforming finite element spaces.

The permitted finite element interpolation estimates give constants $C_u,C_w$ determined by the shape regularity and $\Omega$, independent of $h$, such thatSubstitute these choices in the best approximation bound from part (b):For fixed $u,w$, this constant is independent of $h$. A shape-regular mesh suffices; no uniformity of the individual element sizes is needed, since $h$ is their maximum diameter.

The inclusion of these interpolants is an essential assumption. If the wording “consist of” were read as permitting arbitrary smaller vector subspaces, taking ${\mathcal{U}}_h={\mathcal{W}}_h=\{0\}$ for every $h$ would prevent the claimed convergence for a nonzero solution. The full conforming finite element space interpretation is also needed for the dimension assertion in part (d).

Let $x_1,\dots ,x_N$ be the interior nodes, and let ${\varphi }_1,\dots ,{\varphi }_N$ be the nodal basis of the common conforming finite element space, with ${\varphi }_j(x_i)={\delta }_{ij}$ and zero boundary values. These piecewise-linear hat functions give the expansionsTesting the first equation with ${\varphi }_i$ and then the second equation with ${\varphi }_i$ gives the linear systemHere $S$ is the diffusion stiffness matrix and $M$ the mass matrix. Their entries are symmetric in $i,j$, so both are symmetric matrices.

For $c\in {\mathbb{R}}^N$, put $v_c={\sum }_j{c}_j{\varphi }_j$. The two quadratic forms areIf $c\ne 0$, the linear independence of the basis ensures $v_c\ne 0$, so $c^{\mathsf{T}}Mc>0$. If $c^{\mathsf{T}}Sc=0$, the Poincaré inequality and zero Dirichlet boundary condition force $v_c=0$, contradicting $c\ne 0$. Thus both $S$ and $M$ are positive-definite matrices.

For $X=(p,q)\in {\mathbb{R}}^N\times {\mathbb{R}}^N$, symmetry of $M$ cancels the cross terms:ThereforeThis is the positivity of a quadratic form intended by the last sentence of the question. For $N>0$, $A$ is a nonsymmetric matrix, since $M\ne 0$. Under the convention that “positive-definite matrix” includes symmetry, that last sentence needs this qualification: the symmetric part of a matrix is positive definite. The quadratic lower bound for a nonsymmetric matrix still implies that $A$ is an invertible matrix: $AX=0$ would give $X^{\mathsf{T}}AX=0$, hence $X=0$.

Label the polygon vertices of the nondegenerate triangle $K$ as $a_1,a_2,a_3$, and put $m_{ij}=(a_i+a_j)/2$ for its three edge midpoints. The six monomials $1,x_1,x_2,x_1^2,x_1{x}_2,x_2^2$ form a basis of ${\mathcal{P}}_2(K)$, so its dimension is six. Nondegeneracy matters: $K$ has nonempty interior, and a polynomial vanishing on that interior is the zero polynomial.

Let ${\lambda }_1,{\lambda }_2,{\lambda }_3$ be the affine functions of the barycentric coordinate system, characterized by ${\lambda }_i(a_j)={\delta }_{ij}$ and ${\sum }_i{\lambda }_i=1$. Define six quadratic polynomialsAt $a_k$, $N_i(a_k)={\delta }_{ik}$ and all $N_{ij}(a_k)$ vanish. At $m_{kl}$, exactly ${\lambda }_k$ and ${\lambda }_l$ equal $1/2$, while the third is zero. Thus every $N_i$ vanishes at all three edge midpoints, and $N_{ij}(m_{kl})$ is one for the matching unordered pair and zero otherwise.

These six polynomials have the cardinal property for the six point evaluation functionals. Evaluating any vanishing linear combination at the nodes forces each of its coefficients to vanish. They are therefore six linearly independent members of a six-dimensional vector space, hence a basis.

In particular, every $p\in {\mathcal{P}}_2(K)$ has the unique expansionThis proves both existence and uniqueness for arbitrary six prescribed nodal values, which is precisely an unisolvent finite element. It also identifies the finite element degrees of freedom and shape functions of the quadratic Lagrange triangle explicitly.

Use the Sobolev space with a partial Dirichlet conditionwhere $T$ is the Sobolev trace operator. The restricted Sobolev trace operator is continuous into $L^2({\Gamma }_D)$; its kernel is consequently a closed vector subspace of $H^1(\Omega )$. Thus $V$, with the first-order Sobolev norm, is a Hilbert space. The endpoint convention for ${\Gamma }_D$ and ${\Gamma }_N$ does not affect this definition: finitely many boundary points have zero boundary length. It would be inappropriate to impose pointwise boundary values on arbitrary $H^1$ functions in two dimensions.

For a smooth solution, Green's first identity givesThe boundary term vanishes because $Tv=0$ on ${\Gamma }_D$ and the normal derivative is zero on ${\Gamma }_N$. This suggests the weak formulationHere $a$ is exactly the first-order Sobolev inner product. In particular the Cauchy-Schwarz inequality givesThus $a$ is a bounded bilinear form and a coercive bilinear form, both with constant one. The Lax-Milgram theorem proves existence and uniqueness, together withThe positive reaction term already controls the entire Sobolev norm, so no Poincaré inequality is needed here. The zero Neumann boundary condition is the natural part of the weak formulation; it is not an additional restriction on the test space. For this weak solution, $-\Delta u+u=f$ holds as a distributional identity, and the zero normal derivative on ${\Gamma }_N$ is understood through the corresponding Green's first identity.

For a conforming triangular finite element mesh, takeThis is the quadratic conforming finite element space. Across a common edge of a polygon, the restrictions from neighboring triangles are univariate quadratic polynomials. Their agreement at the two endpoints and the midpoint makes their difference a degree-at-most-two polynomial with three distinct roots of a polynomial, so it vanishes on the entire common edge of a polygon. Thus matching the local finite element degrees of freedom produces the required continuity. Similarly, setting the three nodal values on every boundary edge carrying the Dirichlet boundary condition to zero makes its entire restriction vanish.

Continuity across element interfaces means the piecewise gradients give the global weak gradient, with no jump contribution, so $V_h\subset V$. There are finitely many local finite element degrees of freedom, hence $V_h$ is finite-dimensional. The Galerkin method seeksThe bounded bilinear form and coercive bilinear form estimates from part (b) remain valid on $V_h$. Since a finite-dimensional subspace is closed, the Lax-Milgram theorem gives exactly one $u_h$. Alternatively, the stiffness matrix plus the mass matrix is a positive-definite matrix in any basis of $V_h$. The same existence and uniqueness conclusion holds for any finite-dimensional vector subspace of this full space.

Let $e=u-u_h$. Subtracting the two weak formulations gives Galerkin orthogonality, $a(e,v_h)=0$ for every $v_h\in V_h$. Since $a$ is the first-order Sobolev inner product, for every $v_h\in V_h$ we haveIf $e\ne 0$, division proves $\Vert e{\Vert }_1\le \Vert u-v_h{\Vert }_1$; if $e=0$, the same conclusion is immediate. Thus the Céa lemma has constant one in this case.

Use the quadratic nodal interpolant $I_{h}u$. The assumed $H^3(\Omega )$ regularity gives continuous nodal values by the Sobolev embedding theorem, and its zero Sobolev trace makes the nodal values on ${\Gamma }_D$ zero. The quadratic Lagrange triangle construction then gives $I_{h}u\in V_h$. The allowed quadratic finite element interpolation estimate on a shape-regular mesh iswhere $C_I$ depends on the fixed triangle and shape regularity, but not on $h$. ThereforeThis step uses the inclusion $I_{h}u\in V_h$ explicitly permitted by the question. Arbitrary smaller vector subspaces without the same approximation property would not give this rate.

Put $e=u-u_h\in V$. The Aubin–Nitsche duality argument uses the dual weak formulationBy symmetry of $a$ and the Lax-Milgram theorem, this has a unique solution. Taking $v=z$ and using the Cauchy-Schwarz inequality showsWe must justify $H^2$ elliptic regularity for this dual problem from the result actually supplied. That result concerns $\Delta$, whereas the dual equation is $-\Delta z+z=e$. Put $g=z-e\in L^2(\Omega )$. The given result produces $w\in H^2(\Omega )\cap V$ with $\Delta w=g$, zero normal derivative on ${\Gamma }_N$, andBy Green's first identity, for every $v\in V$,Set $v=w-z$. Then $\nabla (w-z)=0$. On the connected triangle, a Sobolev function with zero weak gradient is constant, and its zero Sobolev trace on ${\Gamma }_D$ forces that constant to be zero. Hence $z=w$, provingThis proves the needed reaction-diffusion regularity from Laplacian regularity without assuming the dual solution already has second derivatives in the L2 space.

Only $H^2$ regularity is available for $z$. Use its continuous piecewise-linear interpolant $L_{h}z$, which belongs to the quadratic conforming finite element space $V_h$ because ${\mathcal{P}}_1(K)\subset {\mathcal{P}}_2(K)$. It has the same zero Dirichlet boundary condition and satisfies the finite element interpolation estimateThus no $H^3$ assumption on the dual solution is necessary. Taking $v=e$ in the dual weak formulation and using Galerkin orthogonality givesIf $\Vert e{\Vert }_0=0$, the desired bound holds. Otherwise divide by $\Vert e{\Vert }_0$ and use part (d):All constants depend only on the fixed triangle and shape regularity. The acute-angle assumption enters through the supplied mixed-boundary elliptic regularity result.

The infinitesimal strain tensor $\epsilon (v)=\frac{1}{2}(\nabla v+(\nabla v{)}^{\mathsf{T}})$ and the divergence are linear in $v$. Consequently the given $a$ is a symmetric bilinear form, andFor a matrix-valued function, $\Vert \cdot {\Vert }_0$ here uses the Frobenius inner product pointwise and then the L2 norm. All expressions are finite for $v\in H^1(\Omega ;{\mathbb{R}}^3)$.

Expand the quadratic term at $\theta v+(1-\theta )w$. The linearity of $\ell$ cancels its contribution to the difference, givingFor $0\le \theta \le 1$, this is nonnegative. Thereforeso $J$ is a convex functional on the whole Sobolev space. No boundary condition is needed for this conclusion. The argument uses nonnegativity of the quadratic term; uniqueness on the zero-boundary Sobolev space will instead use the Korn inequality in part (b).

Write $V=H_0^1(\Omega ;{\mathbb{R}}^3)$. For any $v\in V$ and any real $t$, $u+tv$ is admissible. Expanding the functional givesSince the left side is nonnegative for positive and negative $t$, its derivative at zero vanishes. Thus the first variation gives the weak formulation

To identify its differential equation, define the stress tensorThis is a symmetric matrix field. The Frobenius inner product of a symmetric matrix with a skew-symmetric matrix is zero, so $\sigma (u):\nabla v=\sigma (u):\epsilon (v)$. Since $\mathrm{tr}\epsilon (v)=\mathrm{div}v$, we obtainFor smooth compactly supported test functions, integration by parts in each row givesHence $-\mathrm{div}\sigma (u)=f$ as a distributional identity, that is,Membership in $V$ supplies the zero Sobolev trace. This identifies precisely the required weak solution of the linear elasticity problem, with row-wise divergence.

For uniqueness, use the exact form of the Korn inequality supplied in the paper. Put $E=\epsilon (v)$ and $d=\mathrm{tr}E=\mathrm{div}v$. The deviatoric strain is $D=E-\frac{1}{3}dI$. Because $\mathrm{tr}D=0$ and $I:I=3$, orthogonality for the Frobenius inner product givesIf $C_K$ is the supplied Korn inequality constant and $K$ the permitted Poincaré inequality constant, thenConsequently,If $u_1,u_2$ are two weak solutions, their difference $d_u\in V$ satisfies $a(d_u,v)=0$ for every $v\in V$. Choosing $v=d_u$ and using this coercive bilinear form estimate forces $d_u=0$. The minimizer therefore gives the unique weak solution.

For later use, $|\epsilon (v){|}_F\le |\nabla v{|}_F$ and $|\mathrm{div}v|\le \sqrt{3}|\nabla v{|}_F$ also give the bounded bilinear form and linear functional estimatesTogether with the coercive bilinear form estimate and the Lax-Milgram theorem, these would establish existence independently of the assumed minimizer.

The finite-dimensional vector subspace $V_h\subset V$ is closed, so it is a Hilbert space with the inherited Sobolev norm. The bounded bilinear form, coercive bilinear form, and bounded linear functional estimates from part (b) restrict to $V_h$. Applying the Lax-Milgram theorem proves

For any $v_h\in V_h$, put $d_h=v_h-u_h\in V_h$. Expansion using symmetry of the bilinear form givesThus $u_h$ minimizes $J$ on $V_h$, uniquely by the coercive bilinear form estimate. Similarly, the continuous weak formulation gives the exact identityCombining the two identities provesMore precisely, $J(v_h)-J(u_h)=\frac{1}{2}\Vert v_h-u_h{\Vert }_a^2$ and $J(u_h)-J(u)=\frac{1}{2}\Vert u_h-u{\Vert }_a^2$, with the energy norm justified next. These identities connect the energy excess directly to the squared error.

On $V$, the bilinear form $a$ is symmetric, and part (b) proves $a(v,v)\ge \alpha \Vert v{\Vert }_1^2$. Thus $a$ is an inner product: $a(v,v)\ge 0$, and equality forces $v=0$. Its induced energy norm is $\Vert v{\Vert }_a=\sqrt{a(v,v)}$. For any real scalar $t$,The Cauchy-Schwarz inequality for this inner product follows directly by expanding $a(p-tq,p-tq)\ge 0$ and minimizing over $t$ when $q\ne 0$. Hencewhich proves the triangle inequality. This verifies every norm axiom. In fact, the estimates from part (b) give the equivalent norms

Subtract the discrete weak formulation from the continuous one to get Galerkin orthogonality,For any $v_h\in V_h$, write $u-v_h=(u-u_h)+(u_h-v_h)$. The second summand is in $V_h$, so expansion and Galerkin orthogonality give the Pythagorean theorem in the energy norm:The last term is nonnegative and vanishes exactly when $v_h=u_h$. Choosing $v_h=u_h$ attains the lower bound, provingThus $u_h$ is the unique best approximation; the map $u\mapsto u_h$ is the orthogonal projector onto $V_h$ for the linear elasticity energy norm.