Solving the cubic with Vieta's method
- To: mathgroup at smc.vnet.net
- Subject: [mg71674] Solving the cubic with Vieta's method
- From: "dimitris" <dimmechan at yahoo.com>
- Date: Sun, 26 Nov 2006 03:49:01 -0500 (EST)
Let me consider the Vieta's method for solving the cubic equation.
Off[General::spell1]
$Version
"5.2 for Microsoft Windows (June 20, 2005)"
Here is the standard reduced form of the qubic equation (see
www.mathworld.com)
eq = -q + p*y + y^3 == 0;
Here is the solution provided by Mathematica along with verification.
Solve[eq, y]
Simplify[eq /. %]
{{y -> -(((2/3)^(1/3)*p)/(9*q + Sqrt[3]*Sqrt[4*p^3 + 27*q^2])^(1/3)) +
(9*q + Sqrt[3]*Sqrt[4*p^3 + 27*q^2])^(1/3)/
(2^(1/3)*3^(2/3))}, {y -> ((1 +
I*Sqrt[3])*p)/(2^(2/3)*3^(1/3)*(9*q + Sqrt[3]*Sqrt[4*p^3 +
27*q^2])^(1/3)) -
((1 - I*Sqrt[3])*(9*q + Sqrt[3]*Sqrt[4*p^3 +
27*q^2])^(1/3))/(2*2^(1/3)*3^(2/3))},
{y -> ((1 - I*Sqrt[3])*p)/(2^(2/3)*3^(1/3)*(9*q + Sqrt[3]*Sqrt[4*p^3
+ 27*q^2])^(1/3)) -
((1 + I*Sqrt[3])*(9*q + Sqrt[3]*Sqrt[4*p^3 +
27*q^2])^(1/3))/(2*2^(1/3)*3^(2/3))}}
{True, True, True}
The following inputs follow Vieta's method for solving the cubic (that
is, we will obtain Cardan's formulae for the roots of a reduced cubic
equation).
eq /. y -> w - p/(3*w)
Expand /@ %
(Expand[#1*w^3] & ) /@ %
eq2 = % /. w^(n_) -> v^(n/3)
Expand[v /. FullSimplify[Solve[eq2, v]]]
% /. (a_)*Sqrt[(b_) + (c_)] :> (d*a)*Sqrt[b/d^2 + c/d^2]
% /. d -> 18
({v -> #1} & ) /@ %
solv = % /. Thread[{q/2, q^2/4, p^3/27} -> {R, R^2, Q^3}]
(Solve[w^3 == v /. solv[[#1]], w] & ) /@ Range[2]
solw = % //. {(R + Sqrt[Q^3 + R^2])^(1/3) -> A, (R - Sqrt[Q^3 +
R^2])^(1/3) -> B, -(-1)^(1/3) -> Ï?^2, (-1)^(2/3) -> Ï?}
w /. solw
Outer[Times, %[[1]], %[[2]]]
% //. {A*B -> -Q^3, Ï?^(n_.) -> ((-1)^(2/3))^n}
Position[%, -Q^3]
({solw[[1,#1[[1]]]], solw[[2,#1[[2]]]]} & ) /@ %
solpairs = w /. %
ysols = Apply[Plus, solpairs, {1}]
ysols /. {A -> (R + Sqrt[Q^3 + R^2])^(1/3), B -> (R - Sqrt[Q^3 +
R^2])^(1/3), Ï?^(n_.) -> ((-1)^(2/3))^n}
ysolsrep = % /. {R^(n_.) -> (q/2)^n, Q^(n_.) -> (p/3)^n}
The obtained roots by Cardan's formulae are then y1=A+B, y2=Ï? A+Ï?^2
B, y3=Ï?^2 A+Ï? B}.
Below we give them as they are obtained by above procedure.
ysolsrep
{{y -> (q/2 - Sqrt[p^3/27 + q^2/4])^(1/3) + (q/2 + Sqrt[p^3/27 +
q^2/4])^(1/3)},
{y -> (-(-1)^(1/3))*(q/2 - Sqrt[p^3/27 + q^2/4])^(1/3) +
(-1)^(2/3)*(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3)},
{y -> (-1)^(2/3)*(q/2 - Sqrt[p^3/27 + q^2/4])^(1/3) - (-1)^(1/3)*(q/2
+ Sqrt[p^3/27 + q^2/4])^(1/3)}}
( see
e.g.
http://historical.library.cornell.edu/cgi-bin/cul.math/docviewer?did=01460001&seq=&view=50&frames=0&pagenum=32
http://eqworld.ipmnet.ru/en/solutions/ae/ae0103.pdf )
Substitution of numerical values for {p,q} gives
N[ysolsrep /. {p -> 3, q -> 2}]
{{y -> 1.7142198249542058 + 0.6455631567415321*I}, {y ->
-0.29803581899166054 + 0.5162131809689577*I}, {y -> -1.4161840059625452
- 1.1617763377104895*I}}
which obiously is incorrect.
In contrast
N[Solve[eq, y] /. {p -> 3, q -> 2}]
{{y -> 0.5960716379833217}, {y -> -0.2980358189916608 +
1.8073394944520222*I}, {y -> -0.2980358189916608 -
1.8073394944520222*I}}
So (I believe!) that the choice of branch cuts for qubics within
Mathematica are not in accordance
with Cardan's formula (otherwise I miss something!). This is not a
problem of course but I just mention it.
If we want to use Cardan's formulae for the roots of the cubic within
Mathematica the roots must be given by
ysolsrecor = ({y -> #1} & ) /@ {(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3) -
p/(3*(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3)),
(-1)^(2/3)*(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3) -
p/(3*((-1)^(2/3)*(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3))),
(-(-1)^(1/3))*(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3) -
p/(3*((-(-1)^(1/3))*(q/2 + Sqrt[p^3/27 + q^2/4])^(1/3)))};
Indeed now,
N[ysolsrecor /. {p -> 3, q -> 2}]
{{y -> 0.5960716379833215}, {y -> -0.2980358189916603 +
1.807339494452022*I}, {y -> -0.29803581899166104 -
1.8073394944520218*I}}
and
FullSimplify[eq /. ysolsrecor]
{True, True, True}
Best Regards
Dimitris