Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

1 Answer 1

up vote 4 down vote accepted

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

plot

share|improve this answer
1  
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

Your Answer

 
discard

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.