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