88 views (last 30 days)

Show older comments

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

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.

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. :(

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.

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

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?

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

Start Hunting!