Algorithms for finding the self-intersection of a plane curve?

Does anyone here happen to know about algorithms for finding the self-intersection points of a plane curve? For example, find the self-intersection (X, Y) = (0, 0) in the below plane curve.

Edit: Add an actual example. The scattered points are computed numerically. It is not straightforward to control the spacing between each scattered point. The distribution density of points around the intersection (marked by the red circle) can be improved though. Note also that by eye we can tell that the intersection is not at Fg = 0.

1 Like

detect bifurcation? You monitor the spectrum of the jacobian. BifurcationKit or HomotopyContinuation.jl should be able to do it

3 Likes

I’ve gone through BifurcationKit.jl documentation quickly and found it quite interesting. However, I have difficulties understanding the underlying math. The tutorials are all about ODEs and PDEs. Do you have any tutorials for simple functions (such as y^2 = x^3 + x^2). The self-intersection of the above example happens to be the zero point of this function. But the self-intersection can locate anywhere on the plane. In addition, my actual function is more complicated and contains not only polynomials but natural log functions.

This library may be a good option: IntervalRootFinding

1 Like

The implicit equation provided can be parametrized to: X(t) = t^2 - 1; Y(t) = t^3 - t

Then we can solve the system (X(t), Y(t)) = (X(u), Y(u)), whose non-trivial solution (t != u) is u = 1, which gives the crossing point (0, 0).

Disclaimer: not a mathematician, but the above solution is a physical one, we send two particles at different times on the same trajectory and check their collision.

1 Like

True for this simple example. But my actual function is more complicated and has no explicit analytical expression. So trying to find a numerical approach :slight_smile:

1 Like

Regarding the recently added scatterplot above, would it be possible to extract two sets of points close to the intersection, fit them independently with linear (or other) functions, and then compute the intersection on that basis?

I’ve gone through BifurcationKit.jl documentation quickly and found it quite interesting. However, I have difficulties understanding the underlying math.

It is to find solutions (x,p) with salar p to an equation F(x, p) = 0

The tutorials are all about ODEs and PDEs.

The first part of the ode problems is to find the zeros of a nonlinear function (not just polynomials) as function of a parameter. For example you could search for y as function of x and detect bifurcations.

Do you have any tutorials for simple functions (such as y^2 = x^3 + x^2). The self-intersection of the above example happens to be the zero point of this function. But the self-intersection can locate anywhere on the plane. In addition, my actual function is more complicated and contains not only polynomials but natural log functions.

1 Like

@rveltz is right that at a self intersection there’s a singularity. So if your curve is the solutions of f(x,y) = 0 then the gradient must be orthogonal to the two linearly independent tanget lines and hence

f_x = f_y = 0

at any self intersection. So you could try to solve the nonlinear system

F = \left ( \begin{array}{c} f_x^2 + f_y^2 \\ f \end{array} \right) = \left ( \begin{array}{ll} 0 \\ 0 \end{array} \right)

Any implementation of Newton’s method, like nsol.jl from SIAMFANLEquations.jl would do that. Even if you only have a function to evaluate f and no analytic expression, that’ll be ok with a numerical Jacobian for Newton as long as f is sufficiently smooth and you pay attention to the difference increments for the first (h) and second (\sqrt{h}) derivatives.

You could think of f(x,y) = 0 as an equation where x is the variable and y a parameter. Then if you have a a point on the curve, you could us BIfurcationKit.jl to do continuation and any singularities will be flagged by the software. Even then, you are not guaranteed to find all solutions.

The tough part in your question is to find all solutions. If you had a polynomial system, then you have a computational algebraic geometry problem and Bertini (google it) would work, but there is no Julia interface AFAIK.

2 Likes

This (the numerical solution) should be an easy task for GMT’s module gmtspatial but there is currently an issue that needs to be fixed (gmtspatial -Ii failes and output strange things · Issue #5889 · GenericMappingTools/gmt · GitHub)

In general, it can be framed as a kind of root-finding problem, as several other posters have suggested, and for which there are many numerical algorithms. Formulating the problem more precisely depends on how your curve is parameterized, which you haven’t specified clearly.

2 Likes
using BifurcationKit, LinearAlgebra, Setfield, Plots
const BK = BifurcationKit

F(x, p) = @. -p.y^2 + x^2 + x^3
J(x, p) = ForwardDiff.jacobian(  z->F(z,p), x)
par = (y = 0., )

opts = BK.ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds=0.001, maxSteps = 140, pMin = -1., saveSolEveryStep = 1, detectBifurcation = 3)

br, = continuation(F,J, [-1.], par, (@lens _.y),
	ContinuationPar(opts, ds = -0.001, maxSteps = 800, newtonOptions = NewtonPar(verbose = false, maxIter = 6), plotEveryStep = 10),
	recordFromSolution = (x,p) -> (x = x[1],),
	bothside = true,
	verbosity = 3, plot = true,
	)

plot(br)

The problem is that is gives the fold points as well.

3 Likes

Thanks! This looks great! Will see whether I can formulate my problem to use BifurcationKit.

Here, why do you define another ConinuationPar? Why the first argument is opts? What does it do?

Another good news, I have successfully used this continuation method to draw a binodal line in a ternary phase diagram, which is really awesome! Previously I have to collect all those binodal points by hand which is very cumbersome. It is a more useful application than finding the self-intersection of a plane curve which is only an intermediate step in finding the binodal points. Now, I can simply use BK to find the binodal points for me.

1 Like

Because I typed super fast…sorry. Only one is required

Now, I can simply use BK to find the binodal points for me.

Coool