Solve a system of nonlinear equations with conditions

I have the follwing system of equation which I want to solve . How I will put conditions that all parameters shoud be positive and variables (N,x1,x2,x3,x4,y1,y2) should be equal or greater than zero to not get negative solution?

syms N x1 x2 x3 x4 y1 y2
syms r1 r2 r3 r4;
syms s1 s2 s3 s4 ;
syms epsilon1 epsilon2 omega1 
syms omega2 K11 K22 
syms beta11 beta21 beta31 beta41
syms beta12 beta22 beta32 beta42 gamma12;
syms eta11 eta12 eta13 eta14 ;
syms eta21 eta22 eta23 eta24 ;
syms eta31 eta32 eta33 eta34 ;
syms eta41 eta42 eta43 eta44 ;
syms R c1 c2 c3 c4 ;
syms f11 f12 f21 f22 f31 f32 f41 f42 g12;
F1=R-c1*x1-c2*x2-c3*x3-c4*x4;
F2=x1*(r1*(1-(eta11*x1+eta12*x2+eta13*x3+eta14*x4)/N)-f11*y1-f12*y2)+s1==0;
F3=x2*(r2*(1-(eta21*x1+eta22*x2+eta23*x3+eta24*x4)/N)-f21*y1-f22*y2)+s2==0;
F4=x3*(r3*(1-(eta31*x1+eta32*x2+eta33*x3+eta34*x4)/N)-f31*y1-f32*y2)+s3==0;
F5=x4*(r4*(1-(eta41*x1+eta42*x2+eta43*x3+eta44*x4)/N)-f41*y1-f42*y2)+s4==0;      
F6=y1*(-epsilon1*(1+(y1+omega2*y2)/K11)-g12*y2+beta11*f11*x1+beta21*f21*x2+beta31*f31*x3+beta41*f41*x4)==0;
       
F7=y2*(-epsilon2*(1+(omega1*y1+y2)/K22)+gamma12*g12*y1+beta12*f12*x1+beta22*f22*x2 +beta32*f32*x3+beta42*f42*x4)==0 ;           

eqns=[F2,F3,F4,F5,F6,F7];
diary('SolutionOfEquation.txt')
S = solve(eqns,[x1,x2,x3,x4,y1,y2]);
x1=latex(S.x1)
x2=latex(S.x2)
x3=latex(S.x3)
x4=latex(S.x4)
y1=latex(S.y1)
y2=latex(S.y2)
S.x1
S.x2
S.x3
S.x4
S.y1
S.y2

[Walter Roberson]

In theory,

syms N x1 x2 x3 x4 y1 y2 real
assume(N >= 0 & x1 >= 0 & x2 >= 0 & x3 >= 0 & x4 >= 0 & y1 >= 0 & y2 >= 0)
syms r1 r2 r3 r4 positive
syms s1 s2 s3 s4 positive
syms epsilon1 epsilon2 omega1 positive 
syms omega2 K11 K22 positive
syms beta11 beta21 beta31 beta41 positive
syms beta12 beta22 beta32 beta42 gamma12 positive
syms eta11 eta12 eta13 eta14 positive
syms eta21 eta22 eta23 eta24 positive
syms eta31 eta32 eta33 eta34 positive
syms eta41 eta42 eta43 eta44 positive
syms R c1 c2 c3 c4 positive
syms f11 f12 f21 f22 f31 f32 f41 f42 g12 positive
F1=R-c1*x1-c2*x2-c3*x3-c4*x4;
F2=x1*(r1*(1-(eta11*x1+eta12*x2+eta13*x3+eta14*x4)/N)-f11*y1-f12*y2)+s1==0;
F3=x2*(r2*(1-(eta21*x1+eta22*x2+eta23*x3+eta24*x4)/N)-f21*y1-f22*y2)+s2==0;
F4=x3*(r3*(1-(eta31*x1+eta32*x2+eta33*x3+eta34*x4)/N)-f31*y1-f32*y2)+s3==0;
F5=x4*(r4*(1-(eta41*x1+eta42*x2+eta43*x3+eta44*x4)/N)-f41*y1-f42*y2)+s4==0;      
F6=y1*(-epsilon1*(1+(y1+omega2*y2)/K11)-g12*y2+beta11*f11*x1+beta21*f21*x2+beta31*f31*x3+beta41*f41*x4)==0;
       
F7=y2*(-epsilon2*(1+(omega1*y1+y2)/K22)-gamma12*g12*y1+beta12*f12*x1+beta22*f22*x2 +beta32*f32*x3+beta42*f42*x4)==0 ;           
eqns=[F2,F3,F4,F5,F6,F7];
S = solve(eqns,[x1,x2,x3,x4,y1,y2], 'returnconditions', true);
diary('SolutionOfEquation.txt')
x1=latex(S.x1)
x2=latex(S.x2)
x3=latex(S.x3)
x4=latex(S.x4)
y1=latex(S.y1)
y2=latex(S.y2)
S.x1
S.x2
S.x3
S.x4
S.y1
S.y2

In practice, this is taking a fair time to compute; I suspect it might not finish before I have to leave for an appointment.

[F.O]

Yes, It takes alot of time and the solution are very long. I got 64 points without conditions.

Answer 1

[Walter Roberson]

The calculation had finished by the time I went to have another look.

The solutions were parameterized. I have attached a .mat with the values. (To get the conditions in the form I present them here, simplify() with 'steps', 100)

x1 -> z
x2 -> z1
x3 -> z2
x4 -> z3
y1 -> z4
y2 -> z5

Four independent solutions were created:

Set #1 (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)

0 <= z
0 <= z1
0 <= z2
0 <= z3
z4 == 0
0 <= z5
K22*epsilon2 + epsilon2*z5 == K22*beta12*f12*z + K22*beta22*f22*z1 + K22*beta32*f32*z2 + K22*beta42*f42*z3
s1 == z*(f12*z5 + r1*((eta11*z + eta12*z1 + eta13*z2 + eta14*z3)/N - 1))
s2 == z1*(f22*z5 + r2*((eta21*z + eta22*z1 + eta23*z2 + eta24*z3)/N - 1))
s3 == z2*(f32*z5 + r3*((eta31*z + eta32*z1 + eta33*z2 + eta34*z3)/N - 1))
s4 == z3*(f42*z5 + r4*((eta41*z + eta42*z1 + eta43*z2 + eta44*z3)/N - 1))

Set #2: (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)

0 <= z
0 <= z1
0 <= z2
0 <= z3
0 <= z4
z5 == 0
K11*epsilon1 + epsilon1*z4 == K11*beta11*f11*z + K11*beta21*f21*z1 + K11*beta31*f31*z2 + K11*beta41*f41*z3
s1 == z*(f11*z4 + r1*((eta11*z + eta12*z1 + eta13*z2 + eta14*z3)/N - 1))
s2 == z1*(f21*z4 + r2*((eta21*z + eta22*z1 + eta23*z2 + eta24*z3)/N - 1))
s3 == z2*(f31*z4 + r3*((eta31*z + eta32*z1 + eta33*z2 + eta34*z3)/N - 1))
s4 == z3*(f41*z4 + r4*((eta41*z + eta42*z1 + eta43*z2 + eta44*z3)/N - 1))

Set #3: (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)

N ~= 0
0 <= z
0 <= z1
0 <= z2
0 <= z3
z4 == 0
z5 == 0
N*s1 == r1*z*(eta11*z - N + eta12*z1 + eta13*z2 + eta14*z3)
N*s2 == r2*z1*(eta21*z - N + eta22*z1 + eta23*z2 + eta24*z3)
N*s3 == r3*z2*(eta31*z - N + eta32*z1 + eta33*z2 + eta34*z3)
N*s4 == r4*z3*(eta41*z - N + eta42*z1 + eta43*z2 + eta44*z3)

Set #4: (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)

0 <= z
0 <= z1
0 <= z2
0 <= z3
0 <= z4
0 <= z5
K11*epsilon1 + epsilon1*z4 + K11*g12*z5 + epsilon1*omega2*z5 == K11*beta11*f11*z + K11*beta21*f21*z1 + K11*beta31*f31*z2 + K11*beta41*f41*z3
K22*epsilon2 + epsilon2*z5 + epsilon2*omega1*z4 + K22*g12*gamma12*z4 == K22*beta12*f12*z + K22*beta22*f22*z1 + K22*beta32*f32*z2 + K22*beta42*f42*z3
s1 == z*(f11*z4 + f12*z5 + r1*((eta11*z + eta12*z1 + eta13*z2 + eta14*z3)/N - 1))
s2 == z1*(f21*z4 + f22*z5 + r2*((eta21*z + eta22*z1 + eta23*z2 + eta24*z3)/N - 1))
s3 == z2*(f31*z4 + f32*z5 + r3*((eta31*z + eta32*z1 + eta33*z2 + eta34*z3)/N - 1))
s4 == z3*(f41*z4 + f42*z5 + r4*((eta41*z + eta42*z1 + eta43*z2 + eta44*z3)/N - 1))

Code:

syms N x1 x2 x3 x4 y1 y2 real
assume(N >= 0 & x1 >= 0 & x2 >= 0 & x3 >= 0 & x4 >= 0 & y1 >= 0 & y2 >= 0)
syms r1 r2 r3 r4 positive
syms s1 s2 s3 s4 positive
syms epsilon1 epsilon2 omega1 positive 
syms omega2 K11 K22 positive
syms beta11 beta21 beta31 beta41 positive
syms beta12 beta22 beta32 beta42 gamma12 positive
syms eta11 eta12 eta13 eta14 positive
syms eta21 eta22 eta23 eta24 positive
syms eta31 eta32 eta33 eta34 positive
syms eta41 eta42 eta43 eta44 positive
syms R c1 c2 c3 c4 positive
syms f11 f12 f21 f22 f31 f32 f41 f42 g12 positive
F1=R-c1*x1-c2*x2-c3*x3-c4*x4;
F2=x1*(r1*(1-(eta11*x1+eta12*x2+eta13*x3+eta14*x4)/N)-f11*y1-f12*y2)+s1==0;
F3=x2*(r2*(1-(eta21*x1+eta22*x2+eta23*x3+eta24*x4)/N)-f21*y1-f22*y2)+s2==0;
F4=x3*(r3*(1-(eta31*x1+eta32*x2+eta33*x3+eta34*x4)/N)-f31*y1-f32*y2)+s3==0;
F5=x4*(r4*(1-(eta41*x1+eta42*x2+eta43*x3+eta44*x4)/N)-f41*y1-f42*y2)+s4==0;      
F6=y1*(-epsilon1*(1+(y1+omega2*y2)/K11)-g12*y2+beta11*f11*x1+beta21*f21*x2+beta31*f31*x3+beta41*f41*x4)==0;
       
F7=y2*(-epsilon2*(1+(omega1*y1+y2)/K22)-gamma12*g12*y1+beta12*f12*x1+beta22*f22*x2 +beta32*f32*x3+beta42*f42*x4)==0 ;           
eqns=[F2,F3,F4,F5,F6,F7];
S = solve(eqns,[x1,x2,x3,x4,y1,y2], 'returnconditions', true);

[F.O]

@Walter Roberson I will cheak the solution.

In the meantime do you think that the solution will be different If i consider N as postive constant not variable ? (F1 is not used her )

[Walter Roberson]

If you assume that N is positive constant then the only thing that would change would be that in the third solution, the condition N ~= 0 would not be present, since it would be redundant when you assume N>0

[F.O]

@ Walter Roberson

I downloaded the mat file but didnt know how to get the values .

If I write S.x1 in the command window I get z1 z1 z1 z1

[Walter Roberson]

Yes. S.x1(1) is the set of all values, z such that S.conditions(1) is satisfied.

Another way of saying this is that you can take

subs(S.conditions(1), S.parameters(:),[x1 x2 x3 x4 y1 y2].')

and the result should be rewritten in terms of x1 x2 and so on.

Note that in general you would need to be more careful about the substitutions as the order in the parameters property of the output is not necessarily the same as the order of variables used in your code. You got lucky in this case. The more robust code is a nuisance to type out on my phone though.

and the output should

[F.O]

I am not sure If I understand what you ment there. What I wanted was to get solution in terms of parameters which I cant now.

>> S.x1
 
ans =
 
 z
 z
 z
 z
 
>> subs(S.conditions(1), S.parameters(:),[x1 x2 x3 x4 y1 y2].')
Unrecognized function or variable 'x1'.
 
>> z
Unrecognized function or variable 'z'.
 

[Walter Roberson]

syms x1 x2 x3 x4 y1 y2 real

[F.O]

Ok. That substutes z's with original variables in the conditions. But I want to get solution in terms of parameters like

and others

[Walter Roberson]

Your original question had

S = solve(eqns,[x1,x2,x3,x4,y1,y2], 'returnconditions', true);

which requests that x1, x2, x3, x4, y1, y2 be solved in terms of the other variables, not that the other variables be solved in terms of x1, x2, x3, x4, y1, y2.

You have 6 equations. You can solve for any 6 of the variables in terms of the other ones. Which 6 do you want to solve for ?

[F.O]

What I should get is four solutions for x in terms of parameters not variables and the same for others.

This link shows you the solution I got when not using any conditions on parameters and variables. It is just for x1 in latex format in which I got 64 solution in terms of parametes

https://www.dropbox.com/s/hoxp4bdo0pb9n0s/SolutionInTermsOfParameters.txt?dl=0

[Walter Roberson]

With your code, it is not possible for the readers to know which names you consider to be "parameters" and not "variables".

[F.O]

Well I wrote in the in the introduction that variables are N,x1,x2,x3,x4,y1, y2. I simplified my equation by taking s1,s2,s3,s4= 0 but Matlab can't solve after more than two hours. I even reduced the system by removing F6,F7 and removing terms with y1, y2 in x equations but no answer after half an hour. What ai should do to make the difference clear?

[Walter Roberson]

The answers are huge and take a long time to compute.

The equations can be reduced by solving the first for y1, substituting, solving the second for y2, substituting, solving the third for x4, substituting. After that things get messy. For example the "shortest" solution for x1 by far turns out to involve the root of a polynomial of degree 4, but if you take the explicit forms for those, simplifying them takes many hours.

[F.O]

Did you got the answer? My computer coudnt and got errer messege 'out of memory' . Also if you get the answer then it should get saved in the file determined by diary command.

[Walter Roberson]

Still working on it. It is not easy. Even just for the 4th variable, the partial solutions for one variable alone are over a megabyte.

[F.O]

Without imposing the conditton I got the solution relatviley fast and the size of the text file was 40 MB. I find out today that I have a sign mistake in eq F7 as the sign of gamma12 shoud be + not minus. I will edit my equation. I have nondimentiobalized the system aswell.

[Walter Roberson]

Computation is still going :( Solve for y1, y2, x4; solve for x1 from the 2nd remaining equation; substitute in either the 3rd or 4th solution (the first two are much larger). Now solve x2 from either of the two remaining equations, but that takes a long time and a generous amount of memory....

And once x2 is found and substituted in, there would still be the task of solving for x3, and that is likely to be a zoo.

Once the solutions are found, they are likely to be multiple gigabytes each. I have to ask what you are going to do with those solutions? You cannot publish them in a paper, because they would be page on page on page of combinations of symbols that are essentially meaningless. You aren't going to be able to get away without solving a polynomial of at least degree 4 in symbolic form; if you can read the closed form solution of a quartic and make some kind of intuitive sense of it then your mathematics would be far far superior to my own.

If you are looking for a closed form solution after substitution of particular numeric values, and so were thinking of generating the formulas and then substituting the values in, then be aware that the closed forms will certainly involve numbers with over a thousand digits, and that approximating the numbers can lead to incorrect results. When you have solutions of cubics and quartics and higher degree, the terms of the equations are sometimes very finely balanced, with discontinuities or complex-valued results if you do not calculate to full precision.

[F.O]

I agree that gettng the result may not be very useful. What I am doing here is finding the equilibrum points (solution here) because the system originally is a system of ODE. I may be be able to study the stability of som points that doesnt involve many parameters.

Do you mean if I estimated numerical values for the parameters in equations and solve the system , I will not get correct result either?

[Walter Roberson]

I seem to recall having encounted the idea that the equilibrium points of an ODE can be modeled in terms of eigenvalues of the systems. But you have 35 parameters, and the symbolic eigenvalue of a system with 35 parameters is not even close to being feasible to compute.

[F.O]

Yes, that's correct we need to determine eigenvalues to study stability(Her we linearize the system,find Jacobi and substitute for each point in Jacobi abd find eigenvalues so determine stability) but the good thing is that if we can prove that one eigenvalue is positive for each point then the point is unstable. Is it possible to find eigenvalues symbolically here in Matlab?

[Walter Roberson]

eig() can operate on symbolic expressions, but eig() of a symbolic matrix greater than 4 x 4 is rarely practical (and even 4 x 4 typically takes a couple of hours.) It might be practical for some fairly sparse matrices.

The 35 parameters will ruin any hope of getting a symbolic solution in a reasonable time.

In the case where you can express each of the elements in linear form, then typically it is much much much faster to form a matrix that looks something like

[A11 * x + B11, A12 * x + B12, ... A17 * x + B17
 A21 * x + B21, A22 * x + B22, ... A27 * x + B27
 ....
 A71 * x + B71, ....               A77 * x + B77]
 

but making sure to put in 0 for any entries that are 0. Then take the eigenvalues of that, and then subs() the corresponding expressions for A** and B**.

The eigenvalue step will probably take ... days ... for a 7 x 7, and then the subs() will not exactly be fast. And whatever you do, do not try to simplify() the result !!!

You should try with a 3 x 3 system and see if you can identify the stability of it. 3 x 3 has closed form solutions, but they are large enough to be a nuisance. 4 x 4 has closed form solutions as well, and they are a pain to work with. 5 x 5 and larger is not certain to have closed form solutions unless you can happen to have some special block or diagonal matrix properties.

Anyhow, my point is that if you cannot understand the symbolic solutions for 3 x 3 and then 4 x 4, then 7 x 7 is hopeless.

[F.O]

Thanks for beng patient. I am thinking of reducing my system to three equation( F2 , F6 ,F7,removing variables x2,x3,x4)and study stability and analysis.

[Walter Roberson]

y1 and y2 are easy to solve for, and x1 then just introduces a quadratic and so is not bad to work with. It is the solution for the 4th variable that starts getting costly, and beyond the 4th that is a real problem.

[F.O]

That is why I have changed my mind to start with 3x3 system.

Do you know why matlab can solve the system in case (s1,s2,s3,s4=0 and neglectingF1) when not imposing the condition on variables and parameters and can't do it in case of imposiing condition ?

[Walter Roberson]

Trying to prove that something is positive or find the conditions under which it is positive, can be expensive.