python-3.x - 为什么 Sympy 无法求解我的非线性系统？ Python 解释器将保持执行状态，直到我终止该进程

Question

我需要使用 Sympy 在 Python 中求解非线性方程组。 为此，我编写了下面的代码。但是，当我运行此代码时，Python 仍然很忙，没有返回任何错误消息，此外，也不返回解决方案。

为了进行比较，我在 Matlab 中做了相同的工作，在几秒钟内，程序返回了该系统的两个解决方案。

如何使用 Sympy 解决该系统问题？

问候。

import sympy as sym
import numpy as np

# Variables of the system
S, V, E, A, I, R  = sym.symbols('S, V, E, A, I, R')

# Parameters of the system
N = sym.Symbol("N", positive = True)
mi = sym.Symbol("mi", positive = True)
v = sym.Symbol("v", positive = True)
epsilon = sym.Symbol("epsilon", positive = True)
alpha = sym.Symbol("alpha", positive = True)
gamma_as = sym.Symbol("gamma_as", positive = True)
gamma_s = sym.Symbol("gamma_s", positive = True)
gamma_a = sym.Symbol("gamma_a", positive = True)
lamb = sym.Symbol("lamb", positive = True)
tau = sym.Symbol("tau", positive = True)
beta = sym.Symbol("beta", positive = True)
x = sym.Symbol("x")

# Declaration of the system equations
system = [mi*N - v*S + R - (beta*(A+I)/N)*S - mi*S,\
          v*S - (1-epsilon)*(beta*(A+I)/N)*V - mi*V,\
          (beta*(A+I)/N)*S + (1-epsilon)*(beta*(A+I)/N)*V - sym.exp(-mi*tau)*(beta*(A+I)/N)*S - mi*E,\
          alpha*sym.exp(-mi*tau)*(beta*(A+I)/N)*S - (gamma_as + gamma_a + mi)*A,\
          (1-alpha)*sym.exp(-mi*tau)*(beta*(A+I)/N)*S + gamma_as*A - (gamma_s + mi)*I,\
          gamma_a*A + gamma_s*I - (1+mi)*R]

# Solution
solution_set = sym.nonlinsolve(system, [S, V, E, A, I, R])
pyS, pyV, pyE, pyA, pyI, pyR = solution_set[0]
````

Answer 1

SymPy 通常使用 Groebner 基来求解这样的多项式方程组。为了计算 Groebner 基，SymPy 需要将每个方程识别为给定未知数中的多项式，其系数位于可计算域(“域”)中。您的系数同时涉及 mi 和 exp(-mi*tau)，而 SymPy 的 construct_domain 不喜欢这些系数，因此它放弃构建可计算域，并且使用“EX”域来代替，这非常慢。

解决方案是将 exp(mi*tau) 替换为另一个符号(我将只使用 tau)，然后自己显式计算 Groebner 基:

In [103]: rep = {exp(-mi*tau):tau}

In [104]: system2 = [eq.subs(rep) for eq in system]

In [105]: for eq in system2: pprint(eq)
                        S⋅β⋅(A + I)
N⋅mi + R - S⋅mi - S⋅v - ───────────
                             N     
             V⋅β⋅(1 - ε)⋅(A + I)
S⋅v - V⋅mi - ───────────────────
                      N         
        S⋅β⋅τ⋅(A + I)   S⋅β⋅(A + I)   V⋅β⋅(1 - ε)⋅(A + I)
-E⋅mi - ───────────── + ─────────── + ───────────────────
              N              N                 N         
                     S⋅α⋅β⋅τ⋅(A + I)
-A⋅(γₐ + γₐₛ + mi) + ───────────────
                            N       
                      S⋅β⋅τ⋅(1 - α)⋅(A + I)
A⋅γₐₛ - I⋅(γₛ + mi) + ─────────────────────
                                N          
A⋅γₐ + I⋅γₛ - R⋅(mi + 1)

现在我们可以使用 solve 或 nonlinsolve 但我们自己计算和求解 Groebner 基会更快:

In [106]: %time gb = groebner(system2, [S, V, E, A, I, R])
CPU times: user 3min 1s, sys: 100 ms, total: 3min 1s
Wall time: 3min 1s

Groebner 基础将方程组置于一种几乎已解的形式，称为“有理单变量表示”(RUR)。在这种情况下，它看起来像

S - a*R
V - b*R
E - c*R
A - d*R
I - e*R
R**2 + f*R + g

其中系数 a、b、c、d、e、f、g 是方程中符号参数(α、β 等)的复杂有理函数。从这里我们可以按照以下步骤来求解 Groebner 基:

1. 用 R 求解 S、V、E、A 和 I 的前 5 个线性方程。
2. 求解 R 的最终二次方程，给出两个解 R1 和 R2。
3. 将 R1 和 R2 的解代入 S、V、E、A 和 I 的解。
4. 将其作为两个解决方案元组放在一起。

即:

In [115]: syms = [S, V, E, A, I, R]

In [116]: [lsol] = linsolve(gb[:-1], syms[:-1])

In [117]: R1, R2 = roots(gb[-1], R)

In [118]: sol1 = lsol.subs(R, R1) + (R1,)

In [119]: sol2 = lsol.subs(R, R2) + (R2,)

现在我们有了 nonlinsolve 返回形式的两个解元组。不幸的是，解决方案非常复杂，所以我不会完整地展示它们。您可以通过查看字符串表示的长度来了解其复杂性:

In [122]: print(len(str(sol1)))
794100

In [123]: print(len(str(sol2)))
27850

现在，值得考虑一下您真正想要这些解决方案的目的。也许只是您想用一些显式数值替换符号参数。这里值得注意的是，在尝试以符号方式求解它们之前，首先将这些值代入方程可能会更有效。如果您想了解您的解决方案如何依赖于某些特定参数(例如 mi)，那么您可以用值替换其他所有内容，并更快地获得仅涉及该参数的更简单形式的解决方案:

In [136]: rep = {alpha:1, beta:2, epsilon:3, gamma_as:4, gamma_s:5, gamma_a:6, exp(-mi*tau):7, N:8, v
     ...: :9}

In [137]: system2 = [eq.subs(rep) for eq in system]

In [138]: %time solve(system2, syms)
CPU times: user 3.92 s, sys: 4 ms, total: 3.92 s
Wall time: 3.92 s
Out[138]: 
⎡                              ⎛                                                       ⎛  2        
⎢⎛ 8⋅mi     72              ⎞  ⎜4⋅(mi + 5)⋅(mi + 10)   36⋅(mi + 5)⋅(mi + 10)⋅(mi + 12)⋅⎝mi  + 4⋅mi 
⎢⎜──────, ──────, 0, 0, 0, 0⎟, ⎜────────────────────, ─────────────────────────────────────────────
⎢⎝mi + 9  mi + 9            ⎠  ⎜     7⋅(mi + 9)                  ⎛    4        3         2         
⎣                              ⎝                      7⋅(mi + 9)⋅⎝3⋅mi  + 38⋅mi  + 161⋅mi  + 718⋅mi

    ⎞                                   ⎛  2          ⎞ ⎛    3        2               ⎞            
- 25⎠    24⋅(mi + 1)⋅(mi + 5)⋅(mi + 10)⋅⎝mi  + mi + 50⎠⋅⎝3⋅mi  + 41⋅mi  + 209⋅mi + 787⎠  -4⋅(mi + 1
───────, ──────────────────────────────────────────────────────────────────────────────, ──────────
      ⎞                 ⎛  2            ⎞ ⎛    4        3         2               ⎞                
 + 900⎠     7⋅(mi + 12)⋅⎝mi  + 4⋅mi - 25⎠⋅⎝3⋅mi  + 38⋅mi  + 161⋅mi  + 718⋅mi + 900⎠           (mi +

           ⎛  2          ⎞                ⎛  2          ⎞                  ⎛  2          ⎞ ⎞⎤
)⋅(mi + 5)⋅⎝mi  + mi + 50⎠   -16⋅(mi + 1)⋅⎝mi  + mi + 50⎠   -8⋅(3⋅mi + 25)⋅⎝mi  + mi + 50⎠ ⎟⎥
───────────────────────────, ─────────────────────────────, ───────────────────────────────⎟⎥
     ⎛  2            ⎞                  ⎛  2            ⎞               ⎛  2            ⎞  ⎟⎥
 12)⋅⎝mi  + 4⋅mi - 25⎠        (mi + 12)⋅⎝mi  + 4⋅mi - 25⎠     (mi + 12)⋅⎝mi  + 4⋅mi - 25⎠  ⎠⎦

如果你用值替换所有参数，那么速度会快得多:

In [139]: rep = {alpha:1, beta:2, epsilon:3, gamma_as:4, gamma_s:5, gamma_a:6, exp(-mi*tau):7, N:8, v
     ...: :9, mi:10}

In [140]: system2 = [eq.subs(rep) for eq in system]

In [141]: %time solve(system2, syms)
CPU times: user 112 ms, sys: 0 ns, total: 112 ms
Wall time: 111 ms
Out[141]: 
⎡⎛1200  124200  5224320  -960   -256   -640 ⎞  ⎛80  72            ⎞⎤
⎢⎜────, ──────, ───────, ─────, ─────, ─────⎟, ⎜──, ──, 0, 0, 0, 0⎟⎥
⎣⎝133   55727    67459     23     23     23 ⎠  ⎝19  19            ⎠⎦