Finding the intersection point in a loop curve

I have two waves m and b. Each are scaled monotonically and equivalently by a third parameter x (from zero to unity in this case). I am looking for an efficient way to find the two x values where both m and b are equal.

Imagine that a plot of m versus b. When both m and b are equal at one point along the line, the plot gives a loop as shown in the attached. I am looking for a way to find the intersection point of the loop as the (m,b,x) values.

Before I wrack my head with developing double for-loop iterations over all m with all b, am I missing a simpler way to find what I need?

intersection loop plot

If m and b are truly scaled equivalently, you're looking for a zero in m-b, right, so FindLevel should do it? Or am I missing something?

@Tony -- I don't think it is this easy. I believe that just testing all algebraic combinations such as m + b, m - b, m/b, m*b for zeros or other constant values is not the answer.

The background might help. I have an equation for G = f(x) where G is a Gibb's energy of mixing a binary solution at mole fraction x. At the point where the binary solution separates into two phases (think oil and water), the plot of G versus x has a characteristic shape. See the plot attached. Focus only on the solid blue curve for G(x) (the dashed red is its derivative). The two points indicated by the cursors are the compositions at phase separation. The two points are contained on one line with the same slope m and intercept b.

I have extracted m(x) and b(x) along the entire curve G(x). The two points I want to find are x_lower and x_upper. I imagine I can brute force my way through two for-loops, one in the other, to test all (m,b) combinations to find the two x values where m_lower = m_upper and b_lower = b_upper.

When phase separation occurs, deconstruction from the G(x) curve generates a plot of m(x) versus b(x) that forms a loop. I am wondering if some "intersection" operation exists to find the (m,b) value at the closing point of the loop.

Gibbs energy at phase separation (363.69 KB)

Hi Jeremy,


Would an arc hull approach to the problem work?  Finding the enclosing shape and then test for derivative changes signifying a change from the linear region to the curved region.


Rather than approaching this geometrically, have you thought about using Optimize to minimize the total free energy? You would need some sort of penalty function to satisfy mass balance constraints.

Oh, you're calculating a solvus.

I take it you have functions for G(x) and G'(x), from these you can calculate the intercept b.

Unless you can accurately parameterize G vs b in a way that allows you to solve for an intercept, I think you need nested numerical solvers.

Here is what I would try:

EDIT: my proposed method was incorrect!

I wish I knew how CALPHAD or Thermocalc achieve this.



Perhaps also relevant:

For calculating divariant binary phase loops, I use funcfit to fit a wave {0,0}.

The fit function returns for point 1 the difference in chemical potential of component one in phases A and B, and for point 2 the difference in chemical potential of component two in phases A and B.

To determine the phase A side of the loop, the fit coefficients are an intensive variable (P or T, one or the other is fixed) and the compositional variable (X1) for phase B, and the dependent variable is X1 for phase A. In this way I determine P(x) (or T(x)) that describes one side of the loop.

I use structure-based funcfit because there are a lot of thermodynamic parameters (end-member and mixing parameters) that need to be passed through the fit function to the Gibbs function. Keeping track of everything without structures would be challenging.

Adapting the divariant loop method to a solvus calculation, maybe rewrite the fit function to return the difference in coefficients for the tangents at x1 and x2 (i.e. m2-m1 and b2-b1), and fit the {0,0} wave for x1 and x2 with 0<K0<1 and 0<K1<1 constraints.

@Tony -- Thanks. Yes, I am finding the points on a two phase separation region. While I have an analytical form for the Delta G of mixing, I take the numerical derivative of the G to obtain mu.

I think I may have a clever way to find the intersection point without a for-loop. I will post back once I have tested it.

Could you post the relevant equations for those of us who do not speak "phase separation"...

@AG -- I will post an experiment with the details once I have confirmed my method to find the loop intersection point.

My clever idea was not valid. I attach the experiment. Focus on the plot (m versus b). No loop is visible to report. Click on the non-ideal checkbox. A loop will appear in the (m versus b) graph. The loop intersection is the points 47 and 705 in the (m, b) waves. Notice where the points are on the (G, m) graph. The tie line that connects 47 and 705 have (nearly) the same slope and intercept. The goal is to find the two points on the common tie line.

I was mistaken. At the moment, I am using an analytical derivative equation to obtain mu(x). The analysis to find the loop intersection should not depend on having an analytical expression for mu = dG/dx.

Binary Solid Phase Diagram.pxp (59.17 KB)

Could you do it graphically? Sum up two black and white images (value 0 or 1) and look for the value 2.

Well, your method works if you devise a way to establish the intersection point.

Here is my method using FuncFit:

// this can be called from your update_curves() function
function GetPhaseCompositions(variable bbeta, variable m, variable n)
    make/free wToFit={0,0} 
    Make/free/T/N=2 T_Constraints
    T_Constraints[0] = {"K0 > 0","K0 < 1","K1 > 0","K1 < 1"}
    struct fitFuncStructure s
    make/free coefw = {0.1, 0.9}
    wave s.coefw = coefw
    s.bbeta = bbeta
    s.m = m
    s.n = n
// fit the wave {0,0} to find the X positions of the touching points of the common tangent
    FuncFit zeroes, coefw, wToFit /STRC=s/C=T_Constraints
    print coefw
// if the points are not unique, either there is no solvus
// or the input conditions may lie outside of the solvus

// Non-analytical calculation of the gradient of a tangent to the deltaGmix curve

// Analytical is always better, but this allows for arbitrary mixing
// models that have no obvious analytical function for dG/dx. Since Gmix
// functions are usually well-behaved this simple numerical approximation
// for dG/dx may suffice.
function calc_dDeltamixGdxApprox(variable bbeta, variable m, variable n, variable xx)
    variable dx = 1e-6 
    return (calc_DeltamixG(bbeta, m, n, xx+dx)-calc_DeltamixG(bbeta, m, n, xx-dx))/(2*dx)

structure fitFuncStructure
    Wave coefw
    variable x
    variable bbeta // variables bbeta, m, n as defined in update_curves()
    variable m
    variable n

// for point 0 of the fit wave returns the difference in the gradients of tangents at x=coefw[0] and x=coefw[1]
// for point 1 of the fit wave returns the difference in the intercepts of tangents at x=coefw[0] and x=coefw[1]
// ideally this function should be made robust against returning NaN.
function zeroes(Struct fitFuncStructure &s)
    variable m0 // gradient of tangent at x=coefw[0]
    variable m1 // gradient of tangent at x=coefw[1]
    variable b0 // intercept - value at X=0 - of tangent touching curve at x=coefw[0]
    variable b1 // intercept of tangent at x=coefw[1]
    m0 = calc_dDeltamixGdxApprox(s.bbeta, s.m, s.n, s.coefw[0])
    m1 = calc_dDeltamixGdxApprox(s.bbeta, s.m, s.n, s.coefw[1])
    if (s.x == 0)
        return m1-m0

    b0 = calc_DeltamixG(s.bbeta, s.m, s.n, s.coefw[0]) - m0 * s.coefw[0]
    b1 = calc_DeltamixG(s.bbeta, s.m, s.n, s.coefw[1]) - m1 * s.coefw[1]
    return b1-b0

EDIT: I added some comments to clarify how this works.


@olelytken -- Perhaps. I think whatever image approach is used would have to include a dilation or blur. The (m,b) combinations are nearly the same but not perfectly the same.

@tony -- I will have to test your approach to understand it better. I would especially want to know how to change the approach to handle the case where Delta_{mix}G is defined numerically and mu is correspondingly obtain by numerical differentiation rather than analytically.

Many thanks to you both.