# Lagrange multiplier for boundary conditions in pure Neumann problem

I'm trying to solve $$-u''=\cos(2 \pi x)$$ with boundary conditions $$u'(0)=u'(1)=0$$ and the constraint $$\int_{0}^1 u = 0$$

I have to use linear finite elements, so let's assume that I have $$M$$ degrees of freedom. Since the contributions from the Neumann data on the r.h.s. vanish, the weak formulation is to find a $$u_h \in V_h$$ such that:

$$(\sum_{j=0}^{M-1} U_j\phi_j,\phi_i) = (f,\phi_i)$$ for every $$i=0,\ldots,M-1$$

Now, what happens is that I get the matrix $$A$$ (notice the $$\frac{1}{h}$$, rather than $$\frac{1}{h^2}$$ coming from the integration of shape functions)

$$A = \frac{1}{h} \left[ \begin{array}{ccccccccc} -2 & 1 & & & & \\ 1 & -2 & 1 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & 1 & -2 & 1 \\ & & & & 1 & -2 \\ \end{array} \right]$$

which is known to be singular. As I've read in Larson, Bengzon's book I can add a lagrange multiplier by solving the system:

$$\begin{bmatrix} A & C^T \\ C & 0 \end{bmatrix} \begin{bmatrix} U \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}$$

where $$C$$ is a vector with $$C_{i}=\int_{\Omega} \phi_i dx$$. Since, in my case, the grid is uniform with step-size $$h$$, I have that $$C_{0}=C_{M-1}= \frac{h}{2}$$ and all the other entries are equal to $$h$$

However, my solution is not correct, i.e. the boundary conditions are not fulfilled, as you can see from the picture

Quesiton: What can be wrong? I don't see any possible issue. I assembled the stiffness matrix as if there were no constraint, so I don't see any source of bugs. Any help or hint is highly appreciated.

If you want to reproduce

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

#Set some parameters
M = 81#DoFs
h = 1/(M-1)
A = np.zeros([M,M])
diag = 2*np.diag(np.ones(M))
subdiag = -1*np.diag(np.ones(M-1),-1)
suprdiag = -1*np.diag(np.ones(M-1),+1)

A= (diag + subdiag + suprdiag)/(h)

AA = np.zeros([M+1,M+1]) #matrix with one more dimension
C = np.zeros(M+1)
C[1:-2] = h
C[0]=0.5*h
C[M-1] = 0.5*h
#build matrix to use in the solver, with size M+1
AA[:-1,:-1] = A
AA[-1,:] = C
AA[:,-1] = C

def rhs(x):
return np.cos(2*np.pi*x)

x = np.linspace(0,1,M)
f = np.zeros(M+1)
for i in range(1,M-1): #x_1,...,x_{M-2} are interal dofs

f[M]=0 #constraint

u = np.linalg.solve(AA,f)

plt.plot(x,u[0:-1],'--')

• Your matrix does not seem to be right. The first and last coegficients on the main diagonal should be -1 since there you only have contributions from one side. Nov 3, 2021 at 12:23

I tried implementing the same approach and got the following solution:

Here is my code after pip install scikit-fem==4.0.0:

import numpy as np
from scipy.sparse import bmat, csr_matrix
from skfem import *
from skfem.models import laplace
from skfem.visuals.matplotlib import plot

m = MeshLine().refined(5)
basis = Basis(m, ElementLineP1())

A = laplace.assemble(basis)

@LinearForm
return np.cos(2. * np.pi * w.x[0]) * v

@LinearForm
def unit(v, w):
return v

C = csr_matrix(unit.assemble(basis))

K = bmat([[A, C.T],
[C, None]]).tocsr()
F = np.hstack((f, [0]))

x = solve(K, F)
plot(m, x[:-1])


A[0, 0] = 1/h

It basically does the same exception for the first and the last entries of A as you have already done for f. E.g., the derivative of the last shape function is $$\phi^\prime = \frac{1}{h}$$ and the integral is $$\int_{1-h}^1 (\phi^\prime)^2 dx = \frac{1}{h^2} \cdot h = \frac{1}{h}$$.
• When computing the stiffness matrix if your b.c. were only on the LHS for the equation $f''(x)=0$ wouldn't integration constants intervene and the matrix be $-c_n\delta_{mN}-A_{mn}+\delta_{nN}$, Hence with last line and last column non empty, where the C's are integration constants ? Nov 14, 2021 at 6:00