11. Encountering wrong solutions#
We study finite elements, from the mathematical perspective, as much to uncover their limitations as to prove their capabilities. Indeed, understanding how these methods fail is essential to ensuring their reliability in practice. In a previous notebook, we examined a finite element method that became overtly unstable. We could easily tell from the blow ups there that it would be impossible for the computed solutions to converge, say in the \(L_2\) norm, as meshes became finer. The failure was immediately noticeable.
In this notebook, we explore a more insidious form of failure. Working with a 2D Maxwell boundary value problem, we present a discretization where the computed solutions appear to converge, yet they converge to the wrong solution. To illustrate this phenomenon, we will intentionally bypass Nédélec elements in favor of a plausible Lagrange discretization and end up with a wrong solution. This exercise also highlights why mathematical innovations like Nédélec spaces are not just theoretical curiosities, but indispensable tools in a practical setting.
11.1. Boundary value problem#
The boundary value problem is to find a 2D vector field \(E\) satisfying
on a domain \(\Omega\) with nonconvex corners. Often \(\operatorname{rot} (\operatorname{curl} E)\) is written as \(\nabla \times (\nabla \times E)\), but we continue to use the notation \(\operatorname{rot}\) and \(\operatorname{curl}\), introduced in an early notebook, and used elsewhere too while discussing waveguides. On the domain boundary, \(t\) denotes a unit tangent. Below, we shall select \(\Omega\) as the familiar L-shape.
The boundary source \(g\) is assumed to have an extension \(E_g\) of enough regularity (it is sufficient that \(E_g\) is in \( H(\operatorname{curl})\) and in \(H(\operatorname{div})\)) such that \(E_g \cdot t = g\) on \(\partial\Omega\), so that the problem of finding \(E\) reduces to the problem of finding \(\mathring{E} = E - E_g\) with homogeneous essential boundary conditions. This is the familiar technique of zeroing boundary conditions at the expense of introducing a nonzero volume source.
As we have seen in the prior notebooks with 3D Maxwell examples and 2D waveguide cross sections, boundary value problems akin to (BVP) above arise in important practical models.
11.2. A provably convergent method#
A finite element method with a mathematically rigorous convergence proof is the standard method using Nédélec approximations for \(E\).
We construct, as usual, a weak formulation treating the divergence equation weakly: find \(E\) in \(H(\operatorname{curl})\), \(E = \mathring{E} + E_g,\) where \(\mathring{E} \in \mathring{H}(\operatorname{curl})\) satisfies, together with \(\phi \in \mathring{H}^1\),
Here \((\cdot, \cdot)\) indicates the inner product in \(L_2(\Omega)\) or its Cartesian products. All terms involving \(E_g\) can be moved to the right hand side to see that (MM) is really a problem for \(\mathring{E} \in \mathring{H}(\operatorname{curl})\) and \(\phi\) in \(\mathring{H}^1\) in a standard Galerkin setting where both the test and trial spaces are equal to \(\mathring{H}(\operatorname{curl}) \times \mathring{H}^1\). For discretizing this mixed method (MM) formulation, we use the Nédélec space \(N_{hp}\) (introduced in an early notebook) for approximating \(E\) and the standard Lagrange finite element space \(V_{hp}\) for approximating \(\phi\). But setting the test function \(v = \nabla \psi\) in the first equation of (MM), we see that the exact solution component \(\phi\) must be zero. This zero function can be viewed as a Lagrange multiplier corresponding to the divergence-zero constraint of \(E\).
By verifying the conditions of the standard Brezzi theory of mixed methods, one can prove that the resulting method is stable and convergent.
11.3. An alternative using a simpler space#
Note that (BVP) appears to control both curl and divergence of \(E\). One is tempted to wonder therefore if we could not just work with the Lagrange space? Is there really a need for a special space like the Nédélec space? Digging in, writing \(E = \mathring{E} + E_g\), we observe that the unknown function \(\mathring{E}\) belongs to the space
This results in the following alternate weak formulation: find \(E = \mathring{E} + E_g\) such that \(\mathring{E} \in X\) satisfies
This curl-and-div (CD) formulation appears to be simpler than the prior mixed method (MM) as it does not involve the auxiliary variable \(\phi\). It also appears to be more natural, as it directly controls the curl and divergence of \(E\) in the same equation.
Of course, we need a finite element subspace of \(X\) if we are to compute with (CD). A piecewise polynomial vector field is in \(H(\operatorname{curl})\) if its tangential component is continuous. It is in \(H(\operatorname{div})\) if its normal component is continuous across element interfaces. Hence, it is in \(X\) if all its components are continuous across element interfaces, so the Lagrange finite element space is a conforming finite element subspace of \(X\). These considerations show that the formulation (CD) appears to offer a compelling avenue to compute the solution of (BVP) using the commonly available Lagrange finite element, the simplest finite element.
Unfortunately, we shall see below that widely different approximations are obtained using the two formulations (MM) and (CD). Both of them converge. But only one of them is right.
11.4. The domain#
The problems we intend to exemplify are only visible in a non-convex domain, so we construct an L-shaped domain. We copy over the geometry construction from a prior notebook, but add a few more names to the edges to be able to set boundary conditions later on. We also set a small maxh attribute for the nonconvex corner point (which is the origin in the constructed geometry) to ensure a smaller elements there.
import ngsolve as ng
from ngsolve.webgui import Draw
from ngsolve import Integrate, BilinearForm, LinearForm, dx, Mesh, HCurl, H1
from ngsolve import GridFunction
from ngsolve import sin, cos, acos, x, y, curl, grad, sqrt, CF
from prettytable import PrettyTable
from math import pi
import numpy as np
from netgen.occ import OCCGeometry, WorkPlane, Axes, X, Y, Z, Glue
def MeshLshape(maxh=1/8, pmaxh=0.001):
s1 = WorkPlane(Axes((0,0,0), n=Z, h=X)).Rectangle(1,1).Face() # right box
s2 = WorkPlane(Axes((0,0,0), n=Z, h=-X)).Rectangle(1,1).Face() # top left box
s3 = WorkPlane(Axes((0,0,0), n=Z, h=Y)).Rectangle(1,1).Face() # bottom box
L = s1+s2+s3 # joins shapes, removing internal interfaces
L.edges[Y<=+1e-6].edges[Y>=-1e-6].edges[X>=-1e-6].name = 'midhorizontal'
L.edges[X>=-1e-6].edges[X<=+1e-6].edges[Y<=+1e-6].name = 'midvertical'
L.edges.Min(X).name = 'left'
L.edges.Max(Y).name = 'top'
L.edges.Min(Y).name = 'bottom'
L.edges.Max(X).name = 'right'
p = L.edges['midvertical'].vertices.Max(Y) # point of singularity
p.maxh = pmaxh
geo = OCCGeometry(L, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=maxh))
return mesh
Recall that the branch cut of the polar angle function \(\theta\) posed problems previously, as explained in a prior notebook. Adopting the solution explained there, we obtain the correct \(\theta\) function of the polar coordinate system.
mesh = MeshLshape(maxh=1/8, pmaxh=0.001)
r = sqrt(x*x + y*y)
theta = acos(x/r)
alpha = -(pi/2 + pi/4)
rotatedx = ng.cos(alpha) * x - ng.sin(alpha) * y
rotatedy = ng.sin(alpha) * x + ng.cos(alpha) * y
theta = ng.atan2(rotatedy, rotatedx) - alpha
Draw(theta, mesh, 'theta');
11.5. Choosing an exact solution#
We intentionally choose an exact solution of low regularity to illustrate the above-mentioned discrepancy between the two methods. The exact solution we shall use, in polar coordinates \((r, \theta)\), is
The following observations are pertinent:
The scalar function \(u = r^{2/3} \sin(2\theta/3)\) solves the Laplace equation (and was an exact solution considered in the prior notebook), i.e., \(\operatorname{div}(\nabla u) = 0\). Hence, \(E = \nabla u\) is divergence free.
Since \(E = \nabla u\) also has zero curl, it satisfies both the differential equations of (BVP).
The boundary source in (BVP), namely \(g\), is nonzero only on the edges of \(\Omega\) that do not meet the origin, since \(u=0\) on those edges.
Although \(E\) has a singularity, the data \(g\) is smooth along \(\partial\Omega\), hence it offers a case where no special integration is required to assemble the correct right hand side.
The components of \(E\) are not in \(H^1(\Omega)\).
Here is a plot of the exact solution. The solution becomes infinite because of the strong singularity at the origin (a feature that is not clearly visible in the finite color scale of the next plot, but will be clearer from the subsequent plots).
r = sqrt(x*x + y*y)
rinv = 1.0/ng.IfPos(r-1e-15, r, 1.e-15) # threshold to avoid 0-division
Eexact = CF( ((2/3)*(cos(theta)*sin(2*theta/3) - sin(theta)*cos(2*theta/3))*pow(rinv,1/3),
(2/3)*(cos(theta)*cos(2*theta/3) + sin(theta)*sin(2*theta/3))*pow(rinv,1/3)))
Draw(Eexact, mesh, vectors={'grid_size': 10}, autoscale=0, min=0, max=1);
The \(y\)-component of \(E\) is plotted below shows how it blows up at the singular point. It is zero on one of the edges leading to the singular point.
Draw(Eexact[1], mesh, autoscale=0, min=0, max=1, deformation=True, euler_angles=[-50, 0, -30]);
The \(x\)-component of \(E\), plotted below, is negative and unbounded at the singular point and vanishes at the other edge to the singular point.
Draw(Eexact[0], mesh, autoscale=0, min=-1, max=0);
Keep this picture in mind as we proceed to approximate this electric field by the two methods.
11.6. The Nédélec approximation#
Here is a straightforward implementation of a finite element method for the formulation (MM). Since our test problem has smooth sources, we do not need to implement any additional tricks for accurate integration etc.
def SolveByNedelec(mesh, bc_E, bc_phi, p=1):
""" Given boundary data for E in "bc_E" and for phi in
"bc_phi", solve by the mixed method approach using Nédélec
elements. (The input vector "bc_E" should such that bc_E . t
gives the needed boundary source g.)
"""
N = HCurl(mesh, type1=True, order=p+1-max(1-p,0), dirichlet='[a-z]*')
L = H1(mesh, order=p+1, dirichlet='[a-z]*') # set order so N includes grad L
Xh = ng.FESpace([N, L])
u, phi = Xh.TrialFunction()
v, psi = Xh.TestFunction()
a = BilinearForm(Xh, symmetric=True)
a += (curl(u)*curl(v) + grad(phi)*v + u*grad(psi)) * dx
f = LinearForm(Xh)
Ephi = GridFunction(Xh, 'E_Nedelec')
with ng.TaskManager():
Ephi.components[0].Set(bc_E, ng.BND)
Ephi.components[1].Set(bc_phi, ng.BND)
a.Assemble()
f.Assemble()
r = f.vec.CreateVector()
r.data = f.vec - a.mat * Ephi.vec
Ephi.vec.data += a.mat.Inverse(Xh.FreeDofs()) * r
return Ephi, Xh
To apply the above routine, we need to provide the boundary condition function as input. This is easily provided since we know the exact solution: while Eexact is the exact E-field visualized previously, the exact phi is \(0\).
Ephi, Xh = SolveByNedelec(mesh, Eexact, CF(0), p=4)
(Ex, Ey), phi = Ephi.components
Draw(Ex, mesh, 'Exh', autoscale=0, min=-1, max=0);
This picture of the \(x\)-component of the approximated \(E\)-field clearly looks very similar to the previously displayed exact \(E\)-component. The method, at least visually, appears to be approximating \(E\) correctly. (Below, we will compute convergence rates to be completely sure.)
11.7. The Lagrange approximation#
Now we use Lagrange finite elements to approximate the other formulation (CD). Here we use two copies of the Lagrange space, the first for the \(x\)-component and the second for the \(y\) component of the \(E\)-field. The tangential boundary condition in (BVP) for \(E\cdot t\) involves one or the other of the \(x\) and \(y\)-components of the \(E\)-field. To impose this boundary condition on the separate components, the individual names of boundary segments are useful, as seen in the code below.
def SolveByLagrange(mesh, bc_E, p=1):
""" Solve using the above-mentioned approach using Lagrange
elements (only) for each component of the electric field,
given boundary data in bc_E (so that bc_E . t is the needed g).
"""
# Make Lagrange spaces so that their product has the
# required tangential boundary conditions:
Vx = H1(mesh, order=p, dirichlet='midhorizontal|top|bottom')
Vy = H1(mesh, order=p, dirichlet='right|left|midvertical')
Xh = ng.FESpace([Vx, Vy])
ux, uy = Xh.TrialFunction()
vx, vy = Xh.TestFunction()
u = CF((ux, uy))
v = CF((vx, vy))
# Define the two differential operations required for the form:
def curl2D(w0, w1):
dw0 = grad(w0)
dw1 = grad(w1)
return dw1[0] - dw0[1]
def div2D(w0, w1):
dw0 = grad(w0)
dw1 = grad(w1)
return dw1[1] + dw0[0]
# System:
a = BilinearForm(Xh, symmetric=True)
a += (curl2D(ux, uy) * curl2D(vx, vy) + \
div2D(ux, uy) * div2D(vx, vy)) * dx
f = LinearForm(Xh)
# Solve:
u = GridFunction(Xh, 'E_Lagrange')
u.components[0].Set(bc_E[0], ng.BND)
u.components[1].Set(bc_E[1], ng.BND)
a.Assemble()
f.Assemble()
r = f.vec.CreateVector()
r.data = f.vec - a.mat * u.vec
u.vec.data += a.mat.Inverse(Xh.FreeDofs()) * r
return u, Xh
E, Xh = SolveByLagrange(mesh, Eexact, p=4)
Here are the \(x\) and \(y\) components of the computed \(E\)-approximation:
Draw(E.components[0], mesh, 'ELagrange')
Draw(E.components[1], mesh, 'ELagrange');
This solution doesn’t seem to look anything like the exact solution. In particular, its values near the non-convex corner indicate that it may have completely missed the singularity. Since looks can be deceiving, we proceed to conduct a convergence study.
11.8. Convergence studies#
It is natural to wonder if the visually inaccurate solution obtained by the Lagrange method is a consequence of having used too coarse a mesh. Does the situation improve on finer meshes? To study this, we perform a convergence study: Start with a (coarse) mesh, and solve the same Maxwell problem on successively refined meshes. Each refinement below is obtained by connecting the midpoints of edges of each triangle in the current mesh (as already seen in a prior notebook). We store all the solutions and meshes in lists, so that we can later compute the convergence rates.
11.8.1. Actual and apparent convergence rates#
There are two types of convergence rates that we can compute from the above process of refinements.
The first is the actual convergence rate, which is the rate at which the error with respect to the exact solution goes to zero. This is only computable when we have access to the exact solution.
The second is the apparent convergence rate, also known by various other terms like estimated convergence rate or numerical order of convergence (NOC), etc., which is obtained by treating the finest mesh solution one can manage to compute as the “ground truth”, and then measuring the rate at which solutions from other meshes approach that on the finest mesh.
Specifically, given solutions \(E_i\) on a sequence of meshes of grid size \(h_i= h_0/2^i\), we can estimate the rate of convergence by examining at what rate
Here \(E_{fine}\) is the numerical solution computed on the finest mesh (maximal \(i\)). The rates so computed gives the NOC (Numerical Order of Convergence) of a method. In contrast, when we know the exact solution and are in a position to compute the actual or exact error \(\| E_i - E_{exact} \|_{L^2(\Omega)} \), we see at what rate
Often, in practical problems with no knowledge of the exact solution, we can only compute the apparent convergence rate. Then the engineering standard is to verify correctness of methods using NOC. If a method is converging to the correct solution, then the apparent convergence rate should match the actual convergence rate. In our test problem, since we know the exact solution, we are able to compute both types of convergence rates.
11.8.2. Rates for the Lagrange approximation#
We compute the solutions on successively refined meshes using the formulation (CD) and then compute both the apparent and actual convergence rates.
def SolveByLagrangeSuccessive(hcoarse=1/4, p=1, nrefinements=5):
"""
Starting with a mesh of grid size "hcoarse", solve using the
Lagrange method on successively refined meshes. Store all solutions,
spaces and meshes and return them in lists.
"""
Es = []
Xs = []
meshes = []
mesh = MeshLshape(maxh=hcoarse, pmaxh=hcoarse)
for ref in range(nrefinements):
meshes.append(ng.Mesh(mesh.ngmesh.Copy()))
uh, Xh = SolveByLagrange(mesh, Eexact, p=p)
Es.append(uh)
Xs.append(Xh)
mesh = meshes[-1]
mesh.ngmesh.Refine()
return Es, Xs, meshes
def convergence_study(Es, Xs, meshes):
"""Given solutions on successively refined meshes, return
|| E_i - E_fine|| and ||E_i - E_exact|| for i-th mesh, for all i.
"""
E_diff = []
E_err = []
fine_E = ng.GridFunction(Xs[-1])
with ng.TaskManager():
fine_E.components[0].Set(Es[-1].components[0])
fine_E.components[1].Set(Es[-1].components[1])
for i in range(len(meshes)-1):
diffE0 = fine_E.components[0] - Es[i].components[0]
diffE1 = fine_E.components[1] - Es[i].components[1]
dE = CF((diffE0, diffE1))
E_diff.append(np.sqrt(Integrate(dE*dE, meshes[i])))
dE = CF((Eexact[0] - Es[i].components[0],
Eexact[1] - Es[i].components[1]))
E_err.append(np.sqrt(Integrate(dE*dE, meshes[i])))
return np.array(E_diff), np.array(E_err)
Es, Xs, meshes = SolveByLagrangeSuccessive(hcoarse=1/4, p=1,
nrefinements=6)
E_diff, E_err = convergence_study(Es, Xs, meshes)
Here E_diff indicates the difference between the solution on the \(i\)-th mesh and the solution on the finest mesh.
E_diff
array([0.05933285, 0.03417458, 0.01882264, 0.00941912, 0.00361138])
These numbers certainly look like they are converging. Let’s make a quick function to print the rate of convergence:
def tabrate(name, dat):
col = ['h', name, 'rate']
t = PrettyTable()
t.add_column(col[0], ['1/'+str(2**(2+i)) for i in range(len(dat))])
t.add_column(col[1], ['%.7f'%e for e in dat])
t.add_column(col[2], ['*'] + \
['%1.2f'%r for r in np.log(dat[:-1]/dat[1:])/np.log(2)])
print(t)
tabrate('NOC: ||Eh-Efine||', E_diff)
+------+-------------------+------+
| h | NOC: ||Eh-Efine|| | rate |
+------+-------------------+------+
| 1/4 | 0.0593328 | * |
| 1/8 | 0.0341746 | 0.80 |
| 1/16 | 0.0188226 | 0.86 |
| 1/32 | 0.0094191 | 1.00 |
| 1/64 | 0.0036114 | 1.38 |
+------+-------------------+------+
This observed NOC seems to indicate that the method converges.
The problem is only revealed when we see the exact errors:
tabrate('||Eh-Eexact||', E_err)
+------+---------------+------+
| h | ||Eh-Eexact|| | rate |
+------+---------------+------+
| 1/4 | 0.7450423 | * |
| 1/8 | 0.7207181 | 0.05 |
| 1/16 | 0.7056334 | 0.03 |
| 1/32 | 0.6963176 | 0.02 |
| 1/64 | 0.6905425 | 0.01 |
+------+---------------+------+
E_err
array([0.74504234, 0.72071807, 0.7056334 , 0.69631757, 0.69054246])
These errors do not go to zero, thus showing that the method converges, but to something distant from the exact solution! Such problems are not just theoretical curiosities. It can be a real issue in practice, when the exact solution is not known and one relies on the apparent convergence to verify correctness of methods.
The reason for this spectacular failure in this example can be traced to the fact that while (MM) sets \(E\) in \(\mathring{H}(\operatorname{curl})\), the formulation (CD) sets it in a smaller subspace. If one attempts to construct a theory for (CD), one sees that Lagrange approximations of (CD) solutions can only converge to solutions within \(H^1.\) Since the exact Maxwell solution we constructed is in \(\mathring{H}(\operatorname{curl})\), but not in \(H^1\), the (CD)-based Lagrange method has no chance to get close to the exact solution.
Question for discussion:
Are there rigorous sufficient conditions when (CD) and Lagrange elements are guaranteed to work?
11.8.3. Rates for Nédélec approximation#
We now repeat the same study as above, but now using the mixed formulation (MM) with Nédélec elements, for which we know there is a mathematically rigorous convergence theory.
def SolveByNedelecSuccessive(hcoarse=1/4, p=1, nrefinements=6):
Ephis = []
Xs = []
meshes = []
mesh = MeshLshape(maxh=hcoarse, pmaxh=hcoarse)
for ref in range(nrefinements):
meshes.append(ng.Mesh(mesh.ngmesh.Copy()))
Ephi, Xh = SolveByNedelec(mesh, Eexact, CF(0), p=p)
Ephis.append(Ephi)
Xs.append(Xh)
mesh = meshes[-1]
mesh.ngmesh.Refine()
return Ephis, Xs, meshes
def convergence_study2(Ephis, Xs, meshes):
E_diff = []
E_err = []
fine_Ephi = ng.GridFunction(Xs[-1])
Efine = fine_Ephi.components[0]
with ng.TaskManager():
Efine.Set(Ephis[-1].components[0])
for i in range(len(meshes)-1):
dE = Efine - Ephis[i].components[0]
E_diff.append(sqrt(Integrate(dE*dE, meshes[i])))
dE = Eexact - Ephis[i].components[0]
E_err.append(sqrt(Integrate(dE*dE, meshes[i])))
return np.array(E_diff), np.array(E_err)
Ephis, Xs, meshes = SolveByNedelecSuccessive(hcoarse=1/4, p=0,
nrefinements=6)
Ediff, Eerr = convergence_study2(Ephis, Xs, meshes)
tabrate('NOC: ||Eh-Efine||', Ediff)
tabrate('||Eh-Eexact||', Eerr)
+------+-------------------+------+
| h | NOC: ||Eh-Efine|| | rate |
+------+-------------------+------+
| 1/4 | 0.1704132 | * |
| 1/8 | 0.1094602 | 0.64 |
| 1/16 | 0.0691121 | 0.66 |
| 1/32 | 0.0415876 | 0.73 |
| 1/64 | 0.0222149 | 0.90 |
+------+-------------------+------+
+------+---------------+------+
| h | ||Eh-Eexact|| | rate |
+------+---------------+------+
| 1/4 | 0.1717179 | * |
| 1/8 | 0.1105870 | 0.63 |
| 1/16 | 0.0706518 | 0.65 |
| 1/32 | 0.0449073 | 0.65 |
| 1/64 | 0.0284505 | 0.66 |
+------+---------------+------+
Clearly both the exact errors and the numerically estimated errors with respect to the finest solution are converging to zero. Whatever one’s appreciation for the mathematical elegance of Nédélec elements, their practical value is difficult to dispute, as shown by this example.
11.9. Summary#
We saw a surreptitious form of failure in a finite element method that converges but to something far away from the true solution.
Two formulations for a 2D Maxwell problem were derived, a mixed formulation (MM) and a curl-and-div (CD) formulation. The former sets the \(E\)-field in \(H\)(curl), while the latter places it in a smaller space.
Lagrange elements appeared to give a plausible alternative to the mixed method using Nédélec elements, but it failed on a L-shaped domain.
Because we know the true nature of the singularity on an L-shaped domain, we were able to identify a case where good NOC was not an indicator of convergence to the true solution.
There are true Maxwell solutions in \(H\)(curl) which are not regular enough to be in \(H^1\) and we have \(H\)(curl)-convergent methods that can approximate them.