如何使用Python求解含对数复杂曲线公切线的切点横坐标

问题分析 #

使用SymPy的nonlinsolve求解含对数、幂指项（$x^x$）的双变量非线性公切线方程组无法得到有效结果，核心原因有两点：

• 这类含超越项的方程组不存在闭式解析解，符号求解器本身不具备通用求解能力
• 原代码存在逻辑错误：求解第二条曲线取值时误用了第一条曲线的表达式，fun2 = func(x2)[0]实际是在求单条曲线的自切线，而非两条曲线的公切线

公切线的数学建模逻辑是完全正确的，需要满足两个约束：

1. 两切点处导数值相等：$f'(x_1) = g'(x_2)$
2. 两切点连线斜率等于切线斜率：$f'(x_1) = \frac{f(x_1)-g(x_2)}{x_1-x_2}$

针对这类无解析解的超越方程组，数值求解是唯一稳定高效的方案。

优化方案 #

核心优化点 #

• 利用对数运算法则化简幂指项：$\ln\left((1-x)^{1-x}\cdot x^x\right) = (1-x)\ln(1-x) + x\ln x$，避免先算高次幂再取对数导致的数值溢出问题
• 放弃符号计算链路，全部转为NumPy数值计算，计算速度提升2~3个数量级
• 采用中心差分法计算数值导数，无需手动推导复杂的导函数表达式
• 先在定义域内粗扫残差得到合理初值，再用scipy.optimize.fsolve做精确迭代求解，避免牛顿法初值敏感导致的发散

完整可运行代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# 固定模型参数
a, b, c, d, e, f_coef = -99322.50019502985, -86864.87072433547, -96876.05627516498, -89703.35055202093, -3390.863799999999, -20942.518

# 化简后的第一条曲线数值函数
def f1(x):
    y1 = a + (b - a) * x
    # 化简对数幂指项，提升数值稳定性
    y2 = 12471 * ((1 - x) * np.log(1 - x) + x * np.log(x))
    y3 = 2 * f_coef * x**3 - x**2 * (e + 3 * f_coef) + x * (e + f_coef)
    return y1 + y2 + y3

# 化简后的第二条曲线数值函数
def f2(x):
    y1 = c + (d - c) * x
    y2 = 12471 * ((1 - x) * np.log(1 - x) + x * np.log(x))
    y3 = 2 * f_coef * x**3 - x**2 * (e + 3 * f_coef) + x * (e + f_coef)
    return y1 + y2 + y3

# 中心差分法计算数值导数，精度满足工程需求
def num_deriv(func, x, eps=1e-8):
    return (func(x + eps) - func(x - eps)) / (2 * eps)

# 构造待求解的非线性方程组
def eq_system(vars):
    x1, x2 = vars
    df1 = num_deriv(f1, x1)
    df2 = num_deriv(f2, x2)
    line_slope = (f1(x1) - f2(x2)) / (x1 - x2)
    return [df1 - df2, df1 - line_slope]

# 第一步：在(0,1)定义域内粗扫描，找残差最小的点作为迭代初值
scan_x = np.linspace(0.01, 0.99, 50)
best_init = None
min_res = float("inf")
for x1g in scan_x:
    for x2g in scan_x:
        res = eq_system([x1g, x2g])
        total_res = abs(res[0]) + abs(res[1])
        if total_res < min_res:
            min_res = total_res
            best_init = [x1g, x2g]

# 第二步：fsolve精确迭代求解
x1_sol, x2_sol = fsolve(eq_system, best_init)

# 结果验证
tangent_k = num_deriv(f1, x1_sol)
print(f"f(x)上切点横坐标x1 = {x1_sol:.6f}")
print(f"g(x)上切点横坐标x2 = {x2_sol:.6f}")
print(f"一致性验证：f'(x1)={tangent_k:.4f} | g'(x2)={num_deriv(f2, x2_sol):.4f} | 连线斜率={(f1(x1_sol)-f2(x2_sol))/(x1_sol-x2_sol):.4f}")

# 可视化曲线与公切线
x_plot = np.linspace(0.001, 0.999, 200)
plt.plot(x_plot, f1(x_plot), label="Curve f(x)")
plt.plot(x_plot, f2(x_plot), label="Curve g(x)")
# 绘制公切线
tangent_x = np.linspace(0, 1, 200)
tangent_y = tangent_k * (tangent_x - x1_sol) + f1(x1_sol)
plt.plot(tangent_x, tangent_y, "--", color="red", label="Common tangent")
plt.scatter([x1_sol, x2_sol], [f1(x1_sol), f2(x2_sol)], color="red", zorder=5, s=30)
plt.legend()
plt.ylim(-120000, -70000)
plt.show()

运行说明 #

• 代码依赖numpy、matplotlib、scipy三个第三方库，可通过pip直接安装
• 若需要求解其他曲线的公切线，只需替换f1、f2的函数定义即可，其余求解逻辑通用
• 若曲线定义域不在(0,1)区间，修改粗扫描scan_x的取值范围即可适配

若存在多条公切线，可在粗扫描环节保留多个残差小于阈值的初值，分别代入fsolve求解即可得到所有根。

内容的提问来源于stack exchange，提问作者Sonic Sharma