Incorrect Calculation when using symunit in Symbolic Math Toolbox

clc,clear
syms x a
u = symunit;
a = 1.01325*u.bar
b = -0.049576382*u.bar*u.m^3/u.mol
c = .000003657553277*u.bar*u.m^6/(u.mol)^2
d = -.00000000015676624*u.bar*u.m^9/u.mol^3
eqn = a*x^3+b*x^2+c*x+b == 0
solve(eqn,x,"MaxDegree",3)
vpa(solve(eqn,x,'MaxDegr

Output: 0.016

clc,clear
syms x
u = symunit;
a = 1.01325
b = -0.049576382
c = 3.657553277*10^-6
d = -1.5676624*10^-10
eqn = a*x^3+b*x^2+c*x+b == 0
solve(eqn,x,"MaxDegree",3)
vpa(solve(eqn,x,'MaxDegree',3),10)

Output: 0.382

a = 1.01325
b = -0.049576382
c = 3.657553277*10^-6
d = -1.5676624*10^-10
x = [a b c d];
vpa(roots(x),10)

Output: 0.048

u = symunit;
a = 1.01325*u.bar
b = -0.049576382*u.bar*u.m^3/u.mol
c = .000003657553277*u.bar*u.m^6/(u.mol)^2
d = -.00000000015676624*u.bar*u.m^9/u.mol^3
x = [a b c d];
vpa(roots(x),10)

Output: root(z^3 - ((sym("223272175501577375")*z^2)/sym("4563272322433155072"))*(symunit("m")^3/symunit("mol")) + ((sym("539759593894213625")*z)/sym("149529307461489625399296"))*(symunit("m")^6/symunit("mol")^2) - sym("252691673573204875/1633258782299364015028043776")*(symunit("m")^9/symunit("mol")^3), z, 1)

Wolfram Alpha Answer: 0.048

[Dyuman Joshi]

Where exactly do you get the output 0.016 from the first snippet of code?

Do you mean the sigma2 value? If yes, then you are misunderstanding the solution. The solutions are the elements of the 3x1 vector, where the terms have been grouped. And the information below the solutions show, which terms have been grouped into the respective variables.

Also, there's a mistake in the 2nd piece of code, you have written a*x^3 + b*x^2 + c*x + b == 0, it should be d instead of b. Correcting the mistake leads to the correct output, see below -

clc,clear
syms x a
u = symunit;
a = 1.01325*u.bar;
b = -0.049576382*u.bar*u.m^3/u.mol;
c = .000003657553277*u.bar*u.m^6/(u.mol)^2;
d = -.00000000015676624*u.bar*u.m^9/u.mol^3;
eqn = a*x^3+b*x^2+c*x+b == 0

eqn = 

vpa(solve(eqn,x,"MaxDegree",3), 10)

ans = 

clc,clear
syms x
u = symunit;
a = 1.01325;
b = -0.049576382;
c = 3.657553277*10^-6;
d = -1.5676624*10^-10;
%                     v
eqn = a*x^3+b*x^2+c*x+d == 0;
solve(eqn,x,"MaxDegree",3);
vpa(solve(eqn,x,'MaxDegree',3),10)

ans = 

clear
a = 1.01325;
b = -0.049576382;
c = 3.657553277*10^-6;
d = -1.5676624*10^-10;
x = [a b c d];
vpa(roots(x),10)

ans = 

u = symunit;
a = 1.01325*u.bar;
b = -0.049576382*u.bar*u.m^3/u.mol;
c = .000003657553277*u.bar*u.m^6/(u.mol)^2;
d = -.00000000015676624*u.bar*u.m^9/u.mol^3;
x = [a b c d];
vpa(roots(x),10)

ans = 

[Breyonna]

First of all, thank you SO much. I really appreciate the solution. Do the units negate the calculation? Is it possible to get the numerical output from the code without units into the output with the code with units.

I am basically punching in numbers without understanding the output at this point.

clc,clear
syms x
u = symunit;
a = 1.01325*u.bar;
b = -0.049576382*u.bar*u.m^3/u.mol;
c = .000003657553277*u.bar*u.m^6/(u.mol)^2;
d = -.00000000015676624*u.bar*u.m^9/u.mol^3;
%                     v
eqn = a*x^3+b*x^2+c*x+d == 0;
solve(eqn,x,"MaxDegree",3);
vpa(solve(eqn,x,'MaxDegree',3),10)

I am trying to figure out how to output 0.048 m^3/mol for the above code.

I don't know how to ask it to do that.

Answer 1

[Breyonna]

Got it......somehow. Thanks again.

clc,clear
syms x
u = symunit;
a = 1.01325*u.bar;
b = -0.049576382*u.bar*u.m^3/u.mol;
c = .000003657553277*u.bar*u.m^6/(u.mol)^2;
d = -.00000000015676624*u.bar*u.m^9/u.mol^3;
%                     v
eqn = a*x^3+b*x^2+c*x+d == 0;
vpa(simplify(solve(eqn, x, 'maxdegree', 3)))

ans =