使用 sympy 的符号方程的正确命令序列

Question

注意：我是 sympy 的新手，并试图弄清楚它是如何工作的。

我现在拥有的： 我确实得到了正确的解决方案，但需要 35 - 50 秒。

目标： 通过定义一次符号方程然后将其与不同的变量一起使用来加快计算速度。

设置： 我需要为循环的每次迭代计算一个多项式 G(t)（t = 6 个根）。 （总共 220 次迭代） G(t) 有 6 个其他变量，这些变量在每次迭代中计算出来并且是已知的。 这些变量在每次迭代中都不同。

第一次尝试（慢）： 我只是将每一个都放入一个 python 函数中，在其中我象征性地定义了 Gt 并求解了 t。 它运行了大约 35 - 40 秒。每次迭代都会调用 function_G。

def function_G(f1, f2, a, b, c, d):
    t = sp.symbols('t')
    left = t * ((a * t + b)**2 + f2**2 * (c*t+d)**2)**2 
    right = (a*d-b*c) * (1+ f1**2 * t**2)**2 * (a*t+b) * (c*t+d)
    eq = sp.expand(left - right)
    roots = sp.solveset(Gt, t)
    return roots

• 然后有人给了我一个提示：

您应该只需要（象征性地）求解一次多项式的系数，作为预处理步骤。之后，在处理每次迭代时，您只需计算多项式系数，然后求解根。

• 我要求对此人进行澄清：

所以我定义了函数 g(t)，然后使用 sympy.expand 计算出所有括号/指数，然后 sympy.collect 通过 t 的幂来收集项。最后，我在 collect 的输出上使用了 .coeff 来获取系数以输入 numpy.root。

第二次尝试： 为了遵循建议，我首先以符号方式定义了一个 G(t) 并将其传递给运行循环的函数及其符号参数。因此，函数constructGt() 只被调用一次。

def constructGt():
    t, a, b, c, d, f1, f2 = sp.symbols('t a b c d f1 f2')
    left = t * ((a * t + b)**2 + f2**2 * (c*t+d)**2)**2 
    right = (a*d-b*c) * (1+ f1**2 * t**2)**2 * (a*t+b) * (c*t+d)
    gt = sp.Eq(left - right, 0)
    expanded = sp.expand(gt)
    expanded = sp.collect(expanded, t)

    g_vars = {
      "a": a,
      "b": b,
      "c": c,
      "d": d,
      "f1": f1,
      "f2": f2
    }

    return expanded, g_vars

然后在每次迭代中，我都传递函数及其参数来获取根：

#Variables values:
#a = 0.00011713490404073987 
#b = 0.00020253296124588926 
#c = 4.235688216068313e-07 
#d = 0.012262546040805029 
#f1= -0.012553203944721956 
#f2 = 0.018529776776949003

def function_G(f1_, f2_, a_, b_, c_, d_, Gt, v):
    Gt = Gt.subs([(v['a'], a_), (v['b'], b_), 
                  (v['c'], c_), (v['d'], d_),
                  (v['f1'], f1_), (v['f2'], f2_)])

    roots = sp.solveset(Gt, t)
    return roots 

但它在 56 秒左右变得更慢。

问题： 我不明白我做错了什么？我也不明白这个人是如何使用 .coeff() 然后 np.roots 的结果。

Answer 1

即使您的 f1 和 f2 在变量中是线性的，您也在使用四次多项式并且其根很长。如果这是单变量，那么最好只使用表达式，在已知解的某个常量值处求解它，然后使用该值和相对接近旧常量的新常量并使用 nsolve 获取下一个根目录。如果您对多个解决方案感兴趣，您可能必须使用nsolve 分别“关注”每个根...但我认为您会对整体性能感到更满意。使用real_roots 是另一种选择，尤其是当表达式只是某个变量中的多项式时。

鉴于您正在使用四次方，您应该牢记这一点：一般解决方案是如此冗长和复杂（除了非常特殊的情况），使用一般解决方案并替换值是无效的，因为它们是已知的。但是，求解数值非常容易，而且速度要快得多：

首先创建将要替换的值的符号表达式；使用了无假设的“香草”符号：

t, a, b, c, d, f1, f2 = symbols('t a b c d f1 f2')
left = t * ((a * t + b)**2 + f2**2 * (c*t+d)**2)**2 
right = (a*d-b*c) * (1+ f1**2 * t**2)**2 * (a*t+b) * (c*t+d)
eq = left - right

接下来，定义一个替换字典以替换到表达式中，注意 dict(x=1) 创建 {'x': 1} 并且当它与 subs 一起使用时，将为“x”创建一个普通符号：

reps = dict(
a = 0.00011713490404073987 ,
b = 0.00020253296124588926 ,
c = 4.235688216068313e-07 ,
d = 0.012262546040805029 ,
f1= -0.012553203944721956 ,
f2 = 0.018529776776949003)

计算表达式的实根：

from time import time
t=time();[i.n(3) for i in real_roots(eq.subs(reps))];'%s sec' % round(time()-t)
[-11.5, -1.73, 8.86, 1.06e+8]
'3 sec'

找出表达式的所有 6 个根，但只取实部：

>>> roots(eq.subs(reps))
{-11.4594523988215: 1, -1.73129179415963: 1, 8.85927293271708: 1, 106354884.4365
42: 1, -1.29328524826433 - 10.3034942999005*I: 1, -1.29328524826433 + 10.3034942
999005*I: 1}
>>> [re(i).n(3) for i in _]
[-11.5, -1.73, 8.86, 1.06e+8, -1.29, -1.29]

更改一个或多个值，然后再做一次

reps.update(dict(a=2))
[i.n(3) for i in real_roots(eq.subs(reps))]
[-0.0784, -0.000101, 0.0782, 3.10e+16]

循环更新值：

>>> a = 1
>>> for i in range(3):
...     a += 1
...     reps.update(dict(a=a))
...     a, real_roots(eq.subs(reps))[0].n(3)
...
(2, -0.0784)
(3, -0.0640)
(4, -0.0554)

注意：当使用roots 时，实根将按排序顺序排在最前面，然后虚根将按共轭对出现（否则不按任何给定顺序）。