使用Sympy求解函数公切点二元微分方程仅得平凡解的问题排查

问题描述 #

现有名为F的双势阱函数，函数形态可参考双势阱曲线示例。需要在给定区间内求解函数的两个公切点(x1,F(x1))、(x2,F(x2))，对应约束方程为：

F'(x1)=F'(x2)                     (1)
F'(x1)*(x2-x1)+F(x1)=F(x2)        (2)

代码实现与异常 #

使用sympy.solve编写求解代码如下：

x1 = Symbol('x1')
F1 = F(x1,AA, BB, CC, DD, EE, FF)
x2 = Symbol('x2')
F2 = F(x2,AA, BB, CC, DD, EE, FF)
dF1=diff(F1,x1)
dF2=diff(F2,x2)
print(solve([dF1-dF2, F2-F1-dF1*(x2-x1)], [x1, x2]))

代码运行输出为[(x2, x2)]，仅返回x1=x2的平凡解，无有效非平凡解输出。

待解决疑问 #

• 求解时如何设置x1 != x2的约束，排除无效平凡解？
• 如何获取上述方程的数值求解结果？
• 如何将x1和x2的解限制在指定区间范围内？

正确性验证 #

测试其他求解逻辑均可得到正确结果，测试代码如下：

print(solve([y1+y2,y1-y2], [x1, x2]))
print(solve(dy1, x1))

打印各中间表达式确认正确性，执行输出如下：

-------------------------------------------------------------------
In [1]: x1 = Symbol('x1')
F1 = F(x1,AA, BB, CC, DD, EE, FF)
x2 = Symbol('x2')
F2 = F(x2,AA, BB, CC, DD, EE, FF)
dF1=diff(F1,x1)
dF2=diff(F2,x2)
-------------------------------------------------------------------
In [2]: F1
Out[2]: 𝑥1⋅(1−𝑥1)(−0.0501765666583873𝑥1+0.244249682150038(2𝑥1−1)5−0.535363152864316(2𝑥1−1)4+0.197479965493092(2𝑥1−1)3−0.366799918223659(2𝑥1−1)2−0.203076238788926)
(Out [2] 为函数F的正确表达式)
-------------------------------------------------------------------
In [3]:dF1
Out[3]:𝑥1⋅(1−𝑥1)(−2.93439934578927𝑥1+2.44249682150038(2𝑥1−1)4−4.28290522291453(2𝑥1−1)3+1.18487979295855(2𝑥1−1)2+1.41702310623625)−𝑥1(−0.0501765666583873𝑥1+0.244249682150038(2𝑥1−1)5−0.535363152864316(2𝑥1−1)4+0.197479965493092(2𝑥1−1)3−0.366799918223659(2𝑥1−1)2−0.203076238788926)+(1−𝑥1)(−0.0501765666583873𝑥1+0.244249682150038(2𝑥1−1)5−0.535363152864316(2𝑥1−1)4+0.197479965493092(2𝑥1−1)3−0.366799918223659(2𝑥1−1)2−0.203076238788926)
(Out [3] 为F的导函数，表达式正确)
-------------------------------------------------------------------
In [4]:print(solve([F1+F2,F1-F2], [x1, x2]))  
[(0.0, 0.0), (0.0, 1.00000000000000), (1.00000000000000, 0.0), (1.00000000000000, 1.00000000000000)]
(In [4] 求解结果正确)
-------------------------------------------------------------------
In [5]:print(solve(dF1, x1))
[0.143449671600321, 0.462174698289538, 0.744534434388284, 1.44258261390104, 0.573315370351492 - 0.261496971434863*I, 0.573315370351492 + 0.261496971434863*I]
(In [5] 单变量导函数求根结果正确)
-------------------------------------------------------------------
In [6]:print(solve([dF1-dF2, F2-F1-dF1*(x2-x1)], [x1, x2]))
[(x2, x2)]
(In [6] 为目标方程求解结果，仅返回平凡解，是核心问题)
-------------------------------------------------------------------

内容的提问来源于stack exchange，提问作者Er g Ch