MATLAB Answers

Solving complicated system of symbolic equations

88 views (last 30 days)
Braden Kerr
Braden Kerr on 23 Dec 2020
Answered: J. Alex Lee on 27 Dec 2020
Hello,
I am attempting to solve a system of 3 equations that are all pretty complex and dont simplify down neatly to be used with linsolve.
The three equations at hand are:
% Given
R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
c = [-1; -1; -1];
e = [1; 1; 1];
% Create system of equations
syms x1 x2 x3 u v
x_vector = [x1; x2; x3];
eqn1 = x_vector'*R*x_vector == 1;
eqn2 = -(Q+Q'+2*u*R)^(-1)*(c+v*e) == x_vector;
eqn3 = (-1-e'*((Q+Q'+2*u*R)^-1)*c)/(e'*((Q+Q'+2*u*R)^-1)*e) == v;
We are given R, Q, c and e, but need to find u,v and x where x has components [x1; x2; x3]. I am wondering if this is possible on matlab and if there is another solver I should attempt to use.
If the equations above are confusing, Ive added photos of them typed out normally.
Edit: I wanted to add that I tried to use solve() with my three equations as the three arguments but I get an error that the second argument must be a vector of symbolic variables
  5 Comments
Braden Kerr
Braden Kerr on 23 Dec 2020
To answer your question, it was for a KKT conditions problem and our professor said we could, but it gets messy. I tried, it got messy and so in talking to him he recommended using matlab instead but I havent used matlab for such a system. Im comfortable using it for other nonlinear systems or linear systems but wasnt sure how to code this one.
Thank you, this code helped me figure it out.

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 23 Dec 2020
Edited: John D'Errico on 23 Dec 2020
Let me see....
eqn1 is a scalar, so it adds one piece of information. The same is true of eqn3. eqn2 is three pieces of information. And you have 5 unknowns. So in theory, it seems there are 5 pieces of information, and 5 unknowns. But that is not the complete story.
I would note that (Q + Q' + 2*u*R) will probably be singular for 3 distinct values of u. And that alone will be an issue.
R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
c = [-1; -1; -1];
e = [1; 1; 1];
syms u v
CP = det(Q + Q' + 2*u*R)
u_sing = solve(CP)
u_sing =
-1
- (10*181^(1/2))/59 - 20/59
(10*181^(1/2))/59 - 20/59
So you will need to avoid those places. In fact, I would choose to restrict the search to the 4 intervals defined by the breakpoints:
uintervals = [-inf;sort(double(u_sing));inf]
uintervals =
-Inf
-2.61925831306334
-1
1.94129221136843
Inf
So how would I try to solve this?
Pick ANY value of u in one of the intervals. For example, u = 1. u lies in the interval [-1, 1.9413].
As a function of u, we can solve for v. Everything else is absolutely known.
u = 1;
% As long as we avoid the bad places for u in u_sing, M is well-posed.
M = inv(Q + Q' + 2*u*R);
v = -(1 + e'*M*c)/(e'*M*e);
x = -M*(c + v*e);
% Does x satisfy the requirement that
% x'*R*x == 1?
x'*R*x
0.276636817658454
So, no. But we can view this as a problem of a ONE variable search, where we change u in each of the intervals we have identified. Find the value of u such that x'*R*x == 1, or is as close as possible to that goal.
function eqn1resid = uobj(u,Q,R,e,c)
M = inv(Q + Q' + 2*u*R);
v = -(1 + e'*M*c)/(e'*M*e);
x = -M*(c + v*e);
eqn1resid = abs(x'*R*x - 1);
end
So we wish to MINIMIZE eqn1resid, as a function of u.
[uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),-10,uintervals(2))
uopt =
-2.64230839840025
fval =
5.88144443320893e-05
exitflag =
1
[uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(2),uintervals(3))
uopt =
-1.88901086809246
fval =
3.4982453337884e-05
exitflag =
1
>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(3),uintervals(4))
uopt =
0.244333979030155
fval =
4.84817244350566e-05
exitflag =
1
Finally, in the topmost interval for u, fminbnd does not find any solution.
>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(4),100)
uopt =
1.94134412859712
fval =
0.823090976327026
exitflag =
1
The point being in all of this, is you do not need to do anything fancy. Double precision arithmetic seems to work.
Will you find an algebraic solution? Probably not. But you should be able to find numerical solutions.
  4 Comments
John D'Errico
John D'Errico on 25 Dec 2020
That Walter was able to solve the problem using symbolic methods is a testament to the power of those methods. It also points out there are almost always multiple ways to solve any given problem.
I was lazy yesterday when I wrote my answer though. One of the things I should probably have done was to plot the resulting function I was optimizing, thus plotting x'*R*x-1, as a function of u. That is my fault, because it fails one of the most important rules I tell people to do: PLOT EVERYTHING. And then plot something else. Plot it differently. Think about what you see. That plot forces the analyst to verify they have done something reasonable. Does what you see make sense?
The other thing a plot does is to help you to visualize what is happening. Is there a root on some interval? Are there multiple roots? Is your problem practically solvable?
I'll concede that one thing such an intuitive graphical approach might do is to target the user to the use of numerical solvers, but even with symbolic methods, you should still follow the rule that if you can plot what you are doing, then you should do so, before you try to solve anything.
function res = xRxm1(u,Q,R,e,c)
res = zeros(size(u));
for i = 1:numel(u)
M = inv(Q + Q' + 2*u(i)*R);
v = -(1 + e'*M*c)/(e'*M*e);
x = -M*(c + v*e);
res(i) = x'*R*x - 1;
end
end
I've removed the abs there, so we can look at what it produces more intelligently in terms of finding a zero. I also vectorized the code, so it can handle multiple values for u at once.
R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
c = [-1; -1; -1];
e = [1; 1; 1];
As I said, the intervals of interest for u should be defined by the breaks:
uintervals = [-Inf, -2.61925831306334, -1, 1.94129221136843, Inf];
The inner breaks occur where our matrix was singular, as a function of u. So this method must fail at those values for u. We can verify that fact.
>> xRxm1(-1,Q,R,e,c)
Warning: Matrix is singular to working precision.
> In xRxm1 (line 4)
ans =
NaN
Oh - there is no outer interval to consider. The first interval goes from [-inf,-2.6]. I reduced +/-inf to +/-10 in the code I played with before, because numerical solvers tend not to like infs. So the top interval is [1.94,inf].
Now, plot this function over an interval of interest. When we do so, we will see something interesting. Remember, in the plot here, we are looking for a zero crossing.
fplot(@(u) xRxm1(u,Q,R,e,c),[-.9,1.94])
So it looks like the function has a singularity around x=-0.28 or so. In fact, it almost looks like no zero crossing exists. But if we break the problem down into two further sub-intervals, we see this:
fplot(@(u) xRxm1(u,Q,R,e,c),[-.1,1.94])
yline(0);
fplot(@(u) xRxm1(u,Q,R,e,c),[-.999,-0.4])
yline(0);
So as you see, fzero would have failed to produce a root, if we just gave it the complete interval [-1,1.94]. That is because an objective function that fzero would use will produce end points that have the same sign at each end of that interval. Whereas fminbnd does produce a solution. It will find ONE of those roots.
>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(3),uintervals(4))
uopt =
0.244333979030155
fval =
4.84817244350566e-05
exitflag =
1
And now I'll force it to find the other:
>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(3),-.4)
uopt =
-0.806937125898144
fval =
4.91463840461837e-05
exitflag =
1
I could still use fzero. But I need to be careful about the bracket choice.
[uobj,fval,exitflag] = fzero(@(u) xRxm1(u,Q,R,e,c),[-.999,-.4])
uobj =
-0.806953366404047
fval =
-2.22044604925031e-16
exitflag =
1
Anyway, if we look at each of those intervals, I'd bet that there may actually be TWO such real solutions in most of the intervals I have identified. Except for the top most interval, where none seems to exist at all.
The only other point of interest is to identify why there is a singularity in the function I have been plotting. But really, I need to stop somewhere. :)
As I said at the beginning of thos coment though, I really should have plotted the thing first. That would have made my answer more complete. And vastly longer. :(

Sign in to comment.

More Answers (2)

Walter Roberson
Walter Roberson on 23 Dec 2020
Maple gives a fast response that
RR = roots([3353, 17080, 24623, 6420, - 3300])
u = RR
v = 10059/11420*RR^3 + 3843/1142*RR^2 + 32129/11420*RR + 1413/1142
x1 = -466067/2969200*RR^3 - 26173/296920*RR^2 + 2249503/2969200*RR + 151511/296920
x2 = -55085/118768*RR^3 - 117787/59384*RR^2 - 1787385/831376*RR + 108319/415688
x3 = 17723/28550*RR^3 + 11829/5710*RR^2 + 278251/199850*RR + 9159/39970
There are analytic solutions for RR, as it is a degree 4 polynomial. The numeric values are approximately [0.2443485977827867, -0.8069533664040474, -1.889019590665992, -2.642321360963269] . The analytic solutions are not worth posting as they are too long.
  2 Comments
Walter Roberson
Walter Roberson on 24 Dec 2020
I was surprised how quickly Maple solved it; I had already spent about half an hour trying to get a solution out of MATLAB by analyzing the circumstances under which the ReturnConditions could be true, and Maple responded with a complete solution almost immediately.
Maple does not return in exactly the form I posted: each variable is in terms of RootOf() expressions, occuring multiple times for some of the elements. Pulling those RootOf out into a parameter is easy though:
SOL := solve(EQN);
R := indets(SOL, RootOf)[1];
S := subs(R=RR, SOL);

Sign in to comment.


J. Alex Lee
J. Alex Lee on 27 Dec 2020
I like John's comment about visualizing. To expand on this, and Walter's comment that it is trickier to do this in the complex plane, here is a numerical exploration for complex solutions assuming is the conjugate transpose
clc;
close all;
clear
R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
cnt = [-1; -1; -1];
e = [1; 1; 1];
% real roots by Walter (for visualization)
uRoots = roots([3353, 17080, 24623, 6420, - 3300])
% singular points in u by John (for visualization
uSing = [-1, - (10*181^(1/2))/59 - 20/59, (10*181^(1/2))/59 - 20/59]
% define a grid on the complex plane
[uR,uI] = ndgrid(linspace(-3,2,200),linspace(-0.8,0.8,200));
resMat = nan(size(uR));
% calculate residual for each grid point
for k = 1:numel(uR)
u = uR(k) + 1i*uI(k);
[r,v,x] = resfun1(u , R,Q,cnt,e);
% the residual should be real because of the way it is defined (?)
resMat(k) = real(r);
end
% visualize
figure(1); cla; hold on;
% the actual contour map of residual value
contourf(uR,uI,resMat,linspace(-1,1,35)*0.8,'LineStyle','none')
% draw the curves of solution
contour(uR,uI,resMat,[0,0],'EdgeColor','w','LineWidth',1)
colormap turbo
shading interp
colorbar
set(gca,'DataAspectRatio',[1 1 1])
% plot the singular points on real line
plot(uSing,zeros(size(uSing)),'xw','LineWidth',1,'MarkerSize',6)
% plot the real roots on the real line
plot(uRoots,zeros(size(uRoots)),'ow','LineWidth',1,'MarkerSize',4)
% plot the centers of ostensible circles
cnt = [mean(uRoots(1:2)),mean(uRoots(3:4))]
rad = (uRoots([2,4])'-c)
plot(cnt,zeros(size(cnt)),'+k')
xlabel('real part of $u$','Interpreter','latex')
ylabel('imaginary part of $u$','Interpreter','latex')
title('residual $\left(1-x^T R x\right)$','Interpreter','latex')
function [r,v,x] = resfun1(u , R,Q,c,e)
Z = Q+Q'+2*u*R;
Zi = Z\eye(3);
v = -(1+e'*Zi*c)/(e'*Zi*e);
x = -Zi*(c+v*e);
r = 1 - x'*R*x; % conjugate transpose ensures real result (I think)
end
The solution looks to be circular rings about the mid-points of the 2 consecutive pairs of real solutions found by Walter. The midpoints are the singularities in the residuals found by John.
So there may be a bit nicer structure in the solutions if analyzed more deeply?

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by