【发布时间】: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 时获得数值解并不是我想要的。
【问题讨论】: