CODEWITHSUNDEEP
sympy
user9276239
user9276239

Reputation:

Solve a nonlinear system of equation via sympy

I have a system of ODE and want to find equilibrium points by using nonlinsolve() but when i run it by jubyter or spyder the program keep running without any result.

N,x1,x2,x3,x4,y1,y2,r1,r2,r3,r4,eta1,eta2,eta3,eta4,R,c1,c2,c3,c4,a11,a12,a21,a22,a31,a32,a41,a42,b12,h,h11,h12,h21,h22,h31,h32,h41,h42,s1,s2,s3,s4,epsilon1,epsilon2,omega1,omega2,K11,K22,beta11,beta21,beta31,beta41,beta12,beta22,beta32,beta42,gamma12=sp.symbols('x1,x2,x3,x4,y1,y2,r1,r2,r3,r4,eta1,eta2,eta3,eta4,N,R,c1,c2,c3,c4,a11,a12,a21,a22,a31,a32,a41,a42,b12,h,h11,h12,h21,h22,h31,h32,h41,h42,s1,s2,s3,s4,epsilon1,epsilon2,omega1,omega2,K11,K22,beta11,beta21,beta3a,beta41,beta12,beta22,beta32,beta42,gamma12')


F2=x1*(r1*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a11*y1)/(y1+a11*h11*x1)-(a12*y2)/(y2+a12*h12*x1))+s1

F3=x2*(r2*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a21*y1)/(y1+a21*h21*x2)-(a22*y2)/(y2+a22*h22*x2))+s2

F4=x3*(r3*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a31*y1)/(y1+a31*h31*x3)-(a32*y2)/(y2+a32*h32*x3))+s3

F5=x4*(r4*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a41*y1)/(y1+a42*h41*x4)-(a42*y2)/(y2+a42*h42*x4))+s4

F6=y1*(-epsilon1*(1+(y1+omega2*y2)/K11)-(b12*y2)/(y2+b12*h*y1)\
              +beta11*(a11*x1)/(y1+a11*h11*x1)\
              +beta21*(a21*x2)/(y1+a21*h21*x2)\
              +beta31*(a31*x3)/(y1+a31*h31*x3)\
              +beta41*(a41*x4)/(y1+a41*h41*x4))

F7=y2*(-epsilon2*(1+(omega1*y1+y2)/K22)+gamma12*(b12*y1)/(y2+b12*h*y1)\
              +beta12*(a12*x1)/(y2+a12*h12*x1)\
              +beta22*(a22*x2)/(y2+a22*h22*x2)\
              +beta32*(a32*x3)/(y2+a32*h32*x3)\
              +beta42*(a42*x4)/(y2+a42*h42*x4))   
equ=(F2,F3,F4,F5,F6,F7)
sol=nonlinsolve(equ,x1,x2,x3,x4,y1,y2)   
print(sol)

Is it possible to get the solution in terms of parameters?

Upvotes: 0

Views: 132

Answers (1)

Oscar Benjamin
Oscar Benjamin

Reputation: 14530

Your system of equations can be recast as multivariate polynomial system if you multiply through by all of the denominators. Doing this for example in the case of F2 you would get

In [26]: F2.as_numer_denom()[0]                                                                                                   
Out[26]: 
s₁⋅x₁⋅(a₁₁⋅h₁₁⋅x₂ + y₂)⋅(a₁₂⋅h₁₂⋅x₂ + r₁) + x₂⋅(-a₁₁⋅x₁⋅y₂⋅(a₁₂⋅h₁₂⋅x₂ + r₁) - a₁₂⋅r₁⋅x₁⋅(a₁₁⋅h₁₁⋅x₂ + y₂) + r₂⋅(a₁₁⋅h₁₁⋅x₂ + y₂
)⋅(a₁₂⋅h₁₂⋅x₂ + r₁)⋅(-N⋅y₁ - η₂⋅x₂ - η₃⋅x₃ - η₄⋅x₄ + x₁))

We can see from here that the polynomial is of order 5 since it has terms like x2**3*N*y1 so broadly you have a system of 7 polynomials not of low order. I expect that a general closed form solution will not be possible unless you are lucky.

Upvotes: 1

Related Questions