How to solve 3 simultaneous algebraic equations with a equality constraint.

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;

[Alan Stevens]

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]

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]

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.

Answer 1

[Walter Roberson]

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.

[Yokuna]

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]

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]

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]

Thank you very much. It worked.

[Yokuna]

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]

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]

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]

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

[Yokuna]

Thanks. I understood the problem in my code. I want to learn how to use constraints with my equation and we dont have equation x1+x2+x3==420.

So let the constraints are 50<x1<200, 50<x2<200, 50<x3<200. So it will give me an idea.

[Walter Roberson]

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]

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]

 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]

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]

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]

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]

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]

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]

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]

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]

There does not appear to be any solution for that system of equations.

Answer 2

[Ildeberto de los Santos Ruiz]

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])

[Walter Roberson]

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]

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]

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]

Original post from @Shiv:

If someone could help me to plot x1 vs t from the information.

[Walter Roberson]

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.