scipy 使用python的公切线

eqzww0vc  于 2022-11-10  发布在  Python
关注(0)|答案(2)|浏览(120)

我试图用python找到两条曲线的公切线,但我无法解决它。
这两条曲线的方程很复杂,涉及对数。
在python中有没有一种方法来计算切线的x坐标,该切线通常是两条曲线的公共切线。(x)和g(x)、我想找到x1和x2在一条公切线上的x坐标,其中x1位于f上我尝试f '(x1)= g'(x2)和f '(x1)= f(x1)- f(x2)/(x1- x2)来得到x1和x2,但是由于方程太复杂,我无法使用非线性求解来得到值。

我只想求出公切线的x坐标

有谁能提出一个更好的方法吗?

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()

我在试这个

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])
a7qyws3x

a7qyws3x1#

首先需要做的是简化公式
例如,第一个公式实际上是:

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

这是一个更清晰的二次项,带有对数项。我希望你也可以在对数项上做一些工作,但这是原帖的节选
同样,第二个公式也可以用同样的方法简化,这也是原始发布者的摘录
由此,两个方程都需要对x求导以求出正切,然后使两个公式相等(对于公切线)。
这就彻底解决了问题。
我真的想知道这是一个Python的问题,还是一个纯粹的数学问题。

of1yzvn4

of1yzvn42#

需要注意的是,由于导数是单调的,对于fun1的任何导数值,fun2都有一个解,如果绘制两个导数,这一点很容易看出。
因此,我们需要一个函数,给定一个x1,返回与之匹配的x2。我将使用数值解法,因为系统对于数值解法来说太麻烦了。

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]

我用fmin作为上面的解算器,因为它很好用,而且我知道如何使用它。也许root_scalar能给予更好的结果。
使用上面的函数,我们得到一些导数相等的对(x1, x2)

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)}')

这个结果就是:

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

如果您需要几个x1x2(请注意,在某些情况下,求解器会遇到日志无效的值。请始终检查结果的有效性):

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()

相关问题