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],'--')