5

我使用的是 SymPy 0.7.3 版本，在使用 dsolve 函数时遇到了一些问题。当输入方程的参数太多时，dsolve 似乎有困难。

我尝试解决以下等式：

from sympy import *
p = Function('p')
t, u1, u2, u3, u4, u5 = symbols('t u1 u2 u3 u4 u5')
eq = Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2)))
eq
     Out: Derivative(p(t), t) == (-u3 - u4)*p(t) + exp(-t*(u1 + u2))
%time dsolve(eq)

并得到：

CPU times: user 213.11 s, sys: 0.00 s, total: 213.11 s
Wall time: 213.12 s
       p(t) == (C1 + Piecewise((t*u1/(u1 + u2) + t*u2/(u1 + u2), u3 == u1 + u2 - u4), (-exp(t*u3)*exp(t*u4)/(u1*exp(t*u1)*exp(t*u2) + u2*exp(t*u1)*exp(t*u2) - u3*exp(t*u1)*exp(t*u2) - u4*exp(t*u1)*ex
p(t*u2)), True)))*exp(-t*(u3 + u4))

（耗时 213.12 秒！）

然后我用 u5 替换了 u1+u2：

eq = Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))).subs(u1+u2, u5)
eq
     Out:Derivative(p(t), t) == (-u3 - u4)*p(t) + exp(-t*u5)
%time dsolve(eq)

并得到：

CPU times: user 1.62 s, sys: 0.00 s, total: 1.62 s
Wall time: 1.62 s
        p(t) == (C1 + Piecewise((t, u3 == -u4 + u5), (exp(t*u3)*exp(t*u4)/(u3*exp(t*u5) + u4*exp(t*u5) - u5*exp(t*u5)), True)))*exp(-t*(u3 + u4))

（只有 1.62 秒！）

我尝试过使用不同的提示，但没有帮助......

我还注意到，在更复杂的函数中，dsolve 会崩溃，但是在替换一些常量参数时，它运行得很快。

你知道这种现象的原因是什么吗？有没有办法自动解决这个问题？

4

2 回答 2

2

我注意到 Mathematica 存在同样的问题，限制表达式中符号的数量非常有利。我认为这样做的原因是大多数符号计算工具首先机械地应用通用配方来解决问题，然后尽可能简化结果。

在一般解决方案中，简化器可能很难恢复某些符号仅在给定组合中出现，因此可以将其视为单个符号。所以简化器不必要地处理更大的搜索空间。此外，它必须确保正确处理所有边界情况（u1可能是 < 0、> 0、= 0、复杂、...？）

我发现在将问题交给符号求解器之前尽可能多地进行简化是非常有利的。手动将变量分组在一起（如您的示例中）是一种有用的技术。另一种技术是进行归一化，或者将一个参数设置为 1。例如，在处理多项式时a x^2 + b x + c，大多数时候问题x^2 + B x + C与我们等价。（因为我们确实确定了a != 0，但是忘记告诉求解器，对吗？...）但是对于求解器来说，如果符号变量的数量减少 1，则会产生很大的不同。

在某些时候，这些符号求解器肯定会变得足够聪明，可以在尝试解决问题之前将变量组合在一起。但是，目前看来，仍然需要人为干预。

另一方面，很难想象符号求解器会自动识别更复杂的转换来简化问题，例如从笛卡尔坐标切换到极坐标，或者对变量进行更改，例如l=x+iy, r=x-iy，这些并不明显，但已知对某些问题有利。

更新

看来我们可以为这个问题做的最好的事情就是设置a = u3 + u4和a + b = u1 + u2。除了将参数数量从 4 个减少到 2 个之外，特殊情况现在出现 when b == 0，因此我们可以很容易地指示求解器忽略它： 

from sympy import *
p = Function('p')
t, a = symbols('t a')
b = symbols('b', nonzero=True)
eq = Eq(Derivative(p(t),t), -a*p(t) + exp(-(a + b)*t))
dsolve(eq)
# -> p(t) == (C1 - exp(-b*t)/b)*exp(-a*t)
# (0.75 s)

因此，通过避免特殊情况来帮助求解器将求解时间再次缩短一半（我的系统上的时间与你的相似）。如果特殊情况b == 0实际上是相关的，可以很容易地从 的泰勒级数展开中恢复它exp(-b*t) ~ 1 - b*t。

通常，指定变量为实数、非零、严格大于零等对于避免在特殊情况下使求解器挂起也非常有用。x < 0有时，将x > 0和分开求解实际上更好（避免在不知道 符号的情况下求解器无法进一步简化x == 0的臭名昭著的表达式）。sqrt(x^2)x

于 2014-02-26T13:00:57.820 回答
2

如果您使用"1st_linear_Integral"提示，问题会更清楚一些，它返回将作为未评估的积分执行的操作（我使用1st_linear它是因为这是由 返回的第一个方法classify_ode(eq)，这意味着它是dsolve默认使用的方法）：

In [61]: dsolve(Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))), hint='1st_linear_Integral')
Out[61]:
       ⎛     ⌠                                 ⎞
       ⎜     ⎮                ⌠                ⎟   ⌠
       ⎜     ⎮                ⎮ (u₃ + u₄) dt   ⎟  -⎮ (u₃ + u₄) dt
       ⎜     ⎮  -t⋅u₁  -t⋅u₂  ⌡                ⎟   ⌡
p(t) = ⎜C₁ + ⎮ ℯ     ⋅ℯ     ⋅ℯ               dt⎟⋅ℯ
       ⎝     ⌡                                 ⎠

In [62]: dsolve(Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))).subs(u1+u2, u5), hint='1st_linear_Integral')
Out[62]:
       ⎛     ⌠                          ⎞
       ⎜     ⎮         ⌠                ⎟   ⌠
       ⎜     ⎮         ⎮ (u₃ + u₄) dt   ⎟  -⎮ (u₃ + u₄) dt
       ⎜     ⎮  -t⋅u₅  ⌡                ⎟   ⌡
p(t) = ⎜C₁ + ⎮ ℯ     ⋅ℯ               dt⎟⋅ℯ
       ⎝     ⌡                          ⎠

积分算法对 有问题exp(a*x)*exp(b*x)，而对 没有问题exp((a + b)*x)。基本上，首先调用一个范围更有限的更快算法，然后如果快速算法失败，则调用一个范围更大的较慢算法。快速的可以处理exp((a + b)*x)，但目前不能exp(a*x)*exp(b*x)。

在这种特定情况下，简单的解决方法是使用以下方法将指数合并在一起powsimp：

In [67]: a = dsolve(Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))), hint='1st_linear_Integral')

In [68]: powsimp(a, deep=True)
Out[68]:
       ⎛     ⌠                                  ⎞
       ⎜     ⎮                 ⌠                ⎟   ⌠
       ⎜     ⎮  -t⋅u₁ - t⋅u₂ + ⎮ (u₃ + u₄) dt   ⎟  -⎮ (u₃ + u₄) dt
       ⎜     ⎮                 ⌡                ⎟   ⌡
p(t) = ⎜C₁ + ⎮ ℯ                              dt⎟⋅ℯ
       ⎝     ⌡                                  ⎠

In [69]: %time powsimp(a, deep=True).doit()
CPU times: user 261 ms, sys: 2.11 ms, total: 263 ms
Wall time: 262 ms
Out[69]:
       ⎛     ⎛⎧              t                for u₁ + u₂ - u₃ - u₄ = 0⎞⎞
       ⎜     ⎜⎪                                                        ⎟⎟
       ⎜     ⎜⎪  -t⋅u₁ - t⋅u₂ + t⋅(u₃ + u₄)                            ⎟⎟  -t⋅(u₃ + u₄)
p(t) = ⎜C₁ + ⎜⎨-ℯ                                                      ⎟⎟⋅ℯ
       ⎜     ⎜⎪─────────────────────────────          otherwise        ⎟⎟
       ⎜     ⎜⎪      u₁ + u₂ - u₃ - u₄                                 ⎟⎟
       ⎝     ⎝⎩                                                        ⎠⎠

总的来说，Stefan 的建议可能成立，也可能不成立。理论上，CAS 不应该关心符号常量的复杂程度，因为它们只是常量。现实中是有问题的，因为它结合了常数，然后需要看事物是否取消，等等。此外，像上面这样的微小差异可能会导致实际算法路径的巨大差异。在数学上，两个表达式可以相同，但算法取决于它们在结构上的外观。作为一般经验法则，更简单的表达式往往会做得更好。

如果您发现自己需要用符号替换子表达式很多，您可以使用cse帮助。

假设您的符号是真实的或积极的，通常也有很大帮助，尽管它与这个特定问题无关。

顺便说一句，SymPy 0.7.3 有点老了。最新版本 0.7.5 刚刚发布。

于 2014-03-01T20:15:02.093 回答