如何解决sympy nonlinsolve依赖简单的方程组
我是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,I的解决方案，它会起作用：
# this works
result = sympy.nonlinsolve(eqs,(A,I))
print(result)
给出结果FiniteSet((-AI + A_0,-BI + B_0,-AI - BI + I_0))，但是如果我要就其他变量寻求变量AI,A,I的解决方案，它似乎会死机：
# this hangs
result = sympy.nonlinsolve(eqs,(AI,I))
我在做什么错了？
（正在使用sympy '1.6.2'）
edit：针对@Oscar Benjamin的有用建议，尝试以两种方式解决原始系统：（1）通过将具有线性解的变量代入方程式中，（2）通过简化问题，并假设我们希望我们的解决方案成为函数的7个变量中的3个是已知的。两种方法都给出了大多数错误的答案：nans和非真实的解决方案。有办法解决吗？
import sympy
from sympy import solve,factor,roots,nonlinsolve
A,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,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,k1: 0.2,k3: 0.3,k4: 0.4}]
# all of these give wrong,non-real solutions to the system
orig_roots = [BI1,BI4]
eval_solns(inputs,orig_roots)
# second try: simplify the problem by assuming A0,I0 are given
# substitute them in
eqs2 = [eq.subs({A:As,I:Is}) for eq in eqs2]
eqs2 = [eq.subs({A0: 100,I0: 100}) for eq in eqs2]
eq1,eq2 = eqs2[3:]
print("*\ntrying simpler system with A0,I0 given: ")
print(eqs2)
p = eq1.subs(AI,AI)[0]).as_numer_denom()[0].expand().collect(BI)
BI1,BI)
new_roots = [BI1,BI4]
# solve the simpler system
new_inputs = [{k1: 0.1,{k1: 0.2,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,I0 given:
[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)
在此系统中，所有变量均为正实数，所有解均应为实数。如下所示，用positive=True声明sympy的第一部分没有区别：
A,real=True,positive=True)
即。所有答案仍然不正确。
编辑2：澄清一下，我 am 对解析解决方案感兴趣，但仍然不明白为什么它太复杂而无法导出（然后用于实际值）插入）中的sympy。我想要解析解决方案的原因是，只要给定A0，B0，I0的某些设置，我就可以直接求解k1，k2，k3，k4的特定值。这就是为什么已知A0，B0和I0就能获得解析解的原因。但是，当我给出所有的A0，B0，I0，k1，k2，k3，k4时得到数值解不是我想要的。
解决方法
它看起来很简单，但是可以归结为具有符号系数的三次方，因此，尽管可以通过代数方式找到解，但这并不简单。在这种情况下，nonlinsolve在计算较大的Groebner基础时速度较慢。不过，有一种更快的方法可以解决该问题。
首先求解A，B和I的平凡线性方程：
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,-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₀)
这时，我们在AI和BI中有一个二次多元多项式系统。我们可以为eq2求解AI（因为它在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,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解决方案仅在B0或k4为零的情况下才有效：
In [99]: eq2.subs(BI,BI1)
Out[99]: -B₀⋅k₄
编辑：我在上面的答案旨在显示如何获得一般的符号答案，但这似乎并不是您真正想要的（答案太复杂了，无论如何都无法使用）。由于您想得到特定数字的答案，因此您只需将这些数字代入方程式即可。如果您不寻找分析解决方案，则可以使用nsolve：
In [29]: import sympy
...: from sympy import solve,factor,roots,nonlinsolve
...: A,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,k1: 0.2,k3: 0.3,k4: 0.4}]
In [31]: eqs_sub = [eq.subs(inputs[0]) for eq in eqs]
In [32]: solve(eqs_sub,[A,BI])
Out[32]: [(50.4902894311622,50.4902894311622,0.980578862324382,49.5097105688378,49.5097105688378)]
In [33]: nsolve(eqs_sub,BI],[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,1])
Out[35]:
⎡38.3817627411582⎤
⎢ ⎥
⎢62.4209392826439⎥
⎢ ⎥
⎢0.80270202380213⎥
⎢ ⎥
⎢61.6182372588418⎥
⎢ ⎥
⎣37.5790607173561⎦
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。