0
$\begingroup$

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

enter image description here

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],'--')
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – nicoguaro
    Nov 3 at 12:23
2
$\begingroup$

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

Discrete 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}$.

$\endgroup$
1
  • $\begingroup$ 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 ? $\endgroup$
    – Cretin2
    Nov 14 at 6:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.