# Locate intersection points in a contour plot

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?

• Take a look here: mathematica.stackexchange.com/q/275/12 – Szabolcs Aug 10 '16 at 12:12
• Can be any output produced by ContourPlot and exported through Mma. I guess this would more easily be an image file? – usumdelphini Aug 10 '16 at 12:18
• @Kuba In the case of interest this would be Graphics output then, sorry. – usumdelphini Aug 10 '16 at 12:25
• Sure, here it is - edited the question! – usumdelphini Aug 10 '16 at 12:38
• – Kuba Aug 10 '16 at 12:42

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

• EmptyRegion from 11.0.0, OP's code as copied gives a gap in a blue line that is why the result is empty. – Kuba Aug 10 '16 at 17:56
• 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. – Jason B. Aug 10 '16 at 18:03

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@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],
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],
PlotRange -> 10, Axes -> False, Frame -> True,
Epilog -> {
Red, PointSize[Large],
Point[{curve1[d1] /. sol1, curve1[d1] /. sol2}]
}
] 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}}

• "..., without using any other information on the functions plotted, other than what is seen in the picture." – Kuba Aug 10 '16 at 17:06
• @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. – george2079 Aug 10 '16 at 17:51