Solving a PDE containing DiracDelta

Question

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?

Thanks in advance :)

*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]

Answer 1

The Dirac delta makes the solution discontinuous and Mathematica cannot handle PDEs with discontinuities at the moment.

You will have to construct your own solver in Mathematica to do so. In your case, the location and magnitude of the discontinuity is known analytically, so you can construct derivative matrices that take this into account, and then use the method of lines. For example, see:

https://arxiv.org/abs/1406.4865

Let me know if you need more details.

Answer 2

First of all, I'd like to point out that, in OP's specific example the DiracDelta actually doesn't play a important role i.e. its contribution to the final solution is very weak. In the following answer I'll multiply it with 10^6 for a better illustration.

OK, let's dealing with the DiracDelta. As far as I can tell, the easiest way to handle DiracDelta inside NDSolve is to use a smooth function to approximate it. In OP's case, appro' where

appro = With[{k = 1}, Erf[k #]/2 + 1/2 &];

is already an acceptable approximation:

Plot[appro'[x - σ], {x, -σ, rmax}, PlotRange -> All, Axes -> None]

σ = 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;

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

sol = With[{coef = 10^6}, 
   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[σ] coef  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]} /. 
     DiracDelta -> appro', ρ, {r, σ, rmax}, {t, 0, tmax}, 
    Method -> Union[mol[True, 100], mol[4000, 4]]]];

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

"ScaleFactor"->100 is added because of the reason mentioned here. (Well, but in OP's case, the inconsistency seems not to be that significant, I added it just in case. )

The resulting graph is:

Again, notice I've multiply DiracDelta with 10^6 for a better illustration.

Answer 3

I admire the idea of @b.gatessucks. However, I would like to note that your problem may be has a more straightforward solution. In fact the delta function, such as yours generates a solution with a kink, and fixes a specific boundary conditions in the kink points. Then one can reveal these boundary conditions, and then solve the equation at both side of the kink, where there is no delta function. There the solution is straightforward. I will illustrate this within a much more simple example to only reveal the essence. Let us consider a simple stationary equation:

with the boundary conditions y(x)->0when x turns into the both infinities. Here a and k are constants. This equation is like your one, in the sense that it contains the term delta(x)*y(x).

Now, let us integrate this equation from -b to b, where b>0 and then pass on to the limit b->0:

It is clear that the second term in the left-hand part of the equation vanishes after integration and passing to the limit, if the function is continuous, while the first one represents Limit(y'(b)-y'(-b), b->0)=[y'(0)], where the square brackets denote a jump of a function in the point. This does not vanish, if the point x=0 is a kink point of the continuous function. The term in the right-hand part of the equation also does not vanish. One finds:

[y'(0)]=k y(0).

The latter is the additional boundary condition we were looking for. We can now solve the equation outside of the point x=0 of the kink. The solution is

where A is an arbitrary constant, the sign "-" is valid to the right and "+" - to the left of the kink. Now, our additional boundary condition yields -2a=k which fixes the only value of a at which the problem has the solution, so it is the spectral condition.

I hope this helps. Have fun!