Common tangent using python

Viewed 89

I am trying to find a common tangent to two curves using python but I am not able to solve it.
The equations to the two curves are complicated that involve logarithms.

Is there a way in python to compute the x coordinates of a tangent that is common to both the curves in general. If I have 2 curves f(x) and g(x), I want to find the x-coordinates x1 and x2 on a common tangent where x1 lies on f(x) and x2 on g(x). I am trying f'(x1) = g'(x2) and f'(x1) = f(x1) - f(x2) / (x1 - x2) to get x1 and x2 but I am not able to get values using nonlinsolve as the equations are too complicated.

I want to just find x-coordinates of the common tangent

Can anyone suggest a better way?

import numpy as np
import sympy
from sympy import *
from matplotlib import pyplot as plt


x = symbols('x')
a, b, c, d, e, f = -99322.50019502985, -86864.87072433547, -96876.05627516498, -89703.35055202093, -3390.863799999999, -20942.518


def func(x):
    y1_1 = a - a*x + b*x
    y1_2 = c - c*x + d*x

    c1 = (1 - x) ** (1 - x)
    c2 = (x ** x)
    y2 = 12471 * (sympy.log((c1*c2)))

    y3 = 2*f*x**3 - x**2*(e + 3*f) + x*(e + f)
    eqn1 = y1_1 + y2 + y3
    eqn2 = y1_2 + y2 + y3
    return eqn1, eqn2


val = np.linspace(0, 1)
f1 = sympy.lambdify(x, func(x)[0])(val)
f2 = sympy.lambdify(x, func(x)[1])(val)

plt.plot(val, f1)
plt.plot(val, f2)
plt.show()

I am trying this

x1, x2 = sympy.symbols('x1 x2')

fun1 = func(x1)[0]
fun2 = func(x2)[0]
diff1 = diff(fun1,x1)
diff2 = diff(fun2,x2)
eq1 = diff1 - diff2
eq2 = diff1 - ((fun1 - fun2) / (x1 - x2))

sol = nonlinsolve([eq1, eq2], [x1, x2])
2 Answers

the first thing that needs to be done is to reduce the formulas

for example the first formula is actually this:

formula = x*(1 - x)*(17551.6542 - 41885.036*x) + x*(1 - x)*(41885.036*x - 24333.3818) + 12457.6294706944*x + log((x/(1 - x))**(12000*x)*(1 - x)**12000) - 99322.5001950298
formula = (x-x^2)*(17551.6542 - 41885.036*x) + (x-x^2)*(41885.036*x - 24333.3818) + 12457.6294706944*x + log((x/(1 - x))**(12000*x)*(1 - x)**12000) - 99322.5001950298

# constants
a = 41885.036
b = 17551.6542
c = 24333.3818
d = 12457.6294706944
e = 99322.5001950298
f = 12000

formula = (x-x^2)*(b - a*x) + (x-x^2)*(a*x - c) + d*x + log((x/(1 - x))**(f*x)*(1 - x)**f) - e
formula = (ax^3 -bx^2 + bx - ax^2) + (x-x^2)*(a*x - c) + d*x + log((x/(1 - x))**(f*x)*(1 - x)**f) - e
formula = ax^3 -bx^2 + bx - ax^2 -ax^3 + ax^2 + cx^2 -cx + d*x + log((x/(1 - x))**(f*x)*(1 - x)**f) - e

# collect x terms by power (note how the x^3 tern drops out, so its easier).
formula =  (c-b)*x^2 + (b-c+d)*x  + log((x/(1 - x))**(f*x)*(1 - x)**f) - e

which is much cleaner and is a quadratic with a log term. i expect that you can do some work on the log term too, but this is an excercise for the original poster.

likewise the second formula can be reduced in the same way, which is again an excercise for the original poster.

From this, both equations need to be differentiated with respect to x to find the tangent. Then set both formulas to be equal to each other (for a common tangent).

This would completely solve the question.

I actually wonder if this is a python question at all or actually a pure maths question.....

The important point to note is that, since the derivatives are monotonic, for any value of derivative of fun1, there is a solution for fun2. This can be easily seen if you plot both derivatives.

Thus, we want a function that, given an x1, returns an x2 that matches it. I'll use numerical solution because the system is too cumbersome for numerical solution.

import scipy.optimize

def find_equal_value(f1, f2, x, x1):
    goal = f1.subs(x, x1)
    to_solve = sympy.lambdify(x, (f2 - goal)**2)  # Quadratic functions tend to be better behaved, and the result is the same
    sol = scipy.optimize.fmin(func=to_solve, x0=x1, ftol=1e-8, disp=False)  # The value for f1 is a good starting guess
    return sol[0]

I used fmin as the solver above because it worked and I knew how to use it by heart. Maybe root_scalar can give better results.

Using the function above, let's get some pairs (x1, x2) where the derivatives are equal:

df1 = sympy.diff(func(x)[0])
df2 = sympy.diff(func(x)[1])

x1 = 0.25236537  # Close to the zero derivative
x2 = find_equal_value(df1, df2, x, x1)
print(f'Derivative of f1 in x1: {df1.subs(x, x1)}')
print(f'Derivative of f2 in x2: {df2.subs(x, x2)}')
print(f'Error: {df1.subs(x, x1) - df2.subs(x, x2)}')

This results is:

Derivative of f1 in x1: 0.0000768765858083498
Derivative of f2 in x2: 0.0000681969431752805
Error: 0.00000867964263306931

If you want a x2 for several x1s (beware that in some cases the solver hits a value where the logs are invalid. Always check your result for validity):

x1s = np.linspace(0.2, 0.8, 50)
x2s = [find_equal_value(df1, df2, x, x1) for x1 in x1s]
plt.plot(x1s, x2s); plt.grid(); plt.show()
Related