Solving a 4th order equation in MATLAB

Question

I'm trying to solve a general 4th order equation in one variable in MATLAB. The equation is

$$a x^4 + b x^3 + c x^2 + d x + e = 0$$

I tried this code:

clear;
close all;
clc;

syms a b c d e
syms x

eq1 = a*x^4 + b*x^3 + c*x^2 + d*x + e==0;
x_sol = solve(eq1,x)

But I don't get explicit results. This is what I get:

x_sol=RootOf(a*z^4 + b*z^3 + c*z^2 + d*z + e, z)

I tried in a different way in this code:

clear;
close all;
clc;

syms a b c d e
syms x

coeffs = [a b c d e]; %eq = a*x^4 + b*x^3 + c*x^2 + d*x + e==0;
x_symbolic = roots(coeffs) %solving the equation

a=1; b=0; c=0; d=0; e=-16; %x^4-16=0 (solutions should be [2,-2,2i,-2i]
x_numeric = subs(x_symbolic);
x_numeric = double(x_numeric)

where x_symbolic should give me the symbolic results for the quartic equation. But when I checked the resutls and substituted $[a,b,c,d,e]=[1,0,0,0,-16]$ which means that the equation becomes $x^4-16=0$, I expected to get these 4 results: $[2,-2,2i,-2i]$. but I actually got these four results:

x_numeric =
-1.867e-36 -  1.867e-36i
1.867e-36 -  1.867e-36i
-1.867e-36 +  1.867e-36i
1.867e-36 +  1.867e-36i

Does anybody know why? The goal is to get the correct symbolic solution for the general equation.

Answer 1

Computers: they have limited precision. Once you put floating point numbers in that (ridiculously enormous) root solution, numerical errors get in and screw the result, especially in something like:

https://pastebin.com/YemafsAc (I wanted to put it in here, but the SE doesnt allow me to put more than 30.000 chars in the post).

The solution is correct, but you have errors numerically evaluating the expression.

Is there anything impeding you from doing the following?

syms x
a=1; b=0; c=0; d=0; e=-16; %x^4-16=0 (solutions should be [2,-2,2i,-2i]
coeffs = [a b c d e]; %eq = a*x^4 + b*x^3 + c*x^2 + d*x + e==0;
x_symbolic = roots(coeffs) %solving the equation

Numerically solving these equations is also ridiculously cheap in a computer.

for ii=1:1000
   coeffs=rand(1,5)*20;
   tic;
   x_symbolic = roots(coeffs); %solving the equation
   t(ii)=toc;
end
mean(t)

ans
>> 4.5216e-05
% or a total of 0.045 in my PC.

Compared to the symbolic option:

syms a b c d e
syms x

coeffs = [a b c d e]; %eq = a*x^4 + b*x^3 + c*x^2 + d*x + e==0;
tic
x_symbolic = roots(coeffs); %solving the equation
toc


>> Elapsed time is 1.350836 seconds.

Note that in my PC, you need approximately 100.000 numerical solutions to get to the same time consumption as the symbolic equation solution, and that is only without evaluating the result after. I'll leave as an exercise to the reader the evaluation of the numerical accuracy of roots, but I'll say: its quite accurate!