用python求解多元多项式方程组

Question

编辑：我从中得到的参考方程中包含了几个错误。我已经修好了。现在解决方案可能真的有意义了！

当两层流体在地形上流动时，存在一个数 取决于流动的相对大小 速度和流体中的波速。在

它们被称为“超临界”、“次临界”和“临界”（即 前两个我在这里称之为“额外关键”）。在

以下方程式定义了临界点之间的边界线 （h，U0）参数空间中的超临界行为：

我想消除d_1c（即我不在乎它是什么）并找到 这些方程的解。在

简化因素：

• 我只需要给定的答案
• 我不需要精确的解，只需要一个解的概要 曲线，所以这可以通过解析或数值求解。在
• 我只想在区域（h，U0）=（0，0）到（0.5，1）上绘制。在

我想用热情中可用的模块来解决这个问题 分发（numpy，scipy，sympy），但真的不知道该去哪里 开始。它的消除变量d1c真正混淆 我。在

以下是python中的方程式：

def eq1(h, U0, d1c, d0=0.1):
    f = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0) ** 3) - 1
    return f

def eq2(h, U0, d1c, d0=0.1):
    f = 0.5 * (U0) ** 2 * ((d0 ** 2 / d1c ** 2) - (1 - d0) ** 2 / (1 - d1c - d0) ** 2) + d1c + (h - d_0)
    return f

我期望有一个解决方案，它有许多解决方案分支（不是 总是身体上的，但不用担心），看起来很粗糙 像这样：

我如何着手实施这一点？在

Answer 1

让我先来处理d1c的消去。假设您设法将第一个等式转换成d1c = f(U, h, d0)的形式。然后把它代入第二个方程，在U，h和{}之间有一定的关系。在d0固定的情况下，它为两个变量U和{}定义了一个方程，从中可以找到任何给定的h的{}。这似乎就是你所说的解决方案，基于你最后的草图。 坏消息是，从任何一个方程中得到d1c并不容易。好消息是你不需要。在

fsolve可以取一个方程组，比如两个依赖于两个变量的方程，然后给出解。在本例中：fix h（d0已经修复），并将其作为变量U0和d1c作为变量提供给d0。记录U0的值，对下一个值{}重复上述操作，依此类推。在

请注意，与@duffymo的建议相反，我建议使用fsolve，或者至少从它开始，只有当它耗尽蒸汽时才寻找其他解算器。在

一个可能的警告是，对于给定的h，U0您需要一个开始的猜测，fsolve需要一个开始的猜测，并且没有简单的方法告诉它收敛到解决方案分支之一。如果这是一个问题，请查看brentq解算器。在

另一种方法是观察您可以轻松地从系统中消除U0。这样，您将得到h和d1c的一个方程，为h的每个值求解{}，然后使用原始方程中的任一个来计算给定的d1c和{}。在

使用fsolve的示例：

>>> from scipy.optimize import fsolve
>>> def f(x, p):
...   return x**2 -p
... 
>>> fsolve(f, 0.5, args=(2,))
array([ 1.41421356])
>>> 

这里的args=(2,)是用来告诉fsolve的语法，如果f(x,2)=0，你真正想要解决的是0.5的值。在

Answer 2

您可以使用非线性解算器（如Newton Raphson或BFGS）来解同时存在的非线性方程组。他们对婚姻的起始条件和条件很敏感，所以需要一些照顾。在

Answer 3

半形式化地，您要解决的问题是：给定d0，对h和U0解逻辑公式“存在d1c，使得eq1（h，U0，d1c，d0）=eq2（h，U0，d1c，d0）=0”。在

有一种将公式化简为多项式方程“p（h，U0）=0”的算法，称为“量词消去”，它通常依赖于另一种算法“圆柱代数分解”。不幸的是，这还没有在sympy中实现。在

然而，由于U0很容易被排除，你可以用sympy做一些事情来找到你的答案。从

h, U0, d1c, d0 = symbols('h, U0, d1c, d0')
f1 = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0 * h) ** 3) - 1
f2 = U0**2 / 2 * ((d0 ** 2 / d1c ** 2) + (1 - d0) ** 2 / (1 - d1c - d0 * h)) + d1c + d0 * (h - 1)

现在，从f1中去掉U0，并在f2中插入值（我是“手动”而不是使用solve（）来获得更漂亮的表达式）：

^{pr2}$

f3只依赖于h和d1c。而且，由于它是有理分式，我们只关心它的分子何时变为0，所以我们得到了一个包含两个变量的多项式方程：

p3 = fraction(cancel(f3))

现在，对于一个给定的d0，你应该能够用数值的方法把p3.subs（d0.1）转化成h（d1c），把它插回到U0中，并把（h，U0）作为d1c的函数进行参数化绘图