Locate intersection points in a contour plot

Question

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.

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?

Answer 1

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

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

Answer 2

Ok, so all we have is

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

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
 ]

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

Answer 3

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