用python求解联立多元多项式方程

Question

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

当两层流体流过地形时，存在一个数 取决于流量的相对大小的不同解决方案 流体中的速度和波速。

这些被称为“超临界”、“亚临界”和“临界”（ 前两个我在这里称为“非常关键”）。

以下等式定义了临界值之间的边界线 和 (h, U0) 参数空间中的超临界行为：

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

简化因素：

• 我只需要 given 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

我期待一个具有多个解决方案分支的解决方案（不是 总是身体上的，但不用担心）并且看起来很粗略 像这样：

我该如何实施？

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 并制作 (h, U0) 的参数图为d1c 的函数。

Answer 2

让我先处理消除d1c 的问题。想象一下，您设法将第一个方程按摩成d1c = f(U, h, d0) 的形式。然后将其代入第二个等式，并在U、h 和d0 之间建立一定的关系。固定d0 后，这为两个变量U 和h 定义了一个方程，原则上您可以从中找到U 对于任何给定的h。根据您最后的草图，这似乎就是您所说的解决方案。 坏消息是从你的任何一个方程中得到d1c 并不容易。好消息是您不需要这样做。

fsolve 可以采用方程组，例如两个依赖于两个变量的方程并为您提供解决方案。在这种情况下：修复h（d0 已修复），并提供给您拥有的系统fsolve，将其视为变量U0 和d1c。记录U0 的值，重复h 的下一个值，以此类推。

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

一个可能的警告是，在给定h 的情况下，您期望U0 有多个解决方案：fsolve 需要开始猜测，并且没有简单的方法可以告诉它收敛到解决方案分支之一.如果这是一个问题，请查看brentqsolver。

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

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 值的起始猜测值的语法。

Answer 3

您可以使用 Newton Raphson 或 BFGS 等非线性求解器来求解联立的非线性方程。它们对基质的起始条件和调节很敏感，因此需要小心。