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$
3
  • 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
    Feb 11 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
    Feb 11 at 2:15
  • $\begingroup$ @Ian Ah yes good pick up. Thanks :) $\endgroup$
    – JackT
    Feb 11 at 2:25
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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.