Strong vs Weak solution to one-dimensional elliptic PDE Consider the elliptic PDE
\begin{align}
&-\frac{\mathrm{d}}{\mathrm{d} x} \left(a(x) \frac{\mathrm{d}}{\mathrm{d} x}u(x)\right) = 1, \qquad 0 < x < 1,\\
&u(0) = u(1) = 0.
\end{align}
Here $a \in L^\infty(0,1) \cap C^0(0,1)^C$ is defined as
$$
a(x) = a_1, \text{ if } x \leq 1/2, \quad a(x) = a_2, \text{ if } x > 1/2,
$$
where $a_1$ and $a_2$ are positive real scalars.
I am having troubles with the difference between a strong and a weak solution of this problem. Since $a$ is not continuous, I would expect $u\in H_0^1(0,1) \cap H^2(0,1)^C$, so in particular (since we are in dimension one) I would expect $u \in C^0(0,1) \cap C^1(0,1)^C$. Actually somewhere in between, e.g. Hölder continuous.
I tried to solve the equation "by hand". In particular, I define $u$ piecewise as
$$
u(x) = u_1(x), \text{ if } x \leq 1/2, \quad u(x) = u_2(x), \text{ if } x > 1/2,
$$
Then, one can impose that $u_1$ and $u_2$ solve the equation in the strong sense on each side of $x=1/2$ and obtain
$$
u_1(x) = -\frac{1}{2a_1} x^2 + C_1x + C_2, \quad u_2(x) = -\frac{1}{2a_2} x^2 + D_1x+D_2.
$$
Now we can impose $u_1(0) = 0$ to obtain $C_2 = 0$, and $u_2(1) = 0$ to obtain
$$
D_1 + D_2 = \frac{1}{2a_2}.
$$
We expect the equation to be continuous, so $u_1(1/2) = u_2(1/2)$, which gives the condition
$$
\frac{C_1}{2} - \frac{D_1}{2} - D_2 = \frac{1}{8a_1} - \frac{1}{8a_2}.
$$
We now have three unknowns ($C_1$, $D_1$, and $D_2$), and two conditions. We could impose moreover that $u_1'(1/2) = u_2'(1/2)$ so that the solution $u \in C^1(0,1)$. This yields the condition
$$
C_1 - D_1 = \frac{1}{2a_1} - \frac{1}{2a_2}.
$$
The three linear equations above are solvable and thus define a solution $u$ on the whole domain which is $C^1$.
I tried to compare the solution (obtained by solving the linear equation and fixing $a_1 = 0.5$, and $a_2 = 2$) to a FEM solution on $1000$ elements. The FEM solution converges to the weak solution $u \in H_0^1$ such that
$$
\int_0^1 a(x) u'(x) v'(x) \, \mathrm{d}x = \int_0^1 v(x) \, \mathrm{d}x,
$$
for all $v \in H_0^1(0, 1)$. The result is in the picture below. The FEM solution (in blue) behaves a bit more "similarly" to what I expected. In particular, there is a discontinuity in the derivative at $x = 1/2$ and the solution is "only" Hölder continuous. In particular, the FEM solution is different than the strong solution (in red).
I think there's something wrong with my reasoning that yields the strong solution, but I cannot see why or where. In particular, since the weak solution exists and is unique for this equation, it should be equal to the strong solution when it exists. I somehow trust more the FEM solver in this case.

 A: The assumption $u_1'(1/2)=u_2'(1/2)$ is erroneous. As @Dan Doe mentioned in his answer, the correct compatibility condition is $a_1u_1'(1/2)=a_2u_2'(1/2)$. I will explain why this is the case rigorously.
Since $a(x)$ is smooth away from $x=1/2$ interior regularity theory tells us that $u$ is smooth away from $x=1/2$. Thus, we are justified in saying that $$u(x) = \begin{cases}
u_1(x), &\text{if }0<x<1/2 \\
u_2(x), &\text{if }1/2<x<1 
\end{cases} \tag{$\ast$}$$ where $u_1,u_2$ are smooth functions that solve $-(a_1 u_1')'=1$ in $(0,1/2)$ and $-(a_2 u_2')'=1$ in $(1/2,1)$ respectively. In order to ensure $u \in H^1_0((0,1)$ we must have $u_1(0)=u_2(1)=0$ and $u_1(1/2)=u_2(1/2)$; however, at this point we know nothing about the relationship between $u_1'(1/2)$ and $u_2'(1/2)$.
Now since $u$ is a weak solution of your PDE, $$\int_0^1 \varphi(x) \,dx = \int_0^1 a(x) u'(x) \varphi'(x) \, dx $$  for all $\varphi \in C^\infty_0((0,1))$. By $(\ast)$ and integration by parts, the right hand side of this equality can be written as \begin{align*}
\int_0^{1/2}a_1  u_1'(x) \varphi'(x) \, dx +\int_{1/2}^1a_2 u_2'(x) \varphi'(x) \, dx &=-\int_0^{1/2} (a_1 u_1')' \varphi \, dx -  \int_{1/2}^1 (a_2u_2')' \varphi(x) \, dx \\
&\qquad + \big ( a_1u_1'(1/2)  - a_2u_2'(1/2) \big ) \varphi(1/2)\\
&=\int_0^1  \varphi \, dx  \\
&\qquad + \big ( a_1u_1'(1/2)  - a_2u_2'(1/2) \big ) \varphi(1/2).
\end{align*} Thus, the correct compatibility condition between $u_1'(1/2)$ and $u_2'(1/2)$ is $$a_1u_1'(1/2)=a_2u_2'(1/2).$$

For the particular case $a_1=0.5$ and $a_2=2$, I got $$u_1(x)=-x^2+\frac 7 {10}x, \qquad u_2(x)=-\frac 1 4 x^2+\frac 7 {40} x+\frac 3 {40}.$$ The plot looks like this:  This seems to agree with your numerical solution via FEM.
A: Actually I am not sure if there is any strong solution here. In order to be able to write down $$-\Big(a(x) u'(x)\Big)' = 1 \tag{1}$$ not in an integral sense, you need differentiability of $a(x) u'(x)$ at $x = 0.5$. Your requirement, that $u_1'(x = 0.5) \overset{!}{=} u_2'(x = 0.5)$ does not give you this for discontinuous $a(x)$ at that point.
So I see two approaches to this problem:

*

*Your approach with two piece-wise solutions. Then, if you require $a_1u_1'(x = 0.5) \overset{!}{=} a_2u_2'(x = 0.5)$ you obtain that $\tilde D_1 := \frac{D_1}{a_2} = \frac{C_1}{a_1} =: \tilde C_1$ since you enforce $-0.5 + \tilde C_1 \overset{!}{=} -0.5 + \tilde D_1$.

*If you just integrate $(1)$ over the whole interval, you obtain
$$a(x) u'(x) = -x + \tilde C_1,$$ i.e., again only one constant of integration - the interval-wise ones are then given by division by the respective values of $a_1, a_2$. You can then pursue setting up the solutions on the different intervals.

In any case, you get then again $C_2 = 0, D_2 = \frac{1}{2a_2} - \frac{C_1}{a_2}$. Then, through requiring continuity at $x = 0.5$ in $u$, you have the condition
\begin{align}-\frac{1}{8a_1} + 0.5 \frac{C_1}{a_1}& = -\frac{1}{8a_2} + 0.5 \frac{C_1}{a_2} + \frac{1}{2a_2} - \frac{C_1}{a_2} \\
C_1 &= \frac{1/a_1 + 3 / a_2 }{4 (1/a_1 + 1 /a_2)}\end{align}
and the solution looks for $a_1 =0.5, a_2 = 2$ like this:

Which agrees with the FEM.
