# Math Help - Cubic equations and Matlab

1. ## Cubic equations and Matlab

Anyone have any M-Files for solving Cubic Equations on Matlab? Been searching with no luck.

2. Originally Posted by drhojo
Anyone have any M-Files for solving Cubic Equations on Matlab? Been searching with no luck.
Please be more specific. Do you want to find a solution, do you want
all real solutions, do you want all the solutions including the complex?

Also how is your cubic represented, as an array of coefficients, or
something different?

How have you been searching, third hit on Google for search string "cubic solver matlab"

Code:
function [x,nroot]=cubicfcn(a, b, c, d)
% ----------------------------------------------------------------------
% Usage:       [x,nroot]=cubicfcn(a, b, c, d)
%
% Solve a cubic equation where a, b, c, and d are real.
%   a*x^3 + b*x^2 + c*x + d = 0
%
% Public Variables
%   a, b, c, d  ... coefficients (input)
%   x           ... three (generally) complex solutions (output)
%   nroot       ... number of roots (output)
%
% Instructor: Nam Sun Wang
% ----------------------------------------------------------------------

% Local Variables:
%   y1, y2, y3  ... three transformed solutions
%
% Formula used are given in Tuma, "Engineering Mathematics Handbook", p7
%   (McGraw Hill, 1978).
%   Step 0: If a is 0. use the quadratic formula to avoid dividing by 0.
%   Step 1: Calculate p and q
%           p = ( 3*c/a - (b/a)**2 ) / 3
%           q = ( 2*(b/a)**3 - 9*b*c/a/a + 27*d/a ) / 27
%   Step 2: Calculate discriminant D
%           D = (p/3)**3 + (q/2)**2
%   Step 3: Depending on the sign of D, we follow different strategy.
%           If D<0, thre distinct real roots.
%           If D=0, three real roots of which at least two are equal.
%           If D>0, one real and two complex roots.
%   Step 3a: For D>0 and D=0,
%           Calculate u and v
%           u = cubic_root(-q/2 + sqrt(D))
%           v = cubic_root(-q/2 - sqrt(D))
%           Find the three transformed roots
%           y1 = u + v
%           y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
%           y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
%   Step 3b Alternately, for D<0, a trigonometric formulation is more convenient
%           y1 =  2 * sqrt(|p|/3) * cos(phi/3)
%           y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3)
%           y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3)
%           where phi = acos(-q/2/sqrt(|p|**3/27))
%                 pi  = 3.141592654...
%   Step 4  Finally, find the three roots
%           x = y - b/a/3
% ----------------------------------------------------------------------

%     pi = 3.141592654

% Step 0: If a is 0 use the quadratic formula. -------------------------
if (a == 0.)
[x, nroot] = quadfcn(b, c, d);
return
end

% Cubic equation with 3 roots
nroot = 3;

% Step 1: Calculate p and q --------------------------------------------
p  = c/a - b*b/a/a/3. ;
q  = (2.*b*b*b/a/a/a - 9.*b*c/a/a + 27.*d/a) / 27. ;

% Step 2: Calculate DD (discriminant) ----------------------------------
DD = p*p*p/27. + q*q/4. ;

% Step 3: Branch to different algorithms based on DD -------------------
if (DD < 0.)
%       Step 3b:
%       3 real unequal roots -- use the trigonometric formulation
phi = acos(-q/2./sqrt(abs(p*p*p)/27.));
temp1=2.*sqrt(abs(p)/3.);
y1 =  temp1*cos(phi/3.);
y2 = -temp1*cos((phi+pi)/3.);
y3 = -temp1*cos((phi-pi)/3.);
else
%       Step 3a:
%       1 real root & 2 conjugate complex roots OR 3 real roots (some are equal)
temp1 = -q/2. + sqrt(DD);
temp2 = -q/2. - sqrt(DD);
u = abs(temp1)^(1./3.);
v = abs(temp2)^(1./3.);
if (temp1 < 0.) u=-u; end
if (temp2 < 0.) v=-v; end
y1  = u + v;
y2r = -(u+v)/2.;
y2i =  (u-v)*sqrt(3.)/2.;
end

% Step 4: Final transformation -----------------------------------------
temp1 = b/a/3.;
y1 = y1-temp1;
if (DD < 0.)
y2 = y2-temp1;
y3 = y3-temp1;
else
y2r=y2r-temp1;
end

if (DD < 0.)
x(1) = y1;
x(2) = y2;
x(3) = y3;
elseif (DD == 0.)
x(1) =  y1;
x(2) = y2r;
x(3) = y2r;
else
x(1) = y1;
x(2) = y2r + y2i*i;
x(3) = y2r - y2i*i;
end
RonL

3. Originally Posted by drhojo
Anyone have any M-Files for solving Cubic Equations on Matlab? Been searching with no luck.
One reason you wont find much is that this is done by the Matlab function
roots. See attachment (soory I seem to have cropped off the last root in the
example).

RonL

4. Sorry for not being specific enough, anyway, the example you posted (and which I couldn't find because I was using a different string) is exactly what I was looking for, thanks.