right hand side function is not L-infty bounded

Question

Suppose we want to solve a two dimensional Poisson equation with infinity right hand side function $f$
$$
-\Delta u + u= f
$$
where $f$ is only in L^2, such as $f = r^{-1/6}$. The domain here can be the unit circle.

How to define the function $f$ in this case? As $f$ is infinity large at Origin.

Thanks~

Comments

Here, $(r,\theta)$ denotes the two dimensional polar coordinate. the right hand side function $f$ is in $L^2$. Supplemented with the zero Neumann boundary condition, $\frac{\partial u}{\partial n} = 0$, the problem is well posed. Has anyone solved such kind of problems?

Answer 1

I have two simple workarounds:

Option 1:
Perturbate r to avoid infinity at the origin, for example f= (r+eps)^(-1/6)...

Option 2:
Use Expression:

class F(Expression):
     def eval(self, x, values):
           r = sqrt(x[0]**2+x[1]**2)
           if r == 0.0:
                 values[0] = 1e16 # "infinity"
           else:
                 values[0] = r**(-1/6)

I'm guessing that this might not be what you wish, but it's at least a suggestion.

Comments

Thank you very much for your suggestions. According to your experience, how to select the eps or 1e16 in your second method ? If I want to test the convergence rate of the linear FEM for the above problem, do you think the above trick will affect the accuracy result ?

Answer 2

What is the result if you just naively try $f = r^{-1/6}$ ? As long as f is not evaluated at (or very close to) 0, you might be fine. Here is a related problem.

http://fenicsproject.org/qa/2549/question-on-quadrature-method-for-singular-integrands

Comments

You are right. It seems okay if the singular point is not on the mesh vertex. If I just use $f=r^{-1/6}$ and the Origin is a vertex, I cannot get $u$.

---

From my understanding (note the Novice after my user name) the singular point can occur at a vertex, but not at a Gaussian quadrature point, since this is where f is eventually evaluated. Even if there is a Gaussian quadrature point at r=0, option 2 of the other answer should work. Can you post a small code example to show where the problem lies?

---

please see the following for example

---

from dolfin import *
import math

nx_list = [32]
for nx in nx_list:

mesh = RectangleMesh(-1,-1,1,1,nx,nx)

V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)

F = Expression("pow((1e-10)+x[0]*x[0]+x[1]*x[1],-1.0/6.0)")

F = Expression("pow(x[0]x[0]+x[1]x[1],-1.0/6.0)")

A = assemble(inner(grad(u),grad(v))*dx+u*v*dx)
b = assemble(F*v*dx)
solverp = LinearSolver("lu")
p1 = Function(V)
solverp.solve(A, p1.vector(), b)

print errornorm(F,interpolate(F,V))

plot(mesh)
plot(p1,mesh)

interactive()