使用Python解决多元多项式方程组问题

Question

编辑：我从参考资料中得到的方程式有几个错误。我已经在这里进行了修正。现在解决方案可能更有意义！

当两层流体流过地形时，存在许多不同的解，取决于流速和流体中波速的相对大小。

这些被称为“超临界”，“次临界”和“临界”（我在这里将前两个称为“超临界”）。

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

我想消除d_1c（即我不关心它是什么），并找到在(h，U_0)中这些方程的解。

简化因素：

• 我只需要给定d_0的答案
• 我不需要精确的解决方案，只需要解决方案曲线的概述，因此可以通过解析或数值方法来解决。
• 我只想在区域（h，U0）=（0,0）到（0.5,1）绘制图表。

我希望使用Enthought发行版中可用的模块（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

我希望你能提供一个有多个解决方案分支的解决方案（不一定是物理上的），并且大概看起来像这样： 我该如何实现它？

- aaren

---

我建议您在scicomp.stackexchange.com上提问。我怀疑那里的专业知识更符合您的领域。 - MRocklin

你能具体说明一下你所说的“消除dlc”的意思吗？你是说它应该完全从最终解决方案中取消吗？ - asmeurer

我的意思是我不需要知道d1c。我原本以为我有4个变量的两个方程式（U0，h，d1c，d0），如果我设定其中一个（d0）并消除另一个（d1c），那么我就会剩下两个关于（U0，h）的方程式，从而能够找到表达为U0和h之间关系的解决方案。 - aaren

这个问题已经发布在多个网站上：http://scicomp.stackexchange.com/questions/4843/solving-simultaneous-multivariate-polynomial-equations-with-python - Paul

我已经删除了关于科学计算的问题。我会在 Stack Overflow 上继续发布一段时间。 - aaren

其实，我认为在其他 SE 网站上进行交叉发布是不被鼓励的。如果你想的话，你大概可以把这个问题移到那里。 - asmeurer

Answer 1

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

有一种算法可以将该公式简化为多项式方程“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()函数来获得更漂亮的表达式）：

U2_val = ((f1 + 1)/U0**2)**-1
f3 = f2.subs(U0**2, U2_val)

f3仅依赖于h和d1c。由于它是一个有理分式，我们只关心其分子何时为0，因此我们得到了一个在2个变量中的单一多项式方程：

p3 = fraction(cancel(f3))

现在，对于给定的d0，你应该能够通过数值方法反演p3.subs(d0, .1)来得到h（d1c），将其插入U0中，并作为d1c的函数进行参数绘图（h，U0）。请保留html标记。

- Ronan Lamy

Answer 2

让我先处理消除 d1c。想象一下，您已经成功将第一个方程式变形为 d1c = f(U, h, d0) 的形式。然后，您会将此代入第二个方程式，并得到 U、h 和 d0 之间的某种关系。在固定 d0 的情况下，这定义了一个包含两个变量 U 和 h 的单个方程式，从中您可以原则上找到任何给定 h 的 U。这似乎是您根据最后的草图所称的解决方案。

坏消息是，从您的任一方程式中获得 d1c 不容易。好消息是，您不需要这样做。 fsolve 可以接受一个方程组，例如两个依赖于两个变量的方程式，并给出解决方案。在这种情况下：固定 h（因为 d0 已经固定），并将您拥有的系统作为变量 U0 和 d1c 输入到 fsolve 中。记录 U0 的值，重复下一个 h 的值，以此类推。

请注意，与 @duffymo 的建议相反，我建议使用 fsolve，或者至少从它开始，并在它失去效力时寻找其他解算器。

可能的一个警告是，您期望在给定 h 的情况下有多个 U0 的解： fsolve 需要一个起始猜测，并且没有简单的方法告诉它收敛到其中一个解分支。如果这成为问题，请查看 brentq 解算器。

另一种方法是观察到您可以轻松地从系统中消除 U0。这样，您将得到一个包含 h 和 d1c 的单个方程式，将其解决为每个 h 的 d1c，然后使用您原始方程式之一来计算给定 d1c 和 h 的 U0。

使用 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是x值的起始猜测。在这里，args=(2,)是语法。

- ev-br

Answer 3

您可以使用牛顿拉弗森法或BFGS等非线性求解器来解决同时存在的非线性方程。由于它们对于初始条件和矩阵的条件敏感，因此需要一些仔细的操作。

- duffymo

---

如果你通过分母相乘，它们不是多项式吗？还是它们实际上是非线性多项式？ - aaren

任何高于一次的多项式都是非线性的，这是定义。我看到一个dc1^2项，所以这是二阶的。 - duffymo