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.


2 Answers 2


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$:

enter image description here

Results for $BDM_1 - P_0$:

enter image description here

Results for $[P_2]^2 - P_1$:

enter image description here

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)

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

  def linf(tau, v, w):
      return - 2* np.pi ** 2 * np.sin(np.pi * w.x[0]) * np.sin(np.pi * w.x[1]) *v

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

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

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

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

errors = np.sqrt(np.array(errors))
plt.loglog(hs, errors, 'bo-')
plt.legend(['rate = {:.2f}'.format(np.polyfit(np.log(hs), np.log(errors), 1)[0])])
  • $\begingroup$ 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. $\endgroup$
    – Cherif
    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.

  • $\begingroup$ 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 𝑏. $\endgroup$
    – Cherif
    May 23 at 15:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.