1

我得到了这 6 个方程的非线性系统的意外解。我有 6 个方程、6 个变量和其他 5 个常数（符号）。

import sympy as sy

P1,P2,P,k1,k2,d1,d2,delta,teta,d,H=sy.symbols('P1,P2,P,k1,k2,d1,d2,delta,teta,d,H')

#equilibria
eqV=P+P1+P2
eqM=P1*d*sy.cos(teta)-P2*d*sy.cos(teta)-P*H*sy.sin(teta)

#constitutive
soil1=P1+k1*d1
soil2=P2+k2*d2

#congruents
crot=sy.tan(teta)-(d1-d2)/2/d
cvert=delta-(d1+d2)/2

solution=sy.nonlinsolve((eqM,eqV,crot,cvert,soil1,soil2),[d1,d2,P1,P2,teta,delta])

检查解决方案我没有找到符号常量“H”。对我来说是出乎意料的。

所以问题是：如何使用 SymPy 以正确的方式求解非线性方程组？

4

2 回答 2

0

您偶然发现了nonlinsolve. 它应该作为一个问题报告给 sympy： https ://github.com/sympy/sympy/issues

另一个答案显示了一种解决方法。我将展示另一个基于使用 Groebner 的基础。为此，您需要先将sin/cos变成tan，以便方程都是多项式的：

In [49]: eqM2 = (eqM/cos(teta)).subs(sin(teta), tan(teta)*cos(teta)).cancel()

In [50]: eqM2
Out[50]: -H⋅P⋅tan(teta) + P₁⋅d - P₂⋅d

In [51]: eqs = (eqM2,eqV,crot,cvert,soil1,soil2)

In [52]: syms = [d1,d2,P1,P2,tan(teta),delta]

In [53]: groebner(eqs, syms)
Out[53]: 
             ⎛⎡               2        2                         2        2                     2           2                      2           2            
             ⎜⎢          - H⋅P  - 2⋅P⋅d ⋅k₂                 - H⋅P  - 2⋅P⋅d ⋅k₁               H⋅P ⋅k₁ + 2⋅P⋅d ⋅k₁⋅k₂             H⋅P ⋅k₂ + 2⋅P⋅d ⋅k₁⋅k₂      
GroebnerBasis⎜⎢d₁ + ────────────────────────────, d₂ + ────────────────────────────, P₁ + ────────────────────────────, P₂ + ────────────────────────────, ─
             ⎜⎢                          2                                  2                                  2                                  2         
             ⎝⎣     H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂       H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂       H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂       H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂  H

                                                   2      2         2    ⎤                                                                    ⎞
     P⋅d⋅k₁ - P⋅d⋅k₂                          - H⋅P  - P⋅d ⋅k₁ - P⋅d ⋅k₂ ⎥                                                                    ⎟
─────────────────────────── + tan(teta), δ + ────────────────────────────⎥, d₁, d₂, P₁, P₂, tan(teta), δ, domain=ℤ(d, k₁, k₂, H, P), order=lex⎟
                    2                                             2      ⎥                                                                    ⎟
⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂                  H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂⎦                                                                    ⎠

In [54]: solve(groebner(eqs, syms), syms)
Out[54]: 
⎧           2           2                     2           2                       2        2                        2        2                   2      2   
⎪      - H⋅P ⋅k₁ - 2⋅P⋅d ⋅k₁⋅k₂          - H⋅P ⋅k₂ - 2⋅P⋅d ⋅k₁⋅k₂              H⋅P  + 2⋅P⋅d ⋅k₂                  H⋅P  + 2⋅P⋅d ⋅k₁             H⋅P  + P⋅d ⋅k₁
⎨P₁: ────────────────────────────, P₂: ────────────────────────────, d₁: ────────────────────────────, d₂: ────────────────────────────, δ: ────────────────
⎪                         2                                 2                                 2                                 2                           
⎩    H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂      H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂      H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂      H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂     H⋅P⋅k₁ + H⋅P⋅k₂ 

      2                                              ⎫
 + P⋅d ⋅k₂                     -P⋅d⋅k₁ + P⋅d⋅k₂      ⎪
────────────, tan(teta): ────────────────────────────⎬
     2                                        2      ⎪
+ 4⋅d ⋅k₁⋅k₂             H⋅P⋅k₁ + H⋅P⋅k₂ + 4⋅d ⋅k₁⋅k₂⎭
于 2021-08-17T18:18:13.110 回答
0

在未知数P1、P2、d1、d2中有 4 个线性方程

eqV=P+P1+P2
soil1=P1+k1*d1
soil2=P2+k2*d2
cvert=delta-(d1+d2)/2

你可以解决，只留下未知的增量

sol_P1_P2_d1_d2 = sy.solve((eqV,soil1,soil2, cvert), (P1, P2, d1, d2))

现在，您可以代入其余方程中的第一个

eqM = (P1*d*sy.cos(teta)-P2*d*sy.cos(teta)-P*H*sy.sin(teta)).subs(sol_P1_P2_d1_d2)

并注意它在delta中是线性的，所以你可以再次使用solve

sol_delta = sy.solve((eqM,), delta)

最后

crot = (sy.tan(teta)-(d1-d2)/2/d).subs(sol_P1_P2_d1_d2).subs(sol_delta)
sol_teta = sy.solve((crot,), teta, dict=1)[0]

在这一点上，你可以回替代

sol_delta = sol_delta.subs(sol_teta)
sol_P1_P2_d1_d2 = sol_P1_P2_d1_d2.subs(sol_delta)
于 2021-08-17T17:56:13.740 回答