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 the following Mathematica code

Clear["Global`*"];

V = -(M/Sqrt[b^2 + x^2 + λ*y^2 + (a + Sqrt[h^2 + z^2])^2]);
Vxx = D[V, {x, 2}];
Vyy = D[V, {y, 2}];
Vzz = D[V, {z, 2}];
ρ = 2.325/(4*π*100)*(Vxx + Vyy + Vzz);
ρyz = ρ /. {x -> 0};

M = 9500; a = 3; b = 4; h = 0.15; λ = 1.1;

Syz = ContourPlot[ρyz, {y, -50, 50}, {z, -50, 50}, Contours -> 20, 
ContourStyle -> Black, PlotPoints -> 50, 
RegionFunction -> Function[{y, z}, ρyz < 0], 
PerformanceGoal :> "Speed", FrameLabel -> {"y", "z"}, 
RotateLabel -> False, 
FrameStyle -> Directive[FontSize -> 17, FontFamily -> "Times"], 
Epilog -> Inset[Graphics@Text[Row[{"b = ", b}], 
  BaseStyle -> {17, FontFamily -> "Helvetica", Bold}], 
Scaled[{0.5, 0.95}]], ImageSize -> 550]

which produces this output

enter image description here

This plot shows the contours $\rho(y,z) < 0$. Now I would like to compute the minimum distance at which these contours approach the center (0,0). The distance is defined as $d = \sqrt{y^2 + z^2}$. Inspecting by eye the plot, we see that in this case the minimum distance is approximately $d_{min} \simeq 16$. But of course, this is not sufficient at all. So, my question is how could I compute the minimum distance? Then, I assume it would be easy enough to add a DO loop in order to see the evolution of the minimum distance when a parameter (i.e $b$ or $\lambda$) varies.

Many thanks in advance!

EDIT

Clear["Global`*"];

V = -(M/Sqrt[b^2 + x^2 + λ*y^2 + (a + Sqrt[h^2 + z^2])^2]);
Vxx = D[V, {x, 2}];
Vyy = D[V, {y, 2}];
Vzz = D[V, {z, 2}];
ρ = 2.325/(4*π*100)*(Vxx + Vyy + Vzz);
ρyz = ρ /. {x -> 0};

M = 9500; a = 3; h = 0.15; λ = 1.1;

data = {};
Do[
   sol = FindMinimum[{y^2 + z^2, ρyz == 0}, {y, z}];
   dmin = sol[[1]];
   AppendTo[data, {b, Sqrt[dmin]}],
 {b, 0, 12, 1}

]

share|improve this question

2 Answers 2

up vote 8 down vote accepted

Can use FindMinimum. As below, give the equality constraint that the point lies on the region boundary. Or could use constraint "on or inside" that is, <= instead of ==.

FindMinimum[{y^2 + z^2, \[Rho]yz == 0}, {y, z}]

(* Out[608]= {253.037069551, {y -> -15.728359276, z -> 2.37818923488}} *)
share|improve this answer
    
Oh man, one minute too late. I used NMinimize and had to get rid of the NMinimize::incst :-( Of course +1. –  halirutan Apr 1 '13 at 16:45
    
See my EDIT. I added a DO loop bu the program complaints about tolerance! Why? –  Vaggelis_Z Apr 1 '13 at 16:55
    
And to think I almost posted this as a comment (sound of @halirutan gritting teeth).. –  Daniel Lichtblau Apr 1 '13 at 18:07
    
You shouldn't be scorning @halirutan. He could send you a copy of NKS in retaliation. –  belisarius Apr 2 '13 at 5:17

The solution was already given by Daniel. Just a caveat: Your plot isn't showing the whole truth:

RegionPlot[\[Rho]yz < 0, {y, -50, 50}, {z, -50, 50}, PlotPoints -> 100]

Mathematica graphics

share|improve this answer
    
You are right! I already know this; i just plotted the contours very quick using speed against quality. But again the minimum distance which I want is the same. Am I right? –  Vaggelis_Z Apr 1 '13 at 16:57
    
@Vaggelis_Z Yep, Daniel's answer is right AFAIK –  belisarius Apr 1 '13 at 17:04
    
What do you mean? –  Vaggelis_Z Apr 1 '13 at 17:06
    
@Vaggelis_Z I mean that the FindMinimum[] used by @DanielLichtbau takes care of this issue. –  belisarius Apr 1 '13 at 17:10
    
But when I run the Do loop I get this message: FindMinimum::eit: "The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {0.00193054,11.0354,0}, is returned" –  Vaggelis_Z Apr 1 '13 at 17:12

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.