MATLAB Answers

Solve a system of nonlinear equations with conditions

18 views (last 30 days)
F.O
F.O on 18 Dec 2019
Edited: F.O on 23 Dec 2019
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

  2 Comments

Walter Roberson
Walter Roberson on 18 Dec 2019
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
F.O on 18 Dec 2019
Yes, It takes alot of time and the solution are very long. I got 64 points without conditions.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 18 Dec 2019
Edited: Walter Roberson on 18 Dec 2019
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);

  24 Comments

Walter Roberson
Walter Roberson on 23 Dec 2019
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
F.O on 23 Dec 2019
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
Walter Roberson on 23 Dec 2019
Trying to prove that something is positive or find the conditions under which it is positive, can be expensive.

Sign in to comment.

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!