4
$\begingroup$

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.

enter image description here

$\endgroup$

2 Answers 2

5
$\begingroup$

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: enter image description here This seems to agree with your numerical solution via FEM.

$\endgroup$
4
  • 2
    $\begingroup$ +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. $\endgroup$
    – Ian
    Commented Feb 11, 2022 at 2:13
  • $\begingroup$ 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. $\endgroup$
    – Ian
    Commented Feb 11, 2022 at 2:15
  • $\begingroup$ @Ian Ah yes good pick up. Thanks :) $\endgroup$
    – JackT
    Commented Feb 11, 2022 at 2:25
  • $\begingroup$ @JackT I'm missing the dependency of your solution from $u'(0)$, see Solution $\endgroup$
    – gpmath
    Commented Mar 28, 2023 at 11:47
2
$\begingroup$

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:

enter image description here

Which agrees with the FEM.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .