How to solve 3 simultaneous algebraic equations with a equality constraint.
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
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]
sol = solve(eqn,[x1, x2, x3, t], 'maxdegree', 3)
[sol.x1, sol.x2, sol.x3, sol.t]
vpa(ans,10)
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
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]
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 = [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 23 Jun 2021
Yokuna
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]
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(:))
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
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
]
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(:))
bestx = [X1(bestidx), X2(bestidx), X3(bestidx)]
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])
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
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
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
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
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]
E2 = lhs(eqns)-rhs(eqns)
So at t = 0, 140 is the only possible solution.
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, 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
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:
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 
.
, 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.
And the information given includes three equations plus one constraint equation.
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;
EQ3 = x3 == t*x1+x2+x3;
[x1_12,x2_12] = solve(EQ1,EQ2)
[x1_13,x2_13] = solve(EQ1,EQ3)
[x1_23,x2_23] = solve(EQ2,EQ3)
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)
X1 = subs([x1_12, x1_13, x1_23], t, tsol(2))
X2 = subs([x2_12, x2_13, x2_23], t, tsol(2))
So far, so good, the pairs of equation seem to check out.
X3 = subs(x3, [x1, x2], [X1(1), X2(1)])
... 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.
Categories
Find more on Applications in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
