Symbolic solve with user-specified precision

Hi all,

I have 9 polynomial equations i would love to simplify using

solve(equ,  [beta_bar sigma_bar phi_bar theta_bar phi_pi_bar phi_y_bar g_23],'Real',true,'ReturnConditions', true)

where

equ =

[0 < beta_bar, beta_bar < 1, 0 < sigma_bar, 0 < phi_bar, 0 < phi_pi_bar, 0 < phi_y_bar, g_23^2 - 6835969944064461/9007199254740992 == 0, phi_y_bar - 3907/20000 == 0, phi_pi_bar - 14709/10000 == 0, g_23*theta_bar - (7846842892884325*theta_bar)/9007199254740992 == 0, (11137*g_23)/10000 - (7846842892884325*phi_bar)/9007199254740992 + g_23*phi_bar - 1092378616225659/1125899906842624 == 0, (6250*phi_bar)/6197 - (2659174770252075*theta_bar)/1125899906842624 - (12447*phi_bar*theta_bar)/6197 + phi_bar*theta_bar^2 + (11137*theta_bar^2)/10000 + 55685/49576 == 0, sigma_bar - 11137/10000 == 0, beta_bar*theta_bar - (6197*theta_bar)/6250 == 0, (11137*beta_bar)/10000 - (6197*phi_bar)/6250 + beta_bar*phi_bar - 1243281529372025/1125899906842624 == 0]

so they would look nicer and cleaner, then i can feed them to a optimizer to find each parameter's maximum/minimum. The reason i don't do this simple exercise by hand is i have a lot of similar equations to solve sequentially, so i would love to automate this procedure.

As you can see, there should be infinite solutions and i want parametrized full solutions (so maybe vpasolve is not a good idea?), but matlab threw me a empty solution set. I think it might be because of a precision reason.

Is there any way to get around this ?

Thank you!

[Kylekk]

Even if i run something as simple as

syms x;
equ=[x^2==2, x==1.4142];
equ=vpa(equ,4);
solve(equ,[x],'Real',true,'ReturnConditions', true)

ans = struct with fields:
             x: [0×1 sym]
    parameters: [1×0 sym]
    conditions: [0×1 sym]

I will have an empty solution, what can i do here?

[Torsten]

syms x
equ = x^2 == 2;
sol = solve(equ,x)

sol = 

[Kylekk]

Thank you Torsten, my point is, how to avoid this type of problem in general? I gave the sqrt(2) case just as an example. My actual problem is a little more than dropping a constraint.

[Torsten]

Inequalities are no equations. So you can't specify inequalities in "solve". Maybe you can use "assume" to set conditions on the solution.

The number of equations must be equal to the number of unknowns.

Try to run your code after the necessary changes.

[Kylekk]

I don't think either of the statements is true, e.g.

syms x
solve([x==1, x^2==1, x^3==1, x>0],x,'Real',true,'ReturnConditions', true)

ans = struct with fields:
             x: 1
    parameters: [1×0 sym]
    conditions: symtrue

I believe the problem i have has to do with the numerical precisions, not the setups.

Answer 1

[John D'Errico]

Your belief has preciously little value. Faith in mathematics tends to be a waste of mental energy. Let me expand your equations, looking at them one at a time.

0 < beta_bar

beta_bar < 1

0 < sigma_bar

0 < phi_bar

0 < phi_pi_bar

0 < phi_y_bar

So those are simply lower and upper bounds on the problem. Solve does not understand inequalities, because they never give a unique solution.

syms g_23 phi_y_bar phi_pi_bar theta_bar phi_bar sigma_bar beta_bar
g_23 = solve(g_23^2 - 6835969944064461/9007199254740992 == 0)

g_23 = 

That gave us two potential solutions for g_23. One is positive, one negative, but you never said anything about what g_23 should be.

phi_y_bar = solve(phi_y_bar - 3907/20000 == 0)

phi_y_bar = 

phi_pi_bar = solve(phi_pi_bar - 14709/10000 == 0)

phi_pi_bar = 

Those pair of equations simply assign values to variables. Why would you put lower and upper bounds on something when you then assign it a value?

Next, we see the equation for theta_bar. It has a trivial solution, regardless of the value of g_23.

theta_bar = solve(g_23*theta_bar - (7846842892884325*theta_bar)/9007199254740992 == 0,theta_bar)

theta_bar = 

0

Yep. theta_bar is just ZERO. You can see that, if you factor out theta_bar.

Next, we see an equation for phi_bar. I'll tell vpa to just write out the equation itself. Remember, that you told us that phi_bar MUST be positive.

vpa((11137*g_23)/10000 - (7846842892884325*phi_bar)/9007199254740992 + g_23*phi_bar - 1092378616225659/1125899906842624 == 0,5)

ans = 

Do you see the first case? That tells us that phi_bar is NEGATIVE. Or the second term offers the possibility that phi_bar is EXACTLY zero. But you told us that phi_bar MUST be positive. Neither of the solutions are valid.

The fact is, we can stop right now. There are NO solutions to your problem.

While I could look at the rest of your problem, there is no need to proceed further. Well, I could write them down. This next one again tells us that phi_bar is NEGATIVE. But since we know it is not allowed to be negative...

(6250*phi_bar)/6197 - (2659174770252075*theta_bar)/1125899906842624 - (12447*phi_bar*theta_bar)/6197 + phi_bar*theta_bar^2 + (11137*theta_bar^2)/10000 + 55685/49576 == 0

ans = 

sigma_bar - 11137/10000 == 0

ans = 

beta_bar*theta_bar - (6197*theta_bar)/6250 == 0

ans = 

(11137*beta_bar)/10000 - (6197*phi_bar)/6250 + beta_bar*phi_bar - 1243281529372025/1125899906842624 == 0

ans = 

Again, faith is a waste of mental energy. This is simply not a precision problem, but a question of existence at all.

[Kylekk]

Appreciate your answer, yet I will show you why i was not blindly believing. Again, i want approximate results, and there's infinite of them.

Forget about the bounds for the moment, without them solve() will still give me empty solution.

Let me show you why there should be infinite sols, one equation at a time.

Round all the fractions to 4 digits. I use ~ to indicate approximately equal.

1. g_23 can take two values +/- 0.8712, take 0.8712 for example
2. phi_y_bar = 0.1953,
3. phi_pi_bar=1.4709,
4. 7846842892884325/9007199254740992 ~ 0.8712 ~ g_23, hence the precision point i mentioned.
5. (11137*g_23)/10000 - (7846842892884325*phi_bar)/9007199254740992 + g_23*phi_bar - 1092378616225659/1125899906842624 == 0 gives nothing useful because 1092378616225659/1125899906842624 ~ (11137*g_23)/10000 so this equation always holds!
6. (6250*phi_bar)/6197 - (2659174770252075*theta_bar)/1125899906842624 - (12447*phi_bar*theta_bar)/6197 + phi_bar*theta_bar^2 + (11137*theta_bar^2)/10000 + 55685/49576 == 0 is the equation to be actually solved. For each phi_bar, theta_bar has a quadratic form, hence a continuum of solutions.
7. sigma_bar = 1.1137,
8. beta_bar*theta_bar - (6197*theta_bar)/6250 == 0 gives beta_bar= 0.9915,
9. (11137*beta_bar)/10000 - (6197*phi_bar)/6250 + beta_bar*phi_bar - 1243281529372025/1125899906842624 == 0 always holds if you simply plug in beta_bar.

Now, what bound constraint did i violate above?

It is good to have you spend so much time typing the results, but check your facts before you become personal.

The newly released solve function allows inequalitie and more equations than unknowns, for which the example I've already commented above.

[Torsten]

Use "lsqnonlin" to get the best approximate solution to a system of equations that cannot be satisfied exactly.

"solve" is not useful for this purpose.

[John D'Errico]

Lets go back over it.

Take the positive solution for g_23.

syms g_23 phi_y_bar phi_pi_bar theta_bar phi_bar sigma_bar beta_bar
g_23 = solve(g_23^2 - 6835969944064461/9007199254740992 == 0);
g_23 = g_23(2)

g_23 = 

phi_y_bar = solve(phi_y_bar - 3907/20000 == 0)

phi_y_bar = 

phi_pi_bar = solve(phi_pi_bar - 14709/10000 == 0)

phi_pi_bar = 

No problem so far. And no need to worry yet about precision. Next, we see this equation:

vpa(g_23*theta_bar - (7846842892884325*theta_bar)/9007199254740992 == 0)

ans = 

That tells us that theta bar is either zero, or that the coefficient of theta_bar is zero to within floating point trash. In the latter case, we have no information about theta_bar, so this equation is useless. That either means you have one less equation than you think you have, or it means theta_bar is identically zero.

Now look at the next equation.

(11137*g_23)/10000 - (7846842892884325*phi_bar)/9007199254740992 + g_23*phi_bar - 1092378616225659/1125899906842624 == 0

ans = 

Perhaps this is not so obvious in that form, so let me expand it using VPA.

vpa((11137*g_23)/10000 - (7846842892884325*phi_bar)/9007199254740992 + g_23*phi_bar - 1092378616225659/1125899906842624 == 0,5)

ans = 

Ok, so this equation is also completely useless. We now have two of your 9 equations that offer no information, that have no information content.

At this point, we can stop again, since solve cannot handle underdetermined systems. We can look at the next equation, which has phi_bar and theta_bar.

vpa((6250*phi_bar)/6197 - (2659174770252075*theta_bar)/1125899906842624 - (12447*phi_bar*theta_bar)/6197 + phi_bar*theta_bar^2 + (11137*theta_bar^2)/10000 + 55685/49576 == 0,5)

ans = 

Lets plot that equation. It should be a hyperbola, just looking at the coefficients.

fimplicit((6250*phi_bar)/6197 - (2659174770252075*theta_bar)/1125899906842624 - (12447*phi_bar*theta_bar)/6197 + phi_bar*theta_bar^2 + (11137*theta_bar^2)/10000 + 55685/49576)
xlabel phi_bar
ylabel theta_bar % I think I got the order correct, but not really important.

sigma_bar = solve(sigma_bar - 11137/10000 == 0)

sigma_bar = 

But, next we see this equation. Since theta_bar can take on ANY value, if theta_bar is zero, then the equation is also useless. If theta_bar is non-zero, then it gives a solution for beta_bar.

beta_bar*theta_bar - (6197*theta_bar)/6250 == 0

ans = 

beta_bar = solve(beta_bar*theta_bar - (6197*theta_bar)/6250 == 0,beta_bar)

beta_bar = 

ans =

(11137*beta_bar)/10000 - (6197*phi_bar)/6250 + beta_bar*phi_bar - 1243281529372025/1125899906842624 == 0

ans = 

All of the coefficients, assuming we are willing to round things off, reduce to 0==0. So that equation tells us nothing at all.

In the end, there is nothing for solve to return. Fact. Not belief. And this still is not a precision problem.

[Kylekk]

Thank you Torsten! If i understand correctly you mean to sum the square of the equations and minimize the overall distance between the LHS of equations and zeros? I will try that!

The reason why i was so stuck at solve is it has this nice output structure as [y1,...,yN,parameters,conditions], where optimization based on the free parameters can be extremely fast and accurate. All the other optimization solvers i've tried seem to behave quite poorly when solutions are on a lower dimensional manifold.

[Kylekk]

John D'Errico, again, appreciate your effort, yet you seem not very familiar with the newly released solve function. It allows infinite solutions! (In your words, underdetermined systems!)

So if i have 3 vars and only 1 equation, it will most likely throw one into the parameters structure and let it move between the domain allowed by bounds. Try this in your own computer.

syms x y z
equ=[x^2+y==1, z>1];
sol=solve(equ,[x,y,z],'Real',true,'ReturnConditions', true );
sol.x

ans = 

sol.y

ans = 

sol.z

ans = 

sol.parameters

ans = 

sol.conditions

ans = 

I guess this should put an end to your belief discussion.

Answer 2

[Paul]

Hi @Kylekk

Perhaps I don't fully understand the Question, but I don't think there's a way to force the Symbolic Toolbox to find an "approximate" solution, if that's the desire.

The equations to be solved are:

syms g_23 phi_y_bar phi_pi_bar phi_bar theta_bar sigma_bar beta_bar real
equ(1) = g_23^2 - 6835969944064461/9007199254740992 == 0;
equ(2) = phi_y_bar - 3907/20000 == 0;
equ(3) = phi_pi_bar - 14709/10000 == 0;
equ(4) = g_23*theta_bar - (7846842892884325*theta_bar)/9007199254740992 == 0;
equ(5) = (11137*g_23)/10000 - (7846842892884325*phi_bar)/9007199254740992 + g_23*phi_bar - 1092378616225659/1125899906842624 == 0;
equ(6) = (6250*phi_bar)/6197 - (2659174770252075*theta_bar)/1125899906842624 - (12447*phi_bar*theta_bar)/6197 + phi_bar*theta_bar^2 + (11137*theta_bar^2)/10000 + 55685/49576 == 0;
equ(7) = sigma_bar - 11137/10000 == 0;
equ(8) = beta_bar*theta_bar - (6197*theta_bar)/6250 == 0;
equ(9) = (11137*beta_bar)/10000 - (6197*phi_bar)/6250 + beta_bar*phi_bar - 1243281529372025/1125899906842624 == 0;

The constraints on the variables are

assume([0 < beta_bar, beta_bar < 1, 0 < sigma_bar, 0 < phi_bar, 0 < phi_pi_bar, 0 < phi_y_bar]);

As noted, solve can't find a solution to all nine equations under the stated assumptions.

sol = solve(equ,[beta_bar sigma_bar phi_bar theta_bar phi_pi_bar phi_y_bar g_23],'Real',true,'ReturnConditions',true)

sol = struct with fields:
      beta_bar: [0×1 sym]
     sigma_bar: [0×1 sym]
       phi_bar: [0×1 sym]
     theta_bar: [0×1 sym]
    phi_pi_bar: [0×1 sym]
     phi_y_bar: [0×1 sym]
          g_23: [0×1 sym]
    parameters: [1×0 sym]
    conditions: [0×1 sym]

solve can't find a solution without the assumptions

sol = solve(equ,[beta_bar sigma_bar phi_bar theta_bar phi_pi_bar phi_y_bar g_23],'Real',true,'ReturnConditions',true,'IgnoreProperties',true)

sol = struct with fields:
      beta_bar: [0×1 sym]
     sigma_bar: [0×1 sym]
       phi_bar: [0×1 sym]
     theta_bar: [0×1 sym]
    phi_pi_bar: [0×1 sym]
     phi_y_bar: [0×1 sym]
          g_23: [0×1 sym]
    parameters: [1×0 sym]
    conditions: [0×1 sym]

However, there is a unique solution when considering all but equation 6

sol = solve(equ([1:5 7:9]),[sigma_bar beta_bar phi_bar theta_bar phi_pi_bar phi_y_bar g_23],'Real',true,'ReturnConditions',true)

sol = struct with fields:
     sigma_bar: 11137/10000
      beta_bar: (6723*13671939888128922^(1/2))/24657920000000 + 2540765228537104417/2647623999684608000
       phi_bar: (1262485504*13671939888128922^(1/2))/211445338090507825 - 2196698527725305016/5286133452262695625
     theta_bar: 0
    phi_pi_bar: 14709/10000
     phi_y_bar: 3907/20000
          g_23: 13671939888128922^(1/2)/134217728
    parameters: [1×0 sym]
    conditions: symtrue

However, equation 6 is not close to being satisfied with that solution

vpa(subs(equ(6),[phi_bar theta_bar],[sol.phi_bar sol.theta_bar]))

ans = 

So the equations as written are not consistent.

It sounds like the goal might be an approximate solution? In which case it may be possible to set up some or all of equ as inequalities, like

abs(lhs(equ(1))) < .001

No idea if that's really going to help get to an answer.

[Kylekk]

Thank you Paul! You are right about question, it is about finding an "approximate" solution, but because i know there will be infinitely many, i want them to be parametrized by free parameters. Think about

syms x y z
equ=[x+y+z^2==1, z>0];
sol=solve(equ,[x,y,z],'Real',true,'ReturnConditions', true );
sol.x

ans = 

sol.y

ans = 

u

sol.z

ans = 

v

The (ultimate) goal is to find the upper/lower values of the variables allowed by the equations and inequalities (since some of the vars are free, some are fixed). In the example given, when I rounded it up to 4 digits, they should yield a continuum of solutions to phi_bar and theta_bar, and all the equations/inequalities will still hold. However, as in the example i commented under Torsten's answer, solve cannot seem to tell x==sqrt(2) and x==1.4142 has the same value (even if i set digits(4))

syms x;
equ=[x^2-2==0, x==1.4142];
equ=vpa(equ,4);
solve(equ,[x],'Real',true,'ReturnConditions', true)

ans = struct with fields:
             x: [0×1 sym]
    parameters: [1×0 sym]
    conditions: [0×1 sym]

I did try to used inequalities like

abs(lhs(equ(1))) < .001

ans = 

but then the returned solution is spurious and more messed up than before i solve them, and that is not exactly what i want to have.