【问题标题】:sympy nonlinsolve hangs on simple system of equationssympy nonlinsolve 依赖于简单的方程组
【发布时间】:2020-10-16 13:47:50
【问题描述】:

我是 sympy 的新手,正在尝试用它解决一个相对简单的方程组:

import sympy
A, B, I, AI, BI, A0, B0, I0, k1, k2, k3, k4 = sympy.symbols('A B I AI BI A_0 B_0 I_0 k1 k2 k3 k4', real=True)
eqs = [A0 - AI - A,
B0 - BI - B,
I0 - AI - BI - I,
k1*(A0 - AI)*I - k2*AI,
k3*(B0 - BI)*I - k4*BI]

如果我为A,B,I 寻找解决方案,它会起作用:

# this works
result = sympy.nonlinsolve(eqs, (A,B,I))
print(result)

给出结果FiniteSet((-AI + A_0, -BI + B_0, -AI - BI + I_0)),但如果我就其他变量要求变量AI,BI,A,B,I 的解决方案,它似乎挂起:

# this hangs
result = sympy.nonlinsolve(eqs, (AI,BI,A,B,I))

我在这里做错了什么?

(这里使用的是sympy '1.6.2'

编辑: 为了回应@Oscar Benjamin 的有用建议,我们尝试以两种方式求解原始系统:(1) 通过将具有线性解的变量代入方程, (2) 通过简化问题并假设我们希望我们的解决方案成为函数的 7 个变量中的 3 个是已知的。这两种方法都给出了大多数错误的答案:nans 和非真实的解决方案。有没有办法解决这个问题?

import sympy
from sympy import solve, factor, roots, nonlinsolve
A, B, I, AI, BI, A0, B0, I0, k1, k2, k3, k4 = sympy.symbols('A B I AI BI A_0 B_0 I_0 k1 k2 k3 k4', real=True)
eqs = [A0 - AI - A,
B0 - BI - B,
I0 - AI - BI - I,
k1*(A0 - AI)*I - k2*AI,
k3*(B0 - BI)*I - k4*BI]

print("trying to solve original system:")
# solve linear equations and substitute their solutions in
((As, Bs, Is),) = nonlinsolve(eqs, (A,B,I))
eqs2 = [eq.subs({A:As,B:Bs,I:Is}) for eq in eqs]
eq1, eq2 = eqs2[3:]
p = eq1.subs(AI, solve(eq2, AI)[0]).as_numer_denom()[0].expand().collect(BI)
BI1, BI2, BI3, BI4 = roots(p, BI)

def eval_solns(inputs, roots):
    for vals in inputs:
        for n, root in enumerate(roots):
            print("root %d yields: " %(n+1))
            print(root.subs(vals))

inputs = [{A0: 100, B0: 100, I0: 100, k1: 0.1, k2: 0.1, k3: 0.1, k4: 0.1},
          {A0: 100, B0: 100, I0: 100, k1: 0.2, k2: 0.1, k3: 0.3, k4: 0.4}]
# all of these give wrong, non-real solutions to the system
orig_roots = [BI1, BI2, BI3, BI4]
eval_solns(inputs, orig_roots)

# second try: simplify the problem by assuming A0, B0, I0 are given
# substitute them in
eqs2 = [eq.subs({A:As,B:Bs,I:Is}) for eq in eqs2]
eqs2 = [eq.subs({A0: 100, B0: 100, I0: 100}) for eq in eqs2]
eq1, eq2 = eqs2[3:]     
print("*\ntrying simpler system with A0,B0,I0 given: ")
print(eqs2)
p = eq1.subs(AI, solve(eq2, AI)[0]).as_numer_denom()[0].expand().collect(BI)                                                          
BI1, BI2, BI3, BI4 = roots(p, BI)
new_roots = [BI1, BI2, BI3, BI4]
# solve the simpler system
new_inputs = [{k1: 0.1, k2: 0.1, k3: 0.1, k4: 0.1},
              {k1: 0.2, k2: 0.1, k3: 0.3, k4: 0.4},
              {k1: 0.05, k2: 1.0, k3: 1.0, k4: 1.2}]
# all of these answers are wrong
eval_solns(new_inputs, new_roots)

我得到的输出包括:

...
trying simpler system with A0,B0,I0 given: 
[0, 0, 0, -AI*k2 + k1*(100 - AI)*(-AI - BI + 100), -BI*k4 + k3*(100 - BI)*(-AI - BI + 100)]
root 1 yields: 
100
root 2 yields: 
nan
root 3 yields: 
nan
root 4 yields: 
nan
root 1 yields: 
100
root 2 yields: 
-6.22222222222222 - 33.1666666666667*(9.77105856678793 + 8.49476415953825*I)**(1/3) - 182.875860785409/(9.77105856678793 + 8.49476415953825*I)**(1/3)
root 3 yields: 
-6.22222222222222 - 33.1666666666667*(-1/2 + sqrt(3)*I/2)*(9.77105856678793 + 8.49476415953825*I)**(1/3) - 182.875860785409/((-1/2 + sqrt(3)*I/2)*(9.77105856678793 + 8.49476415953825*I)**(1/3))
root 4 yields: 
-6.22222222222222 - 182.875860785409/((-1/2 - sqrt(3)*I/2)*(9.77105856678793 + 8.49476415953825*I)**(1/3)) - 33.1666666666667*(-1/2 - sqrt(3)*I/2)*(9.77105856678793 + 8.49476415953825*I)**(1/3)
root 1 yields: 
100
root 2 yields: 
104.655319148936 - 100*(-0.00146514175733681 + 0.00423691411709115*I)**(1/3) - 2.718847623359/(-0.00146514175733681 + 0.00423691411709115*I)**(1/3)
root 3 yields: 
104.655319148936 - 100*(-1/2 + sqrt(3)*I/2)*(-0.00146514175733681 + 0.00423691411709115*I)**(1/3) - 2.718847623359/((-1/2 + sqrt(3)*I/2)*(-0.00146514175733681 + 0.00423691411709115*I)**(1/3))
root 4 yields: 
104.655319148936 - 2.718847623359/((-1/2 - sqrt(3)*I/2)*(-0.00146514175733681 + 0.00423691411709115*I)**(1/3)) - 100*(-1/2 - sqrt(3)*I/2)*(-0.00146514175733681 + 0.00423691411709115*I)**(1/3)

在这个系统中,所有变量都是正实数,所有解都应该是实数。在 sympy 中用positive=True 声明第一部分,如下所示,没有任何区别:

A, B, I, AI, BI, A0, B0, I0, k1, k2, k3, k4 = sympy.symbols('A B I AI BI A_0 B_0 I_0 k1 k2 k3 k4', real=True, positive=True)

即。所有的答案仍然不正确。

编辑 2: 澄清一下,我对解析解决方案感兴趣,但仍然不明白为什么它太复杂而无法推导(然后与实际值一起使用插入)同情。我想要解析解的原因是,在给定 A0、B0、I0 的某些设置的情况下,我可以直接求解 k1、k2、k3、k4 的特定值。这就是在已知 A0、B0 和 I0 的情况下可以得到解析解的原因。但是,在给出所有 A0、B0、I0、k1、k2、k3、k4 时获得数值解并不是我想要的。

【问题讨论】:

    标签: python sympy


    【解决方案1】:

    它可能看起来很简单,但这归结为具有符号系数的三次方,因此尽管可以通过代数找到解决方案,但它并不简单。在这种情况下,nonlinsolve 在计算较大的 Groebner 基时会很慢。不过,有一种更快的方法可以找到解决方案。

    首先求解ABI 的平凡线性方程组:

    In [55]: ((As, Bs, Is),) = nonlinsolve(eqs, (A,B,I))                                                                                           
    
    In [56]: As                                                                                                                                    
    Out[56]: -AI + A₀
    
    In [57]: Bs                                                                                                                                    
    Out[57]: -BI + B₀
    
    In [58]: Is                                                                                                                                    
    Out[58]: -AI - BI + I₀
    

    我们现在可以从系统中消除那些未知数和方程:

    In [59]: eqs2 = [eq.subs({A:As,B:Bs,I:Is}) for eq in eqs]                                                                                      
    
    In [60]: eqs2                                                                                                                                  
    Out[60]: [0, 0, 0, -AI⋅k₂ + k₁⋅(-AI + A₀)⋅(-AI - BI + I₀), -BI⋅k₄ + k₃⋅(-BI + B₀)⋅(-AI - BI + I₀)]
    
    In [61]: eq1, eq2 = eqs2[3:]                                                                                                                   
    
    In [62]: eq1                                                                                                                                   
    Out[62]: -AI⋅k₂ + k₁⋅(-AI + A₀)⋅(-AI - BI + I₀)
    
    In [63]: eq2                                                                                                                                   
    Out[63]: -BI⋅k₄ + k₃⋅(-BI + B₀)⋅(-AI - BI + I₀)
    

    此时,我们在AIBI 中有一个二次多元多项式系统。我们可以为AI 求解eq2(因为它在AI 中是线性的)并将其代入eq1 以仅在BI 中得到一个方程:

    In [68]: solve(eq2, AI)                                                                                                                        
    Out[68]: 
    ⎡    2                                            ⎤
    ⎢- BI ⋅k₃ + BI⋅B₀⋅k₃ + BI⋅I₀⋅k₃ + BI⋅k₄ - B₀⋅I₀⋅k₃⎥
    ⎢─────────────────────────────────────────────────⎥
    ⎣                   k₃⋅(BI - B₀)                  ⎦
    
    In [69]: eq1.subs(AI, solve(eq2, AI)[0])                                                                                                       
    Out[69]: 
       ⎛         2                                            ⎞ ⎛               2                                            ⎞      ⎛    2       
       ⎜     - BI ⋅k₃ + BI⋅B₀⋅k₃ + BI⋅I₀⋅k₃ + BI⋅k₄ - B₀⋅I₀⋅k₃⎟ ⎜           - BI ⋅k₃ + BI⋅B₀⋅k₃ + BI⋅I₀⋅k₃ + BI⋅k₄ - B₀⋅I₀⋅k₃⎟   k₂⋅⎝- BI ⋅k₃ + B
    k₁⋅⎜A₀ - ─────────────────────────────────────────────────⎟⋅⎜-BI + I₀ - ─────────────────────────────────────────────────⎟ - ────────────────
       ⎝                        k₃⋅(BI - B₀)                  ⎠ ⎝                              k₃⋅(BI - B₀)                  ⎠                   
    
                                         ⎞
    I⋅B₀⋅k₃ + BI⋅I₀⋅k₃ + BI⋅k₄ - B₀⋅I₀⋅k₃⎠
    ──────────────────────────────────────
         k₃⋅(BI - B₀) 
    

    我们可以将其简化为四次:

    In [73]: p = eq1.subs(AI, solve(eq2, AI)[0]).as_numer_denom()[0].expand().collect(BI)                                                          
    
    In [74]: p                                                                                                                                     
    Out[74]: 
      4 ⎛       2           3⎞     3 ⎛          2                2                3           2              3           2        2   ⎞     2 ⎛  
    BI ⋅⎝- k₁⋅k₃ ⋅k₄ + k₂⋅k₃ ⎠ + BI ⋅⎝- A₀⋅k₁⋅k₃ ⋅k₄ + 2⋅B₀⋅k₁⋅k₃ ⋅k₄ - 3⋅B₀⋅k₂⋅k₃  + I₀⋅k₁⋅k₃ ⋅k₄ - I₀⋅k₂⋅k₃  + k₁⋅k₃⋅k₄  - k₂⋅k₃ ⋅k₄⎠ + BI ⋅⎝2⋅
    
               2        2      2          2      3                2                   3              2             2   ⎞      ⎛       2      2   
    A₀⋅B₀⋅k₁⋅k₃ ⋅k₄ - B₀ ⋅k₁⋅k₃ ⋅k₄ + 3⋅B₀ ⋅k₂⋅k₃  - 2⋅B₀⋅I₀⋅k₁⋅k₃ ⋅k₄ + 3⋅B₀⋅I₀⋅k₂⋅k₃  - B₀⋅k₁⋅k₃⋅k₄  + 2⋅B₀⋅k₂⋅k₃ ⋅k₄⎠ + BI⋅⎝- A₀⋅B₀ ⋅k₁⋅k₃ ⋅k₄
    
         3      3     2         2          2         3     2      2   ⎞     3         3
     - B₀ ⋅k₂⋅k₃  + B₀ ⋅I₀⋅k₁⋅k₃ ⋅k₄ - 3⋅B₀ ⋅I₀⋅k₂⋅k₃  - B₀ ⋅k₂⋅k₃ ⋅k₄⎠ + B₀ ⋅I₀⋅k₂⋅k₃ 
    

    factor 的简单调用表明四个根之一就是BI = B0

    In [84]: factor(p)                                                                                                                             
    Out[84]: 
                  ⎛     2                                  3              3      2     2                   2         2     2                 2   
    -k₃⋅(BI - B₀)⋅⎝A₀⋅BI ⋅k₁⋅k₃⋅k₄ - A₀⋅BI⋅B₀⋅k₁⋅k₃⋅k₄ + BI ⋅k₁⋅k₃⋅k₄ - BI ⋅k₂⋅k₃  - BI ⋅B₀⋅k₁⋅k₃⋅k₄ + 2⋅BI ⋅B₀⋅k₂⋅k₃  - BI ⋅I₀⋅k₁⋅k₃⋅k₄ + BI ⋅I₀
    
          2     2      2     2                 2      2                                       2                      2         2⎞
    ⋅k₂⋅k₃  - BI ⋅k₁⋅k₄  + BI ⋅k₂⋅k₃⋅k₄ - BI⋅B₀ ⋅k₂⋅k₃  + BI⋅B₀⋅I₀⋅k₁⋅k₃⋅k₄ - 2⋅BI⋅B₀⋅I₀⋅k₂⋅k₃  - BI⋅B₀⋅k₂⋅k₃⋅k₄ + B₀ ⋅I₀⋅k₂⋅k₃ ⎠
    

    也许该解决方案对您来说已经足够了。其他三个要复杂得多。您可以通过以下方式获得它们:

    In [85]: BI1, BI2, BI3, BI4 = roots(p, BI)                                                                                                     
    
    In [86]: BI1                                                                                                                                   
    Out[86]: B₀
    

    我不会显示BI2 等的输出,因为它太复杂了。并非所有这些根都必然是原始系统的解决方案,因为要进行转换才能得到它们,因此您必须向后工作以找出哪些是。

    我认为这意味着BI = B0 解决方案仅在B0k4 为零时才有效,例如:

    In [99]: eq2.subs(BI, BI1)                                                                                                                     
    Out[99]: -B₀⋅k₄
    

    编辑:我上面的答案旨在展示如何获得一般的象征性答案,但这似乎不是您真正想要的(答案太复杂而没有多大用处)。由于您想获得特定数字的答案,您只需将这些数字代入方程式即可。如果您不是在寻找解析解决方案,您可以只使用 nsolve:

    In [29]: import sympy 
        ...: from sympy import solve, factor, roots, nonlinsolve 
        ...: A, B, I, AI, BI, A0, B0, I0, k1, k2, k3, k4 = sympy.symbols('A B I AI BI A_0 B_0 I_0 k1 k2 k3 k4', positive=True) 
        ...: eqs = [A0 - AI - A, 
        ...: B0 - BI - B, 
        ...: I0 - AI - BI - I, 
        ...: k1*(A0 - AI)*I - k2*AI, 
        ...: k3*(B0 - BI)*I - k4*BI]                                                                                                               
    
    In [30]: inputs = [{A0: 100, B0: 100, I0: 100, k1: 0.1, k2: 0.1, k3: 0.1, k4: 0.1}, 
        ...:           {A0: 100, B0: 100, I0: 100, k1: 0.2, k2: 0.1, k3: 0.3, k4: 0.4}]                                                            
    
    In [31]: eqs_sub = [eq.subs(inputs[0]) for eq in eqs]                                                                                          
    
    In [32]: solve(eqs_sub, [A, B, I, AI, BI])                                                                                                     
    Out[32]: [(50.4902894311622, 50.4902894311622, 0.980578862324382, 49.5097105688378, 49.5097105688378)]
    
    In [33]: nsolve(eqs_sub, [A, B, I, AI, BI], [1, 1, 1, 1, 1])                                                                                   
    Out[33]: 
    ⎡50.4902894311622 ⎤
    ⎢                 ⎥
    ⎢50.4902894311622 ⎥
    ⎢                 ⎥
    ⎢0.980578862324382⎥
    ⎢                 ⎥
    ⎢49.5097105688378 ⎥
    ⎢                 ⎥
    ⎣49.5097105688378 ⎦
    
    In [34]: eqs_sub = [eq.subs(inputs[1]) for eq in eqs]                                                                                          
    
    In [35]: nsolve(eqs_sub, [A, B, I, AI, BI], [1, 1, 1, 1, 1])                                                                                   
    Out[35]: 
    ⎡38.3817627411582⎤
    ⎢                ⎥
    ⎢62.4209392826439⎥
    ⎢                ⎥
    ⎢0.80270202380213⎥
    ⎢                ⎥
    ⎢61.6182372588418⎥
    ⎢                ⎥
    ⎣37.5790607173561⎦
    

    【讨论】:

    • 谢谢,但我对剩余变量(A0、B0、I0、k1、k2、k3、k4)的 A、B、I、AI、BI 表达式感兴趣,所以不幸的是,任何使用 A、B、I、AI 或 BI 的表达式都是没有用的。
    • 我已经(主要)向您展示了如何获得您想要的东西。然而最终的结果是相当复杂的。看看BI2。根据您的要求,这给出了 BI 的剩余变量。虽然这是一个复杂的表达式,但其他表达式也是如此。这就是我没有给出最终答案的原因。
    • 谢谢。你能多说一点,为什么当 nonlinsolve 使用的自动化程序挂起时这些手替换似乎有效?我不这样做
    • BI2 和其余的根也不正确。 BI2.subs({A0: 100, B0: 100, I0: 100, k1: 1.0, k2: 1.0, k3: 1.0, k4: 1.0}) 产生 nan。所有正实数都应该在这里给出正实数的答案。不知道出了什么问题。
    • nonlinsolve 函数在这种情况下很慢,因为正如我所说,它正在计算一个大的 Groebner 基。
    猜你喜欢
    • 2022-08-17
    • 2013-03-11
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2013-03-20
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多