首页 文章

用python求解同时多元多项式方程

提问于
浏览 次
6

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

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

critical-flow

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

以下等式定义了（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

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

critical-regime-diagram

我该如何实现呢？

python math scipy polynomial-math sympy

3 回答

  • 3

    半正式地，您要解决的问题如下：给定d0，求解逻辑公式“存在d1c使得eq1（h，U0，d1c，d0）= eq2（h，U0，d1c，d0）=对于h和U0，为0“ .

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

    但是，由于U0很容易被淘汰，因此您可以通过同情来找到答案 . 从...开始

    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 .

    回复于
  • 0

    让我先讨论消除 d1c . 想象一下，你设法按下第一个等式，形成 d1c = f(U, h, d0) . 然后你将它替换为第二个等式，并在 U ， h 和 d0 之间有一定的关系 . 在 d0 固定的情况下，这为两个变量 U 和 h 定义了一个等式，原则上，对于任何给定的 h ，您可以从中找到 U . 根据您的上一个草图，这似乎是您所谓的解决方案 . 坏消息是从你的任何一个方程中得到 d1c 并不容易 . 好消息是你不需要 .

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

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

    一个可能的警告是，你希望有一个以上的解决方案，因为 h ： fsolve 需要一个初始猜测，并且没有简单的方法告诉它收敛到其中一个解决方案分支 . 如果结果证明是个问题，请查看 brentq 求解器 .

    另一种方法是观察您可以轻松地从系统中消除 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 的值的起始猜测，你实际想要解决的问题 .

    回复于
  • 4

    您可以使用Newton Raphson或BFGS等非线性求解器来求解同时的非线性方程 . 他们对初始条件和对基础的调节敏感，因此需要一些护理 .

    回复于

相关问题