Lagrange multiplier for boundary conditions in pure Neumann problem

Question

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[i]=(integrate.quad(lambda y: (1/h)*(y-x[i-1])*rhs(y),x[i-1],x[i])[0] + integrate.quad(lambda y: (1/h)*(x[i+1]-y)*rhs(y),x[i],x[i+1])[0])

f[0] = integrate.quad(lambda y: (1/h)*(-y+x[1])*rhs(y),x[0],x[1])[0]
f[M-1]=integrate.quad(lambda y: (1/h)*(+y-x[-1])*rhs(y),x[M-1],x[M-2])[0]
f[M]=0 #constraint

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

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

Answer 1

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
def load(v, w):
  return np.cos(2. * np.pi * w.x[0]) * v

f = load.assemble(basis)

@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])

In your code, try adding the following lines:

A[0, 0] = 1/h
A[M-1, M-1] = 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}$.