I'm studying this problem:
$$ \begin{equation} \left (\mathcal P\right ) \quad \begin{cases}\label{P} -u''(x) = f(x), \quad \text{sur} \quad I = ]0,1[\\ u'(0) = 1 \\ u(1) = 0 \end{cases} \end{equation} $$
where
$$ f(x) = H\left(x-\frac{1}{2} \right), \quad \text{H = Heaviside function} $$
$$ \begin{align*} -u'(x) & = (x-1/2)H\left(x-\frac{1}{2}\right) + c_1, \quad u'(0) = 1 \implies c_1 = 1 \\ & = \frac{1}{2}(2x-1)H\left(x-\frac{1}{2}\right) + 1 \end{align*} $$
$$ \int (2x-1)H(x-{1}{2}) \, dx = \left[ \frac{1}{2}(2x-1)^2H\left(x-\frac{1}{2}\right) \right] - \int (2x-1)H\left(x-\frac{1}{2}\right) \, dx $$
$$ \implies \int (2x-1)H\left(x-\frac{1}{2}\right) \, dx = \frac{1}{4}(2x-1)^2H\left(x-\frac{1}{2}\right) $$
$$ u(x) = -\frac{1}{8}(2x-1)^2H\left(x-\frac{1}{2}\right) + x + c_2, \quad u(1) = 0 \implies c_2 = -\frac{7}{8} $$
the exact solution is:
$$ u(x) = -\frac{1}{8}(2x-1)^2H\left(x-\frac{1}{2}\right) + x -\frac{7}{8}, \quad u\left(\frac{1}{2}\right) = \frac{-3}{8} $$
We know that we need to search the solution in H^1(I) with the boundary conditions. How can i show that it's also in H^2(I)?
Also, at the end with the P1 finite element aproximation, i have this matrix:
$$ A = \begin{pmatrix} \dfrac{2}{h} & -\dfrac{1}{h} & 0 & \dots & \dots & 0\\ -\dfrac{1}{h} & \dfrac{2}{h} & -\dfrac{1}{h} & 0 & \dots & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots & \vdots\\ 0 & \ddots & \ddots & \ddots & \vdots & 0\\ 0 & \dots & \ddots & -\dfrac{1}{h} & \dfrac{2}{h} & -\dfrac{1}{h} \\ 0 & \dots & \dots & 0 & -\dfrac{1}{h} & \dfrac{2}{h} \end{pmatrix} $$
which i easily implemented with in Fortran. I'm wondering how can i implement the source term?
The variational formulation of this problem is $$ \begin {cases} u\in V, \forall v\in V\\ \int_{0}^{1} u'(x)v'(x) \, dx = \int_{0}^{1} f(x)v(x) \, dx - v(0) \end {cases} $$
For the $H^2$-regularity, since the space $V=\{v\in H^1(0,1):v(1)=0\}$, and $C^\infty_0(0,1)\subset V$, from the weak formulation, we immediately get that
$$\int_0^1u'(x)v'(x)=-\int_0^1(-f)v\quad\forall v\in C^\infty_0(0,1),$$
and so $-f$ is the second weak derivative of $u$. Since $-f\in L^2(0,1)$, it follows that $u\in H^2(0,1)$.
For the right-hand side implementation, I assume you have a uniform finite element grid, and that the middle gridpoint lies on the point $x=1/2$, and that the finite element space is piecewise linear. I'm also assuming that $H$ is zero on $[0,1/2)$ and $1$ otherwise.
You then just have to calculate $L(v_i):=\int_0^1 fv_i\,dx-v_i(0)=\int_{1/2}^1v_i\,dx-v_i(0)$ for each linear basis function $v_i$. Interestingly for $v_0$ (the basis function that is $1$ at $x=0$ and $0$ at the first internal gridpoint $x_1=h$, we have $L(v_0) = -v_0(0)=-1$. For all the $v_i$'s that are supported only in $[h,1/2]$, we have $L(v_i)=0$. For the $v_i$ that is $1$ at $x=1/2$, we have $L(v_i) = \int_{1/2}^{1/2+h}v_i = \sqrt{h^2+1}$, then for all the rest of the $v_i$'s $L(v_i)=\int_{1/2}^{1}v_i = 2\sqrt{h^2+1}$ (since we integrate over the full support of $v_i$).