# Lagrange multiplier dominates total residual

I’m trying to solve a PDE Lf=0 (for some operator L) with the mass constraint \int f dx=1. I suspect there is a non-trivial solution f\not\equiv 0 but I know for sure that f\equiv 0 is also a solution for my model. I wrote the below implementation with a Lagrange multiplier, but when I test the residual against f it shows that the PDE residual (excluding the Lagrange multiplier part) is of the same order (10e-1) as that of the Lagrange multiplier, whereas the total residual is 10e-14 after using solve. How can I ensure that both residuals are small, i.e. the PDE residual and the total residual are both small?

from dolfin import *
import matplotlib.pyplot as plt

# Sub domain for Dirichlet boundary condition
class NofluxBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((x < DOLFIN_EPS or x > (1.0 - DOLFIN_EPS)) \
and on_boundary)

# Sub domain for Periodic boundary conditional
class PeriodicBoundaryC(SubDomain):

# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x < DOLFIN_EPS and x > -DOLFIN_EPS and on_boundary)

# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y = x - 1.0
y = x

class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x < DOLFIN_EPS and x > -DOLFIN_EPS and on_boundary)

# Create mesh and finite element
mesh = UnitSquareMesh(32, 32)
pbc = PeriodicBoundaryC()
nofluxb = NofluxBoundary()
pb = PeriodicBoundary()
Q = FiniteElement("CG", triangle, 1)
LM = FiniteElement("Real", triangle, 0)
M = FunctionSpace(mesh, MixedElement([Q, LM]), constrained_domain=pbc)

# Create domain and boundaries
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
pb.mark(boundaries, 1)
nofluxb.mark(boundaries, 2)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

drc = Expression('-sin(pi*x)',degree=2)

# Create form
u = Function(M)
f, lm = split(u)
wt = TestFunction(M)
v, lmt = split(wt)

r = Expression('x',degree=2)
cospsi = Expression('cos(2*pi*x)',degree=2)
sinpsi = Expression('sin(2*pi*x)',degree=2)

+ lm*v*dx + lmt*f*dx - lmt*1*dx

solve(B == 0, u)

f,lm = u.split()

# Check residuals testing against f itself