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.

Let's consider the following simple case:

V = 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
h = 0.3;
rad = 6*h;

S1 = ContourPlot[V, {x, -2, 2}, {y, -2, 2}, Contours -> {h}, 
ContourShading -> False, ContourStyle -> {Darker[Green], Thick}, 
PlotPoints -> 100, PerformanceGoal :> "Speed", AspectRatio -> 1];
S2 = ContourPlot[x^2 + y^2, {x, -2, 2}, {y, -2, 2}, Contours -> {rad},
ContourShading -> False, ContourStyle -> {Blue, Thick, Dashed}, 
PlotPoints -> 100, PerformanceGoal :> "Speed", AspectRatio -> 1];
P0 = Show[{S1, S2}, FrameLabel -> {"x", "y"}, RotateLabel -> False, 
FrameStyle -> Directive[FontSize -> 20, FontFamily -> "Helvetica"], 
PlotRange -> 1.5, ImageSize -> 550]

which gives

chans

We see, that the blue dashed circle intersects the green contour at the three openings (exits). I would like to draw the symmetrical parts of the circle inside these exits. By the term symmetrical, I mean with reverse curvature. I added by hand and with red color in the plot the symmetrical parts.

Any ideas how to obtain this?

EDIT

A copy of the above code (in version 9) can be found here:code.

share|improve this question
    
@Anon See my edit where I provide a copy of the code in version 9. –  Vaggelis_Z Nov 10 '13 at 11:04
    
You're aware that plotting the contour of $x^2+y^2$ at $\mathrm{rad}$ results in a circle of radius $\sqrt{\mathrm{rad}}$, I hope? –  Rahul Nov 10 '13 at 11:17
    
@Anon This file was created in version 9.0.0 and works fine in v 6.0.0, v 7.0.1 and v 8.0.0 and with some modifications (removing the PerformanceGoal option) also in v 5.2.0! So, it should work in v 9.0.1. –  Vaggelis_Z Nov 10 '13 at 11:17
    
@RahulNarain Perfectly aware! –  Vaggelis_Z Nov 10 '13 at 11:19
    
Sorry, I had to restart the kernel. I have already given the question a +1 because I thought it was well written btw. –  Pickett Nov 10 '13 at 11:20

6 Answers 6

First we need to find the intersection points. It's more convenient to find them as angles around the circle, rather than $(x,y)$ coordinates, otherwise we might not be sure what order they're in.

r = Sqrt[rad];
angles = θ /. Solve[{(V /. {x -> r Cos[θ], y -> r Sin[θ]}) == h, 0 <= θ < 2 π}, θ]

{1.32755, 1.81404, 3.42195, 3.90843, 5.51634, 6.00283}

By inspection we see that the arcs we want are from the first to the second, then the third to the fourth, and the fifth to the sixth point. With a little geometry, we can define a function to draw the reflected arc between any two points on a circle.

reflectedArc[{x_, y_}, r_, {θ1_, θ2_}] :=
 Circle[{x, y} + r {Cos[θ1] + Cos[θ2], Sin[θ1] + Sin[θ2]}, r, {π + θ1, π + θ2}]

We apply this onto our pairs of points and produce the final figure.

arcs = Graphics[{Red, reflectedArc[{0, 0}, r, #] & /@ Partition[angles, 2]}];
Show[{S1, S2, arcs}, FrameLabel -> {"x", "y"}, RotateLabel -> False, 
 FrameStyle -> Directive[FontSize -> 20, FontFamily -> "Helvetica"], 
 PlotRange -> 1.5, ImageSize -> 550]

enter image description here

(By the way, I must point out that you could have defined S2 much more easily as Graphics[{Blue, Thick, Dashed, Circle[{0, 0}, Sqrt[rad]]}].)

share|improve this answer

This solution is different from the other solutions in that it does not use the fact that the blue curve is a circle. Instead, we search the interior of the blue curve for points which will lie on the boundary when reflected.

I recast your problem statement as

f1[x_, y_] := 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
f2[x_, y_] := x^2 + y^2;
h = 0.3; rad = 6*h;

S1 = ContourPlot[f1 @@ {x, y}, {x, -2, 2}, {y, -2, 2}, 
   Contours -> {h}, ContourShading -> False, 
   ContourStyle -> {Darker[Green], Thick}, PlotPoints -> 100, 
   PerformanceGoal :> "Speed", AspectRatio -> 1];
S2 = ContourPlot[f2 @@ {x, y}, {x, -2, 2}, {y, -2, 2}, 
   Contours -> {rad}, ContourShading -> False, 
   ContourStyle -> {Blue, Thick, Dashed}, PlotPoints -> 100, 
   PerformanceGoal :> "Speed", AspectRatio -> 1];

Then I found the points of intersection, and lines to be reflected across. This is a weak point in the solution, as it depends on manually identifying the lines to reflect across.

pts = {x, y} /. Solve[{f2[x, y] == rad, f1[x, y] == h}, {x, y}, Reals]; // Quiet
lines = Partition[pts, 2];

Define reflections geometrically (Projection is built-in, that was new to me)

Rejection[u_, v_] := u - Projection[u, v];
Reflection[{a_, b_}, c_] := a + Projection[c - a, b - a] - Rejection[c - a, b - a]

Then use my Reflection and the RegionFunction option

S3 = ContourPlot[f2 @@ Reflection[#, {x, y}], {x, -2, 2}, {y, -2, 2}, 
   RegionFunction -> Function[{x, y}, f2[x, y] < rad], 
   Contours -> {rad}, ContourShading -> False, 
   ContourStyle -> {Red, Thick, Dashed}, PlotPoints -> 100, 
   PerformanceGoal :> "Speed", AspectRatio -> 1] & /@ lines;

And finally combine plots

P0 = Show[{S1, S2, S3}, FrameLabel -> {"x", "y"}, 
  RotateLabel -> False, 
  FrameStyle -> Directive[FontSize -> 20, FontFamily -> "Helvetica"], 
  PlotRange -> 1.5, ImageSize -> 550]

Here is a picture with your function, as well as one with f2[x_,y_]:=x^4+y^6 Reflected Contours

share|improve this answer
    
Neat approach. FYI, ReflectionTransform and ReflectionMatrix are built in, too. –  Michael E2 Nov 10 '13 at 15:33
    
Thanks @MichaelE2, I haven't played with those enough to be comfortable with them. What I don't like is that they use the normals to the lines/surfaces which often means some more effort on my part to find the normals before I can use them. –  Timothy Wofford Nov 10 '13 at 15:42
    
Conveniently Cross[{a, b}] makes a vector normal to the 2D vector {a, b}. –  Michael E2 Nov 10 '13 at 15:51
    
Very nice! Computing S3 is a little slow, but you can make it faster by // Evaluateing the f2 @@ Reflection[#, {x, y}]. –  Rahul Nov 10 '13 at 19:37
    
If f1[x_, y_] := 1/2*(x^2 + y^2) - x*y^2; the method does not work and the reflected arcs are wrong! any suggestions? –  Vaggelis_Z Nov 18 '13 at 17:16

Let me join this party.

I propose a purely graphical approach (it looks like quite compact):

f1[x_, y_] := 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
f2[x_, y_] := x^2 + y^2;
h = 0.3; rad = 6 h;

Show[Normal@ContourPlot[f2[x, y] == rad, {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
    ContourStyle -> {Red,Thick}, RegionFunction -> (f1[#1, #2] < h &), MaxRecursion -> 4] /.
     Line[{p1_, pp___, p2_}] :> GeometricTransformation[Line[{p1, pp, p2}], 
       ReflectionTransform[Cross[p2 - p1], p2]], 
 ContourPlot[f1[x, y] == h, {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
  ContourStyle -> {Darker[Green], Thick}]]

Here Normal converts GraphicsComplex to Lines. Then I reflect this lines around ending points.

enter image description here

share|improve this answer
    
Rats! Beat me!. :) +1. –  Michael E2 Nov 10 '13 at 15:50
    
@MichaelE2 Thank you for tip with Cross! Now it looks a bit better. –  ybeltukov Nov 10 '13 at 15:55
    
I guess this is the only approach that both (i) works for things other than circles, and (ii) doesn't require manual inspection that the intersection points are in the right order. Nice! –  Rahul Nov 10 '13 at 19:43
    
Probably the most general solution..+1 –  PlatoManiac Nov 10 '13 at 19:55
    
@ybeltukov very instructive, compact. Learned a lot. –  ubpdqn Nov 11 '13 at 2:40

We can try as the following.

V = 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
h = 0.3; rad = 6*h; center = {0, 0};
pts = {x, y} /.NSolve[{x^2 + y^2 == 6*h, V == h}, {x, y},Reals];(*find intersections*)
midpts = (Mean /@ Partition[pts, 2]);(*Find mid points the black ones*)
dist =Norm /@ midpts;(*Distance from center*)
Off[Divide::infy];(*equations of straight line joining center and mid points*)
eqs = (y - #[[2]] == (x - #[[1]]) First@Ratios@(center - #) & /@ midpts) /. 
   a_ == ComplexInfinity -> (x == 0) /.ComplexInfinity == b_ -> (y == 0);
(*pattern check for straight lines x=0 and y=0*)
Off[NSolve::ratnz];(*Find equation of the reverse circles*)res = 
 MapThread[Total[({x, y}-First@(Select[{x, y} /.Chop@NSolve[{#1, Norm[#2 - {x, y}] == #3},
   {x, y},Reals], # != center &]))^2] == rad &, {eqs, midpts,dist}];

We find collinear point on the line joining the center and the mid point (in black) lying at the same distance from the mid points as the distance from the midpoint to the center. In res we have the equations of the reverse circles centered at those points. You can ContourPlot these circle equations now on a right x-domain to get the right segment. Red and black points are drawn on you last Show command using Epilog.

reverse = 
  MapThread[
   ContourPlot[#1, Evaluate@Flatten@{x, #2[[1]]}, {y, -2, 2}, 
     Contours -> {rad}, ContourShading -> False, 
     ContourStyle -> {Red, Dashed}, PlotPoints -> 100, 
     PerformanceGoal :> "Speed", AspectRatio -> 1] &, {res,Transpose /@ Partition[pts, 2]}];

enter image description here

share|improve this answer

This is not particularly clever but works:

V = 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
cp = ContourPlot[{V == 0.3, x^2 + y^2 == 1.8}, {x, -2, 2}, {y, -2, 2}];
sol = NSolve[{V == 0.3, x^2 + y^2 == 1.8}, {x, y}];
roots = {x, y} /. sol
rem[u_] := Module[{f}, f[x_] := (f[Reverse@x] = Sequence[]; x); f /@ u];
ans = rem[Nearest[roots, #, 2] & /@ roots];
cent[{{p_, q_}, {r_, s_}}] := 
 First@Cases[{x, y} /. 
    Chop[NSolve[{(x - p)^2 + (y - q)^2 == 
        1.8, (x - r)^2 + (y - s)^2 == 1.8}, {x, y}]], Except[{0, 0}]];
vec = {#, cent[#]} & /@ ans;

arcs = {#[[2]], 
    Arg[Complex @@@ (Map[Function[x, x - #[[2]]], #[[1]]])]} & /@ vec;
Show[cp, Graphics[Circle[#[[1]], Sqrt[1.8], #[[2]]] & /@ arcs]]

Essentially, just find the intersections, second circle centre and arcs then compose. Styling could be altered.

enter image description here

share|improve this answer

This is like Rahul's solution only imagine that we didn't realize that we could use trigonometry to find the origin of the reflected arc analytically, and we also didn't realize polar coordinates were more suited for the task hehe.

(* Find the points of intersection. *)
sol = Solve[{
    V == h,
    x^2 + y^2 == rad
    }, {x, y}, Reals];

(* Find three sets with two points in each. Group by how close they \
are to each other. *)
pts = FindClusters[{x, y} /. sol, 3];

(* Find the origin of the other circle that intersects each set of \
coordinates *)
origos = {x, y} /. First@FindInstance[{
                Norm[{x, y} - #[[1]]] == Norm[{x, y} - #[[2]]],
                Norm[{x, y} - #[[1]]] == Sqrt[rad],
                Norm[{x, y}] > Sqrt[rad]
                }, {x, y}, Reals] & /@ pts;

(* Create data list: {origin,intersection 1,intersection 2} *)

data = MapThread[{#, #2[[1]], #2[[2]]} &, {origos, pts}];

(* Generate the circle segments *)
circles = Circle[
          #[[1]], Norm[#[[1]] - #[[2]]], {
       ArcTan @@ (#[[2]] - #[[1]]),
       ArcTan @@ (#[[3]] - #[[1]])
      }
          ] & /@ data;

(* Draw. *)
Show[P0, Graphics[{Red, circles}] ]

plot

share|improve this answer

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.