cubic equation in J
Author Message cubic equation in J

The following script gives a closed form for the solution of a cubic equation
with real coefficients. This is the method published by Vieta in 1615. The
basis of the solution is due originally to Scipio del Ferro, who died in
1526, but who never published his solution. It was rediscovered by Tartaglia
in 1535, who also kept the method to himself. He divulged it under an oath of
secrecy to Cardano, who published it in 1545, crediting Tartaglia. After this
Cardano and Tartaglia debated bitterly and exchanged insults. Cardano's
method didn't handle the case where three real roots were obtained as the
sums of complex numbers. Bombelli in 1572 developed a method for working
freely with complex numbers, and so was able to solve the cubic equation in
greater generality. Vieta converted the solution to the trigonometric one
given below.

NB. Cubic equation: see Numerical Recipes in C, Press et al., Section 5.5, p.
157,
NB. or Handbook of Mathematical Functions by Abramowitz and Stegun, Section
3.8.2, p. 17.
NB. For quartic equation, see Abramowitz and Stegun, op. cit., Section 3.8.3,
pp 17-18.

NB. Uses arccos and cos from 'trig' utility script

NB. 'cubic' applies to a 4-item list of real number coefficients.
NB. The leading item is the constant term, then the linear term, etc.
NB. If last item is not 1, apply 'N' to argument before applying 'cubic'.

NB. For example, cubic _30 31 _10 1 is 2 5 3
NB. and cubic N _60 62 _20 2 is also 2 5 3
NB. and Cubic _60 62 _20 2 is also 2 5 3.

NB. Note that 2+5+3 is 10 (-a2, the quadratic term),
NB. (2*5)+(2*3)+(5*3) is 31 (a1, the linear term),
NB. and 2*5*3 is 30 (-a0, the constant term).
NB. This gives a simple way of checking the result.

N =: % {: NB. divide by tail to normalize coefficients with 1 as
high-order term
a0 =: {. NB. constant term
a1 =: 1&{ NB. linear term
a2 =: 2&{ NB. quadratic term

cubic =: (_2: * ([: %: Q) * [: cos (0p1 2p1 4p1"_ + theta) % 3:) - a2 % 3:

A simpler, more general, and much more efficient closed form is also
available in J, namely the polynomial primitive, which solves polynomials of
arbitrary order and with complex coefficients.

p.

Eugene McDonnell
Palo Alto

Wed, 04 Aug 1999 03:00:00 GMT

 Page 1 of 1 [ 1 post ]

Relevant Pages