Regularization of a discontinuous source term in an elliptic pde

Question

Suppose I'm solving $$\frac{d}{dx}\left(K(x)\frac{du}{dx}\right)=f \text{ in }\Omega,$$ $$u=g \text{ on } \partial\Omega$$where $K(x)$ is smooth and $$ f(x) = \left\{ \begin{array}{ll} 0 & \quad x \neq x_0 \\ 1 & \quad x =x_0 \end{array} \right. $$

I tried solving this problem by a finite difference method with f as given, but I noticed some strange behavior in the solution as the number of discrete nodes gets larger. I'm pretty sure the strange behavior is due to the discontinuity at $x_0$.

So, I propose to replace $f$ by a function $f_\epsilon$ such that $f_\epsilon=0$ everywhere except in the small region $B_\epsilon=\{x\text{ s. t. }|x-x_0|<\epsilon\}$. In $B_\epsilon$, $f_\epsilon$ is a polynomial such that

$$\begin{array}{ll} f_\epsilon(x_0)=1 & \quad f_\epsilon'(x_0)=0 \\ f_\epsilon(x_0)=0 & \quad f_\epsilon'(x_0-\epsilon)=0 \\ f_\epsilon(x_0+\epsilon)=0 & \quad f_\epsilon'(x_0+\epsilon)=0 \end{array}$$

and $f_\epsilon=0$ throughout the rest of the domain. This gives me a polynomial of degree 5 and ensures $f_\epsilon\in C^1(\Omega)$.

However, I suspect that since finite difference methods seek classical solutions, we probably need more differentiability; that is, $f_\epsilon\in C^2(\Omega)$. It is not immediately obvious to me what other conditions I can/should impose on $f_\epsilon$ to ensure this extra regularity.

So, I pose the following questions:

1. Does $f_\epsilon$ need to be in $C^2(\Omega)$?
2. What other conditions can I impose on my polynomial to ensure this extra regularity?
3. Are other (non-polynomial) functions better suited to approximate the discontinuous source term?
4. To observe the true behavior of the original discontinuous source term, should I solve a sequence of regularized problems with $f_\epsilon$ such that $\epsilon\rightarrow 0$?

Answer 1

A more typical regularization would be to use

$$ f_\epsilon (x) = e^{1-\frac{1}{1-((x-x_0)/\epsilon)^2}},\, {\rm for }\, |x-x_0|<\epsilon $$ and zero otherwise.

This bump function is $C^\infty$, and you can control its width and height arbitrarily.

With regard to your questions, $f$ would typically need to be less regular than your solution $u$, so that's not your problem. And, yes, you would want to solve a sequence of problems, but once you have no grid points inside $|x-x_0| < \epsilon$ besides the one at $x_0$, there's no point in shrinking $\epsilon$ any more. You'd need to refine the mesh at that point.

Finally, it's kind of an odd forcing function. Are you sure you didn't mean to use a delta function?

Answer 2

That can't be the problem you want to solve. Given your right hand side has nonzero values only at $x=x_0$ but is finite there, you have that $\|f\|_{L^2}=0$. On the other hand, if you had zero boundary values, you have the stability estimate $\|\bar u\|_{H^1}\le C\|f\|_{L^2}$ where $\bar u$ is the solution with zero boundary values; in other words, it is zero. The solution to your problem is then the same as to a problem with a zero right hand side (i.e., in 1d, a linear function connecting the boundary values at the left and right).

That you don't converge to this solution is an artifact resulting from the fact that the finite difference method wants to sample point values of the right hand side, but that these are not defined for $L^2$ functions. You can't resolve this by mollification because every appropriately mollified function of zero $L^2$ norm also needs to have zero $L^2$ norm -- i.e., the only useful mollification would be the zero function.