2
$\begingroup$

I have the following output from ContourPlot, and I would like to be able to find what the intersection points of the two curves are, without using any other information on the functions plotted, other than what is seen in the picture.

enter image description here

The image was produced using the following code:

r[x_, y_] := Sqrt[x^2 + y^2];
Θ[x_, y_] := ArcTan[x, y];

f[x_, y_, t_] := 
 1/2 - BesselK[0, r[x, y]]/(
  2 BesselK[0, 1]) + ϵ1[t] BesselK[0, r[x, y]]/
   BesselK[0, 1] + 
  D1[t] BesselK[1, r[x, y]]/(2 BesselK[1, 1])Cos[Θ[x, y] + ζ1[t]]

g[x_, y_, t_] := 
 Sin[Θ[x, y] + ζ2[t]]/
  r[x, y] + ϵ2/2 r[x, y]^(-1/2) Cos[Θ[x, y]/2]

ϵ1[t_] = -0.5 + 10^-6 t;
ζ2[t_] = π/3 + 10^-4 t;
ζ1[t_] = π/3 + 10^-6 t;
ϵ2 = 0.5;
D1[t_] = 0.5 - 10^-3 t;

 ContourPlot[{g[x, y, 0] == 0, f[x, y, 0] == 0}, {x, -10, 10}, {y, -10, 10}]

How do I go about doing this?

$\endgroup$
5
$\begingroup$

If plot is defined as your ContourPlot, then you can get the intersection points by extracting the Line objects from the plot and using RegionIntersection

plot // Normal // Cases[#, Line[_], Infinity] & // Apply@RegionIntersection

(* Point[{{-1.02539, 1.31753}, {0.379181, -1.26021}}] *)

and verify the results

Show[plot, Graphics[{Red, PointSize@Large, %}]]

Mathematica graphics

Version 11

For some reason version 11 will produce a disconnected line for the first contour, and so when you extract the lines you have more than 2, and there is no point where they all intersect and so the result is an EmptyRegion. A workaround is to find the points where any two lines intersect.

plot // Normal // Cases[#, Line[_], Infinity] & // Subsets[#, {2}] & //
   Map@Apply@RegionIntersection // DeleteCases[_EmptyRegion]
(* {Point[{{0.379181, -1.26021}}], Point[{{-1.031, 1.31378}}]} *)
$\endgroup$
  • $\begingroup$ EmptyRegion[2] from 11.0.0, OP's code as copied gives a gap in a blue line that is why the result is empty. $\endgroup$ – Kuba Aug 10 '16 at 17:56
  • $\begingroup$ Looks like ContourPlot produces a single line for g[x, y, 0] == 0 in version 10, but in version 11 it's two unconnected lines. $\endgroup$ – Jason B. Aug 10 '16 at 18:03
4
$\begingroup$

Ok, so all we have is

plot = ContourPlot[{g[x, y, 0] == 0, f[x, y, 0] == 0}, {x, -10, 10}, {y, -10, 10}]

Mathematica graphics

In order to view all of the information we have about the curves we can use FullForm[plot]. We see that we have all of the points that make up the lines, and that they are wrapped by Line. We can extract the points like this:

lines = Cases[Normal[plot], _Line, Infinity];
{pts1, pts2} = lines /. Line -> Identity;

We now have the points corresponding to each of the lines. In order to find the intersection between the lines it would be suitable to linearly interpolate between the points. In order to interpolate with Interpolation we need to introduce an independent variable, and I will introduce the distance along the curve as that independent variable.

createInterpolation[pts_] := Interpolation[Transpose[{
    Prepend[0]@Accumulate@Map[Norm, Differences@pts],
    pts
    }], InterpolationOrder -> 1]
curve1 = createInterpolation[pts1];
curve2 = createInterpolation[pts2];

We can plot the interpolated functions to see that we got it right:

Show[
 ParametricPlot[curve1[d], {d, 0, 23.9}],
 ParametricPlot[curve2[d], {d, 0, 9.48}, 
  PlotStyle -> ColorData[97][2]],
 PlotRange -> 10, Axes -> False, Frame -> True
 ]

Mathematica graphics

Now that we have functions that are continuous over the relevant domain, we can use NMinimize to find the intersections.

sol1 = Last@NMinimize[{
    Norm[curve1[d1] - curve2[d2]]^2,
    0 < d1 < 23.9 && 0 < d2 < 9.49
    }, {d1, d2}]

{d1 -> 12.2012, d2 -> 4.63393}

To find the other intersection we can change the constraints to exclude the intersection that has already been found:

sol2 = Last@NMinimize[{
    Norm[curve1[d1] - curve2[d2]]^2,
    0 < d1 < 10 && 0 < d2 < 9.48
    }, {d1, d2}]

{d1 -> 8.97608, d2 -> 9.48}

Finally we can plot the solutions to gauge whether they are correct:

Show[
 ParametricPlot[curve1[d], {d, 0, 23.9}],
 ParametricPlot[curve2[d], {d, 0, 9.48}, 
  PlotStyle -> ColorData[97][2]],
 PlotRange -> 10, Axes -> False, Frame -> True,
 Epilog -> {
   Red, PointSize[Large], 
   Point[{curve1[d1] /. sol1, curve1[d1] /. sol2}]
   }
 ]

Mathematica graphics

$\endgroup$
1
$\begingroup$

If you need precision results I'd use this:

FindMinimum[f[x, y, 0]^2 + g[x, y, 0]^2, {{x, -1}, {y, 1}}]
FindMinimum[f[x, y, 0]^2 + g[x, y, 0]^2, {{x, 0}, {y, -1}}]

where the initial guesses are taken from the contour plot.

{3.51393*10^-25, {x -> -1.02897, y -> 1.31148}}

{5.69013*10^-19, {x -> 0.378089, y -> -1.26037}}

$\endgroup$
  • 1
    $\begingroup$ "..., without using any other information on the functions plotted, other than what is seen in the picture." $\endgroup$ – Kuba Aug 10 '16 at 17:06
  • $\begingroup$ @kuba, fair enough, but someone else might encounter the question and not have that restriction. The results taken from the plot are off by ~0.005, which is obviously good enough if you just want to put a marker on the plot. $\endgroup$ – george2079 Aug 10 '16 at 17:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.