Solving complicated system of symbolic equations
22 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Braden Kerr
el 23 de Dic. de 2020
Respondida: J. Alex Lee
el 27 de Dic. de 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/469097/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/469102/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/469107/image.png)
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 comentarios
J. Alex Lee
el 23 de Dic. de 2020
Editada: J. Alex Lee
el 24 de Dic. de 2020
Any reason to expect analytical solutions with symbolics?
Here are 3 more real solutions (obtained numerically), according to Walter's answer it is the complete set
x1 = -0.177337903314
x2 = 0.370352455572
x3 = 0.806985447758
u = -1.889019590663
v = 1.993495950432
x1 = 0.687844023858
x2 = -0.389942909133
x3 = 0.702098885193
u = 0.244348597712
v = 2.138523547078
x1 = -0.076001268894
x2 = 0.947579094374
x3 = 0.128422174517
u = -0.806953366403
v = 0.695475571311
Respuesta aceptada
John D'Errico
el 23 de Dic. de 2020
Editada: John D'Errico
el 23 de Dic. de 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 comentarios
J. Alex Lee
el 25 de Dic. de 2020
Editada: J. Alex Lee
el 25 de Dic. de 2020
Good points. I guess still, depending on the minimization algorithm, you wouldn't necessarily be able to find all roots within a particular interval if there are multiple ones. I guess case in point, this method missed the root u = -0.80695 in the 3rd interval.
I'm fascinated that the system is reducible to a polynomial in u, as Walter showed with Maple, which must be the better answer for this specific problem...but numerically (if we didn't know the equation for u was polynomial with determinable coefficients) I guess one can scan over some reasonable range of u to look for cross-overs, especially if some more analysis of the equations and known values could narrow the range from [-Inf,Inf] to something finite.
Any way to do better than the [-Inf,Inf] outer bracket?
John D'Errico
el 25 de Dic. de 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])
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/470712/image.jpeg)
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/470717/image.jpeg)
fplot(@(u) xRxm1(u,Q,R,e,c),[-.999,-0.4])
yline(0);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/470722/image.jpeg)
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. :(
Más respuestas (2)
Walter Roberson
el 23 de Dic. de 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 comentarios
J. Alex Lee
el 24 de Dic. de 2020
so, maple just spat out that the problem reduces a 4th order polynomial in u?
and also, that catches a mistake of one of my invalid solutions since you show there should only be 4 real ones
Walter Roberson
el 24 de Dic. de 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);
J. Alex Lee
el 27 de Dic. de 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/471507/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/471512/image.png)
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?
0 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!