Analytic solution of quatric polynomial does not add up

88 次查看（过去 30 天）
thejinn1314
thejinn1314 2025-8-29，17:09
评论： Steven Lord 2025-8-31，13:46
Dear friends,
I currently ran into a math-issue when trying to get the analytical roots of a polynomial of order 4.
What I want to solve is:
The solution shall run on a microcontroller with dynamically changing parameters [a,d,e]. All parameters are real. I am searching for the real results, imaginary solutions are to be prevented. When applying Matlab symbolic, The parameters [a,d,e] are simplified and represent longer sub-terms. In accordance to the article Wiki: quartic function, I will receive quite a long solution, that is futile to reproduce here. I used the following short code to generate it:
syms x a d e
y= a*x^4+d*x+e
sol=solve(y==0,x,'MaxDegree',4,'ReturnConditions',true)
To check the solution, I tested it on an easier sub-problem: When d=0 is chosen follows:
When chosing a>0, e<0 this clearly yields some real results, besides the imaginary ones. Now the problem: the above given long solution cannot reproduce this result. Even if I force assumptions that should solve it, follows something that does not look equal:
syms x a d e
assume(a>=0)
assume(d>=0)
assume(e<=0)
y= a*x^4+d*x+e
sol=solve(y==0,x,'MaxDegree',4,'ReturnConditions',true)
testresult=simplify(subs(sol.x, d,0))
latex(testresult)
However, if I change the equation before solving the quartic function, the solution of it looks as on paper:
y= a*x^4+e
assume(a>=0)
assume(d>=0)
assume(e<=0)
sol=solve(y==0,x,'MaxDegree',4,'ReturnConditions',true)
testresult2=sol.x
latex(testresult2)
Numerical testing via inserting numbers in the two solutions according to the assumptions a>0, e<0 confirms the difference of the results. Chosing a=1 and e=-2 yields:
1) the general solution of the quartic polynomial (inserting d=0 after solving the quartic polynomial):
vpares1=vpa(subs(testresult,{a,e},{1,-2}))
vpares1 =
0
0
0
0
2) the solution of the simplified quartic polynomial (inserting d=0 before solving):
vpares2=vpa(subs(testresult2,{a,e},{1,-2}))
vpares2 =
1.1892071150027210667174999705605
-1.1892071150027210667174999705605
-1.1892071150027210667174999705605i
1.1892071150027210667174999705605i
Visibly, the trivial solution is not permissible given the solution of case 2). Interestingly, if the requirement is switched from {a>0,e<0} to {a<0, e>0} the solutions match again.
Similar difficulties have shown, when I hand-coded the general solution from the mentioned Wikipedia article. As I need the solution for , I now have trust-issues with the long solution. Therefore, I ask for your help or comment, why the "general solution" does not end in the same result as the simplified problem.
Thank you in advance for your help!
  1 个评论
Steven Lord
Steven Lord 2025-8-31，13:46
The solution shall run on a microcontroller with dynamically changing parameters [a,d,e]. All parameters are real. I am searching for the real results, imaginary solutions are to be prevented. When applying Matlab symbolic,
Since you're planning to target a microcontroller, you probably don't want to use Symbolic Math Toolbox. While you could solve the system in MATLAB to generate a symbolic solution then generate C code from that symbolic solution with the ccode function, I'd probably use a different approach, one that would work even if you later wanted to modify your code to solve a quintic equation or higher.
The roots function in MATLAB supports the "C/C++ Code Generation" extended capability, so you could generate code from it for use on your microcontroller. Yes, the notes on that documentation page for that extended capability state that the values returned are always returned as a complex array. But your code could check which root(s) have a small magnitude imag part (that function also supports that extended capability) and select the real part of that root or roots.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2025-8-30，20:09
Ran in:
Hi thejinn1314,
simplify-ing the general solution with IgnoreAnalyticConstraints might be useful, though perhaps could cause some unease
syms x
syms a d e real
y = a*x^4+d*x+e;
sol = solve(y==0,x,'MaxDegree',4,'ReturnConditions',true);
sol.conditions
ans = 
[num,den] = numden(sol.x);
Each solution has the same denominator
isAlways(den == den(1))
ans = 4×1 logical array
1 1 1 1
so we only need to work with the first entry of den. Simplify it and num with IgnoreAnalyticConstraints = true
den = simplify(den(1),'IgnoreAnalyticConstraints',true);
num = simplify(num,'IgnoreAnalyticConstraints',true)
num = 
Form the solution vector
sol = num/den;
Verify
simplify(subs(y,x,sol),100)
ans = 
ans = 
Solution for the special case of d = 0
sold0a = simplify((subs(sol,d,0)),100)
sold0a = 
I tried to simplify that further, but didn't get very far.
Solutin to the original equation setting d = 0 from the outset.
sold0b = solve(a*x^4 + e);
sold0b = sold0b([2,4,3,1])
sold0b = 
Compare some cases with a > 0 and a < 0
vpa([subs(sold0a,[a,e],[2,2]), subs(sold0b,[a,e],[2,2])])
ans = 
vpa([subs(sold0a,[a,e],[-2,2]), subs(sold0b,[a,e],[-2,2])])
ans = 
  2 个评论
thejinn1314
thejinn1314 2025-8-31，12:16
Hey Paul,
Thank you for your answer. I especially like the idea to treat nominator and denominator individually.
Indeed I can reproduce the result for the three cases. First for the generalized solution after inserting d=0:
syms x
syms a d e
y = a*x^4+d*x+e;
sol = solve(y==0,x,'MaxDegree',4,'ReturnConditions',true);
sol.conditions
[num,den] = numden(sol.x)
den = simplify(den(1),'IgnoreAnalyticConstraints',true)
num = simplify(num,'IgnoreAnalyticConstraints',true)
sol = num/den
sold0a = simplify((subs(sol,d,0)),100)
res11=vpa([subs(sold0a,[a,e],[1,-2])])
res12=vpa([subs(sold0a,[a,e],[-1,2])])
outputs for {a>0,e<0}, here [a=1,e=-2]:
res11 =
- 1.1892071150027210667174999705605 - 9.1835496157991211560057541970488e-41i
1.8367099231598242312011508394098e-40 - 1.1892071150027210667174999705605i
- 1.8367099231598242312011508394098e-40 + 1.1892071150027210667174999705605i
1.1892071150027210667174999705605 + 9.1835496157991211560057541970488e-41i
and for {a<0,e>0}, here [a=-1,e=2]:
res12 =
- 1.1892071150027210667174999705605 - 6.4284847310593848092040279379342e-40i
- 6.4284847310593848092040279379342e-40 + 1.1892071150027210667174999705605i
6.4284847310593848092040279379342e-40 - 1.1892071150027210667174999705605i
1.1892071150027210667174999705605 + 6.4284847310593848092040279379342e-40i
Which fits to the solution of the reduced problem:
y= a*x^4+e
sol=solve(y==0,x,'MaxDegree',4,'ReturnConditions',true)
result2=sol.x
res2=vpa(subs(result2,{a,e},{1,-2}))
outputs:
res2 =
1.1892071150027210667174999705605
-1.1892071150027210667174999705605
-1.1892071150027210667174999705605i
1.1892071150027210667174999705605i
This was not possible with my solution given in the question above. Therefore, it fulfills my requirements given by this test and for the question.
However, as you mentioned too, the condition "IgnoreAnalyticConstraints" could also lead to unforseen effects. What confuses me is, that if I turn it off, the result is not applicable again:
syms x
syms a d e
y = a*x^4+d*x+e;
sol = solve(y==0,x,'MaxDegree',4,'ReturnConditions',true);
sol.conditions
[num,den] = numden(sol.x)
den = simplify(den(1))
num = simplify(num)
sol = num/den
sold0a = simplify((subs(sol,d,0)),100)
res11=vpa([subs(sold0a,[a,e],[1,-2])])
res12=vpa([subs(sold0a,[a,e],[-1,2])])
outputs for {a>0,e<0} and {a<0, e>0} only an error for division by zero:
Error using symengine
Division by zero.
Error in sym/subs>mupadsubs (line 160)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);
Error in sym/subs (line 145)
G = mupadsubs(F,X,Y);
Error in script250829 (line 82)
res11=vpa([subs(sold0a,[a,e],[1,-2])])
only {a>0,e>0} will lead to a result.
So whatever analytic constraint is ignored, it is doing something good.
This of course leaves some remaining uncertainty, but is a sufficient answer for now. Thank you a lot for helping out!
Paul
Paul 2025-8-31，13:15
Yes, if you turn off IgnoreAnalyticConstraints, then the case with d = 0 and e and a with opposite signs actually leads to a 0/0 situation. I tried to deal with that a few different ways, but the best I could do was as above. I agree that there is some remaining uncertainty; maybe you can continue to work the problem and get the solution on more solid ground.

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2025-8-29，18:09
编辑：John D'Errico 2025-8-30，0:01
Ran in:
One of the things you need to learn about symbolic algebra, is two expressions may not look the same, but that does not mean they are not the same.
syms x a d e
assume(a>=0)
assume(d>=0)
assume(e<=0)
y= a*x^4+d*x+e
y = 
x_d = solve(y==0,x,'MaxDegree',4,'ReturnConditions',true)
x_d = struct with fields:
x: [4×1 sym] parameters: [1×0 sym] conditions: [4×1 sym]
Now substitute d==0,
x_d_0 = subs(x_d,d,0);
x_d_0 = simplify(x_d_0.x)
x_d_0 = 
Now, substitute d==0, into y. And then solve.
x_0 = solve(subs(y,d,0)==0,'MaxDegree',4,'ReturnConditions',true);
x_0 = x_0.x
x_0 = 
Now, they don't look the same. But this does not mean anything. First, the two sets of solutions are not in the same order. So I'll shuffle one of the two sets so they correspond.
simplify(x_0 - x_d_0([4,1,3,2]))
ans = 
And, yes, I know they look different. But that does not mean they are different. To make it a bit more clear, I'll insert some numbers in for a and e.
vpa(subs(x_0 - x_d_0([4,1,3,2]),[a,e],[1,2]))
ans = 
And there we see the two vectors are the same, with only 40 digit floating point slop remaining. Do you see? These are not different solutions. They just appear that way. If you wear a mask at a halloween party, are you a different person? Of course not.
The ultimate problem is not that they are different, but that complex numbers inside radicals can often appear wildly differnt, yet still be identical. The quartic solution is a complete mess, with layers of radicals. And the result is it produces something that only looks strange and different. It still represents the same numbers in the end. Both roads still led to Rome in the end.
A second thing you should learn is if you can simplify a problem first, that often makes the symbolic solution much simpler. You can see exactly that in this case, where when you substitute d==0 into y, it made the problem reduce nicely, so that solve did not need to involve the full blown quartic solution. We can easily give other examples, where no solution at all can be found to a fully general problem, but if you simplify the probem, the solution is now trivial.
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2025-8-30，20:22
编辑：John D'Errico 2025-8-30，21:43
Ran in:
As I try to follow a moving target...
I'm sorry. I substituted the wrong sign for e in my example. But that is not relevant. First of all, since you care only about real roots, I would ask why not just use a numerical solver? For example, fzero would trivially solve for the root.
format long g
yfun = @(x,a,e) a*x.^4 + e;
[xdoub,ydoub,exitflag] = fzero(@(x) yfun(x,-1,2),[eps,1000])
xdoub =
1.18920711500272
ydoub =
2.22044604925031e-16
exitflag =
1
Even better, you could just as simply have used nthroot, instead of getting bogged down in a mess of stuff..
e = 2;a = -1;
nthroot(-e/a,4)
ans =
1.18920711500272
and be done with it.
But most likely, this is not your real problem. It never is. Or you are going through all of this to try to understand what is happening. But rather than jumping to the conclusion that the code is unreliable, it is the case here that you don't understand how to make the code work properly.
Looking at what you show, it seems you want to solve the case where a=-1, and e=2. Note that this might throw a curve at the solver, which does not really understand this is the same problem as solving for a=1, and e=-2. The problem is, that factor of -1 on both coefficient will introduce complex numbers in different places. But even so, you just need to look CAREFULLY at what you are doing, and at the results.
syms x real positive
syms y a d e
assume(a>=0)
assume(d>=0)
assume(e<=0)
Note that often, assumptions are not used terribly intelligently by solve. And since you want to see real positive solutions for x, I posed that when I defined x. This may help things, at least in one of the cases.
y = a*x^4+d*x+e;
sol=solve(y==0,x,'MaxDegree',4,'ReturnConditions',true)
sol = struct with fields:
x: [4×1 sym] parameters: [1×0 sym] conditions: [4×1 sym]
result1 = simplify(subs(sol.x, d,0))
result1 = 
sol2 = solve(subs(y,d,0)==0,x,'MaxDegree',4,'ReturnConditions',true)
sol2 = struct with fields:
x: (-e)^(1/4)/a^(1/4) parameters: [1×0 sym] conditions: a ~= 0 & e ~= 0
result2 = simplify(sol2.x)
result2 = 
In result1, MATLAB currently does not know which solution is real and positive, because it still does not know anything about a and e. In result2, because things were much simpler, it was able to decide there is only one root you will care about.
subs(result2,[a,e],[-1,2])
ans = 
result1vpa = vpa(subs(result1,[a,e],[-1,2]))
result1vpa = 
itol = 1e-30;
result1vpa(real(result1vpa) > 0 & abs(imag(result1vpa)) < itol)
ans = 
And of course, exactly ONE of the roots in result1 is both real and positive. I wonder if it might be the 4th root of 2?
So, yes, the solution is indeed perfectly reliable! You are doing a lot of stuff, I think not really understanding what is happening, but then immediately jumping to the conclusion there must be a bug in the software, that it is not reliable. Do you see what I am saying? When I get a result that confuses me, the very first thing I do is CAREFULLY go back through what I did, and I assume I made a mistake. Then I look for my mistake, long before I will ever assume the code is unreliable. (In this case, for exmple, it seems it helps to tell MATLAB that x must be real and positive.) In the end though, code is no better than what you pass into the tools, and how you use them. The pertinent factor here is how you use the code, and not if the code is reliable.
thejinn1314
thejinn1314 2025-8-31，11:15
First of all, you got all emotional, when interpreting my text as "Matlab code does not work" although this was never my implication nor intention. If you would have read my answer "carefully" - yeah caps-lock is really professional - you would have seen, that it was always my goal to understand the results and how to apply them to my application and not to blame Matlab for wrong code. That is what I mean by "the solution is unreliable [for my application]" and its following line in the comment "please point out my mistake". So please excuse me for causing this misconception about my intention. I would never first assume that a nearly 500 year old solution is faulty. What I want is guidance on understanding the application of the result or a different formulation that fits my application.
However, given that I formulated politely and on point, I think your answer is unnecessary sassy and not focused on the problem. Can numerical solvers do the job? For sure, but not for a real-time application in my microcontroller (as mentioned from the beginning). Also: it does not matter for the problem stated here. You then focus completly about the real solution, although the problem I point out is the discrepancy for the cases of {a>0, e<0} and {a<0, e>0} and then only show the case {a<0,e>0} that I have shown before as "working", too. You say that "You are doing a lot of stuff, I think not really understanding what is happening [...]" - well thank you for pointing this out, maybe that is why I ask the question here?(!). Then please enlighten me, but about the actual problem. As you struggle to see my actual problem, here the summary of it again:
Why does for the aforementioned result derived from general solution after inserting d=0:
for the case {a>0, e<0}, here e.g. [a=1,e=-2] produce x=0:
vpares1=vpa(subs(result1,{a,e},{3,-4}))
vpares1 =
0
0
0
0
instead of the correct
vpares2 =
1.1892071150027210667174999705605
-1.1892071150027210667174999705605
-1.1892071150027210667174999705605i
1.1892071150027210667174999705605i
derived from the simplified problem (case 2) or the case {d=0, a<0, e>0}, e.g. [a=1,e=-2] as shown by my first answer to your post and by you in your answer?
This discrepancy probalby lies in an condition or assumption to the result that I struggle to see. However, obviously I need to obey it get reliable results for cases of d that are different from zero.
However, please refrain from further wasting your time on this problem, if my way of asking is upsetting you and preventing you from writing polite and supporting answers.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by