solve non linear system of equations in sympy

1
votes

I get a not-expected solution of this 6 equations' non linear system. I have 6 equations, 6 variables and other 5 constants (symbolic).

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])

inspecting solution I don't find symbolic constant 'H'. For me is unexpected.

So question is: How to solve in a correct way a system of nonlinear equations using SymPy?

2
sympynonlinear-equation
Or, to make it even simpler, if you divide eqM by cos(θ) and introduce a new variable, τ=tan(θ) in the 2 equations eqM and crot, then you have a linear system in d₁,d₂,P₁,P₂,τ,δ - gboffi
{P1: (-H*P**2*k1 - 2*P*d**2*k1*k2)/(H*P*k1 + H*P*k2 + 4*d**2*k1*k2), P2: (-H*P**2*k2 - 2*P*d**2*k1*k2)/(H*P*k1 + H*P*k2 + 4*d**2*k1*k2), d1: (H*P**2 + 2*P*d**2*k2)/(H*P*k1 + H*P*k2 + 4*d**2*k1*k2), d2: (H*P**2 + 2*P*d**2*k1)/(H*P*k1 + H*P*k2 + 4*d**2*k1*k2), delta: (H*P**2 + P*d**2*k1 + P*d**2*k2)/(H*P*k1 + H*P*k2 + 4*d**2*k1*k2), tau: (-P*d*k1 + P*d*k2)/(H*P*k1 + H*P*k2 + 4*d**2*k1*k2)} - gboffi
thank you, actually I didn't think was a linear system, but it is a starting point, k1 and k2 should become function of d1 and d2... - Filippo Forlani

2 Answers

0
votes

You have 4 linear equations in the unknowns P1, P2, d1, d2

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

that you can solve, leaving unknown only delta

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

Now, you can substitute in the first of the remaining equations

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

and note that it's linear in delta, so again you can use solve

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

and finally

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]

at this point you can back substitute

sol_delta = sol_delta.subs(sol_teta)
sol_P1_P2_d1_d2 = sol_P1_P2_d1_d2.subs(sol_delta)
0
votes

You've stumbled across a bug in nonlinsolve. It should be reported as an issue to sympy: https://github.com/sympy/sympy/issues

The other answer shows one workaround for this. I'll show another based on using the Groebner basis. For this you need to turn the sin/cos into tan first so that the equations are all polnomial:

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₂⎭