Convergence rates (2D Poisson equation)
In this notebook we numerically verify the theoretical converge rates of the finite element discretization of an elliptic problem.
Specifically, for a domain with boundary , we consider the the boundary value problem (BVP):
Here, denotes the part of the boundary where we prescribe Dirichlet boundary conditions, and denotes the part of the boundary where we prescribe Neumann boundary conditions. denotes the unit normal of pointing outside .
The coefficient , , are chosen such that the analytical solution is .
To obtain the weak form, we define the functional spaces and .
Then, the weak formulation of the boundary value problem reads
Find :
Finally, to obtain the finite element discretization we introduce a uniform triangulation (mesh) of the domain and finite dimensional subspace . The space consists of globally continuous functions that are piecewise polynomial with degree on each element of the mesh , that is where denotes the space of polynomial function with degree .
By letting and , the finite element discretization of the BVP reads:
Find such that
Finite element error estimates
Assuming that the analytical solution is regular enough (i.e. ), the following error estimates hold

Energy norm:

norm:
where denote the mesh size.
Step 1. Import modules
from __future__ import print_function, division, absolute_import
import dolfin as dl
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import logging
logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)
dl.set_log_active(False)
Step 2. Finite element solution of the BVP
Here we implement a function solveBVP
that computes the finite element discretizions and solves the discretized problem.
As input, it takes the number of mesh elements n
in each direction and the polynomial degree degree
of the finite element space. As output, it returns the errors in the and energy norm.
def solveBVP(n, degree):
# 1. Mesh and finite element space
mesh = dl.UnitSquareMesh(n,n)
Vh = dl.FunctionSpace(mesh, 'CG', degree)
# 2. Define the right hand side f and g
f = dl.Constant(0.)
g = dl.Expression("DOLFIN_PI*exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])", degree=degree+2)
# 3. Define the Dirichlet boundary condition
u_bc = dl.Expression("sin(DOLFIN_PI*x[0])", degree = degree+2)
def dirichlet_boundary(x,on_boundary):
return (x[1] < dl.DOLFIN_EPS or x[0] < dl.DOLFIN_EPS or x[0] > 1.0  dl.DOLFIN_EPS) and on_boundary
bcs = [dl.DirichletBC(Vh, u_bc, dirichlet_boundary)]
# 4. Define the weak form
uh = dl.TrialFunction(Vh)
vh = dl.TestFunction(Vh)
A_form = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
b_form = f*vh*dl.dx + g*vh*dl.ds
# 5. Assemble and solve the finite element discrete problem
A, b = dl.assemble_system(A_form, b_form, bcs)
uh = dl.Function(Vh)
dl.solve(A, uh.vector(), b, "cg", "petsc_amg")
# 6. Compute error norms
u_ex = dl.Expression("exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])", degree = degree+2, domain=mesh)
err_L2 = np.sqrt( dl.assemble((uhu_ex)**2*dl.dx) )
grad_u_ex = dl.Expression( ("DOLFIN_PI*exp(DOLFIN_PI*x[1])*cos(DOLFIN_PI*x[0])",
"DOLFIN_PI*exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])"), degree = degree+2, domain=mesh )
err_energy = np.sqrt(dl.assemble(dl.inner(dl.grad(uh)  grad_u_ex, dl.grad(uh)  grad_u_ex)*dl.dx))
return err_L2, err_energy
Step 3. Generate the converges plots
The function make_convergence_plot
generates the converges plots.
It takes as input a numpy.array n
that contains a sequence of number of mesh elements and the polynomial degree degree
of the finite element space.
It plots the energy norm of the error (in red) and the norm of the error (in blue) in a loglog scale. The axis is the mesh size .
The slope of the lines in the loglog scale represents the order of converge and dotted lines represents the expected convergence rate.
def make_convergence_plot(n, degree):
errs_L2 = np.zeros(n.shape)
errs_En = np.zeros(n.shape)
for i in np.arange(n.shape[0]):
print(i, ": Solving problem on a mesh with", n[i], " elements.")
eL2, eE = solveBVP(n[i], degree)
errs_L2[i] = eL2
errs_En[i] = eE
h = np.ones(n.shape)/n
plt.loglog(h, errs_L2, '*b', label='L2')
plt.loglog(h, errs_En, '*r', label='Energy')
plt.loglog(h, .7*np.power(h,degree+1)*errs_L2[0]/np.power(h[0],degree+1), 'b', label = 'order {0}'.format(degree+1))
plt.loglog(h, .7*np.power(h,degree)*errs_En[0]/np.power(h[0],degree), 'r', label = 'order {0}'.format(degree))
plt.legend()
plt.grid()
plt.show()
Converges rate of piecewise linear finite element (k=1)
degree = 1
n = np.power(2, np.arange(2,7)) # n = [2^2, 2^3, 2^4 2^5, 2^6]
make_convergence_plot(n, degree)
0 : Solving problem on a mesh with 4 elements.
1 : Solving problem on a mesh with 8 elements.
2 : Solving problem on a mesh with 16 elements.
3 : Solving problem on a mesh with 32 elements.
4 : Solving problem on a mesh with 64 elements.
Converges rate of piecewise quadratic finite element (k=2)
degree = 2
make_convergence_plot(n, degree)
0 : Solving problem on a mesh with 4 elements.
1 : Solving problem on a mesh with 8 elements.
2 : Solving problem on a mesh with 16 elements.
3 : Solving problem on a mesh with 32 elements.
4 : Solving problem on a mesh with 64 elements.
Copyright © 2018, The University of Texas at Austin & University of California, Merced. All Rights reserved. See file COPYRIGHT for details.
This file is part of the hIPPYlib library. For more information and source code availability see https://hippylib.github.io.
hIPPYlib is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License (as published by the Free Software Foundation) version 2.0 dated June 1991.