\(\DeclareMathOperator{\opdiv}{div}\)

49. Finite Element Error Analysis

Some lucky accident happens in the analysis for mixed methods for second order problems:

Let \(I_h\) be the Raviart Thomas interpolation operator, and \(P_h\) the \(L_2\)-projection.

We start the error analysis using the discrete \(\inf-\sup\)-stability, as usual for mixed methods:

\[\begin{eqnarray*} \| \sigma_h - I_h \sigma \|_{H(\opdiv)} + \| u_h - P_h u \|_{L_2} & \leq & \sup_{\tau_h, v_h} \frac{B(\sigma_h-I_h\sigma, u_h - P_h u; \tau_h, v_h)}{\|\tau_h\|_{H(\opdiv)} + \|v_h \|_{L_2}} \\ & = & \sup_{\tau_h, v_h} \frac{B(\sigma-I_h\sigma, u - P_h u; \tau_h, v_h)}{\|\tau_h\|_{H(\opdiv)} + \|v_h \|_{L_2}} \end{eqnarray*}\]

Now, the famous commuting diagram comes into play:

\[\begin{split} \DeclareMathOperator{\Div}{div} \begin{array}{ccc} H(\Div) & \stackrel{\opdiv}{\longrightarrow} & L^2 \\[8pt] \downarrow I_h & & \downarrow P_h \\[8pt] RT_k & \stackrel{\Div}{\longrightarrow} & P_k \: \\[3ex] \end{array} \end{split}\]

the interpolation operators commute with the \(\opdiv\)-operator:

\[ \opdiv I_h = P_h \opdiv \]

The term in the numerator above is

\[\begin{eqnarray*} && a(\sigma - I_h \sigma, \tau_h) + b(\sigma - I_h \sigma, q_h) + b(\tau_h, u - P_h u) \\ & = & \int (\sigma - I_h \sigma) \tau_h + \int \opdiv(\sigma - I_h \sigma) q_h + \int \opdiv \tau_h (u-P_h u) \end{eqnarray*}\]

Using the commuting diagram property we get:

\[ \int (\sigma - I_h \sigma) \tau_h + \int \underbrace{(I - P_h) \opdiv \sigma}_{\in V_h^\bot} \;\underbrace{ q_h}_{\in V_h} + \int \underbrace{ \opdiv \tau_h}_{\in V_h} \, \underbrace{ (u-P_h u) }_{\in V_h^\bot} \]

Thanks to orthogonality, the second and third term vanish !

Thus, we get the error estimate

\[ \| \sigma_h - I_h \sigma \|_{H(\opdiv)} + \| u_h - P_h u \|_{L_2} \preceq \sup_{\tau_h} \frac{ \int (\sigma - I_h \sigma) \tau_h} { \| \tau_h \|_{H(\opdiv)}} \leq \| \sigma - I_h \sigma \|_{L_2} \]

Using the triangle inequality for \(\sigma\), we get

\[ \| \sigma_h - \sigma \|_{L_2} + \| u_h - P_h u \|_{L_2} \preceq \| \sigma - I_h \sigma \|_{L_2} \]

The flux error is as good as we can interpolate into the flux space. Since the finite element space for \(u_h\) is of lower order, the error \(u - u_h\) is in general of lower order. But, the filtered error \(\| u_h - P_h u_h \|\) has the better order.

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

order=1
Sigma = HDiv(mesh,order=order)
V = L2(mesh, order=order-1)
X = Sigma * V

sigma,u = X.TrialFunction()
tau,v = X.TestFunction()

a = BilinearForm(X)
a += (sigma*tau + div(sigma)*v + div(tau)*u)*dx
a.Assemble()

f = LinearForm(10*v*dx).Assemble()

gfu = GridFunction(X)
gfu.vec.data = a.mat.Inverse() * f.vec

sol_sigma = gfu.components[0]
sol_u = gfu.components[1]
print ("The flux")
Draw (sol_sigma, mesh, "sigma");

The flux

The scalar part of the solution:

Draw (sol_u, mesh, "u");

49.1. Local post-processing

Since \(\nabla u = \lambda^{-1} \sigma\), we can reconstruct a better approximation \(\widetilde u\) by small, element-wise problems:

\[ \widetilde u = \operatorname{arg}\min_{v_h \in P^{k+1} \atop \int_T v_h = \int_T u_h} \| \lambda \nabla v_h -  \sigma \|_{L_2, \lambda^{-1}}^2 \]

This optimization problems can be written as a mixed variational problem:

Find: \(\widetilde u \in P^{k+1,dc}\) and \(p_h \in P^0\):

\[\begin{split} \begin{array}{ccccll} \sum_T \int_T \lambda \nabla \widetilde u \nabla \widetilde v & + & \int_{\Omega} \widetilde v_h p_h & = & \sum_T \int_T \sigma_h \nabla \widetilde v_h & \forall \, \widetilde v_h \\ \int_{\Omega} \widetilde u_h q_h & & & = & \int_{\Omega} u_h q_h & \forall q_h \end{array} \end{split}\]

This requires to solve decoupled problems on every element, what is cheap.

V2 = L2(mesh, order=order+1)
Q2 = L2(mesh, order=0)
X2 = V2*Q2

u2,p2 = X2.TrialFunction()
v2,q2 = X2.TestFunction()
a2 = BilinearForm(X2)
a2 += (grad(u2)*grad(v2) + u2*q2 + v2*p2)*dx
f2 = LinearForm(X2)
f2 += (sol_sigma*grad(v2) + sol_u * q2)*dx

a2.Assemble()
f2.Assemble()
gfu2 = GridFunction(X2)
gfu2.vec.data = a2.mat.Inverse() * f2.vec

Draw(gfu2.components[0], mesh, "upost");