6

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

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

临界流

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

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

eq1

eq2

我想消除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

我期待一个具有许多解决方案分支的解决方案（并不总是物理的，但不要担心）并且看起来大致如下：

关键制度图

我该如何实施呢？

4

3 回答 3

4

半正式地，您要解决的问题如下：给定 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。

于 2012-12-11T21:51:17.000 回答
3

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

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

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

U0一个可能的警告是，您期望对于given有多个解决方案h：fsolve需要一个初始猜测，并且没有简单的方法可以告诉它收敛到一个解决方案分支。如果这被证明是一个问题，请查看brentq求解器。

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

使用示例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。

于 2012-12-11T19:21:00.687 回答
0

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

于 2012-12-11T15:59:22.110 回答