# Taylor-Hood elements for Darcy's equation

I would like to know if Taylor-Hood elements $$P_2$$-$$P_1$$ form a stable pair for the mixed approximation of Darcy's equation ( or Poisson's equation) with Dirichlet B.C. In the literature I only find them only in the context of the Stokes problem.

I found some evidence that this method is not optimal, see the introductions in the following papers:

I tried doing some experimental numerics to check the convergence rates. I'm working with the mixed formulation \begin{aligned} (\sigma, \tau) + (u, \mathrm{div}\, \tau) & = 0 \\ (\mathrm{div}\, \sigma, v) &=-(f,v) \end{aligned} and $$f=2 \pi^2 \sin(\pi x) \sin(\pi y)$$. Then I'm evaluating $$\| \sigma - \sigma_h \|_{0}$$ for different elements as a function of the mesh parameter $$h$$.

Results for $$RT_0 - P_0$$: Results for $$BDM_1 - P_0$$: Results for $$[P_2]^2 - P_1$$: While you would expect $$\| \sigma - \sigma_h \|_0 = O(h^2)$$ for such basis functions, I'm instead observing $$\| \sigma - \sigma_h \|_0 = O(h^{1.5})$$ which might be caused by the missing (uniform?) compatibility condition. Despite this, I'm observing a quadratic convergence rate for $$\|u - u_h\|_0$$. However, I believe the energy norm is $$\| \sigma - \sigma_h \|_{div} + \| u - u_h\|_0$$ and for the first term I'm observing $$\| \sigma - \sigma_h\|_{div} = O(h^{0.5})$$ for $$[P_2]^2-P_1$$ while for the other elements ($$RT_0$$ and $$BDM_1$$) I'm observing $$O(h)$$.

So I guess it depends on what you mean by "stable". There is inf-sup condition and there is ellipticity on the kernel. To get optimal rates you need these to hold uniformly with constants independent of the mesh. I think this is not the case here since we observe suboptimal rates.

Edit: Here is the source code after pip install scikit-fem[all]==6.0.0:

import numpy as np
import matplotlib.pyplot as plt
from skfem import *
from skfem.helpers import dot, div

hs, errors = [], []

for refs in [3, 4, 5]:
m = MeshTri.init_sqsymmetric().refined(refs)
e = ElementVector(ElementTriP2()) * ElementTriP1()
#e = ElementTriRT0() * ElementTriP0()
#e = ElementTriBDM1() * ElementTriP0()
basis = Basis(m, e)

@BilinearForm
def bilinf(sigma, u, tau, v, w):
return dot(sigma, tau) + div(sigma) * v + div(tau) * u

@LinearForm
def linf(tau, v, w):
return - 2* np.pi ** 2 * np.sin(np.pi * w.x) * np.sin(np.pi * w.x) *v

A = asm(bilinf, basis)
b = asm(linf, basis)
x = solve(A, b)
(sigma, sigmabasis), (u, ubasis) = basis.split(x)

@Functional
def error(w):
return (w['u'] - np.sin(np.pi * w.x) * np.sin(np.pi * w.x)) ** 2

@Functional
def derror(w):
return (div(w['sigma']) + 2 * np.pi ** 2 * np.sin(np.pi * w.x) * np.sin(np.pi * w.x)) ** 2

hs.append(m.param())
#errors.append(error.assemble(ubasis, u=u))
errors.append(derror.assemble(sigmabasis, sigma=sigma))

errors = np.sqrt(np.array(errors))
plt.title(str(e.elems.__class__))
plt.loglog(hs, errors, 'bo-')
plt.legend(['rate = {:.2f}'.format(np.polyfit(np.log(hs), np.log(errors), 1))])

• Thanks a lot for the detailed explanation. Apparently the uniform ellipticity is not required for the well posedness of the discrete problem nor for the convergence. May 23 at 15:09

Yes, they are. I don't remember a reference to the literature, but the proof goes along the same line as for the Stokes case.

• For Stokes problem, the uniform ellipticity of 𝑎 holds naturally for $H^1(\Omega)^2.$ In Poisson problem the $H(div)$ uniform ellipticity holds only on the kernel of 𝑏. May 23 at 15:02