Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Plotting a function -
  • Next by Date: RE: Area of ellipse between major axis and ray through focus, given angle
  • Previous by thread: Re: Word permutations - frustrated by lists.
  • Next by thread: Integration (change of variable)