You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to solve 3 simultaneous algebraic equations with a equality constraint.
6 views (last 30 days)
Show older comments
If someone could help me to plot x1 vs t from the information. Please if someone could give any idea.
%Initial conditions
x1=140; x2=140; x3=140;
%Equations
x1 =t*x1+x2+t*x3;
x2 = 2*t*x1+t*x2+x3;
x3 = t*x1+x2+x3;
% Equality constraint
x1+x2+x3=420;
3 Comments
Alan Stevens
on 23 Jun 2021
As you give "initial conditions" and want a plot of x1 as a function of t, should those equatios be the following differential equations (you have written algebraic equations):
dx1/dt =t*x1+x2+t*x3;
dx2/dt = 2*t*x1+t*x2+x3;
dx3/dt = t*x1+x2+x3;
MSM Farhan
on 27 Jul 2021
N1=976:1:10000;
Q1=976:1:10000;
N=reshape (N1, 95, 95);
Q=reshape (q1, 95, 95);
Z0=4
Z1=1
f0=9.0336
f1=-0.3275
D1=10
D2=10
d1=10
d2=10
H1=10
H2=10
n1=10
K1=10
c11=N*Z1-f1/4-(3/4)*cos(2*q)*(2*f0+ f1)-2*(sin(q).^2-cos(q).^2)*D1+4*(cos(q).^2)*
(cos (q).^2-3*(sin(q).^2))*d1+H1*sin(q)+H2*cos (q)-2*n1+4*K1*sin(2*q)
c22=N*Z1-f1/4-(3/4 )*cos(2*q)*(2*f0+f1)-2*(sin(q).^2-cos(q).^2)*D2+4*(cos(q).^2)*(cos
(q).^2-3*(sin(q).^2))*d2+H1*sin(q)+H2*cos (q)-2*n1+4*K1*sin(2*q)
c21=-N*Z1+f1/4-(3/4 )*cos(2*q)*f1+2*n1
a1=-(3/4 )*sin(2*q)*(f0+f1)+D1*sin(2*q)+2*d1*(cos(q).^2)*sin(2*q)-H1*cos(q)+H2*sin(q)- 2*K1*cos(2*q)
a2=-(3/4 )*sin(2*q)*(f0+ f1)+D2*sin(2*q)+2*d2*(cos(q).^2 )*sin(2*q)-H1*cos(q)+H2*sin(q)- 2*K1*cos(2*q)
C11=(c21+c22)/(2*(c11*c22-c21.^2))
C21=(c11+c21)/(2*(c21.^2-c11*c22))
e1=(a2-a1)*C11
e2=(a2-a1)*C21
E0=-N*Z0+f0/4-N*Z1-f1/4+(3/4)*cos(2*q)*(f0+f1)-(cos(q).^2)*(D1+D2)-(cos(q).^4)*(d1+d2)-2*(H1 *sin(q)+H2*cos(q)-n1+K1*sin(2*q)
E11=-(3/4)*sin(2*q)*(f0+f1)*(e1+e2)+sin(2*q)*(D1*e1+D2*e2)+2*(cos(q).^2)*sin(2*q)* (d1*e1+d2*e2)-H1*cos (q)*(e1 + e2)+H2 *sin(q)*(e1+e2)-2*K1*(e1+e2)
E2=(1/2)*(N*Z1-f1/4)*((e1-e2).^2)-(3/8)*cos(2*q)*(2*f0*(e1.^2+e2.^2)+f1*((e1+e2).^2))-(sin(q).^2-cos(q).^2)*(D1 *e1.^2+D2*e2.^2)+2*(cos(q).^2*(cos(q).^2-3*(sin(q).^2))* (d1*e1.^2+d2*e2.^2)+(1/2)*H1*sin(q)*(e1.^2+e2.^2)+(1/2)*H2*cos(q)*(e1.^2+e2.^2)-n1* ((e1-e2).^2)+2*K1*sin(2*q)*(e1.^2+e2.^2)
E3=(1/8)*sin(2*q)*(4*f0*(e1.^3+e2.^3)+f1*((e1+e2).^3))-(4/3)*cos(q)*(sin(q)*(D1*e1.^3+ D2*e2.^3)-4*cos(q)*sin(q)*((5/3)cos(q).^2-(sin(q).^2)*(d1*e1.^3+d2*e2.^3)+(1/6)* H1*cos(q)*(e1.^3+e2.^3)-(1/6)*H2*sin (q)*(e1.^3+e2.^3)+(4/3)*K1*cos(2*q)*(e1.^3+ e2.^3)
E4=-(1/24)*(N*Z1-f1/4)*((e1-e2).^4)+(1/32)*cos(2*q)*(8* f0*(e1.^4 +e2.^4)+f1*((e1+e2).^4))-(1/3)*cos(2*q)*(D1*e1.^4+D2*e2.^4)-(1/3)*cos(2*q)*(5*cos(q).^2-3*sin(q).^2)*(d1*e1.^4+d2*e2.^4)- (1/24)*H1*sin(q)*(e1.^4+e2.^4)-(1/24)*H2*cos(q)*(e1.^4 +e2.^4)+(1/12)*n1*((e1-e2).^4)-(2/3)*K1*sin (2*q)*(e1.^4+e2.^4)
R=E0+E1+E2+E3+E4;
Figure (1)
h=surf (R);
h=xlabel (‘J/\omega‘);
h=ylabel (‘angle \theta (radians)‘);
h=Zlabel (‘E(\theta) /\omega‘);
Walter Roberson
on 27 Jul 2021
That does not appear to be related? please open a new question, and when you do please be more clear what you are asking for.
Accepted Answer
Walter Roberson
on 23 Jun 2021
You cannot usefully plot x1 vs t. Your system defines three specific sets of points, two of which are complex-valued
syms x1 x2 x3 t
eqn = [x1 == t*x1+x2+t*x3;
x2 == 2*t*x1+t*x2+x3;
x3 == t*x1+x2+x3;
x1+x2+x3==420]
eqn =

sol = solve(eqn,[x1, x2, x3, t], 'maxdegree', 3)
sol = struct with fields:
x1: [3×1 sym]
x2: [3×1 sym]
x3: [3×1 sym]
t: [3×1 sym]
[sol.x1, sol.x2, sol.x3, sol.t]
ans =

vpa(ans,10)
ans =

so the only real-valued solution is x1 about -235, x2 about 732, x3 about -76, and t about 3.1 .
You might perhaps be expecting all-positive results, but look at your equations:
x3 = t*x1+x2+x3;
x3 appears with coefficient 1 on both sides, so you can subtract it from both sides, leading to
0 == t*x1 + x2
and if t and x1 and x2 are all positive, then that equation cannot be satisfied. It can potentially be satisfied if t and x2 are both 0
If you substitute t = 0 into your first three equations, you can come out with a consistent solution only if x1 = x2 = x3 = 0. However, that does not satisfied the constraint. This establishes that there is no consistent solution for arbitrary times.
20 Comments
Yokuna
on 23 Jun 2021
Edited: Yokuna
on 23 Jun 2021
The initial conditions are x1=140; x2=140; x3=140;
The original equations are:-
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x1+x2+x3==420
I want to plot x1 vs t. Is it possible for these equations.
Walter Roberson
on 23 Jun 2021
Can we assume real numbers? It is not clear that we can, since you have expressions such as
abs((18*x2)/125 - (21*x3)/100 + 91/100)^2
but the abs() is redundant there unless x2 or x3 can be complex valued.
Walter Roberson
on 23 Jun 2021
Edited: Walter Roberson
on 23 Jun 2021
syms x1 x2 x3 t
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x1+x2+x3==420]
eqns =

E2 = lhs(eqns)-rhs(eqns)
E2 =

F = matlabFunction(E2, 'vars', {[x1,x2,x3], t})
F = function_handle with value:
@(in1,t)[in1(:,1)-t.*(abs(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1).^2.*sign(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1).*1.0e+1+sqrt(abs(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1)).*sign(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1).*1.0e+1)+t.*(abs(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).^2.*sign(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).*1.0e+1+sqrt(abs(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2)).*sign(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).*1.0e+1)-1.4e+2;in1(:,2)+t.*(abs(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1).^2.*sign(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1).*1.0e+1+sqrt(abs(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1)).*sign(in1(:,1).*(-2.4e+1./1.25e+2)+in1(:,2).*(1.8e+1./1.25e+2)+1.11e+2./5.0e+1).*1.0e+1)+t.*(abs(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2).^2.*sign(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2).*1.0e+1+sqrt(abs(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2)).*sign(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2).*1.0e+1)-t.*(abs(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).^2.*sign(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).*1.0e+1+sqrt(abs(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2)).*sign(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).*1.0e+1).*2.0-1.4e+2;in1(:,3)-t.*(abs(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2).^2.*sign(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2).*1.0e+1+sqrt(abs(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2)).*sign(in1(:,2).*(1.8e+1./1.25e+2)-in1(:,3).*(2.1e+1./1.0e+2)+9.1e+1./1.0e+2).*1.0e+1)+t.*(abs(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).^2.*sign(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).*1.0e+1+sqrt(abs(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2)).*sign(in1(:,1).*(2.4e+1./1.25e+2)-in1(:,2).*(3.6e+1./1.25e+2)+in1(:,3).*(2.1e+1./1.0e+2)-3.13e+2./1.0e+2).*1.0e+1)-1.4e+2;in1(:,1)+in1(:,2)+in1(:,3)-4.2e+2]
T = linspace(0,3,500);
nT = length(T);
sols = zeros(nT, 3);
x0 = [140, 140, 140];
options = optimoptions(@fsolve, 'Algorithm', 'levenberg-marquardt', 'display', 'none');
have_warned = false;
for tidx = 1 : nT
[thissol, ~, exitflag] = fsolve(@(x) F(x,T(tidx)), x0, options);
sols(tidx, :) = thissol;
x0 = thissol;
if exitflag <= 0 & ~have_warned
warning('solution failure code %d starting at time = %g', exitflag, T(tidx));
have_warned = true;
end
end
plot(T, sols);
legend({'x1', 'x2', 'x3'});

The values seem to be nearly steady-state after that.
Yokuna
on 26 Jun 2021
If x1 has a bound for example
50<x1<100
Then i tried it by giving it in the equtions defined above but the E2 defined in above code is causing the problem as it is LHS-RHS. Could you help me resolve it?
Walter Roberson
on 26 Jun 2021
syms x1 x2 x3 t
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x1+x2+x3==420]
eqns =

x3_partial = solve(eqns(4), x3);
e2 = subs(eqns(1:3), x3, x3_partial);
ERR = sum((lhs(e2)-rhs(e2)).^2);
ERRT = subs(ERR, t, 1/10);
ERRF = matlabFunction(ERRT, 'vars', {x1, x2});
[X1, X2] = ndgrid(linspace(50,100), linspace(-200,200));
Z = ERRF(X1, X2);
surf(X1, X2, Z); xlabel('x1'); ylabel('x2'); zlabel('error')

So obviously we can rule out negative x2 as solutions.
[X1, X2] = ndgrid(linspace(50,100), linspace(140,190));
Z = ERRF(X1, X2);
surf(X1, X2, Z); xlabel('x1'); ylabel('x2'); zlabel('error')

min(Z(:))
ans = 1.3673e+04
This tells us that for t = 1/10 that if you constrain x1 to be between 50 and 100, that you cannot get a solution over the reals.
Yokuna
on 26 Jun 2021
Yes, got it!
If x1+x2+x3==420 is not present. If we want inequality constraint on 50<x1<100. I wrote this inequality in my eqns. But it is giving the previous results only.
Walter Roberson
on 26 Jun 2021
format long g
syms x1 x2 x3 t
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
]
eqns =

err = sum((lhs(eqns)-rhs(eqns)).^2);
T = 1/100;
errf = matlabFunction(subs(err,t,T), 'vars', [x1, x2, x3]);
[X1, X2, X3] = ndgrid(linspace(50,100,100), linspace(130,150,100), linspace(120,140,100));
ERR = errf(X1, X2, X3);
[besterr, bestidx] = min(ERR(:))
besterr =
1615.17056667144
bestidx =
725300
bestx = [X1(bestidx), X2(bestidx), X3(bestidx)]
bestx = 1×3
100 140.505050505051 134.545454545455
Notice the besterr is appreciably above 0. This tells us that the three equations cannot be simultaneously satisfied.
You can extend X2 and X3 into negatives and up to (say) 500, but it doesn't help.
errf2 = @(x) errf(x(1),x(2),x(3));
[betterx, bettererr] = fmincon(errf2, bestx, [], [], [], [], [50, -inf, -inf], [100, inf, inf])
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
betterx = 1×3
99.9999998947007 140.47716303407 134.488125323314
bettererr =
1615.16508546622
Walter Roberson
on 26 Jun 2021
Create vectors
lb = [50, 50, 50]
ub = [200, 200, 200]
Now if you are taking the ndgrid() approach to create a grid of points to evaluate the function at, then replace
[X1, X2, X3] = ndgrid(linspace(50,100,100), linspace(130,150,100), linspace(120,140,100));
with
[X1, X2, X3] = ndgrid(linspace(lb(1),ub(1),100), linspace(lb(2),ub(2),100), linspace(lb(3),ub(3),100));
and if you are using fmincon then
[betterx, bettererr] = fmincon(errf, x0, [], [], [], [], lb, ub)
It is also possible to use constraints with solve(), but you have to be careful with interpreting the results, as solve() will typically just tell you one "representative" value out of the possible values. For example if the internal logic of solve was able to calculate that 1/10 < x1 < 22/7 then solve() would hide that information and would instead report back that x1 is 1.
To use constraints with vpasolve() use
vpasolve(EQUATIONS, LIST_OF_VARIABLES, [lb(:), ub(:)])
LIST_OF_VARIABLES is needed so that it knows which variable corresponds to which bound. The [lb(:), ub(:)] you are constructing is an array with two columns, with the first column being lower bound and the second column being upper bound, and the rows corresponding to variables.
Yokuna
on 27 Jun 2021
I want to plot x1 vs t. I used fmincon but not getting how to proceed after this stage.
close all
clear all
clc
x(1)=140; x(2)=140; x(3)=140; x0=[140;140;140];
D=420;
k1=10; k2=10; mu=0.5; nu=2;lamda2=1;n=3;
alpha(1)=51; alpha(2)=31; alpha(3)=78; beta(1)=1.22; beta(2)=3.44; beta(3)=2.53; gama(1)=0.096; gama(2)=0.072; gama(3)=0.105;
sigma=0.144;
syms x1 x2 x3 t
x1 =t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140;
x2 = 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140;
x3 = t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140;
lb = [50, 50, 50 ]
ub = [200, 200, 200 ]
[X1, X2, X3] = ndgrid(linspace(lb(1),ub(1),100), linspace(lb(2),ub(2),100), linspace(lb(3),ub(3),100 ));
[betterx, bettererr] = fmincon(@myfun,x0,[],[],[],[],lb,ub,@nonlcon)
%% Function I want to mimimize
function f = myfun(x)
f=gama(1)*x1^2+beta(1)*x1+alpha(1)+gama(2)*x2^2+beta(2)*x2+alpha(2)+gama(3)*x3^2+beta(3)*x3+alpha(3);
end
%% Equality constraints
function [c,ceq]=nonlcon(x1,x2,x3)
c = [];
ceq(1)=x1 - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - 140;
ceq(2)=x2 + t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - 140;
ceq(3)=x3 - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) + t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - 140;
end
Walter Roberson
on 28 Jun 2021
f=gama(1)*x1^2+beta(1)*x1+alpha(1)+gama(2)*x2^2+beta(2)*x2+alpha(2)+gama(3)*x3^2+beta(3)*x3+alpha(3);
You use gama in the function but you do not pass it in. The only way that can work is if gama is a shared variable in nested functions, but you did not code that way.
Yokuna
on 28 Jun 2021
Changed my function also but getting the error of having "Not enough input arguments."
function f = myfun(x1,x2,x3,t)
alpha(1)=51; alpha(2)=31; alpha(3)=78; beta(1)=1.22; beta(2)=3.44; beta(3)=2.53; gama(1)=0.096; gama(2)=0.072; gama(3)=0.105;
f=gama(1)*x1.^2+beta(1)*x1+alpha(1)+gama(2)*x2.^2+beta(2)*x2+alpha(2)+gama(3)*x3.^2+beta(3)*x3+alpha(3);
end
function [c,ceq]=nonlcon(x1,x2,x3,t)
c = [];
ceq(1)=x1 - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - 140;
ceq(2)=x2 + t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - 140;
ceq(3)=x3 - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) + t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - 140;
end
Walter Roberson
on 28 Jun 2021
Your new problem cannot be solved. Your constructs rely upon t, but your code does not define t.
x1 =t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140;
I already showed you that you need to transform those into equations.
Walter Roberson
on 28 Jun 2021
function f = myfun(x1,x2,x3,t)
Remember I showed you
errf = matlabFunction(subs(err,t,T), 'vars', [x1, x2, x3]);
This was to generate a function handle for an anonymous function using symbolic equations as a start. You should stil be doing that, even if the equations changed.
Also notice where I did
errf2 = @(x) errf(x(1),x(2),x(3));
[betterx, bettererr] = fmincon(errf2, bestx, [], [], [], [], [50, -inf, -inf], [100, inf, inf])
That errf2 is important: it wraps a function that expects multiple variables (like your myfun does) inside a function that expects only a single variable that is a vector. It takes the vector of values and distributes them as separate inputs into the function that expects independent variables. This is why I did not get errors about not enough input arguments.
Likewise, the nonlinear constraints function that you tell fmincon() to use, will be passed a vector of values. If your actual function uses distinct parameters, then you need the same kind of wrapper to separate out the values and pass them as distinct parameters.
Yokuna
on 1 Jul 2021
In the code you gave, it seems independent of initial condition x0. whatever we take x0, it always takes [140,140,140]. Why is it so? I tried with [170,110,140] which even does not hold equation no. 4 in eqns.
syms x1 x2 x3 t
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x1+x2+x3==420]
E2 = lhs(eqns)-rhs(eqns)
F = matlabFunction(E2, 'vars', {[x1,x2,x3], t})
T = linspace(0,3,500);
nT = length(T);
sols = zeros(nT, 3);
x0 = [170, 110, 140];
options = optimoptions(@fsolve, 'Algorithm', 'levenberg-marquardt', 'display', 'none');
have_warned = false;
for tidx = 1 : nT
[thissol, ~, exitflag] = fsolve(@(x) F(x,T(tidx)), x0, options);
sols(tidx, :) = thissol;
x0 = thissol;
if exitflag <= 0 & ~have_warned
warning('solution failure code %d starting at time = %g', exitflag, T(tidx));
have_warned = true;
end
end
plot(T, sols);
legend({'x1', 'x2', 'x3'});
Walter Roberson
on 1 Jul 2021
For any given t value, there is only one real-valued solution. The initial guess, x0, just helps guide it toward the (one) solution.
I suspect you are talking about t = 0. So let's see:
t = 0
t = 0
syms x1 x2 x3
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x1+x2+x3==420]
eqns =

E2 = lhs(eqns)-rhs(eqns)
E2 =

So at t = 0, 140 is the only possible solution.
Yokuna
on 17 Jul 2021
Edited: Yokuna
on 17 Jul 2021
Hi, I want to solve 5 equations instead of 4, but fsolve gives error. Could you please help?
syms x1 x2 x3 x_L t
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x_L==0.002*x1*x1
x1+x2+x3==420+x_L]
E2 = lhs(eqns)-rhs(eqns)
F = matlabFunction(E2, 'vars', {[x1,x2,x3,x_L], t});
T = linspace(0,3,500);
nT = length(T);
sols = zeros(nT, 3);
x0 = [140, 140, 140];
options = optimoptions(@fsolve, 'Algorithm', 'levenberg-marquardt', 'display', 'none');
have_warned = false;
for tidx = 1 : nT
[thissol, ~, exitflag] = fsolve(@(x) F(x,T(tidx)), x0, options);
sols(tidx, :) = thissol;
x0 = thissol;
if exitflag <= 0 & ~have_warned
warning('solution failure code %d starting at time = %g', exitflag, T(tidx));
have_warned = true;
end
end
plot(T, sols);
legend({'x1', 'x2', 'x3'});
Walter Roberson
on 17 Jul 2021
syms x1 x2 x3 x_L t
eqns = [
x1 == t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x2 == 2*t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) - t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^2*sign((18*x2)/125 - (24*x1)/125 + 111/50) + 10*abs((18*x2)/125 - (24*x1)/125 + 111/50)^(1/2)*sign((18*x2)/125 - (24*x1)/125 + 111/50)) + 140
x3 == t*(10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^2*sign((18*x2)/125 - (21*x3)/100 + 91/100) + 10*abs((18*x2)/125 - (21*x3)/100 + 91/100)^(1/2)*sign((18*x2)/125 - (21*x3)/100 + 91/100)) - t*(10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^2*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100) + 10*abs((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)^(1/2)*sign((24*x1)/125 - (36*x2)/125 + (21*x3)/100 - 313/100)) + 140
x_L==0.002*x1*x1
x1+x2+x3==420+x_L]
eqns =

E2 = lhs(eqns)-rhs(eqns)
E2 =

F = matlabFunction(E2, 'vars', {[x1,x2,x3,x_L], t});
T = linspace(0,3,500);
nT = length(T);
sols = zeros(nT, 4);
x0 = [140, 140, 140, 140^2*0.002];
options = optimoptions(@fsolve, 'Algorithm', 'levenberg-marquardt', 'display', 'none');
have_warned = false;
have_warned2 = false;
for tidx = 1 : nT
[thissol, ~, exitflag, output] = fsolve(@(x) F(x,T(tidx)), x0, options);
sols(tidx, :) = thissol;
x0 = thissol;
if exitflag <= 0
if ~have_warned
warning('solution failure code %d starting at time = %g', exitflag, T(tidx));
have_warned = true;
disp(output)
end
elseif ~have_warned2
warning('Solutions resumed at time = %g', T(tidx));
have_warned2 = true;
end
end
Warning: solution failure code -2 starting at time = 0
iterations: 3
funcCount: 20
stepsize: 0.0071
cgiterations: []
firstorderopt: 2.1854e-04
algorithm: 'levenberg-marquardt'
message: '↵No solution found.↵↵fsolve stopped because the last step was ineffective. However, the vector of function↵values is not near zero, as measured by the value of the function tolerance. ↵↵<stopping criteria details>↵↵fsolve stopped because the sum of squared function values, r, is changing by less ↵ than options.FunctionTolerance = 1.000000e-06 relative to its initial value.↵However, r = 3.670370e+02, exceeds sqrt(options.FunctionTolerance) = 1.000000e-03.↵↵'
plot(T, sols);
legend({'x1', 'x2', 'x3', 'x_L'});

Despite the graph, the messages tell you there was no solution.
Walter Roberson
on 18 Jul 2021
There does not appear to be any solution for that system of equations.
More Answers (1)
ILDEBERTO DE LOS SANTOS RUIZ
on 21 Jul 2021
You only need to express
in terms of t and plot that relationship:

syms x1 x2 x3 t
x3 = solve(x1+x2+x3 == 420,x3)
EQ1 = x1 == t*x1+x2+t*x3;
EQ2 = x2 == 2*t*x1+t*x2+x3;
[x1,x2] = solve(EQ1,EQ2)
ezplot(x1,[0,2])

5 Comments
Walter Roberson
on 21 Jul 2021
you missed x3 = t*x1+x2+x3;
There are four equations in four variables (including t). There is one solution that is all real and two that involve complex values. Because you have the same number of equations and variables and there are no redundant equations, you cannot end up with t as a free parameter.
ILDEBERTO DE LOS SANTOS RUIZ
on 21 Jul 2021
Edited: ILDEBERTO DE LOS SANTOS RUIZ
on 21 Jul 2021
The original post (from Shiv) does not ask to find a unique solution for
, but only to plot
versus t. Due to the constraint,
is not a true variable, since it results from an affine combination of
and
.





Walter Roberson
on 21 Jul 2021
The original post was created thinking that there would be a general solution to the system, that for each t value there would be an x1, x2, x3 value that would satisfy all of the equations and the constraint as well.
Your solution satisfies the constraint, and the first two equations, but not the third equation. You "pick and choose" which of the equations to use; you would have gotten a different answer if you had picked the second and third equations instead of the first and second.
ILDEBERTO DE LOS SANTOS RUIZ
on 21 Jul 2021
Edited: ILDEBERTO DE LOS SANTOS RUIZ
on 21 Jul 2021
If someone could help me to plot x1 vs t from the information.
Walter Roberson
on 21 Jul 2021
And the information given includes three equations plus one constraint equation.
syms x1 x2 x3 t
x3 = solve(x1+x2+x3 == 420,x3)
x3 = 

EQ1 = x1 == t*x1+x2+t*x3;
EQ2 = x2 == 2*t*x1+t*x2+x3;
EQ3 = x3 == t*x1+x2+x3;
[x1_12,x2_12] = solve(EQ1,EQ2)
x1_12 =

x2_12 =

[x1_13,x2_13] = solve(EQ1,EQ3)
x1_13 =

x2_13 =

[x1_23,x2_23] = solve(EQ2,EQ3)
x1_23 =

x2_23 =

ezplot(x1_12,[-2,5])
hold on
ezplot(x1_13,[-2 5])
ezplot(x1_23,[-2 5])
hold off
legend({'EQ1,EQ2', 'EQ1,EQ3', 'EQ2,EQ3'}, 'location', 'southwest');

Three very different pairwise answers. It looks like there might be a common answer near t = 3; let us see:
ezplot(x1_12,[2.5,3.5])
hold on
ezplot(x1_13,[2.5,3.5])
ezplot(x1_23,[2.5,3.5])
hold off
legend({'EQ1,EQ2', 'EQ1,EQ3', 'EQ2,EQ3'}, 'location', 'southwest');

tsol = vpasolve(x1_12 == x1_13, 3)
tsol =

X1 = subs([x1_12, x1_13, x1_23], t, tsol(2))
X1 = 

X2 = subs([x2_12, x2_13, x2_23], t, tsol(2))
X2 = 

So far, so good, the pairs of equation seem to check out.
X3 = subs(x3, [x1, x2], [X1(1), X2(1)])
X3 = 

... which is the solution from the first row of solutions I posted in https://www.mathworks.com/matlabcentral/answers/862790-how-to-solve-3-simultaneous-algebraic-equations-with-a-equality-constraint#answer_731385
That is, there is only one real-valued solution to all of the equations simultaneously.
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)