# Symbolic solve with user-specified precision

2 views (last 30 days)
Kylekk on 27 Jul 2022
Edited: Kylekk on 27 Jul 2022
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 on 27 Jul 2022
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.

John D'Errico on 27 Jul 2022
Edited: John D'Errico on 27 Jul 2022
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 on 27 Jul 2022
Edited: Kylekk on 27 Jul 2022
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.

Paul on 27 Jul 2022
Edited: Paul on 27 Jul 2022
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.
##### 1 CommentShowHide None
Kylekk on 27 Jul 2022
Edited: Kylekk on 27 Jul 2022
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.

R2022a

### Community Treasure Hunt

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

Start Hunting!