Solve a system of nonlinear equations with conditions

4 ビュー (過去 30 日間)
F.O
F.O 2019 年 12 月 18 日
編集済み: F.O 2019 年 12 月 23 日
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 件のコメント
Walter Roberson
Walter Roberson 2019 年 12 月 18 日
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 2019 年 12 月 18 日
編集済み: F.O 2019 年 12 月 18 日
Yes, It takes alot of time and the solution are very long. I got 64 points without conditions.

サインインしてコメントする。

採用された回答

Walter Roberson
Walter Roberson 2019 年 12 月 18 日
編集済み: Walter Roberson 2019 年 12 月 18 日
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 件のコメント
F.O
F.O 2019 年 12 月 23 日
編集済み: F.O 2019 年 12 月 23 日
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 2019 年 12 月 23 日
Trying to prove that something is positive or find the conditions under which it is positive, can be expensive.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeEquation Solving についてさらに検索

製品


リリース

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by