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 have a density function $\rho(r,z)$ and I want to calculate the minimum distance $d_{min} = \sqrt{r^2 + z^2}$ from the center (0,0) in which $\rho$ becomes negative. The easiest way is to find where $\rho = 0$, since for all greater distances is always $\rho < 0$. In my system, there is a variable parameter $\delta$ which takes values in the interval [0, 1]. So, I want to monitor the evolution of $d_{min}$ as a function of $\delta$. The corresponding Mathematica code is the following:

Clear["Global`*"];
Off[FindMinimum::eit];

Vg = v0^2/2*Log[r^2 + (1 + δ)*(a + Sqrt[h^2 + b*z^2])^2 + cb^2]; 
Vn = (-G*Mn)/Sqrt[r^2 + z^2 + cn^2]; 
Vh = (-G*Mh)/Sqrt[r^2 + z^2 + ch^2];
Vt = Vg + Vn + Vh;
Vr = D[Vt, r];
Vz = D[Vt, z];
Vrr = D[Vr, r];
Vzz = D[Vz, z];
ρ = 2.325/(4*π*100)*(Vrr + Vr/r + Vzz);

G = 1;
v0 = 25; a = 3; h = 0.12; b = 1.5; cb = 2.5;
Mn = 400; cn = 0.25;
Mh = 10000; ch = 25;

data = {};
Do[
  sol = FindMinimum[{r^2 + z^2, ρ == 0}, {r, z}];
  dmin = sol[[1]];
  AppendTo[data, {δ, Sqrt[dmin]}],
  {δ, 0, 1, 0.05}
  ]

S1 = ListPlot[data, Joined -> True, Frame -> True, 
FrameLabel -> {"δ", Subscript["d", "min"]}, 
FrameStyle -> Directive[FontSize -> 17, FontFamily -> "Times"], 
PlotStyle -> Black, ImageSize -> 600, PlotRange -> {0, 200}]

Here is the plot (S1) I get using the values obtained from the FindMinimum inside the Do loop

plot

We see, that there are two horizontal areas regarding the evolution of $d_{min}$, the first one when $0 \leq \delta \leq 0.3$ and the second when $0.7 \leq \delta \leq 1$. I am quite sure that for some unknown to me reason these areas are in fact wrong. In particular, when $\delta = 0$ density is always positive so, $d_{min}$ should be infinite not almost zero as indicated in the plot. In fact, the plot should be like this

plot2

I suspect, that there might be some issue with the precision of the calculations. Any ideas?

share|improve this question
add comment

1 Answer

Maybe a starting point :

rho[r_, z_, \[Delta]_] = 2.325/(4*\[Pi]*100)*(Vrr + Vr/r + Vzz)

tol = 10^-3;
sol = {#, FindMinimum[{r^2 + z^2, Abs[rho[r, z, #]] <= tol}, {r, z}]} & /@ Range[0, 1, 0.05] 

ListLinePlot[{#[[1]], Sqrt[#[[2, 1]]]} & /@ sol]

plot

share|improve this answer
    
A good starting point indeed but not exactly what I want. I need to know where $\rho = 0$ not where $\rho \leq 0.001$. In fact, if I replace in my original code the $\rho = 0$ inside \FindMinimum` with $\rho = 10^{-16}$ everything is OK. –  Vaggelis_Z Jul 8 '13 at 8:11
add comment

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.