We have our PDE with Neumann boundary conditions, with $f \in L^2([0,1])$
$$
-\Delta u = f \quad \in [0,1], \quad \partial_x u(0) = \partial_x u(1) = 0
$$
With compatibility conditons $\int f dx = 0 $ and some uniqueness condition (say $\int u = 0 $ ).
The weak formulation of this problem is to find $u \in H^1$ satisfying
$$ \int_0^1 \dot u\ \dot v\ dx = \int_0^1 f v dx, \quad \forall v \in H^1$$
Typically, one solves numerically this problem employing the P1 finite elements: Take a spatial discretization $0 = x_0 < x_1 = 1/N < …. < x_n = 1 $, and then we replace the space $H^1$ by
$$ V_h = \{ v \in C([0,1])\ :\ v\big|_{[x_i, x_{i+1}]} \text{ is affine } \} $$
Then we take as ansatz $u$, our seeked solution, to belong in $V_h$, and we test with functions $v$ in $V_h$ too. This leads to the standard system
$$
\sum_{i=0}^N u_i \int_0^1 \dot \varphi_i(x) \dot \varphi_j(x) dx = \int_0^1 f(x) \varphi_j(x) dx, \quad j =0,…,N
$$
where $\varphi_i$ are the hat functions that satisfies $\varphi_i(x_j) = \delta_{ij}$.
My question is: Does the obtained discrete solution $u$ satisfies Neumann boundary conditions?
To my understanding, Neumann boundary conditions in the space $V_h$ would correspond to $u_0 = u_1, u_{N-1} = u_N$.
Now, checking the system of equations at $j=0$ gives
$$ u_0 \int_0^1 \dot \varphi_0 (x) \dot \varphi_0(x) + u_1 \int_0^1 \dot \varphi_1 (x) \dot \varphi_0(x) = \int_0^1 f(x) \varphi_0(x) dx
$$
replacing the values of the hat function gives
$$ u_0/N – u_1/N = \int_0^{1/N} f(x) (1 – Nx) dx $$
Fair enough, if I have not made any mistake so far, it seems that the solution doesn't satisfies Neumann boundary condition. But, does refining the mesh gives a solution with coefficients $|u_0 – u_1| \rightarrow 0$ ?
I considered $ f(x) =1 $ for $x <1/2$ and $f(x) = -1$ for $x > 1/2$, then on the equation associated to $j=0$ we have
$$ \frac{u_0}{N} – \frac{u_1}{N} = \int_0^{1/N} (1-Nx) dx = \frac{1}{N} – N \frac{1}{2N^2} = \frac{1}{2N} $$
Which would imply $u_0 – u_1 = 1/2$ nomatter how refined is the problem.
My intuition tells me that this cannot be possible, being FEM the standard to numerically solve elliptic PDE's. What am I doing wrong?
Best Answer
I've found my mistake, I misscalculated the equation associated to j=0, it is
$$ Nu_0 - Nu_1 = \int_0^{1/N} f(x) (1 - Nx) dx $$
From which we can see that indeed, the Neumann boundary conditions are not satisfied by the discrete solution, but when taking the limit of the the discretization width to zero ($\frac{1}{N} \rightarrow 0 $) then the values of $|u_0 - u_1| \rightarrow 0$.