common tangent of double-well function

Hello,

I have a complex problem and I am looking for a way to solve it using existing functions.

Let say that I have a double-well function y=f(x). Both y and x are stored in arrays. What I would like to do is to find a common tangent of two points near the local minima (red line). From a mathematical point of view, there are only two-point existing that can have a common tangent (red line). Others crossing function with more than 2 points (secants), they are touching only one point (unique tangent) or they are inside part of the function (common secants).

Those two points are not necessarily local minimus in extremal cases this can be a descending function in most of its range and have only one minimum but still, there is only one line that touches only two-point without crossing others.

Is there any function that can give me both coordinates or just y or x values of these two and only two points. I was thinking of a function that will create those lines and check if only two points are touched. I will appreciate any idea, tip link, or anything that can help build such a piece of code.

[Matt J]

From a mathematical point of view, there are only two-point existing that can have a common tangent (red line).

That would be true if you had f() in functional form, but since you only have samples x,y there is an infinite spectrum of different functions passing through the points, and a correspondingly infinite spectrum of tangents. You need to define some sort of interpolation model to connect the points continuously.

Answer 1

[Matt J]

This might give you an acceptable discrete approximation, but note my comment.

p=polyfit([0,1,2,3,0.5, 2.5],[0 0 0 0 -1, -10],4); %fake data
x=linspace(-1,4,10000);
y=polyval(p,x);

%% Calculation begins here %%%

k=convhull(x,y);

[xk,yk]=deal(x(k),y(k));
j = diff(xk)>0;
[xk,yk]=deal(xk(j),yk(j));

interDist=vecnorm(  diff([xk(:),yk(:)],1,1) ,2,2);

[~,imax]=max(interDist);

x1=xk(imax); y1=yk(imax);      %the two tangent points (x1,y1) and (x2,y2)
x2=xk(imax+1); y2=yk(imax+1);

%% Plot the results %%%
close all
hold on
plot(x,polyval(p,x))
 pl=polyfit([x1,x2],[y1,y2],1);
plot(x,polyval(pl,x))
xlim([0,3])
hold off

[Mateusz Brzezinski]

I was thinking about sth like that.

I know that without fitting it may be impossible to do.

Isn't better to make two separate fittings near to the region that I expect it to be? Or one function may be enough? Hard to tell right? Also, I was thinking of what order of the polynomial use.

[Mateusz Brzezinski]

From an analytical point of view, this function doesn't have a certain equation. All inputs change for each x iteration so it's almost impossible to work out an equation - that is why I only have numerical approximated values.

[Matt J]

I'm not following you. There are no curve fitting steps in my calculation of the two tangent points. (I use polyfit only to generate example x,y data).

[Mateusz Brzezinski]

I used your code and I put an extreme case.

In this case that should be no common tangent so x1=x2

(why? - 2nd derivative of f(x) has only one root - so 1st derivative of f(x) has no local maxima and minima and function becomes single-well function). And after reducing the amount of data to fitting it can give a quite ok approximation. x1=0.452... and x2=0.459... so ok-ish. Do you think It can be improved if I divide the function into two parts (only data near local minima) and use polyfit on two new "different" functions?

Anyway thanks a lot it works better than whatever I tried so far.

[Mateusz Brzezinski]

" I'm not following you. There are no curve fitting steps in my calculation of the two tangent points. (I use polyfit only to generate example x,y data)."

I used polyfit to fit it. I simply replaced your fake data with my x and y.

[Matt J]

I used your code and I put an extreme case ... In this case that should be no common tangent so x1=x2

I don't understand what the result should be in the "extreme case". You mean when the function is convex, you want a horizontal tangent line through the global minimum? That's an easy modification.

k=convhull(x,y);

[xk,yk]=deal(x(k),y(k));
j = diff(xk)>0;
[xk,yk]=deal(xk(j),yk(j));

if numel(xk)==numel(x) %graph is convex

    [ymin,imin]=min(y);
    [x1,x2,y1,y2]=deal(x(imin), x(imin), ymin, ymin);
    
else
    
interDist=vecnorm(  diff([xk(:),yk(:)],1,1) ,2,2);
[~,imax]=max(interDist);

x1=xk(imax); y1=yk(imax);      %the two tangent points (x1,y1) and (x2,y2)
x2=xk(imax+1); y2=yk(imax+1);

end

[Mateusz Brzezinski]

Ok, let me tell you more precisely what I am working on. I am working on a modified phase diagram of two polymers. Binodal decomposition is calculated as a list of two points (those which have one common tanged) for the desire temperature range. f(x) depends on T so for each T I construct a corresponding double-well function. The extreme case means that I reached a peak of binodal decomposition so x1=x2 - and double-well becomes single-well.

Normally if you use Flory Huggins model the f(x) can be given as equation - I have a modified version and I can not develop one certain equation (purer numerical approach). That is why my function needs to be fitted somehow. As better fitting is as better results I will get. So I was just wondering if for most cases (that have two local minima) I can use two separate regional fittings in regions close to those minima. In your first graph let say the first x range is 0.2 to 0.8 and the second rage is 2.2 up to 2.8. How I can change this code to do this? So.. I take my x and y data points I am reducing the vector length x=0:0.1:1 so xnewleft = x(3:9, :) and xnewright=x(23:29, :) and I am doing the same for y ynewleft = y(3:9, :) and ynewright=y(23:29, :) and I am using two fittings p1=polyfit(xnewleft,ynewleft,n) and p2=polyfit(xnewright,ynewright,n) wher n = 2. This way i created two separated 2nd order polynomial function (parabols) that I would like to find common tangent line of. How I can change your code to do this?

[Mateusz Brzezinski]

Long story short how can I modify your code to find a common tangent of two parabolas like those?

[Matt J]

Combine the two parabolas into a single double-well curve,

y=min(yleft,yright);

Then you can apply my code to the resulting x,y. But I do not see why you think it is necessary or better to introduce quadratic approximations. My original method will work directly on the original y=f(x) and its accuracy is limited only be how finely you sample x. Conversely, if you make polynomial fits, your accuracy will be fundamentally limited based on how quadratically-shaped f() actually is.

[Mateusz Brzezinski]

It's hard to give one always working fitting. I am just trying to work out the best idea to do this. In my data sample, a calculated x1 value is shifted (to what is physically possible) when this curve is starting to flatten - look on a graph that I called "extremal case". Binodla points start to go inside spinodal what is fundamentally wrong by the physical definition - and yes this is only due to approximation and stepping.

[Haiying Yang]

Hi Matt,

I have a similar question about finding the common tangents of two local minima of a symmetric triple well function. The function is shown below:

The function is:

The following figures shows the common tangents I want to find out (drawn in the red line). I want to find out the common tangents near the specifc two local minima.

I tried your code using my function, however, it can only give me the two common tangents shown in the red dots in the following figure, which are not what I want. How can I get the two common tangents I want here? How to revise the code?@Matt J

Thank you very much.

Best,

[Matt J]

I would try narrowing the range of x.

x=linspace(-0.25,2,10000);

[Haiying Yang]

Thank you Matt, it works.

Answer 2

[John D'Errico]

This is actually not a difficult problem to solve, IF you think about what a tangent line means.

For example, consider the polynomial

syms x
y(x) = expand((x-1)*(x+1)*(x+2)*(x-2)*(x-2.5)*(x+2.5))

y(x) = 

(I don't do homework problems for people. Sorry. So this is my own function. Feel free to look at what I wrote and to use it. But that would force you to think about what I wrote, and how to modify it to solve a different problem.)

fplot(y,[-4,4]),ylim([-40,30])

The goal is to find lines that are doubly tangent to the curve. How would we define that in terms of mathematics? First, we need the derivative of the function y.

dy = diff(y,x)

dy(x) = 

What do points of common tangency have in common? The slope MUST be the same at both points. Can we write that in mathematics? Assume the x xoordinates of the points of common tangency are at the locations u and v.

syms u v
eq(1) = dy(u) == dy(v)

eq = 

That just says the slope of the curve at x=u mst be the same as the slope at x=v. But we need more than that. We can define a line by two points along the line. Think of them as (u,y(u)) and (v,y(v)). If the tngent line connects tose two points, then we must have this simple relation (THINK ABOUT WHAT IT MEANS!)

eq(2) = y(u) + (v - u)*dy(u) == y(v)

eq = 

Now we can just try using solve on the problem.

uvsol = solve(eq)

uvsol = struct with fields:
    u: [13×1 sym]
    v: [13×1 sym]

And that must fail (if you look at it, because we are effectively trying to solve an implicit polynomial of degree 13), so we must convert those results to floating point numbers to resolve them. Actually, we should be happy that solve finds these solutions.

uv = vpa([uvsol.u,uvsol.v]);

Again, 13 solutions were found overall. We want to exclude the cases where we have u==v. We also want to worry only about the cases where v > u, since we seem to have duplicates. That is, if we found the solution (u,v), then the solution (v,u) should be excluded, as it would be the same line. We can resolve both issues by including only those sets where u < v (thus excluding those where u>=v).

uv(uv(:,1) >= uv(:,2),:) = []

uv = 

That leaves 6 distinct solutions. We can plot them simply enough, as

UV = double(uv')

UV = 2×6

-1.4931 -2.2958 -2.3373 -2.2754 -1.4469 0.1251 1.4931 2.2958 -0.1251 1.4469 2.2754 2.3373

yfun = matlabFunction(y);
yUV = yfun(UV);
fplot(y,[-3,3]),ylim([-30,20])
hold on
plot(UV,yUV,'o-')
grid on

Those are the 6 lines of double tangency for this specific function. (At least for MY choice of function.) Change the function, and the solutions will be different. There may be different numbers of doubly tangent lines.

A nice feature of this solution is these lines are indeed the exact solutions to the problem, and they effectively encompass all possible solutions. (Exact to within floating point trash, that is.)

[Haiying Yang]

Hi John,

Thank you very much for your detailed explaination of the code and how you did. I got the common tangents for my function. The following is the figure.

Best,

[将实]

Hello，

im so sorry to ask this question because i'm beginner to use matlab. I'm trying to find the tangent line on my equation of Flory Huggins FH= x*log(x)/N+(1-x)*log(1-x)+chi*x*(1-x) so that we have N=2 , chi=2 for example, but i dont know why ur code doesn't work why i use to this function to find the tangente line. Could you explain it please?

[John D'Errico]

@将实

Can I explain what? Note that not all equations you will write down have a simple analytical solution.

Here is your function:

N = 2;
chi = 2;
syms x
F(x) = x*log(x)/N+(1-x)*log(1-x)+chi*x*(1-x);
fplot(F,[0,1])

Note that it is undefined for x < 0 and x > 1, so I restricted the plot to that interval.

Any solution must lie in that interval. Is there a line of common tangency? Well, yes.

syms u v
dF(x) = diff(F,x)

dF(x) = 

A line of common tangency must satisfy a few simple relations...

uv = solve(dF(u) == dF(v), F(v) == F(u) + dF(u)*(v - u),[u,v])

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

uv = struct with fields:
    u: 0
    v: 0

Solve found no solution. But these are NOT simple polynomial equations. So solve found nothing useful. Instead, use a numerical solver, vpasolve. It will need starting values.

uv = vpasolve(dF(u) == dF(v), F(v) == F(u) + dF(u)*(v - u),[u,v],[.1,.9])

uv = struct with fields:
    u: 0.059566504583593630872615139319929
    v: 0.8490591390394413942116043378796

fplot(F,[0,1],'b')
hold on
grid on
plot(double([uv.u,uv.v]),double([F(uv.u),F(uv.v)]),'or-')
hold off

[将实]

Hi John,

Thank you for such detailed explanation. This is very helpful to me. But I still want to ask why the solution I found using ''fsolve''' without using toolbox is that u and v are very close, that is to say, it can find the solutions but the solution it finds is not what I want.

[John D'Errico]

Think about it. Suppose in my example, I had chosen u==v. Is that a solution to that set of equations?

Now look at the starting values I chose. They were far apart. AND they were close to the solution I knew to exist.

Do you understand that a nonlinear system of equations can have multiple sets of solutions? That is not at all unusual. Some of those solutions may be trivial ones. In fact, I would expect this system of equations will have infinitely many pairs of trivial solutions, all where u==v. So if you started the solver with u near v, then you will get u==v, probably to within the convergence tolerance.

When you use a nonlinear solver, it is imperative that you choose starting values that are at least reasonable, else you will likely get garbage for a solution. That is possibly why your attempt failed.