Solving a PDE containing DiracDelta

I want to get the answer from a PDE:

\begin{align*} \frac{\partial \rho(r,t)}{\partial t}&=Dr^{-2}\frac{\partial}{\partial r}r^2h(r)e^{-U(r)}\frac{\partial}{\partial r}e^{U(r)}\rho(r,t)-\left(\frac{\lambda}{g(\sigma)}\frac{\delta(r-\sigma)}{4\pi\sigma^2}+\frac{3\lambda}{4\pi r^6}\right)\rho(r,t)\\ U(r)&=-\frac4r\text{; }\sigma=D=\lambda=h(r)=1\\ \rho(r,t=0)&=g(r)=e^{-U(r)}\text{; }\lim_{r\to\infty}\rho(r,t)=1\text{; }\left.4\pi\sigma^2 Dh(\sigma)\frac{\partial}{\partial t}e^{U(r)}\rho(r,t)\right|_{r=\sigma}=\lambda e^{U(\sigma)}\rho(\sigma,t) \end{align*}

The PDE has one DiracDelta[] function. When I tried to solve the equation with NDSolve[], it gave me

NDSolve::ndnum: Encountered non-numerical value for a derivative at t==0.


So I wanted to know the reason of message.
Using r-sigma instead of the delta function, DiracDelta[r-sigma], I could get an answer from NDSolve without the above error message. However, it was not the answer of interest.

How can I solve the problem with DiracDelta[] function?

*Codes start from here :

σ = 1
rc = -4
λ = 1
sink[r_, λ_] = (3 λ)/(4 π r^6)
U[r_] = rc/r
g[r_] = Exp[-U[r]]
Dif = 1
h[r_] = 1
r0 = 1;
tmax = 0.01;
rmax = 1000 Sqrt[6 Dif tmax] + r0;

sol = NDSolve[{D[ρ[r, t], {t, 1}] ==
Dif r^-2 D[
r^2 h[r] Exp[-U[r]] D[Exp[U[r]] ρ[r, t], {r, 1}], {r,
1}] - (λ/g[σ] DiracDelta[r - σ]/(
4 π σ^2) + sink[r, λ]) ρ[r,
t], ρ[r, 0] == Exp[-(U[r] - rc/rmax)], ρ[rmax, t] ==
1, 4 π Dif h[σ] (
D[Exp[U[r]] ρ[r, t], {r, 1}] /.
r -> σ) == λ Exp[U[σ]] ρ[σ,
t]}, ρ, {r, σ, rmax}, {t, 0, tmax},
PrecisionGoal -> 4, StartingStepSize -> 0.0001,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 100, "MaxPoints" -> 200}}]

ρn = First[ρ /. sol]
Plot3D[ρn[r, t], {r, σ, rmax/100}, {t, 0, tmax}, PlotRange -> All]

-
Your code works as it is for me. There are warnings though, like inconsistency of boundary and initial conditions. – b.gatessucks Sep 27 '12 at 7:47
@b.gatessucks I tried this code with MatheMatica7. As your advice, it can be solved by MatheMatica8. Thank you. – Jaehoon Kim Sep 27 '12 at 8:49

It helps if you change the DiracDelta[r- σ] with D[HeavisideTheta[r - σ], r] :

sol = NDSolve[{D[ρ[r, t], {t, 1}] == Dif r^-2 D[r^2 h[r] Exp[-U[r]] D[Exp[U[r]] ρ[r, t], {r, 1}], {r,
1}] - (λ/g[σ] D[HeavisideTheta[r - σ], r]/(4 π σ^2) + sink[r, λ]) ρ[r, t], ρ[r, 0] == Exp[-(U[r] - rc/rmax)], ρ[rmax, t] == 1, 4 π Dif h[σ] (D[Exp[U[r]] ρ[r, t], {r, 1}] /.
r -> σ) == λ Exp[U[σ]] ρ[σ,
t]}, ρ, {r, σ, rmax}, {t, 0, tmax}]


-
smart but why? isn't it the same? – chris Sep 27 '12 at 8:31
@b.gatessucks, Thank you. It works well. :) – Jaehoon Kim Sep 27 '12 at 8:47
@b.gatessucks is it fair to say this is not really an answer? I.e. changing the dirac delta makes no difference? – chris Sep 27 '12 at 13:22
@chris It is fair, after checking I realized the reduced amount of warnings is due to removing the PrecisionGoal and Method parts, not the substitution. – b.gatessucks Sep 27 '12 at 14:11