Google Answers Logo
View Question
 
Q: for mathtalk-ga only: Need more help with cubic equations ( Answered 5 out of 5 stars,   2 Comments )
Question  
Subject: for mathtalk-ga only: Need more help with cubic equations
Category: Science > Math
Asked by: donkan-ga
List Price: $5.00
Posted: 12 Dec 2004 02:42 PST
Expires: 11 Jan 2005 02:42 PST
Question ID: 441538
mathtalk-ga, I tried to boil down your answer at
<http://answers.google.com/answers/threadview?id=433886> to something
I thought I could use to write a Python script:

x^3 + a x^2 + b x + c = 0 (I assume a,b,c to be real)

M  =  b - a^2/3

N  =  -(2(a^3)/27 - ab/3 + c)

z is a root of z^3 + Mz = N  if and only if
  
z = w - M/(3w) and w^3 is a root of (w^3)^2 - N(w^3) - (M^3)/27 = 0

I will need to compute the 3 cube roots of w^3 (which I can now do),
and then solve for z. There will be 3 identical pairs of z's.
Discarding one member of each pair, what remain will be the 3 roots of
the cubic equation.

Is my understanding correct? If not, please correct it.

Thanks,

donkan

Request for Question Clarification by mathtalk-ga on 12 Dec 2004 07:43 PST
Hi, donkan-ga:

You did leave out the part where I said:

  x is a root of x^3 + a x^2 + b x + c = 0
  
  if and only if
  
  x = z - a/3 and z is a root of z^3 + Mz = N

Maybe for the sake of your Python program, it makes sense to describe
the outer and inner problems as subroutines.

The outer routine takes three "real" numbers a,b,c and (we'll assume)
prints three "complex" roots of x^3 + a x^2 + b x + c = 0.

To do so requires a call to an inner routine, which takes two "real"
numbers M,N and returns three "complex" roots of z^3 + Mz = N.

All the outer routine needs to do is compute M,N (from a,b,c) as
inputs to the inner routine, and then subtract a/3 from the roots
returned by the inner routine to produce the roots for the "original"
problem.

If this is the sort of clarification you need, I'll post as an Answer
detailed comments on how to write the inner routine, which solves z^3
+ Mz = N.

regards, mathtalk-ga

Clarification of Question by donkan-ga on 12 Dec 2004 07:56 PST
mathtalk-ga,

Yes, please.

donkan

Clarification of Question by donkan-ga on 12 Dec 2004 08:37 PST
mathtalk,

You didn't mention w^3. Don't I have to compute w at some point?

donkan
Answer  
Subject: Re: for mathtalk-ga only: Need more help with cubic equations
Answered By: mathtalk-ga on 17 Dec 2004 06:07 PST
Rated:5 out of 5 stars
 
Hi, donkan-ga:

Below is an implementation of the analytical formulas for solving the
cubic equation, using three functions:

  solveCubic (module name)
    solveCubicMonic(a,b,c)
    solveCubicReduced(m,n)
    cubeRoots(x)

Each function returns a triple of complex numbers.  Take a look and
see if you have any questions.  Since Python treats "whitespace" as
significant in some ways, look carefully to see that the posting here
in this textbox hasn't "broken" the code in some way (I'll check it
also).

regards, mathtalk-ga

>>>>>>>>>>> begin solveCubic.py >>>>>>>>>>>>>>>>>>

# -*- coding: UTF-8 -*-

import math, cmath

def solveCubicMonic(a,b,c):
    """Solve a monic real cubic polynomial.
    
    Input three real coefficients a,b,c and return
    three complex roots of Uł + aU˛ + bU + c = 0.
    """
    m = b - (a*a)/3
    n = (b + m + m)*a/9 - c
    x,y,z = solveCubicReduced(m,n)
    x,y,z = x - a/3, y - a/3, z - a/3
    return (x,y,z)

def solveCubicReduced(m,n):
    """Solve a reduced monic real cubic polynomial.

    Input two real coefficients m,n and return
    three complex roots of uł + mu = n.
    """
    if m == 0:
        return cubeRoots(n)
    if n == 0:
        return (0,sqrt(-m + 0j),-sqrt(-m + 0j))
    m = m/3 + 0j
    n = n/2
    w = n + cmath.sqrt(n*n + m*m*m)
    x,y,z = cubeRoots(w)
    return (x - m/x, y - m/y, z - m/z)

def cubeRoots(x):
    """Find three complex cube roots of real or complex x."""
    z = x + 0j
    s = 1
    if z.real < 0:
        z,s = -z,-s
    t = math.atan2(z.imag, z.real)/3
    a = s * (abs(z)**(1./3))
    w = complex(-0.5,math.sqrt(0.75))
    w_ = w.conjugate()
    u = a * complex(math.cos(t),math.sin(t))
    return (u, u*w, u*w_)

Clarification of Answer by mathtalk-ga on 17 Dec 2004 06:13 PST
Okay, I found one mistake.  In the middle of the program:

        return (0,sqrt(-m + 0j),-sqrt(-m + 0j))

should qualify function sqrt with module name cmath:

        return (0,cmath.sqrt(-m + 0j),-cmath.sqrt(-m + 0j))

as it is a couple of lines below in the same function.

regards, mathtalk-ga
donkan-ga rated this answer:5 out of 5 stars and gave an additional tip of: $10.00
Wonderful, mathtalk-ga. You appear to have learned a lot of Python to do this.
I'd like to give you 10 stars to go along with my tip. 

Thank you!!

BTW When I ran what you gave me, IDLE suggested a different coding,   
         "# -*- coding: cp1252 -*-".
You can see your work at <http://www.rcblue.com/Python/solveCubicEquation.py>

donkan

Comments  
Subject: Re: for mathtalk-ga only: Need more help with cubic equations
From: mathtalk-ga on 12 Dec 2004 11:52 PST
 
The interested reader may want to check out these links:

[python -- the home page]
http://www.python.org/

The Python programming language was created at CWI in Amsterdam in
1990, and has become a very popular interpretive platform for rapid
prototyping because of its object-oriented and extensible design. 
"Python" is a humorous homage to Monty Python.  The most recent major
release is Python 2.4 (Nov. 30, 2004), which is provided for Windows
platforms starting with this release as an .msi (Microsoft installer)
file.

[Dive Into Python]
http://diveintopython.org/

This site hosts an online and up-to-date copy of the book Dive Into
Python by Mark Pilgrim.  The text includes recommendations on choosing
a Python implementation.

[ActivePython -- freeware Python]
http://www.activestate.com/Products/ActivePython/

Windows users may be interested in Mark Hammond's PythonWin IDE, which
is bundled with the Windows version available here (and at
SourceForge).

[A Brief Introduction to Python for Scientific Computing]
http://volsung.gotdns.org/gem-website/python_intro.html

A concise introductory overview with links to related sites, from the
point of view of using Python for scientific computing/numerical
methods.

[Python Complex Number Arithmetic]
http://faculty.ed.umuc.edu/~swalsh/Python/PCNA.html

As an interpreted language it is often unnecessary to declare
variables, but this raises the question of how the interpreter should
know whether real or complex arithmetic is meant.  An example here
illustrates writing complex numbers using "j" as the imaginary square
root of -1 to implicitly signal that a complex value is intended, i.e.
4 + 0j.  Here "j" functions as a modifier rather than as an
independent name, sort of like using ".0" to force an integer to be
treated as a floating point value.

[5.8 cmath -- Mathematical functions for complex numbers]
http://docs.python.org/lib/module-cmath.html

The "base" module cmath in Python does not have "convenience" function
cbrt for complex arguments, though it does implement sqrt for complex
numbers (and base module math implements cbrt for real ones).  I'll
discuss a workaround in the Answer posted above.

regards, mathtalk-ga
Subject: Re: for mathtalk-ga only: Need more help with cubic equations
From: donkan-ga on 12 Dec 2004 16:02 PST
 
And here's a Python function, written by Tim Peters, for computing the
3 cube roots of a complex number:
<http://www.rcblue.com/Python/croots.txt>

donkan

Important Disclaimer: Answers and comments provided on Google Answers are general information, and are not intended to substitute for informed professional medical, psychiatric, psychological, tax, legal, investment, accounting, or other professional advice. Google does not endorse, and expressly disclaims liability for any product, manufacturer, distributor, service or service provider mentioned or any opinion expressed in answers or comments. Please read carefully the Google Answers Terms of Service.

If you feel that you have found inappropriate content, please let us know by emailing us at answers-support@google.com with the question ID listed above. Thank you.
Search Google Answers for
Google Answers  


Google Home - Answers FAQ - Terms of Service - Privacy Policy