# 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. 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 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.

• +1. I'll point out that there is a natural physical interpretation of this condition. One can interpret this BVP as an equilibrium heat equation where you have a rod in perfect thermal contact with some cold material on the ends and surrounded in a warmer heat bath elsewhere. This compatibility condition requires that all the heat that enters/exits $(0,1/2)$ at $1/2$ exits/enters $(1/2,1)$ there as well. So it is really a conservation law.
– Ian
Feb 11 at 2:13
• That said there is a typo at the beginning. It is correct to have $u_1=u_2$ at $1/2$, it should be that $u_1'=u_2'$ at $1/2$ is erroneous.
– Ian
Feb 11 at 2:15
• @Ian Ah yes good pick up. Thanks :) Feb 11 at 2:25

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:

1. 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$$.
2. 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.