Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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$?
share|improve this question

2 Answers 2

up vote 4 down vote accepted

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?

share|improve this answer
The forcing function is finite in magnitude, but its support is very small in comparison to the domain... so small that I treat it as if it were a discrete point source. –  Paul Feb 14 '13 at 21:17
Cool. Is its actual form worse to compute than the bump I gave? Maybe you should just use it with a broader base of support. –  Bill Barth Feb 14 '13 at 21:29
That can work... Thanks! Incidentally... if I wanted to extend this function into two dimensions, I should use $$f_\epsilon=e^{1-\frac{1}{1-\left((x-x_0)/ \epsilon\right)^2-\left((y-y_0)/\epsilon\right)^2}}$$, right? –  Paul Feb 15 '13 at 0:58
The best way to extend to multiple dimensions is to use the product $f_\epsilon(x)f_\epsilon(y)$ (which I don't think is equivalent to what you wrote). –  Bill Barth Feb 15 '13 at 20:02

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.

share|improve this answer
This is very intriguing... Initially (with the discontinuous source terms as defined above), my finite difference solution was piecewise linear function with a sharp edge at $x_0$. As I refine my mesh, the solution actually did converge to a linear function connecting the boundary values as you stated. I believed this to be a mistake, since the effect of the source term diminished as the mesh was refined. –  Paul Feb 16 '13 at 14:19
Perhaps I should start with box function like $$f_\epsilon=\left\{\begin{array}{ll} 0 &\quad x<x_0-\epsilon\\ 1 &\quad x_0-\epsilon\leq x\leq x_0+\epsilon \\ 0 &\quad x>x_0+\epsilon \end{array}\right.$$ and regularize it instead... –  Paul Feb 16 '13 at 14:27
Ultimately you need to let yourself be guided by the physics of your problem. The physical interpretation of my statement is essentially this: You have a force density $f(x)$ that is finite and only acts on a set of size zero. Consequently, the total force that acts is zero, and the system isn't going to show any reaction to it. You'll have to think about what you want the equation to describe first, not how to regularize a formulation that doesn't make any physical sense. –  Wolfgang Bangerth Feb 16 '13 at 20:39
I think @Paul's response to my question about the nature of the forcing function is the most illuminating. He's been approximating the function implied by the physics with the form he originally gave above. The "real" function has finite amplitude and tiny, but finite support. I think that he should just use the form of this function with much larger support (so that he gets a few samples in the non-zero part of the function) and then refine the mesh and support together until it converges. –  Bill Barth Feb 16 '13 at 20:57
Yes, that would make sense: Regularize the function in such a way that its mass remains constant. –  Wolfgang Bangerth Feb 17 '13 at 18:31

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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